Package com.nr.ode

Source Code of com.nr.ode.Stochsim

package com.nr.ode;

import static com.nr.NRUtil.*;
import static java.lang.Math.*;
import com.nr.la.NRsparseCol;
import com.nr.ran.Ran;

/**
* stochastic simulation of ODEs
* Copyright (C) Numerical Recipes Software 1986-2007
* Java translation Copyright (C) Huang Wen Hui 2012
*
* @author hwh
*
*/
public class Stochsim {
  public double[] s;
  public double[] a;
  public double[][] instate, outstate;
  public NRsparseCol[] outchg, depend;
  public int[] pr;
  public double t, asum;
  private Ran ran;

  // begin user section
  static final int mm=3;
  static final int nn=4;

  public double k0,k1,k2;

  public Stochsim(final double[] sinit) {
    this(sinit, 1);
  }
 
  public Stochsim(final double[] sinit, final int seed) {
    s = sinit;
    a = new double[mm];
    outchg = new NRsparseCol[mm];
    depend = new NRsparseCol[mm];
    for(int i=0;i<mm;i++) {
      outchg[i] = new NRsparseCol();
      depend[i] = new NRsparseCol();
    }
    pr = new int[mm];
    t =0;
    asum=0;
    ran = new Ran(seed);
   
    int i,j,k,d;
    describereactions();
    sparmatfill(outchg,outstate);
    double[][] dep = new double[mm][mm];
    for (i=0;i<mm;i++) for (j=0;j<mm;j++) {
      d = 0;
      for (k=0;k<nn;k++)
        d = d | ((int)instate[k][i] & (int)outstate[k][j]);
        //d = d || ((int)instate[k][i] && (int)outstate[k][j]);
      dep[i][j] = d;
    }
    sparmatfill(depend,dep);
    for (i=0;i<mm;i++) {
      pr[i] = i;
      if(i==0)
        a[i] = rate0();
      else if(i==1)
        a[i] = rate1();
      else
        a[i] = rate2();
      asum += a[i];
    }
  }
 
  private double rate0() {return k0*s[0]*s[1];}
  private double rate1() {return k1*s[1]*s[2];}
  private double rate2() {return k2*s[2];}
 
  public void describereactions () {
    k0 = 0.01;
    k1 = .1;
    k2 = 1.;
    double indat[] = {
      1.,0.,0.,
      1.,1.,0.,
      0.,1.,1.,
      0.,0.,0.
    };
    instate = buildMatrix(nn,mm,indat);
    double outdat[] = {
      -1.,0.,0.,
      1.,-1.,0.,
      0.,1.,-1.,
      0.,0.,1.
    };
    outstate = buildMatrix(nn,mm,outdat);
  }
  // end user section


  public double step() {
    int i,n,m,k=0;
    double tau,atarg,sum,anew;
    if (asum == 0.) {t *= 2.; return t;}
    tau = -log(ran.doub())/asum;
    atarg = ran.doub()*asum;
    sum = a[pr[0]];
    while (sum < atarg) sum += a[pr[++k]];
    m = pr[k];
    if (k > 0){
      swap(pr,k,k-1);
    }
    if (k == mm-1) asum = sum;
    n = outchg[m].nvals;
    for (i=0;i<n;i++) {
      k = outchg[m].row_ind[i];
      s[k] += outchg[m].val[i];
    }
    n = depend[m].nvals;
    for (i=0;i<n;i++) {
      k = depend[m].row_ind[i];
      if(k==0)
        anew = rate0();
      else if(k==1)
        anew = rate1();
      else
        anew = rate2();
      asum += (anew - a[k]);
      a[k] = anew;
    }
    if (t*asum < 0.1)
      for (asum=0.,i=0;i<mm;i++) asum += a[i];   
    return (t += tau);
  }
 
  public static void sparmatfill(final NRsparseCol[] sparmat, final double[][] fullmat) {
    int n,m,nz,nn=fullmat.length,mm=fullmat[0].length;
    if (sparmat.length != mm) throw new IllegalArgumentException("bad sizes");
    for (m=0;m<mm;m++) {
      for (nz=n=0;n<nn;n++) if (fullmat[n][m]!=0) nz++;
      sparmat[m].resize(nn,nz);
      for (nz=n=0;n<nn;n++) if (fullmat[n][m]!=0) {
        sparmat[m].row_ind[nz] = n;
        sparmat[m].val[nz++] = fullmat[n][m];
      }
    }
  }
}
TOP

Related Classes of com.nr.ode.Stochsim

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.