MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: matlab code to mma?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg52202] Re: matlab code to mma?
  • From: sean_incali at yahoo.com (sean kim)
  • Date: Sun, 14 Nov 2004 04:31:20 -0500 (EST)
  • References: <cmvf41$smp$1@smc.vnet.net> <cn4ln3$16e$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Hi Daniel. 

Thanks for your help with this. 

But the calpropensity function has two function definitions in it. and
Mathematica seems to be having problems with the inner function having the
argument v0_ (as it is used in out function as an argument)

How do I fix this? 

thanks all for any insights. 


sean 


PS. I can now run the matlab file ( it took me a while to figure out )




D Herring <dherring at at.uiuc.dot.edu> wrote in message news:<cn4ln3$16e$1 at smc.vnet.net>...
> Hi all, (message also sent directly to Sean)
> 
> I think the following is a rather faithful, line-by-line translation of 
> Sean's Matlab code to Mathematica code.  It is, however, completely 
> untested.
> 
> Later,
> Daniel
> 
> 
> (* Note: Matlab doesn't automatically check the number or type of 
> arguments being passed to a function.  This is left to the function's 
> author. At run time, 'nargin' contains actual number of arguments 
> passed; this is often used to provide default arguments. *)
> 
> gillespie::usage="gillespie[tf_, v0_]\n
>   runs the Gillespie algorithm for constitutive expression up to
>   a time tf. the value of v0 may be specified but has a default value
>   of 0.01";
> 
> gillespie[tf_,v0_:0.01]:=
> Block[{calpropensities,},
> 	(* Local function *)
> 	calpropensities[s_,v0_]:=
> 	Block[{v1,d0,d1,D,M,N,a},
> 		(* rates *)
> 		v1=0.04;
> 		d0=Log[2]/180;
> 		d1=Log[2]/3600;
> 
> 		(* conversion *)
> 		D=s[[1]];
> 		M=s[[2]];
> 		N=s[[3]];
> 
> 		(* propensities *)
> 		a=Table[0, {4}];
> 		a[[1]]=v0*D;
> 		a[[2]]=d0*M;
> 		a[[3]]=v1*M;
> 		a[[4]]=d1*N;
> 
> 		Return[a];
> 	];
> 
> 	(* Initial amounts *)
> 	D=1;
> 	M=0;
> 	N=0;
> 	t=0;
> 
> 	(* Initialize other variables *)
> 	j=2;
> 	nr=4; (* number of reactions *)
> 
> 	(* Matlab-specific note:
> 		create some large result matrix
> 		(this is not essential but speeds up the program) *)
> 	result=Table[0,{1000},{4}];
> 
> 	(* store data for first time point *)
> 	result[[1,All]]={t, D, M, N};
> 
> 	While[t<tf,
> 		(* conversion *)
> 		s={D,M,N};
> 
> 		(* store data *)
> 		result[[j,All]]={t,D,M,N};
> 		j=j+1;
> 
> 		(* calculate propensities *)
> 		a=calpropensities[s,v0];
> 		a0=Total[a];
> 
> 		(* generate two random numbers *)
> 		r=Table[Random[],{2}];
> 
> 		(* calculate time of next reaction *)
> 		dt=Log[1/r[[1]]]/a0;
> 
> 		(* calculate which reaction is next *)
> 		Do[
> 			If[Total[Take[a,mu]] > r[[2]]*a0,
> 				Break[];
> 			];
> 			,{mu, nr}
> 		];
> 
> 		(* carry out reaction *)
> 		t=t+dt;
> 		Which[
> 			mu==1, M=M+1,
> 			mu==1, M=M-1,
> 			mu==3, N=N+1,
> 			mu==4, N=N-1
> 		];
> 	];
> 
> 	(* store data for last time point *)
> 	result[[j,All]]={t,D,M,N};
> 
> 	(* remove any excess zeros from result *)
> 	result=Drop[result, j-Length[result]];
> ];
> 
> Steven M. Christensen wrote:
> > [Please send comments on this directly to Sean - moderator]
> > 
> >  i have received some codes written in matlab, and i would like to
> >  convert it to Mathematica.
>  ....
> >  function result= gillespie(tf, v0)
> >  % function result= gillespie(tf, v0)
> >  %
> >  % runs the Gillespie algorithm for constitutive expression up to 
> >  % a time tf. the value of v0 may be specified but has a default value
> >  % of 0.01
> >  
> >  if nargin == 1
> >      v0= 0.01;
> >  end;
> >  
> >  % initial amounts
> >  D= 1;
> >  M= 0;
> >  N= 0;
> >  t= 0; 
> >  % initialize other variables
> >  j= 2;
> >  nr= 4;              % number of reactions
> >  
> >  % create some large result matrix (this is not essential but
> >  % speeds up the program)
> >  result= zeros(100000, 4);
> >  
> >  % store data for first time point
> >  result(1,:)= [t D M N];
> >  
> >  while t < tf
> >      
> >      % conversion
> >      s(1)= D;
> >      s(2)= M;
> >      s(3)= N;
> >      
> >      % store data
> >      result(j,:)= [t D M N];
> >      j= j+1;
> >          
> >      % calculate propensities
> >      a= calpropensities(s, v0);
> >      a0= sum(a);
> >      
> >      % generate two random numbers
> >      r= rand(1,2);
> >      
> >      % calculate time of next reaction
> >      dt= log(1/r(1))/a0;
> >      
> >      % calculate which reaction is next
> >      for mu= 1:nr
> >          if sum(a(1:mu)) > r(2)*a0
> >              break;
> >          end;
> >      end;
> >      
> >      % carry out reaction
> >      t= t+dt;
> >      if mu == 1
> >          M= M+1;
> >      elseif mu == 2
> >          M= M-1;
> >      elseif mu == 3
> >          N= N+1;
> >      elseif mu == 4
> >          N= N-1;
> >      end;
> >      
> >  end;
> >  
> >  % store data for last time point
> >  result(j,:)= [t D M N];
> >  
> >  % remove any excess zeros from result
> >  result(j+1:end,:)= [];
> >  
> >  
> >  
> > 
> > 
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> > 
> >  function a= calpropensities(s, v0)
> >  
> >  % rates
> >  v1= 0.04;
> >  d0= log(2)/180; 
> >  d1= log(2)/3600;
> >  
> >  % conversion
> >  D= s(1);
> >  M= s(2);
> >  N= s(3);
> >  
> >  % propensities
> >  a(1)= v0*D;
> >  a(2)= d0*M;
> >  a(3)= v1*M;
> >  a(4)= d1*N;


  • Prev by Date: Re: equal distribution of last digits base ten in the primes by b-normality
  • Next by Date: Re: Counting Runs
  • Previous by thread: Re: matlab code to mma?
  • Next by thread: scroll problem with mouse