       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/180;
d1=Log/3600;

(* conversion *)
D=s[];
M=s[];
N=s[];

(* propensities *)
a=Table[0, {4}];
a[]=v0*D;
a[]=d0*M;
a[]=v1*M;
a[]=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[]]/a0;

(* calculate which reaction is next *)
Do[
If[Total[Take[a,mu]] > r[]*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:
>
>  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