SteadyStateApproximation/sinko_model_rates.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
133
134
135
136
137
138
139
140
# -*- coding: utf-8 -*-
#Created on Oct 14, 2015
#@author: Inom Mirzaev


"""
    Model rates and parameters used for generation of existence regions of the
    Sinko-Streifer population model (see Sinko, J. W. and Streifer, W. (1967),
    Ecology, 48(6):910-918. )
"""


from __future__ import division
from functools import partial


import numpy as np


"""
    Number of CPUs used for computation of existence and stability regions.
    For faster computation number of CPUs should be set to the number of cores available on
    your machine.
"""

ncpus = 2


# Minimum and maximum floc sizes
x0 = 0
x1 = 1


#Renewal rate
def renewal(x , a):
    #Should return a vector
    return a*(x + 1)


#Growth rate
def growth(x ,  b ):
    #Should return a vector
    return b*(x+1)


#Removal rate
def removal( x , c ):

     #Should return a vector
     return c * ( x + 1 - x )


"""
    Initialization of intervals for generation of existence and stability regions.
    For smoother plots more discretization points should be given.
"""

#Interval initialization renewal function

# a minimum
amin = 0

# a maximum
amax = 1

#Number of discretization  points in a interval
apart = 40

# b minimum
bmin = 0

# b maximum
bmax = 1

#Number of discretization points in b interval
bpart = 40

# c minimum
cmin = 0

# c maximum
cmax = 1

# Number of discretization points in c interval
cpart = 40


a_ = np.linspace( amax , amin , apart , endpoint=False ).tolist()
a_.sort()

b_ = np.linspace( bmax , bmin , bpart , endpoint=False ).tolist()
b_.sort()


c_ = np.linspace( cmax , cmin , cpart , endpoint=False ).tolist()
c_.sort()


grid_x, grid_y, grid_z = np.meshgrid(a_, b_, c_, indexing='ij')


#Initializes uniform partition of (x0, x1) and approximate operator F_n
def sinko_initialization(N , a, b, c, x1=x1 , x0=x0):

    #delta x
    dx = ( x1 - x0 ) / N

    #Uniform partition into smaller frames
    nu = x0 + np.arange(N+1) * dx


    #Initialize renewal function with parameter a
    renew = partial( renewal , a=a)

    #Initialize growth function with paramter b
    grow = partial( growth , b=b )

    #Initialize removal function with paramter c
    rem = partial( removal , c=c )

    #Removal operator
    Remov_mat =  - np.diag( rem( nu[range( 1 , N + 1 )] ) )

    #Growth operator
    Growth_mat = np.zeros( ( N , N ) )

    #Renewal operator
    Renew_mat = np.zeros( ( N , N ) )
    Renew_mat[0,:] = Renew_mat[0,:] + renew( nu[range( 1 , N+1 ) ] )

    for jj in range(N-1):
        Growth_mat[jj,jj] = -grow( nu[jj+1] ) / dx
        Growth_mat[jj+1,jj] = grow( nu[jj+1] ) / dx

    Growth_mat[N-1, N-1] = -grow( nu[N] ) / dx

    return ( Renew_mat , Growth_mat , Remov_mat , nu , N , dx)