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;