matlab to MMA, need some help converting this code.

From: sean kim (sean_incali_at_yahoo.com)
Date: 11/10/04 

Date: 10 Nov 2004 09:37:47 -0800

does anyone know both mathematica AND matlab?

if so can any of you help me out with this?

i have received some codes written in matlab, and i would like to
convert it to mma.

anyone up for the challenge *and* have the time and energy? It will be
nice to see line by line explanation of the conversion.

following is the code which I'm having problems understanding.

thanks for any insights

sean

-------

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;

Relevant Pages

  • Re: matlab to MMA, need some help converting this code.
    ... I don't much about matlab. ... >nice to see line by line explanation of the conversion. ... >% store data for first time point ... > % calculate which reaction is next ...
    (sci.math.symbolic)
  • Re: C to MATLAB from scratch.Where to start?
    ... My question is about the reverse conversion from .m to C/C++, ... Start from scratch. ... Matlab doesn't handle trivial program elemnts like for loops very well. ... m code to some extent is ambiguious, the automatic C compiler ...
    (comp.soft-sys.matlab)
  • Converting between java-matlab
    ... the running of my algorithm. ... In matlab, I can use two solutions to add the next element ... which requires no further matlab conversion. ... scans the vector but we all know how matlab cycles are NOT ...
    (comp.soft-sys.matlab)
  • Re: Converting sequential MATLAB code to a parallel form
    ... one of the high level goals of this conversion experiment is to have the code ... When you run with a matlabpool, we attempt to synchronise the MATLAB path ... between the desktop client and the workers. ... There are several ways that you can proceed - either you can place your M-code ...
    (comp.soft-sys.matlab)
  • Comparing CSVs in from 2 files
    ... I'm really new to Matlab so don't laugh if this question is ... so trivial that it makes you sick. ... matlab that compares 2 x,y,z csv data files to make sure they are ... If one of the values is wrong after going through the conversion I ...
    (comp.soft-sys.matlab)
Loading