Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87486004
deltaNotch3.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Oct 12, 23:21
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Oct 14, 23:21 (2 d)
Engine
blob
Format
Raw Data
Handle
21604180
Attached To
R9411 tisue modeling
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
Log In to Comment