Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/runDeltaNotch.cpp
b7305f7b5a8bmaster

Last Change
Raw File

runDeltaNotch.cpp
View Options

#include <iostream>
#include <fstream>
#include "matrix.h"
#include "mex.h"
//#include "lattice.hpp"
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "deltaNotch4.hpp"
#define RECORD_TRAJ
using namespace std;
using namespace boost::numeric::odeint;
typedef runge_kutta_dopri5< double , double , double , double , vector_space_algebra > stepper_type;
typedef runge_kutta_cash_karp54< state_type , double > error_stepper_type;
void mexFunction(int nlhs,mxArray* plhs[],int nrhs, const mxArray *prhs[]){
const mxArray *cells = prhs[0];
const mxArray *bonds = prhs[1];
const mxArray *nbonds = prhs[2];
const mxArray *par = prhs[3];
const mxArray *verts = prhs[4];
const mxArray *inhib = prhs[5];
const mxArray *deltaprod = prhs[6];
double *params = mxGetPr(par);
uint ncells = mxGetN(cells);
int maxbonds = mxGetM(cells);
int totbonds = mxGetN(bonds);
int nbparams = mxGetN(par);
int nverts = mxGetN(verts);
uint steps =0;
// cout<<"c running..."<<params[0]<<endl;
Lattice lat(mxGetPr(cells), mxGetPr(bonds), ncells, mxGetPr(nbonds),totbonds,maxbonds, mxGetPr(verts), nverts);
srand(time(NULL));
DeltaNotch2 dncom;
dncom.setLattice(&lat);
dncom.initState();
dncom.writeState();
state_type& y0 = dncom.getState();
#ifdef RECORD_TRAJ
vector< state_type >states;
vector< double > times;
vector<uint> cidx;
lat.forAllCellsAndBoundaries(DeltaNotch2::collectCellIdx,(void *)&cidx,NULL,NULL);
push_back_sel_state_and_time obs(states,times,cidx);
//push_back_state_and_time obs(states,times);
#else
KeepLast& obs =dncom.observeLast();
#endif
/*
if(nbparams>=1)dncom.setDiffusion(params[0]);
if(nbparams>=2)dncom.setNotchProduction(params[1]);
if(nbparams>=3)dncom.setDeltaProduction(params[2]);
if(nbparams>=3)dncom.setNotchCleavageRate(params[3]);
// if(nbparams>=4)dncom.setCisNotchCleavageRate(params[4]);
if(nbparams>=5)dncom.setTransBindingRate(params[5]);
if(nbparams>=6)dncom.setTransUnbindingRate(params[6]);
if(nbparams>=7)dncom.setCisBindingRate(params[7]);
if(nbparams>=8)dncom.setCisUnbindingRate(params[8]);
*/
// uint steps = integrate_const(dncom,y,0.0,2.0,2.0);
steps = integrate_adaptive(make_controlled<error_stepper_type>(0.001,0.01),dncom,y0,0.0,2500.0,1.0,obs);
cout<<endl<<"integration done in "<<steps<<" steps"<<endl;
state_type& y = obs.getState();
dncom.writeState();
cout<<"n out args: "<<nlhs<<endl;
plhs[0] = mxCreateDoubleMatrix(ncells,1,mxREAL);
plhs[1] = mxCreateDoubleMatrix(ncells,1,mxREAL);
double *delta = mxGetPr(plhs[0]);
double *notch = mxGetPr(plhs[1]);
// if(mxGetN(delta)*mxGetM(delta)<n){
// mxFree(plhs[0]);//needed?
// plhs[0] = mxCreateNumericMatrix(n, 1, mxDOUBLE_CLASS, mxREAL);
// delta = mxGetPr(plhs[0]);
// }
// if(mxGetN(notch)*mxGetM(notch)<n){
// mxFree(plhs[1]); //needed ?
// plhs[1] = mxCreateNumericMatrix(n, 1, mxDOUBLE_CLASS, mxREAL);
// notch = mxGetPr(plhs[1]);
// }
for(uint i=0;i<ncells;i++){
uint idx = lat.getCell(i).stateVectorIndex();
cout<<idx<<" ";
delta[i] = dncom.getDeltaProduction(y0[idx]);
}
#ifdef RECORD_TRAJ
ofstream fout("traj.txt");
for(uint i=0; i<=steps; i+=5 ){
for(uint j=0;j<states[i].size();j++){
fout<< states[i][j]<<" ";
}
fout<<endl;
}
fout.close();
#endif
cout<<"quitting c..."<<endl;
}
c4science · Help