CnrLwlss

import numpy as np
import math
import matplotlib.pyplot as plt
import webbrowser

def makePlot(xvals,yvals,fname="harmonograph.png",fopen=False,w=50,h=50,width=0.1):
    maxs=[np.max(np.abs(xvals)),np.max(np.abs(yvals))]
    farthest=max(maxs)

    drw=plt.figure(1)

CnrLwlss

#install.packages(c("mvtnorm","fields"))
library(fields)
library(mvtnorm)

grmf=function(nrow,ncol,sigsq=1,V=1){
	pos=expand.grid(row=1:nrow,col=1:ncol)
	mat=matrix(nrow=nrow,ncol=ncol,0)

	dists=rdist(pos)
	#dists[row(dists)==col(dists)]=1

CnrLwlss

count_DDY835=c(10,13,9,7,0,15,7,0,7,0,0,10,9,5,0,10,12,6,9,8,0,10,10,0,11,0,0,9,0,0,0,0)
count_DDY836=c(15,11,8,17,0,10,7,0,7,0,0,12,8,8,0,9,12,13,10,12,0,13,10,0,8,0,0,10,0,0,0,0)
gtype=c("wt","nmd2D","exo1D","rad24D","cdc13D","nmd2_exo1D","nmd2D_rad24D","nmd2D_cdc13D",
"exo1D_rad24D","exo1D_cdc13D","rad24D_cdc13D","nmd2D_exo1D_rad24D","nmd2D_exo1D_cdc13D","nmd2D_rad24D_cdc13D",
"exo1D_rad24D_cdc13D","nmd2D_exo1D_rad24D_cdc13D","rif1_wt","rif1_nmd2D","rif1_exo1D","rif1_rad24D","rif1_cdc13D",
"rif1_nmd2_exo1D","rif1_nmd2D_rad24D","rif1_nmd2D_cdc13D","rif1_exo1D_rad24D","rif1_exo1D_cdc13D","rif1_rad24D_cdc13D","rif1_nmd2D_exo1D_rad24D",
"rif1_nmd2D_exo1D_cdc13D","rif1_nmd2D_rad24D_cdc13D","rif1_exo1D_rad24D_cdc13D","rif1_nmd2D_exo1D_rad24D_cdc13D")

df=data.frame(gene=gtype,count_DDY835=count_DDY835,count_DDY836=count_DDY836,stringsAsFactors=FALSE)

CnrLwlss

# Problem with bedpe file format as generated by bedtools
# https://code.google.com/p/bedtools/issues/detail?id=152
#SRR364781 -> wt_sample
#SRR364782 -> wt_replicate
#SRR364783 -> pol32_sample
#SRR364784 -> pol32_replicate

#source("https://bioconductor.org/biocLite.R")
#biocLite("ShortRead")
#biocLite("org.Sc.sgd.db")

CnrLwlss

double MCMC_base_truncate_low(double truncate, struct_data *D,struct_para *D_para,struct_priors *D_priors,double *accept,double *h,double para,double (*foo)(struct struct_data *D,struct struct_para *D_para,struct struct_priors *D_priors,double,int,int),int l, int m){
  double logu,logaprob,can;
   //can=rnorm(para,*h); /*can=para+gsl_ran_gaussian(RNG,*h);*/
   // if(can<=(truncate)){
   //   can=para;
   // }

   can=truncate-1.0;
   while(can<=truncate){
	 can=rnorm(para,*h);  

CnrLwlss

simDSLogistic=function(K,r,N0){
  # Unusually, for this model, we know the number of events a priori
  eventNo=K-N0
  # So we can just generate all required random numbers (quickly) in one go
  unifs=runif(eventNo)
  # Every event produces one cell and consumes one unit of nutrients
  clist=(N0+1):K
  # Simulate time between events by generating 
  # exponential random numbers using the inversion method
  dts=-log(1-unifs)/(r*clist*(1-clist/K))

CnrLwlss

import numpy as np

def simDSLogistic(K,r,N0):
    '''Discrete stochastic logistic model simulation.  Carrying capacity K,
intrinsic growth rate r, initial population size (inoculum) N0.  Returns an
array with a (continuous) time column and a (discrete) population size column.'''
    # Unusually, for this model, we know the number of events a priori
    eventNo=K-N0
    # So we can just generate all required random numbers (quickly) in one go
    unifs=np.random.uniform(size=eventNo)

CnrLwlss

library(DEoptim)
library(rgenoud)
library(snow)

# CirclePacking.R from https://gist.github.com/4571810
source("CirclePacking.R")

# Number of circles and dimensions of bounding rectangle 
N=100; W=10; H=10
root="N64_100"

CnrLwlss

# Simulating n steps of a biased random walk
# starting from x0 and with decreasing steps 
# occurring with probability theta
randomWalk=function(n=100,x0=0,theta=0.5){
	x=x0
	res=array(x0,dim=n+1)
	unifs=runif(n)
	for(i in 0:(n-1)){
		if (unifs[i+1]<=theta){x=x-1}else{x=x+1}
		res[i+2]=x	

CnrLwlss

simCellsHybrid=function(K,r,N0,NSwitch,detpts=100){
  # Every event produces one cell and consumes one unit of nutrients
  if(NSwitch>N0){
    # Unusually, for this model, we know the number of events a priori
    eventNo=NSwitch-N0
    # So we can just generate all required random numbers (quickly) in one go
    unifs=runif(eventNo)
    clist=(N0+1):NSwitch
    # Time between events
    dts=-log(1-unifs)/(r*clist*(1-clist/K))