Page MenuHomec4science
Files F92038193

runDeltaNotchCell.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 16, 20:20

runDeltaNotchCell.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 "deltaNotchCell.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[]){
bool with_scale = (nrhs>6);
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 *state = prhs[5];
const mxArray *scale = with_scale?prhs[6]:NULL;
double *params = mxGetPr(par);
uint ncells = mxGetN(cells);
int maxbonds = mxGetM(cells);
int totbonds = mxGetN(bonds);
int nbparams = mxGetN(par);
uint nverts = mxGetN(verts);
double *scale_pt = with_scale?mxGetPr(scale):NULL;
uint steps =0;
// cout<<"c running..."<<params[0]<<" "<<ncells<<endl;
Lattice lat(mxGetPr(cells), mxGetPr(bonds), ncells, mxGetPr(nbonds),totbonds,maxbonds, mxGetPr(verts), nverts,scale_pt);
// cout<<"conversion done"<<endl;
srand(time(NULL));
DeltaNotchCell dncom;
// cout<<"setting lattice"<<endl;
dncom.setLattice(&lat);
dncom.initState(mxGetPr(state));
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
double duration = params[0];
state_type dy ;
steps = integrate_adaptive(make_controlled<error_stepper_type>(0.001,0.01),dncom,y0,0.0,duration,1.0,obs);
// cout<<endl<<"integration done in "<<steps<<" steps"<<endl;
state_type& y = obs.getState();
// dncom.writeState();
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]);
notch[i] = y0[idx];
delta[i] = y0[idx+1];
}
#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;
}

Event Timeline

c4science · Help