SteadyStateApproximation/sinko_existence_region.py at master · MathBioCU/SteadyStateApproximation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
# -*- coding: utf-8 -*-
#Created on Feb 18, 2016
#@author: Inom Mirzaev

"""
    This code generates existence region for the steady states of the famous Sinko-Streifer model.
    Infinitesimal generator G is apporixmated by an n-by-n matrix G_n. Consequently,
    steady states of G is approximated by zeros of the matrix G_n.

    Computed regions are plotted in 3D and saved in 'images' folder.

    The model rates should be specified in the 'sinko_model_rates.py' file.

    The program has been written in parallel. Therefore, for the faster
    computation, the parameter 'ncpus' in 'sinko_model_rates' should
    be set to maximum number of cores available.
"""

from __future__ import division
from sinko_model_rates import *


import multiprocessing as mp
import numpy as np
import time , os


start = time.time()


#Given a singular matrix this function returns nullspace of that matrix
def null(a, rtol=1e-5):

    u, s, v = np.linalg.svd(a)
    rank = (s > rtol*s[0]).sum()

    return v[rank:].T.copy()


x = np.ravel( grid_x )
y = np.ravel( grid_y )
z = np.ravel( grid_z )

myarray = np.array([x, y, z]).T


#Initialize approximate matrice
Renewal_mat , Growth_mat , Removal_mat , nu , N , dx  = sinko_initialization( 100 , 1 , 1,  1 )


def region_plots( nn , Renewal_mat=Renewal_mat,
                       Growth_mat=Growth_mat,
                       Removal_mat=Removal_mat,
                       myarray = myarray):

    An =  myarray[nn , 0]*Renewal_mat + myarray[nn , 1]*Growth_mat  + myarray[nn, 2]*Removal_mat

    pos_sol = 0
    eigs = 0

    #if dimension of the nullspace is nonzero. Return positive steady state exists
    if np.sum( np.abs( null(An) ) ) > 0:
        pos_sol = 1
        eigs = np.max(  np.real ( np.linalg.eig(  An )[0] ) )


    return ( myarray[nn , 0] , myarray[nn , 1] , myarray[nn, 2] , pos_sol , eigs )


if __name__ == '__main__':

    #Number of CPUs to be used
    pool = mp.Pool( processes = ncpus )
    ey_nana = range( len( myarray) )
    result = pool.map( region_plots , ey_nana )

    #The output is saved in the data_files folder
    output = np.asarray(result)
    output = output[ np.nonzero( output[: , 3 ] ) ]
    fname = 'sinko_data'
    np.save( os.path.join( 'data_files' , fname ) , output )


from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import matplotlib.tri as mtri

fname = 'sinko_data.npy'
output=np.load( os.path.join( 'data_files' , fname ) )
output = output[ np.nonzero( output[: , 3 ] ) ]

points = output[ :, 0:3]
values = output[ : , -1 ]


eigs = griddata( points , values , ( grid_x , grid_y , grid_z ) )

out = np.array([np.ravel(grid_x) , np.ravel(grid_y)  , np.ravel(grid_z) , np.ravel(eigs)] ).T
mypts = out[ np.nonzero( np.isnan(out[:, 3] )==False )[0] ]


"""
    Plots the existence region for Sinko-Streifer population model
"""
plt.close('all')

fig = plt.figure(0)
ax = fig.add_subplot(111, projection='3d')

tri = mtri.Triangulation( mypts[ : , 0 ] , mypts[ : , 1] )

ax.plot_trisurf(mypts[:, 0] , mypts[:, 1] , mypts[:, 2] ,
                triangles = tri.triangles , edgecolor='none', cmap='jet', shade=True)


ax.view_init( azim=100 , elev=26  )
ax.set_xlabel( '$a$' , fontsize=20 )
ax.set_ylabel( '$b$' , fontsize=20 )
ax.set_zlabel( '$c$' , fontsize=20 )
ax.set_xlim( amin , amax )
ax.set_ylim( bmin , bmax )
ax.set_zlim( cmin , cmax )

plt.savefig( os.path.join( 'images' , 'sinko_exist_region.png' ) , dpi=400)


end = time.time()


print "Elapsed time", round( end - start   , 2 ) ,  "seconds "