Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

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

Search the Archive

Re: matlab code to mma?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg52165] Re: matlab code to mma?
  • From: D Herring <dherring at at.uiuc.dot.edu>
  • Date: Sat, 13 Nov 2004 04:40:15 -0500 (EST)
  • References: <cmvf41$smp$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

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: Why these definitions are so slow
  • Previous by thread: Re: TraditionalForm of expressions WITHOUT evaluating.
  • Next by thread: Re: matlab code to mma?