Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/testDN4.cpp
b7305f7b5a8bmaster

Last Change
Raw File

testDN4.cpp
View Options

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include "deltaNotch4.hpp"
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 collectCellIdx(Cell *c,void *arg){
vector<uint> *v = (vector<uint> *)arg;
v->push_back(c->stateVectorIndex());
}
void collectECDIdx(Lattice *lat, vector<uint>& v){
Cell& ce = lat->getCell(1);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
v.push_back(b->extracellMatrix()->stateVectorIndex());
}
}
void collectFreeDelta(Lattice *lat, vector<uint>& v){
Cell& ce = lat->getCell(1);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
v.push_back(b->stateVectorIndex());
}
}
void collectOppFreeDelta(Lattice *lat, vector<uint>& v){
Cell& ce = lat->getCell(1);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
v.push_back(b->opposite()->stateVectorIndex());
}
}
int main(int argc, char *argv[]){
//Lattice lat(TWO_CELLS);
Lattice lat(SEVEN_CELLS);
DeltaNotch4 dn;
dn.setDefaultParameters();
dn.setLattice(&lat);
dn.initState();
dn.deltaUp(0);
vector<state_type> x_vec;
vector<double> times;
for(int i=1;i<argc;i++){
char *p;
p=strtok(argv[i],"=");
if(!strcmp(p,"bn")){
dn.setNotchProduction(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"bd")){
dn.setDeltaProduction(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"ktp")){
dn.setTransBindingRate(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"kcp")){
dn.setCisBindingRate(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"ks")){
dn.setNotchCleavageRate(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"diff")){
dn.setDiffusion(atof(strtok(NULL,"=")));
}
if(!strcmp(p,"gd")){
dn.setNotchDegradation(atof(strtok(NULL,"=")));
}
}
state_type& y0 = dn.getState();
cout<<"dimension "<<dn.getStateDimension()<<endl;
vector<uint> idx;
lat.forAllCellsAndBoundaries(collectCellIdx,(void *)&idx,NULL,NULL);
//collectECDIdx(&lat,idx);
// collectFreeDelta(&lat,idx);
// collectOppFreeDelta(&lat,idx);
uint steps = integrate_adaptive(make_controlled<error_stepper_type>(0.001,0.01),dn,y0,0.0,10.0,2.0, push_back_sel_state_and_time(x_vec,times,idx));
cout<<"done"<<endl;
// for(uint i=0; i<=steps; i++ ){
// cout << times[i] <<" ";
// for(uint j=0;j<x_vec[i].size();j++){
// cout<< x_vec[i][j]<<" ";
// }
// cout<<endl;
// }
for(uint i=0; i<=steps; i++ ){
cout << times[i] <<" ";
for(uint j=0;j<x_vec[i].size();j++){
cout<< x_vec[i][j]<<" ";
}
cout<<endl;
}
/*
Vec2 v1(1,2);
Vec2 v2(5,4);
v2-v1;
Vec2 v3 = v2-v1;
cout<<v3[0]<<" "<<v3[1]<<endl;
*/
return 0;
};
c4science · Help