1""" 2PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net) 3 4Copyright (C) 2004-2020 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved, 5 6Brett G. Olivier (bgoli@users.sourceforge.net) 7Triple-J Group for Molecular Cell Physiology 8Stellenbosch University, South Africa. 9 10Permission to use, modify, and distribute this software is given under the 11terms of the PySceS (BSD style) license. See LICENSE.txt that came with 12this distribution for specifics. 13 14NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK. 15Brett G. Olivier 16""" 17from __future__ import division, print_function 18from __future__ import absolute_import 19from __future__ import unicode_literals 20 21from pysces.version import __version__ 22__doc__ = '''PySCeS unittest module''' 23 24import os,sys,shutil 25import unittest 26import pysces # get rid of me 27from time import sleep 28from numpy.testing import assert_array_equal, assert_equal 29from numpy.testing import assert_almost_equal, assert_array_almost_equal 30from scipy import array, zeros, ones, logspace 31from scipy.optimize import fsolve 32from scipy.integrate import odeint 33 34def CopyTestModels(dirIn=os.path.join(pysces.install_dir,'pscmodels'),dirOut=pysces.model_dir,overwrite=0): 35 """ 36 CopyTestModels(dirIn=os.path.join(pysces_install,'pscmodels'),dirOut=pysces_model,overwrite=0) 37 38 'Copy all PySCeS test model files in dirIn to dirOut, defaults to: 39 in: pysces/pscmodels 40 out: pysces.model_dir 41 42 Arguments: 43 ========= 44 dirIn [default=os.path.join(pysces_install,'pscmodels')]: target directory 45 dirOut [default=pysces_model]: destination directory 46 overwrite [default=0]: automaitcally (1) overwrite target files 47 48 """ 49 if not os.path.exists(dirOut): 50 os.makedirs(dirOut) 51 52 if os.path.exists(dirIn): 53 if os.path.exists(dirOut): 54 print('src : ' + dirIn) 55 print('dest: ' + dirOut) 56 flist = os.listdir(dirIn) 57 if overwrite == 0: 58 flist2 = os.listdir(dirOut) 59 maxlen = 0 60 for x in flist: 61 if len(x) > maxlen: maxlen = len(x) 62 if len(flist) != 0: 63 for File in flist: 64 if File[:12] == 'pysces_test_' and File[-4:] == '.psc': 65 if overwrite == 0: 66 try: 67 a = flist2.index(File) 68 #print File + (maxlen-len(File))*'.' + ' skipped' 69 except: 70 shutil.copy(os.path.join(dirIn,File),os.path.join(dirOut,File)) 71 print(File + (maxlen-len(File))*'.' + ' ok') 72 else: 73 shutil.copy(os.path.join(dirIn,File),os.path.join(dirOut,File)) 74 print(File + (maxlen-len(File))*'.' + ' ok') 75 else: 76 print('Empty directory?') 77 else: 78 print(dirOut + ' does not exist') 79 else: 80 print(dirIn + ' does not exist') 81 82 83class PyscesTest: 84 '''PySCeS test suite: takes a test level as an argument''' 85 __version__ = __version__ 86 model_dir = os.path.join(pysces.model_dir,'tests') 87 88 def __init__(self,lvl=2,std2file=0): 89 #copy models from server to local model store 90 print('\nCopying pysces_test models if necessary ...') 91 CopyTestModels() 92 print('done.') 93 self.basic_runner = unittest.TextTestRunner() 94 self.BasicTest = unittest.makeSuite(PyscesBasicTest,'test') 95 self.ExtendedTest = unittest.makeSuite(PyscesExtendedTest,'test') 96 self.ExternalTest = unittest.makeSuite(PyscesExternalTest,'test') 97 if lvl > 3: 98 lvl = 3 99 100 class NullWriter: 101 def __init__(self,dres=0): 102 if dres: 103 self.Fout = open(os.path.join(os.getcwd(),'pysces_test_results.txt'),'w') 104 self.fof = 1 105 else: 106 self.fof = 0 107 def write(self,s): 108 if self.fof: 109 self.Fout.write(s) 110 111 def close(self): 112 if self.fof: 113 self.Fout.flush() 114 self.Fout.close() 115 116 self.__dfi__ = NullWriter(std2file) 117 tmpSTD = sys.stdout 118 for x in range(lvl): 119 print('\nLevel ' + str(x+1) + ' tests') 120 sys.stdout = self.__dfi__ 121 getattr(self, 'lvl_'+str(x+1))() 122 sys.stdout = tmpSTD 123 sleep(1) 124 self.__dfi__.close() 125 126 def lvl_1(self): 127 self.basic_runner.run(self.BasicTest) 128 129 def lvl_2(self): 130 self.basic_runner.run(self.ExtendedTest) 131 132 def lvl_3(self): 133 self.basic_runner.run(self.ExternalTest) 134 135class PyscesBasicTest(unittest.TestCase): 136 '''Basic test class, tests low-level numerical algorithms''' 137 __version__ = __version__ 138 model_dir = os.path.join(pysces.model_dir,'tests') 139 140 MathArrayTest = pysces.PyscesStoich.MathArrayFunc() 141 142 def test_swaprow_d(self): 143 arr = array([[1,2,3],[4,5,6],[7,8,9]],'d') 144 karr = array([[7,8,9],[4,5,6],[1,2,3]],'d') 145 arr = self.MathArrayTest.SwapRowd(arr,0,2) 146 assert_array_equal(arr,karr) 147 148 def test_swapcol_d(self): 149 arr = array([[1,2,3],[4,5,6],[7,8,9]],'d') 150 karr = array([[3,2,1],[6,5,4],[9,8,7]],'d') 151 arr = self.MathArrayTest.SwapCold(arr,0,2) 152 assert_array_equal(arr,karr) 153 154 def test_swaprow_z(self): 155 arr = array([[1,2,3],[4,5,6],[7,8,9]],'D') 156 karr = array([[7,8,9],[4,5,6],[1,2,3]],'D') 157 arr = self.MathArrayTest.SwapRowz(arr,0,2) 158 assert_array_equal(arr,karr) 159 160 def test_swapcol_z(self): 161 arr = array([[1,2,3],[4,5,6],[7,8,9]],'D') 162 karr = array([[3,2,1],[6,5,4],[9,8,7]],'D') 163 arr = self.MathArrayTest.SwapColz(arr,0,2) 164 assert_array_equal(arr,karr) 165 166 def test_matrixfloatfix(self): 167 arr = array([[1.0e-15,2,3],[4,-0.0,6],[7,8,1.0e-15]]) 168 karr = array([[0.0,2,3],[4,0.0,6],[7,8,0.0]]) 169 self.MathArrayTest.MatrixFloatFix(arr,val=1.0e-13) 170 assert_array_equal(arr,karr) 171 172 def test_FSOLVE_linear(self): 173 ds = zeros((3),'d') 174 s = zeros((3),'d') 175 def linear1_ode(s): 176 ds[0] = 100. - 6.*s[0] + s[1] 177 ds[1] = 5.*s[0] - 4.*s[1] + s[2] 178 ds[2] = 3.*s[1] - 3.*s[2] + 1. 179 return ds 180 181 res = fsolve(linear1_ode,[1.,1.,1.]) 182 known = [23.1025641, 38.61538462, 38.94871795] 183 assert_array_almost_equal(res,known) 184 185 ## def test_FSOLVE_moiety(self): 186 ## """This set of ODE's is numerically unstable and might fail without error""" 187 ## ds = zeros((3),'d') 188 ## s = zeros((3),'d') 189 ## def moibranch_ode(s): 190 ## s[2] = 1. - s[1] 191 ## ds[0] = s[1]/(2. + 2.*s[1]) - 2.*s[0]/(1. + s[0]) 192 ## ds[1] = 10.*s[2]/(2. + 2.*s[2]) - s[1]/(2. + 2.*s[1]) 193 ## ds[2] = s[1]/(2. + 2.*s[1]) - 10.*s[2]/(2. + 2.*s[2]) 194 ## return ds 195 ## res = fsolve(moibranch_ode,[0.01,0.01,0.01],maxfev=1000,xtol=1.e-6) 196 ## known = [0.1385856, 0.94882133, 0.05117867] 197 ## assert_array_almost_equal(res,known) 198 199 200 def test_ODEINT_linear(self): 201 ds = zeros((3),'d') 202 s = zeros((3),'d') 203 def linear1_ode(s,t): 204 ds[0] = 100. - 6.*s[0] + s[1] 205 ds[1] = 5.*s[0] - 4.*s[1] + s[2] 206 ds[2] = 3.*s[1] - 3.*s[2] + 1. 207 return ds 208 209 sim_res = odeint(linear1_ode,[1.,1.,1.],array([0,1,2,3,4,5,6,7,8,9,10])) 210 known = array([[ 1. , 1. , 1. ],\ 211 [ 20.70698909, 27.55450322, 20.9819434 ],\ 212 [ 22.4418045 , 35.47695146, 33.60751644],\ 213 [ 22.91227529, 37.71013683, 37.40388317],\ 214 [ 23.04761842, 38.35410015, 38.50268374],\ 215 [ 23.08668173, 38.53997688, 38.81993704],\ 216 [ 23.09809306, 38.59337738, 38.91160251],\ 217 [ 23.10124276, 38.60909619, 38.93798616],\ 218 [ 23.10217397, 38.61358525, 38.94561234],\ 219 [ 23.10245088, 38.61486396, 38.94781836],\ 220 [ 23.10253812, 38.61526501, 38.94851136]]) 221 222 assert_array_almost_equal(sim_res,known,3) 223 224 def test_ODEINT_moiety(self): 225 """This set of ODE's is numerically unstable and might fail without error""" 226 ds = zeros((3),'d') 227 s = zeros((3),'d') 228 def moibranch_ode(s,t): 229 s[2] = 1. - s[1] 230 ds[0] = s[1]/(2. + 2.*s[1]) - 2.*s[0]/(1. + s[0]) 231 ds[1] = 10.*s[2]/(2. + 2.*s[2]) - s[1]/(2. + 2.*s[1]) 232 ds[2] = s[1]/(2. + 2.*s[1]) - 10.*s[2]/(2. + 2.*s[2]) 233 return ds 234 235 sim_res = odeint(moibranch_ode,[0.7,1.,0.3],array([0,1,2,3,4,5,6,7,8,9,10])) 236 known = array([[ 0.7 , 1. , 0.3 ],\ 237 [ 0.43521812, 0.94587168, 0.05412832],\ 238 [ 0.21547096, 0.94879365, 0.05120635],\ 239 [ 0.15590709, 0.94882109, 0.05117891],\ 240 [ 0.14233313, 0.9488214 , 0.0511786 ],\ 241 [ 0.13938884, 0.94882111, 0.05117889],\ 242 [ 0.13875741, 0.94882111, 0.05117889],\ 243 [ 0.1386223 , 0.94882133, 0.05117867],\ 244 [ 0.1385936 , 0.94882134, 0.05117866],\ 245 [ 0.13858742, 0.94882133, 0.05117867],\ 246 [ 0.13858601, 0.94882133, 0.05117867]]) 247 248 assert_array_almost_equal(sim_res[-5:,:],known[-5:,:],3) 249 250class PyscesExtendedTest(unittest.TestCase): 251 '''Extended test class, tests modelling related methods''' 252 __version__ = __version__ 253 model_dir = os.path.join(pysces.model_dir,'tests') 254 255 def test_statemetab_linear1(self): 256 lin = pysces.model('pysces_test_linear1.psc', self.model_dir) 257 lin.doLoad() 258 lin.hybrd_mesg = 0 259 lin.State() 260 linmet = array([23.1025641, 38.61538462, 38.94871795],'d') 261 assert_array_almost_equal(lin.state_species,linmet) 262 263 def test_statemetab_branch1(self): 264 bra = pysces.model('pysces_test_branch1.psc', self.model_dir) 265 bra.doLoad() 266 bra.hybrd_mesg = 0 267 bra.State() 268 bramet = array([4.8583996,1.88547254,1.49124431,1.49124431],'d') 269 assert_array_almost_equal(bra.state_species,bramet) 270 271 def test_statemetab_moiety1(self): 272 moi = pysces.model('pysces_test_moiety1.psc', self.model_dir) 273 moi.doLoad() 274 moi.hybrd_mesg = 0 275 moi.State() 276 moimet = array([3.6886875,16.25569882,7.3113125,4.39229787,41.02504596,2.60770213,0.42718994,2.57281006,2.44791155,17.0012171],'d') 277 assert_array_almost_equal(moi.state_species,moimet) 278 279 def test_stateflux_linear1(self): 280 lin = pysces.model('pysces_test_linear1.psc', self.model_dir) 281 lin.doLoad() 282 lin.hybrd_mesg = 0 283 lin.State() 284 linflux = array([76.8974359, 76.8974359, 76.8974359, 76.8974359],'d') 285 assert_array_almost_equal(lin.state_flux,linflux) 286 287 def test_stateflux_branch1(self): 288 bra = pysces.model('pysces_test_branch1.psc', self.model_dir) 289 bra.doLoad() 290 bra.hybrd_mesg = 0 291 bra.State() 292 braflux = array([2.42139889,2.42139889,1.21069945,1.21069945,1.21069945,1.21069945],'d') 293 assert_array_almost_equal(bra.state_flux,braflux) 294 295 def test_stateflux_moiety1(self): 296 moi = pysces.model('pysces_test_moiety1.psc', self.model_dir) 297 moi.doLoad() 298 moi.hybrd_mesg = 0 299 moi.State() 300 moiflux = array([250.01825652,250.01825652,250.01825652,250.01825652,250.01825652,250.01825652,250.01825652],'d') 301 assert_array_almost_equal(moi.state_flux,moiflux) 302 303 def test_elas_linear1(self): 304 lin = pysces.model('pysces_test_linear1.psc', self.model_dir) 305 lin.doLoad() 306 lin.hybrd_mesg = 0 307 lin.State() 308 lin.EvalEvar() 309 line = [-0.30043348,0.,0., 1.50216739,-0.50216739,0., 0.,1.50650217,-0.50650217, 0.,0.,1.01300433] 310 line_new = [] 311 for x in range(lin.elas_var.shape[0]): 312 for y in range(lin.elas_var.shape[1]): 313 line_new.append(round(lin.elas_var[x,y],2)) 314 for x in range(len(line)): 315 line[x] = round(line[x],2) 316 317 assert_array_almost_equal(line,line_new) 318 319 def test_elas_branch1(self): 320 bra = pysces.model('pysces_test_branch1.psc', self.model_dir) 321 bra.doLoad() 322 bra.hybrd_mesg = 0 323 bra.State() 324 bra.EvalEvar() 325 brae = [-0.66930781448548349, 0.0, 0.0, 0.0, 0.78845903183743249,\ 326-0.52920041809870422, 0.0, 0.0, 0.0, 0.95441603882382564,\ 327-0.60578219088834051, 0.0, 0.0, 0.95441603882382564, 0.0,\ 328-0.60578219088834051, 0.0, 0.0, 0.9421058792599355, 0.0, 0.0,\ 3290.0, 0.0, 0.9421058792599355] 330 brae_new = [] 331 for x in range(bra.elas_var.shape[0]): 332 for y in range(bra.elas_var.shape[1]): 333 brae_new.append(round(bra.elas_var[x,y],2)) 334 for x in range(len(brae)): 335 brae[x] = round(brae[x],2) 336 assert_array_almost_equal(brae,brae_new) 337 338 def test_elas_moiety1(self): 339 moi = pysces.model('pysces_test_moiety1.psc', self.model_dir) 340 moi.doLoad() 341 moi.hybrd_mesg = 0 342 moi.State() 343 moi.EvalEvar() 344 moie = [1.4753672613878632,-0.47536726138786306,-0.47536726138786312,\ 3450.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.4278931518740323,0.0,\ 3461.4278931518740325,-0.42789315187403271,-0.42789315187403271,\ 3470.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0514524437327983,0.0,\ 3481.0514524437327983,-0.051452443732798468,0.0,0.0,0.0,0.0,\ 3490.0,-0.04300468618304179,0.0,1.0430046861830418,0.0,0.0,\ 350-0.043004686183041797,0.0,-0.073768363069393092,0.0,\ 3511.0737683630693931,0.0,0.0,0.0,0.0,0.0,1.0737683630693931,\ 3520.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.029048874655981143,\ 3531.0290488746559812,0.0,-0.029048874655981143,0.0,0.0,0.0,\ 3540.0,0.0,0.0,0.0,0.0,0.0, 1.0199985395848277] 355 moie_new = [] 356 for x in range(moi.elas_var.shape[0]): 357 for y in range(moi.elas_var.shape[1]): 358 moie_new.append(round(moi.elas_var[x,y],2)) 359 for x in range(len(moie)): 360 moie[x] = round(moie[x],2) 361 assert_array_almost_equal(moie,moie_new) 362 363 def test_cc_linear1(self): 364 lin = pysces.model('pysces_test_linear1.psc', self.model_dir) 365 lin.doLoad() 366 lin.hybrd_mesg = 0 367 lin.State() 368 lin.EvalEvar() 369 lin.EvalCC() 370 lincc=[0.02564102564102570300,0.05128205128205126600,\ 3710.15384615384615383000,0.76923076923076916000,0.02564102564102571000,0.05128205128205128000,\ 3720.15384615384615385000,0.76923076923076938000,0.02564102564102571000,0.05128205128205128000,\ 3730.15384615384615385000,0.76923076923076938000,0.02564102564102570300,0.05128205128205126600,\ 3740.15384615384615383000,0.76923076923076916000,-0.30636428644396779000,-0.61272857288793536000,\ 3750.15318214322198387000,0.76591071610991934000,-0.08534676570192656400,-0.17069353140385327000,\ 376-0.51208059421155983000,0.76812089131733974000,-0.96185074526088343000,0.05062372343478333000,\ 3770.15187117030435002000,0.75935585152175000000] 378 lincc_new = [] 379 for x in range(lin.cc_all.shape[0]): 380 for y in range(lin.cc_all.shape[1]): 381 lincc_new.append(round(lin.cc_all[x,y],2)) 382 for x in range(len(lincc)): 383 lincc[x] = round(lincc[x],2) 384 assert_array_almost_equal(lincc,lincc_new) 385 386 def test_cc_branch1(self): 387 bra = pysces.model('pysces_test_branch1.psc', self.model_dir) 388 bra.doLoad() 389 bra.hybrd_mesg = 0 390 bra.State() 391 bra.EvalEvar() 392 bra.EvalCC() 393 bracc=[0.25338970963029361000,-0.13797075277724413000,\ 394-0.21457061516904127000,0.32372624345386269000,0.39406892242342095000,0.38135649243870806000,\ 395-0.13797075277724413000,0.25338970963029356000,0.39406892242342095000,0.32372624345386258000,\ 396-0.21457061516904130000,0.38135649243870812000,-0.13797075277724413000,0.25338970963029356000,\ 3970.39406892242342095000,0.32372624345386258000,-0.21457061516904130000,0.38135649243870812000,\ 3980.05770947842652475500,0.05770947842652472700,0.08974915362718985400,0.32372624345386269000,\ 3990.08974915362718984000,0.38135649243870817000,0.25338970963029361000,-0.13797075277724413000,\ 400-0.21457061516904127000,0.32372624345386269000,0.39406892242342095000,0.38135649243870806000,\ 4010.05770947842652474100,0.05770947842652471300,0.08974915362718984000,0.32372624345386264000,\ 4020.08974915362718982600,0.38135649243870806000,-0.23751396180748979000,-0.23751396180748979000,\ 403-0.36937913195668803000,0.55728841856739109000,-0.36937913195668803000,0.65649776896096446000,\ 404-0.08622262758262820600,-0.08622262758262820600,-0.13409249329650016000,-0.48367318660804903000,\ 405-0.13409249329650016000,0.92430342836630586000,-0.79249085140642661000,-0.14644930661681688000,\ 406-0.22775636995026044000,0.34361981023636490000,0.41828517483934929000,0.40479154289778968000,\ 407-0.14644930661681685000,-0.79249085140642661000,0.41828517483934929000,0.34361981023636484000,\ 408-0.22775636995026050000,0.40479154289778979000] 409 bracc_new = [] 410 for x in range(bra.cc_all.shape[0]): 411 for y in range(bra.cc_all.shape[1]): 412 bracc_new.append(round(bra.cc_all[x,y],2)) 413 for x in range(len(bracc)): 414 bracc[x] = round(bracc[x],2) 415 #os.sys.stderr.write(str(bracc[x]) + ' --> ' + str(bracc_new[x])+'\n') 416 assert_array_almost_equal(bracc,bracc_new) 417 418 def test_cc_moiety1(self): 419 moi = pysces.model('pysces_test_moiety1.psc', self.model_dir) 420 moi.doLoad() 421 moi.hybrd_mesg = 0 422 moi.State() 423 moi.EvalEvar() 424 moi.EvalCC() 425 moicc = moicc=[0.01114694087227875000,0.07381787576415140000,\ 4260.18139104475836643000,0.04685617079784056700,0.25208717198140801000,0.04329618737852313600,\ 4270.39140460844743269000,0.01114694087227862500,0.07381787576415056700,0.18139104475836440000,\ 4280.04685617079784004000,0.25208717198140518000,0.04329618737852265100,0.39140460844742830000,\ 4290.01114694087227874000,0.07381787576415133100,0.18139104475836623000,0.04685617079784051800,\ 4300.25208717198140773000,0.04329618737852309500,0.39140460844743230000,0.01114694087227875000,\ 4310.07381787576415140000,0.18139104475836643000,0.04685617079784056700,0.25208717198140801000,\ 4320.04329618737852313600,0.39140460844743269000,0.01114694087227875000,0.07381787576415140000,\ 4330.18139104475836643000,0.04685617079784056700,0.25208717198140801000,0.04329618737852313600,\ 4340.39140460844743269000,0.01114694087227876100,0.07381787576415146900,0.18139104475836659000,\ 4350.04685617079784060900,0.25208717198140823000,0.04329618737852317800,0.39140460844743302000,\ 4360.01114694087227875200,0.07381787576415141400,0.18139104475836645000,0.04685617079784057400,\ 4370.25208717198140806000,0.04329618737852314300,0.39140460844743280000,-0.96946517151898726000,\ 4380.07237057005414701500,0.17783461222620814000,0.04593748812318156800,0.24714464011293155000,\ 4390.04244730330314599300,0.38373055769937375000,-0.00493043463469603270,-0.03265059135931840800,\ 440-0.08023158100034956400,0.14913606725290945000,\ 4410.02650472034900546900,0.11529508789995771000,-0.17312326850750914000,-0.00258941341724495350,\ 442-0.01714775382110165000,-0.04213679882643606200,0.25950964810000904000,-0.07785637143685617000,\ 443-0.02885675558547331700,-0.09092255501289707400,-0.00651179467430067560,-0.04312275948862360300,\ 444-0.10596460973080281000,-0.02003169777867189200,0.40783108111031047000,-0.00355036168876400080,\ 445-0.22864985774914801000,-0.07520202175595115700,-0.49800690277268156000,1.11329060303376930000,\ 4460.28758054022385904000,1.54718927875470390000,0.26573108181778565000,-2.64058257930148920000,\ 4470.08511194652599511600,-0.37976718223865702000,-0.93319355560031703000,-0.24105862936563618000,\ 448-1.29690043219020780000,-0.22274375836758617000,2.98855161123640300000,-0.01413200623441426100,\ 4490.06305662607990611400,0.15494766227242857000,0.04002542759392844300,0.21533763168639128000,\ 4500.03698442240380646300,-0.49621976380204558000,0.01332315621836589900,0.08822932693215754200,\ 4510.21680398717628285000,-0.25121007495788511000,0.32322678367369212000,-0.85819164176559659000,\ 4520.46781846272298522000,0.01096817827331699600,0.07263406439629339900,0.17848209108569948000,\ 4530.03374050370795465800,-0.68693259335574453000,0.00598007183654026410,0.38512768405594056000,\ 4540.00513245176028354240,0.03398840011328187900,0.08351894906745090100,-0.51437161070192172000,\ 4550.15431837495286327000,0.05719670138973560700,0.18021673341830685000] 456 moicc_new = [] 457 for x in range(moi.cc_all.shape[0]): 458 for y in range(moi.cc_all.shape[1]): 459 moicc_new.append(round(moi.cc_all[x,y],2)) 460 for x in range(len(moicc)): 461 moicc[x] = round(moicc[x],2) 462 #os.sys.stderr.write(str(moicc[x]) + ' --> ' + str(moicc_new[x])+'\n') 463 assert_array_almost_equal(moicc,moicc_new) 464 465 466 467class PyscesExternalTest(unittest.TestCase): 468 '''Extended test class, tests external/add-in numerical algorithms''' 469 __version__ = __version__ 470 model_dir = os.path.join(pysces.model_dir,'tests') 471 472 def test_PITCON1(self): 473 import scipy 474 print(''' 475 C PCPRB1.FOR The Freudenstein-Roth function. 476 C Limit points in the first variable occur at: 477 C (14.28309, -1.741377, 0.2585779) 478 C (61.66936, 1.983801, -0.6638797) 479 ''') 480 def fx(X): 481 FX[0] = X[0] - ((X[1] - 5.0)*X[1] + 2.0) * X[1] - 13.0 + 34.0 * (X[2]-1.0) 482 FX[1] = X[0] + ((X[1] + 1.0)*X[1] - 14.0) * X[1] - 29.0 + 10.0 *(X[2] - 1.0) 483 return FX 484 def fjac(s): 485 return 486 487 parm = 3; iwork = scipy.zeros((30 + parm),'i'); rwork = scipy.zeros((30 + (6*parm)*parm),'d') 488 ipar = scipy.zeros((parm),'i'); fpar = scipy.zeros((parm),'d'); xr = scipy.zeros((parm),'d') #output array 489 FX = scipy.zeros((3),'d') #function array 'v' 490 491 iwork[0] = 0;iwork[1] = 2;iwork[2] = 0;iwork[3] = 0;iwork[4] = 0;iwork[5] = 1;iwork[6] = 0 492 iwork[7] = 6;iwork[8] = 2;iwork[9] = 0;iwork[10] = 0;iwork[11] = 0;iwork[12] = 30 493 iwork[13] = len(iwork);iwork[14] = 30 + 4*parm;iwork[15] = len(rwork);iwork[16] = 20 494 495 rwork[0] = 0.00001;rwork[1] = 0.00001;rwork[2] = 0.01;rwork[3] = 20.0;rwork[4] = 0.3 496 rwork[5] = 1.0;rwork[6] = 1.0;rwork[7] = 0.0;rwork[19] = 3.0 497 498 xr[0] = 15.0;xr[1] = -2.0;xr[2] = 0.0 499 500 output = [];limits = [];iterations = 10 501 502 for run in range(iterations): 503 ierror,iwork,rwork,xrout = pysces.pitcon.pitcon1(fjac,fpar,fx,ipar,iwork,rwork,xr) 504 if iwork[0] == 4: 505 print('\nLimit point in run = ' + repr(run + 1)) 506 print(xrout) 507 limits.append(xrout) 508 output.append(xrout) 509 510 output = scipy.array(output) 511 limit0 = [14.28309125, -1.74137688, 0.25857788] 512 limit1 = [61.66936258, 1.98380112, -0.66387974] 513 out6 = [4.48781323e+001, 4.87659925e-001, 5.95322940e-002] 514 output6 = [] 515 for x in range(3): 516 output6.append(output[6,x]) 517 518 assert_array_almost_equal(limits[0],limit0,4) 519 assert_array_almost_equal(limits[1],limit1,4) 520 assert_array_almost_equal(output6,out6,4) 521 522 #def test_PITCON2(self): # unreliable 523 #mod = pysces.model('pysces_test_pitcon.psc', self.model_dir) 524 #mod.doLoad() 525 #mod.V2 = mod.V3 = 100.0 526 #mod.A05 = 10.0 527 #mod.pitcon_iter = 5 528 #mod.pitcon_par_space = logspace(0,1,10) 529 #res = mod.PITCON('P') 530 ## pysces.plt.plot2D(res, 0, style='.'); pysces.plt.logxy() 531 532 #os.sys.stderr.write('\n\"Iteration terminates after NITMAX and Newton method fails to converge\" warnings can be ignored\n') 533 #benchmark = array([ 534 #[8.39949877e-01,1.61211501e+01,4.84068424e+00,7.12011640e+01,7.12011640e+01,7.12011640e+01], 535 #[5.40984421e-01,1.57720807e+01,4.31351073e+00,7.27545617e+01,7.27545617e+01,7.27545617e+01], 536 #[-5.14362451e-02,1.45210385e+01,3.06361620e+00,7.64860963e+01,7.64860963e+01,7.64860963e+01], 537 #[-4.66310864e-01,1.31962735e+01,1.95980785e+00,8.04668684e+01,8.04668684e+01,8.04668683e+01], 538 #[1.13190191e+00,1.61011691e+01,5.26326909e+00,6.96411063e+01,6.96411063e+01,6.96411063e+01], 539 #[8.30873009e-01,1.61149804e+01,4.82588194e+00,7.12478479e+01,7.12478479e+01,7.12478479e+01], 540 #[4.39606946e-01,1.56001789e+01,4.11842961e+00,7.33077027e+01,7.33077027e+01,7.33077027e+01], 541 #[-3.07339668e-01,1.37808508e+01,2.41799659e+00,7.87206889e+01,7.87206889e+01,7.87206889e+01], 542 #[1.64608357e+00,1.37635551e+01,5.46388538e+00,6.53427535e+01,6.53427535e+01,6.53427535e+01], 543 #[1.57039736e+00,1.45248830e+01,5.53382469e+00,6.63454792e+01,6.63454792e+01,6.63454792e+01], 544 #[1.36033858e+00,1.56423070e+01,5.48938893e+00,6.81979734e+01,6.81979734e+01,6.81979734e+01], 545 #[3.90033764e-01,1.55086998e+01,4.02030887e+00,7.35869382e+01,7.35869382e+01,7.35869382e+01], 546 #[1.90300355e+00,4.61218050e-01,8.11483528e-01,1.67232557e+01,1.67232557e+01,1.67232557e+01], 547 #[1.65604156e+00,1.14527548e+00,1.32711044e+00,2.91604808e+01,2.91604810e+01,2.91604810e+01], 548 #[1.60792245e+00,3.36846985e+00,2.47903648e+00,4.55723011e+01,4.55723011e+01,4.55723011e+01], 549 #[1.76260920e+00,1.02529273e+01,4.77285351e+00,6.09994737e+01,6.09994737e+01,6.09994737e+01], 550 #[2.48410121e+00,1.76577346e-01,5.36816396e-01,7.17264824e+00,7.17264824e+00,7.17264824e+00], 551 #[2.04216203e+00,3.36094901e-01,6.94802141e-01,1.31279273e+01,1.31279273e+01,1.31279273e+01], 552 #[1.68318189e+00,9.87364532e-01,1.22120154e+00,2.69666997e+01,2.69666997e+01,2.69666997e+01], 553 #[1.60195518e+00,3.09604637e+00,2.35831199e+00,4.43144773e+01,4.43144773e+01,4.43144773e+01], 554 #[3.29451248e+00,1.10989590e-01,5.07996700e-01,3.71775348e+00,3.71775351e+00,3.71775348e+00], 555 #[2.53835638e+00,1.67389490e-01,5.28506741e-01,6.75388116e+00,6.75388117e+00,6.75388117e+00], 556 #[2.13066473e+00,2.84747678e-01,6.44229065e-01,1.14218470e+01,1.14218470e+01,1.14218470e+01], 557 #[1.85930898e+00,5.18944892e-01,8.62231652e-01,1.81726016e+01,1.81726016e+01,1.81726016e+01], 558 #[4.93995542e+00,1.06865480e-01,6.50716299e-01,2.37791780e+00,2.37791781e+00,2.37791780e+00], 559 #[5.51998987e+00,1.12930767e-01,7.15671948e-01,2.26203167e+00,2.26203167e+00,2.26203167e+00], 560 #[7.20197168e+00,1.35245848e-01,9.14818799e-01,2.13476038e+00,2.13476038e+00,2.13476038e+00], 561 #[1.21464994e+01,2.10315787e-01,1.52485159e+00,2.11433598e+00,2.11433598e+00,2.11433598e+00], 562 #[6.29277796e+00,1.22633783e-01,8.05857247e-01,2.18036060e+00,2.18036060e+00,2.18036060e+00], 563 #[6.71463987e+00,1.28373921e-01,8.56137906e-01,2.15469266e+00,2.15469266e+00,2.15469266e+00], 564 #[7.96638113e+00,1.46372121e-01,1.00775931e+00,2.11668562e+00,2.11668562e+00,2.11668562e+00], 565 #[1.16768618e+01,2.02994936e-01,1.46632082e+00,2.11150946e+00,2.11150946e+00,2.11150946e+00], 566 #[8.04040339e+00,1.47466151e-01,1.01680363e+00,2.11553054e+00,2.11553054e+00,2.11553055e+00], 567 #[8.46158253e+00,1.53734006e-01,1.06838323e+00,2.11040547e+00,2.11040547e+00,2.11040547e+00], 568 #[9.21463978e+00,1.65086727e-01,1.16102019e+00,2.10586653e+00,2.10586653e+00,2.10586653e+00], 569 #[1.07956247e+01,1.89324600e-01,1.35672050e+00,2.10728985e+00,2.10728985e+00,2.10728985e+00], 570 #[1.02976931e+01,1.81645169e-01,1.29494069e+00,2.10576586e+00,2.10576586e+00,2.10576586e+00], 571 #[1.07186771e+01,1.88135560e-01,1.34716589e+00,2.10700665e+00,2.10700665e+00,2.10700665e+00], 572 #[1.13140312e+01,1.97355266e-01,1.42115741e+00,2.10957637e+00,2.10957637e+00,2.10957637e+00], 573 #[1.21560089e+01,2.10464256e-01,1.52603757e+00,2.11439651e+00,2.11439651e+00,2.11439651e+00]]) 574 #assert_array_almost_equal(benchmark[:5,:],res[:5,:],3) 575 #assert_array_almost_equal(benchmark[-5:,:],res[-5:,:],1) 576 577 578 579 580if __name__ == '__main__': 581 unittest.main() 582 583