Diffusion tisue modeling (master) Browse

runMechanicalForces2.cpp
b7305f7b5a8bmaster

runMechanicalForces2.cpp
View Options

			#define WITH_OBSERVER_MICHA
			#include <signal.h>
			#include <iostream>
			#include <fstream>
			#include "matrix.h"
			#include "mex.h"
			#include <boost/array.hpp>
			#include <boost/range/algorithm.hpp>
			#include <boost/numeric/odeint.hpp>

			#include "collierDoubleSignaling.hpp"
			// includes the cpp files containing template class methods
			#include "collierSignaling.cpp"
			#include "signalAndMechanicalForces.cpp"


			using namespace std;
			using namespace boost::numeric::odeint;


			void INThandler(int sig)
			{
			printf("Ctrl-C pressed\n");
			exit(0);
			}

			//SignalAndMechanicalForces <CollierSignaling <HillCoupling> > forces;
			SignalAndMechanicalForces <CellInteraction > forces;
			//SignalAndMechanicalForces <CollierDoubleSignaling> forces;

			//bool stop(const state_type &x){return forces.findSmallCell(x)!=-1;}
			bool stop(const state_type &x){return 0;}

			//typedef runge_kutta_dopri5< double , double , double , double , vector_space_algebra > stepper_type;
			typedef runge_kutta_cash_karp54< state_type , double > error_stepper_type;
			//typedef runge_kutta_fehlberg78 <state_type> error_stepper_type;
			//typedef rosenbrock4_controller< state_type , double > error_stepper_type;
			//typedef runge_kutta4< state_type> 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 = prhs[6];
			const mxArray *frozen_cells = prhs[7];
			const mxArray *bc = prhs[8];

			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 = mxGetPr(scale);

			uint steps =0;

			bool with_delamination= false;

			signal(SIGINT, INThandler);

			// 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;
			lat.freeze(mxGetPr(frozen_cells));
			// cout<<"freezing done"<<endl;
			if(*mxGetPr(bc) > 0.5){
			lat.setPeriodicBC();
			}

			srand(time(NULL));

			// MechanicalForces forces;
			forces.setLattice(&lat);
			//forces.setLattice(&lat,mxGetPr(state));

			// cout<<"lattice set"<<endl;
			double duration = params[0];

			// state_type& y = obs.getState();
			state_type y;
			//forces.initState(y);
			cout<<"initialization..."<<endl;
			forces.initState(y,mxGetPr(state));
			//forces.noMechanics();
			cout<<"initialization done "<<y.size()<<endl;
			KeepLast obs;
			obs.init(y.size());
			// steps = integrate(forces,y,0.0,duration,0.001,obs);
			auto stepper = make_controlled<error_stepper_type>(0.001,0.001);
			// steps = integrate_adaptive(stepper,forces,y,0.0,duration,0.01,obs);
			uint nmstep=10;
			double curtime=0.0;

			for (uint msteps=0;msteps<nmstep;msteps++){
			double mstep_dur = duration/nmstep*(msteps+1);
			while(mstep_dur-curtime>0.001){
			// cout<<"stop condition"<<endl;
			stop(y);
			cout<<"integrating..."<<endl;
			// if(msteps)
			auto iter= boost::find_if(make_adaptive_range(stepper,forces,y,curtime,mstep_dur,0.01,obs),stop);
			// forces.fillSignalVector(mxGetPr(plhs[4]));
			cout<<"integration done"<<endl;
			int smallcell = forces.findSmallCell(y);
			forces.fromState(y);
			if(smallcell >= 0 && with_delamination){
			cout<<"delaminating..." <<smallcell<<endl;
			forces.delaminate(smallcell);
			// lat.delaminate(smallcell);
			cout<<"T1s..."<<endl;
			// forces.findT1Transitions();// not very efficient...
			cout<<"reinit..."<<endl;
			forces.updateState(y);
			cout<<"cont"<<endl;
			}
			// cout<<"initializing..."<<endl;
			curtime = obs.getTime();
			cout<<"time: "<<curtime<<endl;
			}

			}
			// steps = integrate_adaptive(make_controlled<error_stepper_type>(0.001,0.001),forces,y,0.0,duration,0.01,obs);
			// stepper_type stepper;
			//steps = integrate_const(stepper,forces,y,0.0,duration,0.0001,obs);
			// cout<<"done in "<<iter<< " steps"<<endl;
			// y = obs.getState();
			// cout<<"preparing output... "<<totbonds<<" "<<lat.getNbBoundaries()<<endl;

			uint newmaxbonds= lat.maxNbBoundariesInCell();
			forces.updateSignal(y);

			// fill up return values
			plhs[0] = mxCreateDoubleMatrix(2*nverts,1,mxREAL);// vertices
			double *ver0= mxGetPr(plhs[0]);
			for(uint i=0;i<2*nverts;i++){ver0[i]= y[i];}

			plhs[1] = mxCreateDoubleMatrix(4,totbonds,mxREAL);// boundaries
			lat.getBoundaryMatrix(mxGetPr(plhs[1]));

			plhs[2] = mxCreateDoubleMatrix(newmaxbonds,lat.getNbCells(),mxREAL);// cells
			lat.getCellMatrix(mxGetPr(plhs[2]),newmaxbonds);

			plhs[3] = mxCreateDoubleMatrix(2,lat.getNbCells(),mxREAL);//cell centers
			lat.getCellCenterMatrix(mxGetPr(plhs[3]));

			uint nso = forces.getSignalingOutputNumber(); // signals
			for(uint sigo=0;sigo<nso;sigo++){
			pair<int,int> si = forces.getSignalingOutputSize(sigo);
			plhs[sigo+4] = mxCreateDoubleMatrix(si.first,si.second,mxREAL);
			forces.fillSignalingVector(mxGetPr(plhs[4+sigo]),sigo);
			}
			}

tismod/src/runMechanicalForces2.cppb7305f7b5a8bmaster

runMechanicalForces2.cppView Options

tismod/src/runMechanicalForces2.cpp
b7305f7b5a8bmaster

runMechanicalForces2.cpp
View Options