Page MenuHomec4science
Files F120433188

deltaNotch3.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Jul 4, 08:31

deltaNotch3.cpp
View Options

#include "deltaNotch3.hpp"
#include "Matrix.h"
void DeltaNotch3::operator() (const state_type &y , state_type &dydt , const double t)
{
uint nbcells = lattice-> getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
uint nbond = ce.getNbBoundaries();
double dsi = 0;
uint sind = ce.stateVectorIndex();
double si = y[sind];
double bd = getDeltaProduction(si/nbond);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint ind = b->stateVectorIndex(); // index for cell state
Boundary *bo= b->opposite();
double ddi,dni,dnidi,dnidj;
double dj,nj,njdi;
// cell state
double di=y[ind];
double ni=y[ind+1];
double nidi=y[ind+2];
double nidj=y[ind+3];
// values from the other side of the boundary if(bo!=NULL){
if(bo!=NULL){
uint nind =bo->stateVectorIndex();// index for neigbouring cell
dj = y[nind];
nj = y[nind+1];
njdi = y[nind+3];
}
else{// no neigboring cell
dj=nj=njdi=0.0;
}
//values from previous and next boundaries
double np = y[b->prev()->stateVectorIndex()+1];
double nn = y[b->next()->stateVectorIndex()+1];
//// delta
ddi= bd -gd*di; // production and degr.
ddi += ktm*njdi-ktp*di*nj; // trans binding
ddi += ks*njdi+kcs*nidi; //release after cleaving (cis and trans)
ddi += -kcp*ni*di+kcm*nidi; //// cis binding
//// notch
dni = bn -gn*ni;
dni += ktm*nidj -ktp*dj*ni; // trans binding
dni += - kcp*ni*di+kcm*nidi; // cis binding
dni += diffu*(np+nn-2*ni); //diffusion
//// cis complex
dnidi = kcp*ni*di - (kcm+ki+kcs)*nidi;
//// trans complex
dnidj = -ktm*nidj + ktp*dj*ni - ks*nidj;
//// cleaved notch
dsi += ks*nidj + kcs*nidi;
////
dydt[ind] = ddi;
dydt[ind+1] = dni;
dydt[ind+2] = dnidi;
dydt[ind+3] = dnidj;
}
dsi += -gs*si;
dydt[sind] = dsi;
}
}
// typedef double T;
// bool DeltaNotchAxis::InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
// using namespace boost::numeric::ublas;
// typedef permutation_matrix<std::size_t> pmatrix;
// // create a working copy of the input
// matrix<T> A(input);
// // create a permutation matrix for the LU-factorization
// pmatrix pm(A.size1());
// // perform LU-factorization
// int res = lu_factorize(A,pm);
// if( res != 0 )
// return false;
// // create identity matrix of "inverse"
// inverse.assign(ublas::identity_matrix<T>(A.size1()));
// // backsubstitute to get the inverse
// lu_substitute(A, pm, inverse);
// return true;
// }
Vec2 DeltaNotchAxis::cellAxis(Cell& c, const state_type& state){
Vec2 axis(0,1);
Vector y;
Matrix X(c.getNbBoundaries(),5),Xt,Xi;
BoundaryIterator it(&c);
uint cnt=0;
/*
for(Boundary *b=it.first();!it.isLast();b=it.next()){
//cout<<"b "<<b->center()<<endl;
// cout<<"c "<<c.center()<<endl;
Vec2 dir = b->center()-c.center();
double x1=dir[0];
double x2= dir[1];
X(cnt,1)=x1*x1;
X(cnt,2)=x2*x2;
X(cnt,3)=x1*x2;
X(cnt,4)=x1;
X(cnt,5)=x2;
y(cnt) = state[b->stateVectorIndex()+3];
cnt++;
}
*/
X.Transpose(Xt);
Xi = (Xt*X).Inverse();
if(Matrix::IsInverseOk()){
Vector beta = (Xi*Xt)*y;
Matrix hess(2,2,false);
Vector eval;// eigenvalues
Matrix evec; //eigenvectors
hess(0,0)=beta[0];hess(0,1)=hess(1,0)=beta[2];hess(2,2)=beta[1];
hess.TriEigen(eval,evec);
axis = evec.GetColumn(0);
axis.Normalize();
double diff = MAX(0,MAX(0,eval[0])-MAX(0,eval[1]));
diff = diff*diff;
double fac = diff/(1+diff);
axis *= fac;
}
else{
axis[0]=axis[1]=0.0;
}
return axis;
}
// matrix<double> X;
// matrix<double> Xi;
// vector<double> y;
// vector<double> b;
// matrix<double> M;
// bool inverted = InvertMatrix(prod(trans(X),X),Xi);
// beta = prod(prod(Xi,trans(X)),y);
// B
void DeltaNotchAxis::operator() (const state_type &y , state_type &dydt , const double t){
uint nbcells = lattice-> getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
uint nbond = ce.getNbBoundaries();
double dsi = 0;
uint sind = ce.stateVectorIndex();
double si = y[sind];
double bd = getDeltaProduction(si/nbond);
Vec2 axis = cellAxis(ce,y);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint ind = b->stateVectorIndex(); // index for cell state
Boundary *bo= b->opposite();
double ddi,dni,dnidi,dnidj;
double dj,nj,njdi;
// cell state
double di=y[ind];
double ni=y[ind+1];
double nidi=y[ind+2];
double nidj=y[ind+3];
// values from the other side of the boundary
if(bo!=NULL){
uint nind =bo->stateVectorIndex();// index for neigbouring cell
dj = y[nind];
nj = y[nind+1];
njdi = y[nind+3];
}
else{// no neigboring cell
dj=nj=njdi=0.0;
}
//values from previous and next boundaries
double np = y[b->prev()->stateVectorIndex()+1];
double nn = y[b->next()->stateVectorIndex()+1];
// axial modulation of delta
Vec2 vv(0,0);
//b->center()-ce.center() ;
double bdb = bd*(1-axis.Norm2()*axis.cosDoubleAngle(
vv));
//b->center()-ce.center()));
// double bdb = bd;
//// delta
ddi= bdb -gd*di; // production and degr.
ddi += ktm*njdi-ktp*di*nj; // trans binding
ddi += ks*njdi; //release after cleaving
ddi += -kcp*ni*di+kcm*nidi; //// cis binding
//// notch
dni = bn -gn*ni;
dni += ktm*nidj -ktp*dj*ni; // trans binding
dni += - kcp*ni*di+kcm*nidi; // cis binding
dni += diffu*(np+nn-2*ni); //diffusion
//// cis complex
dnidi = kcp*ni*di - kcm*nidi - ki*nidi;
//// trans complex
dnidj = -ktm*nidj + ktp*dj*ni - ks*nidj;
//// cleaved notch
dsi += ks*nidj;
////
dydt[ind] = ddi;
dydt[ind+1] = dni;
dydt[ind+2] = dnidi;
dydt[ind+3] = dnidj;
}
dsi += -gs*si;
dydt[sind] = dsi;
}
}

Event Timeline

c4science · Help