1import tempfile 2import shutil 3import os 4 5import numpy as np 6from numpy import pi 7from numpy.testing import (assert_array_almost_equal, 8 assert_equal, assert_warns) 9import pytest 10from pytest import raises as assert_raises 11 12from scipy.odr import (Data, Model, ODR, RealData, OdrStop, OdrWarning, 13 multilinear, exponential, unilinear, quadratic, 14 polynomial) 15 16 17class TestODR: 18 19 # Bad Data for 'x' 20 21 def test_bad_data(self): 22 assert_raises(ValueError, Data, 2, 1) 23 assert_raises(ValueError, RealData, 2, 1) 24 25 # Empty Data for 'x' 26 def empty_data_func(self, B, x): 27 return B[0]*x + B[1] 28 29 def test_empty_data(self): 30 beta0 = [0.02, 0.0] 31 linear = Model(self.empty_data_func) 32 33 empty_dat = Data([], []) 34 assert_warns(OdrWarning, ODR, 35 empty_dat, linear, beta0=beta0) 36 37 empty_dat = RealData([], []) 38 assert_warns(OdrWarning, ODR, 39 empty_dat, linear, beta0=beta0) 40 41 # Explicit Example 42 43 def explicit_fcn(self, B, x): 44 ret = B[0] + B[1] * np.power(np.exp(B[2]*x) - 1.0, 2) 45 return ret 46 47 def explicit_fjd(self, B, x): 48 eBx = np.exp(B[2]*x) 49 ret = B[1] * 2.0 * (eBx-1.0) * B[2] * eBx 50 return ret 51 52 def explicit_fjb(self, B, x): 53 eBx = np.exp(B[2]*x) 54 res = np.vstack([np.ones(x.shape[-1]), 55 np.power(eBx-1.0, 2), 56 B[1]*2.0*(eBx-1.0)*eBx*x]) 57 return res 58 59 def test_explicit(self): 60 explicit_mod = Model( 61 self.explicit_fcn, 62 fjacb=self.explicit_fjb, 63 fjacd=self.explicit_fjd, 64 meta=dict(name='Sample Explicit Model', 65 ref='ODRPACK UG, pg. 39'), 66 ) 67 explicit_dat = Data([0.,0.,5.,7.,7.5,10.,16.,26.,30.,34.,34.5,100.], 68 [1265.,1263.6,1258.,1254.,1253.,1249.8,1237.,1218.,1220.6, 69 1213.8,1215.5,1212.]) 70 explicit_odr = ODR(explicit_dat, explicit_mod, beta0=[1500.0, -50.0, -0.1], 71 ifixx=[0,0,1,1,1,1,1,1,1,1,1,0]) 72 explicit_odr.set_job(deriv=2) 73 explicit_odr.set_iprint(init=0, iter=0, final=0) 74 75 out = explicit_odr.run() 76 assert_array_almost_equal( 77 out.beta, 78 np.array([1.2646548050648876e+03, -5.4018409956678255e+01, 79 -8.7849712165253724e-02]), 80 ) 81 assert_array_almost_equal( 82 out.sd_beta, 83 np.array([1.0349270280543437, 1.583997785262061, 0.0063321988657267]), 84 ) 85 assert_array_almost_equal( 86 out.cov_beta, 87 np.array([[4.4949592379003039e-01, -3.7421976890364739e-01, 88 -8.0978217468468912e-04], 89 [-3.7421976890364739e-01, 1.0529686462751804e+00, 90 -1.9453521827942002e-03], 91 [-8.0978217468468912e-04, -1.9453521827942002e-03, 92 1.6827336938454476e-05]]), 93 ) 94 95 # Implicit Example 96 97 def implicit_fcn(self, B, x): 98 return (B[2]*np.power(x[0]-B[0], 2) + 99 2.0*B[3]*(x[0]-B[0])*(x[1]-B[1]) + 100 B[4]*np.power(x[1]-B[1], 2) - 1.0) 101 102 def test_implicit(self): 103 implicit_mod = Model( 104 self.implicit_fcn, 105 implicit=1, 106 meta=dict(name='Sample Implicit Model', 107 ref='ODRPACK UG, pg. 49'), 108 ) 109 implicit_dat = Data([ 110 [0.5,1.2,1.6,1.86,2.12,2.36,2.44,2.36,2.06,1.74,1.34,0.9,-0.28, 111 -0.78,-1.36,-1.9,-2.5,-2.88,-3.18,-3.44], 112 [-0.12,-0.6,-1.,-1.4,-2.54,-3.36,-4.,-4.75,-5.25,-5.64,-5.97,-6.32, 113 -6.44,-6.44,-6.41,-6.25,-5.88,-5.5,-5.24,-4.86]], 114 1, 115 ) 116 implicit_odr = ODR(implicit_dat, implicit_mod, 117 beta0=[-1.0, -3.0, 0.09, 0.02, 0.08]) 118 119 out = implicit_odr.run() 120 assert_array_almost_equal( 121 out.beta, 122 np.array([-0.9993809167281279, -2.9310484652026476, 0.0875730502693354, 123 0.0162299708984738, 0.0797537982976416]), 124 ) 125 assert_array_almost_equal( 126 out.sd_beta, 127 np.array([0.1113840353364371, 0.1097673310686467, 0.0041060738314314, 128 0.0027500347539902, 0.0034962501532468]), 129 ) 130 assert_array_almost_equal( 131 out.cov_beta, 132 np.array([[2.1089274602333052e+00, -1.9437686411979040e+00, 133 7.0263550868344446e-02, -4.7175267373474862e-02, 134 5.2515575927380355e-02], 135 [-1.9437686411979040e+00, 2.0481509222414456e+00, 136 -6.1600515853057307e-02, 4.6268827806232933e-02, 137 -5.8822307501391467e-02], 138 [7.0263550868344446e-02, -6.1600515853057307e-02, 139 2.8659542561579308e-03, -1.4628662260014491e-03, 140 1.4528860663055824e-03], 141 [-4.7175267373474862e-02, 4.6268827806232933e-02, 142 -1.4628662260014491e-03, 1.2855592885514335e-03, 143 -1.2692942951415293e-03], 144 [5.2515575927380355e-02, -5.8822307501391467e-02, 145 1.4528860663055824e-03, -1.2692942951415293e-03, 146 2.0778813389755596e-03]]), 147 ) 148 149 # Multi-variable Example 150 151 def multi_fcn(self, B, x): 152 if (x < 0.0).any(): 153 raise OdrStop 154 theta = pi*B[3]/2. 155 ctheta = np.cos(theta) 156 stheta = np.sin(theta) 157 omega = np.power(2.*pi*x*np.exp(-B[2]), B[3]) 158 phi = np.arctan2((omega*stheta), (1.0 + omega*ctheta)) 159 r = (B[0] - B[1]) * np.power(np.sqrt(np.power(1.0 + omega*ctheta, 2) + 160 np.power(omega*stheta, 2)), -B[4]) 161 ret = np.vstack([B[1] + r*np.cos(B[4]*phi), 162 r*np.sin(B[4]*phi)]) 163 return ret 164 165 def test_multi(self): 166 multi_mod = Model( 167 self.multi_fcn, 168 meta=dict(name='Sample Multi-Response Model', 169 ref='ODRPACK UG, pg. 56'), 170 ) 171 172 multi_x = np.array([30.0, 50.0, 70.0, 100.0, 150.0, 200.0, 300.0, 500.0, 173 700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, 10000.0, 174 15000.0, 20000.0, 30000.0, 50000.0, 70000.0, 100000.0, 150000.0]) 175 multi_y = np.array([ 176 [4.22, 4.167, 4.132, 4.038, 4.019, 3.956, 3.884, 3.784, 3.713, 177 3.633, 3.54, 3.433, 3.358, 3.258, 3.193, 3.128, 3.059, 2.984, 178 2.934, 2.876, 2.838, 2.798, 2.759], 179 [0.136, 0.167, 0.188, 0.212, 0.236, 0.257, 0.276, 0.297, 0.309, 180 0.311, 0.314, 0.311, 0.305, 0.289, 0.277, 0.255, 0.24, 0.218, 181 0.202, 0.182, 0.168, 0.153, 0.139], 182 ]) 183 n = len(multi_x) 184 multi_we = np.zeros((2, 2, n), dtype=float) 185 multi_ifixx = np.ones(n, dtype=int) 186 multi_delta = np.zeros(n, dtype=float) 187 188 multi_we[0,0,:] = 559.6 189 multi_we[1,0,:] = multi_we[0,1,:] = -1634.0 190 multi_we[1,1,:] = 8397.0 191 192 for i in range(n): 193 if multi_x[i] < 100.0: 194 multi_ifixx[i] = 0 195 elif multi_x[i] <= 150.0: 196 pass # defaults are fine 197 elif multi_x[i] <= 1000.0: 198 multi_delta[i] = 25.0 199 elif multi_x[i] <= 10000.0: 200 multi_delta[i] = 560.0 201 elif multi_x[i] <= 100000.0: 202 multi_delta[i] = 9500.0 203 else: 204 multi_delta[i] = 144000.0 205 if multi_x[i] == 100.0 or multi_x[i] == 150.0: 206 multi_we[:,:,i] = 0.0 207 208 multi_dat = Data(multi_x, multi_y, wd=1e-4/np.power(multi_x, 2), 209 we=multi_we) 210 multi_odr = ODR(multi_dat, multi_mod, beta0=[4.,2.,7.,.4,.5], 211 delta0=multi_delta, ifixx=multi_ifixx) 212 multi_odr.set_job(deriv=1, del_init=1) 213 214 out = multi_odr.run() 215 assert_array_almost_equal( 216 out.beta, 217 np.array([4.3799880305938963, 2.4333057577497703, 8.0028845899503978, 218 0.5101147161764654, 0.5173902330489161]), 219 ) 220 assert_array_almost_equal( 221 out.sd_beta, 222 np.array([0.0130625231081944, 0.0130499785273277, 0.1167085962217757, 223 0.0132642749596149, 0.0288529201353984]), 224 ) 225 assert_array_almost_equal( 226 out.cov_beta, 227 np.array([[0.0064918418231375, 0.0036159705923791, 0.0438637051470406, 228 -0.0058700836512467, 0.011281212888768], 229 [0.0036159705923791, 0.0064793789429006, 0.0517610978353126, 230 -0.0051181304940204, 0.0130726943624117], 231 [0.0438637051470406, 0.0517610978353126, 0.5182263323095322, 232 -0.0563083340093696, 0.1269490939468611], 233 [-0.0058700836512467, -0.0051181304940204, -0.0563083340093696, 234 0.0066939246261263, -0.0140184391377962], 235 [0.011281212888768, 0.0130726943624117, 0.1269490939468611, 236 -0.0140184391377962, 0.0316733013820852]]), 237 ) 238 239 # Pearson's Data 240 # K. Pearson, Philosophical Magazine, 2, 559 (1901) 241 242 def pearson_fcn(self, B, x): 243 return B[0] + B[1]*x 244 245 def test_pearson(self): 246 p_x = np.array([0.,.9,1.8,2.6,3.3,4.4,5.2,6.1,6.5,7.4]) 247 p_y = np.array([5.9,5.4,4.4,4.6,3.5,3.7,2.8,2.8,2.4,1.5]) 248 p_sx = np.array([.03,.03,.04,.035,.07,.11,.13,.22,.74,1.]) 249 p_sy = np.array([1.,.74,.5,.35,.22,.22,.12,.12,.1,.04]) 250 251 p_dat = RealData(p_x, p_y, sx=p_sx, sy=p_sy) 252 253 # Reverse the data to test invariance of results 254 pr_dat = RealData(p_y, p_x, sx=p_sy, sy=p_sx) 255 256 p_mod = Model(self.pearson_fcn, meta=dict(name='Uni-linear Fit')) 257 258 p_odr = ODR(p_dat, p_mod, beta0=[1.,1.]) 259 pr_odr = ODR(pr_dat, p_mod, beta0=[1.,1.]) 260 261 out = p_odr.run() 262 assert_array_almost_equal( 263 out.beta, 264 np.array([5.4767400299231674, -0.4796082367610305]), 265 ) 266 assert_array_almost_equal( 267 out.sd_beta, 268 np.array([0.3590121690702467, 0.0706291186037444]), 269 ) 270 assert_array_almost_equal( 271 out.cov_beta, 272 np.array([[0.0854275622946333, -0.0161807025443155], 273 [-0.0161807025443155, 0.003306337993922]]), 274 ) 275 276 rout = pr_odr.run() 277 assert_array_almost_equal( 278 rout.beta, 279 np.array([11.4192022410781231, -2.0850374506165474]), 280 ) 281 assert_array_almost_equal( 282 rout.sd_beta, 283 np.array([0.9820231665657161, 0.3070515616198911]), 284 ) 285 assert_array_almost_equal( 286 rout.cov_beta, 287 np.array([[0.6391799462548782, -0.1955657291119177], 288 [-0.1955657291119177, 0.0624888159223392]]), 289 ) 290 291 # Lorentz Peak 292 # The data is taken from one of the undergraduate physics labs I performed. 293 294 def lorentz(self, beta, x): 295 return (beta[0]*beta[1]*beta[2] / np.sqrt(np.power(x*x - 296 beta[2]*beta[2], 2.0) + np.power(beta[1]*x, 2.0))) 297 298 def test_lorentz(self): 299 l_sy = np.array([.29]*18) 300 l_sx = np.array([.000972971,.000948268,.000707632,.000706679, 301 .000706074, .000703918,.000698955,.000456856, 302 .000455207,.000662717,.000654619,.000652694, 303 .000000859202,.00106589,.00106378,.00125483, .00140818,.00241839]) 304 305 l_dat = RealData( 306 [3.9094, 3.85945, 3.84976, 3.84716, 3.84551, 3.83964, 3.82608, 307 3.78847, 3.78163, 3.72558, 3.70274, 3.6973, 3.67373, 3.65982, 308 3.6562, 3.62498, 3.55525, 3.41886], 309 [652, 910.5, 984, 1000, 1007.5, 1053, 1160.5, 1409.5, 1430, 1122, 310 957.5, 920, 777.5, 709.5, 698, 578.5, 418.5, 275.5], 311 sx=l_sx, 312 sy=l_sy, 313 ) 314 l_mod = Model(self.lorentz, meta=dict(name='Lorentz Peak')) 315 l_odr = ODR(l_dat, l_mod, beta0=(1000., .1, 3.8)) 316 317 out = l_odr.run() 318 assert_array_almost_equal( 319 out.beta, 320 np.array([1.4306780846149925e+03, 1.3390509034538309e-01, 321 3.7798193600109009e+00]), 322 ) 323 assert_array_almost_equal( 324 out.sd_beta, 325 np.array([7.3621186811330963e-01, 3.5068899941471650e-04, 326 2.4451209281408992e-04]), 327 ) 328 assert_array_almost_equal( 329 out.cov_beta, 330 np.array([[2.4714409064597873e-01, -6.9067261911110836e-05, 331 -3.1236953270424990e-05], 332 [-6.9067261911110836e-05, 5.6077531517333009e-08, 333 3.6133261832722601e-08], 334 [-3.1236953270424990e-05, 3.6133261832722601e-08, 335 2.7261220025171730e-08]]), 336 ) 337 338 def test_ticket_1253(self): 339 def linear(c, x): 340 return c[0]*x+c[1] 341 342 c = [2.0, 3.0] 343 x = np.linspace(0, 10) 344 y = linear(c, x) 345 346 model = Model(linear) 347 data = Data(x, y, wd=1.0, we=1.0) 348 job = ODR(data, model, beta0=[1.0, 1.0]) 349 result = job.run() 350 assert_equal(result.info, 2) 351 352 # Verify fix for gh-9140 353 354 def test_ifixx(self): 355 x1 = [-2.01, -0.99, -0.001, 1.02, 1.98] 356 x2 = [3.98, 1.01, 0.001, 0.998, 4.01] 357 fix = np.vstack((np.zeros_like(x1, dtype=int), np.ones_like(x2, dtype=int))) 358 data = Data(np.vstack((x1, x2)), y=1, fix=fix) 359 model = Model(lambda beta, x: x[1, :] - beta[0] * x[0, :]**2., implicit=True) 360 361 odr1 = ODR(data, model, beta0=np.array([1.])) 362 sol1 = odr1.run() 363 odr2 = ODR(data, model, beta0=np.array([1.]), ifixx=fix) 364 sol2 = odr2.run() 365 assert_equal(sol1.beta, sol2.beta) 366 367 # verify bugfix for #11800 in #11802 368 def test_ticket_11800(self): 369 # parameters 370 beta_true = np.array([1.0, 2.3, 1.1, -1.0, 1.3, 0.5]) 371 nr_measurements = 10 372 373 std_dev_x = 0.01 374 x_error = np.array([[0.00063445, 0.00515731, 0.00162719, 0.01022866, 375 -0.01624845, 0.00482652, 0.00275988, -0.00714734, -0.00929201, -0.00687301], 376 [-0.00831623, -0.00821211, -0.00203459, 0.00938266, -0.00701829, 377 0.0032169, 0.00259194, -0.00581017, -0.0030283, 0.01014164]]) 378 379 std_dev_y = 0.05 380 y_error = np.array([[0.05275304, 0.04519563, -0.07524086, 0.03575642, 381 0.04745194, 0.03806645, 0.07061601, -0.00753604, -0.02592543, -0.02394929], 382 [0.03632366, 0.06642266, 0.08373122, 0.03988822, -0.0092536, 383 -0.03750469, -0.03198903, 0.01642066, 0.01293648, -0.05627085]]) 384 385 beta_solution = np.array([ 386 2.62920235756665876536e+00, -1.26608484996299608838e+02, 1.29703572775403074502e+02, 387 -1.88560985401185465804e+00, 7.83834160771274923718e+01, -7.64124076838087091801e+01]) 388 389 # model's function and Jacobians 390 def func(beta, x): 391 y0 = beta[0] + beta[1] * x[0, :] + beta[2] * x[1, :] 392 y1 = beta[3] + beta[4] * x[0, :] + beta[5] * x[1, :] 393 394 return np.vstack((y0, y1)) 395 396 def df_dbeta_odr(beta, x): 397 nr_meas = np.shape(x)[1] 398 zeros = np.zeros(nr_meas) 399 ones = np.ones(nr_meas) 400 401 dy0 = np.array([ones, x[0, :], x[1, :], zeros, zeros, zeros]) 402 dy1 = np.array([zeros, zeros, zeros, ones, x[0, :], x[1, :]]) 403 404 return np.stack((dy0, dy1)) 405 406 def df_dx_odr(beta, x): 407 nr_meas = np.shape(x)[1] 408 ones = np.ones(nr_meas) 409 410 dy0 = np.array([beta[1] * ones, beta[2] * ones]) 411 dy1 = np.array([beta[4] * ones, beta[5] * ones]) 412 return np.stack((dy0, dy1)) 413 414 # do measurements with errors in independent and dependent variables 415 x0_true = np.linspace(1, 10, nr_measurements) 416 x1_true = np.linspace(1, 10, nr_measurements) 417 x_true = np.array([x0_true, x1_true]) 418 419 y_true = func(beta_true, x_true) 420 421 x_meas = x_true + x_error 422 y_meas = y_true + y_error 423 424 # estimate model's parameters 425 model_f = Model(func, fjacb=df_dbeta_odr, fjacd=df_dx_odr) 426 427 data = RealData(x_meas, y_meas, sx=std_dev_x, sy=std_dev_y) 428 429 odr_obj = ODR(data, model_f, beta0=0.9 * beta_true, maxit=100) 430 #odr_obj.set_iprint(init=2, iter=0, iter_step=1, final=1) 431 odr_obj.set_job(deriv=3) 432 433 odr_out = odr_obj.run() 434 435 # check results 436 assert_equal(odr_out.info, 1) 437 assert_array_almost_equal(odr_out.beta, beta_solution) 438 439 def test_multilinear_model(self): 440 x = np.linspace(0.0, 5.0) 441 y = 10.0 + 5.0 * x 442 data = Data(x, y) 443 odr_obj = ODR(data, multilinear) 444 output = odr_obj.run() 445 assert_array_almost_equal(output.beta, [10.0, 5.0]) 446 447 def test_exponential_model(self): 448 x = np.linspace(0.0, 5.0) 449 y = -10.0 + np.exp(0.5*x) 450 data = Data(x, y) 451 odr_obj = ODR(data, exponential) 452 output = odr_obj.run() 453 assert_array_almost_equal(output.beta, [-10.0, 0.5]) 454 455 def test_polynomial_model(self): 456 x = np.linspace(0.0, 5.0) 457 y = 1.0 + 2.0 * x + 3.0 * x ** 2 + 4.0 * x ** 3 458 poly_model = polynomial(3) 459 data = Data(x, y) 460 odr_obj = ODR(data, poly_model) 461 output = odr_obj.run() 462 assert_array_almost_equal(output.beta, [1.0, 2.0, 3.0, 4.0]) 463 464 def test_unilinear_model(self): 465 x = np.linspace(0.0, 5.0) 466 y = 1.0 * x + 2.0 467 data = Data(x, y) 468 odr_obj = ODR(data, unilinear) 469 output = odr_obj.run() 470 assert_array_almost_equal(output.beta, [1.0, 2.0]) 471 472 def test_quadratic_model(self): 473 x = np.linspace(0.0, 5.0) 474 y = 1.0 * x ** 2 + 2.0 * x + 3.0 475 data = Data(x, y) 476 odr_obj = ODR(data, quadratic) 477 output = odr_obj.run() 478 assert_array_almost_equal(output.beta, [1.0, 2.0, 3.0]) 479 480 def test_work_ind(self): 481 482 def func(par, x): 483 b0, b1 = par 484 return b0 + b1 * x 485 486 # generate some data 487 n_data = 4 488 x = np.arange(n_data) 489 y = np.where(x % 2, x + 0.1, x - 0.1) 490 x_err = np.full(n_data, 0.1) 491 y_err = np.full(n_data, 0.1) 492 493 # do the fitting 494 linear_model = Model(func) 495 real_data = RealData(x, y, sx=x_err, sy=y_err) 496 odr_obj = ODR(real_data, linear_model, beta0=[0.4, 0.4]) 497 odr_obj.set_job(fit_type=0) 498 out = odr_obj.run() 499 500 sd_ind = out.work_ind['sd'] 501 assert_array_almost_equal(out.sd_beta, 502 out.work[sd_ind:sd_ind + len(out.sd_beta)]) 503 504 @pytest.mark.skipif(True, reason="Fortran I/O prone to crashing so better " 505 "not to run this test, see gh-13127") 506 def test_output_file_overwrite(self): 507 """ 508 Verify fix for gh-1892 509 """ 510 def func(b, x): 511 return b[0] + b[1] * x 512 513 p = Model(func) 514 data = Data(np.arange(10), 12 * np.arange(10)) 515 tmp_dir = tempfile.mkdtemp() 516 error_file_path = os.path.join(tmp_dir, "error.dat") 517 report_file_path = os.path.join(tmp_dir, "report.dat") 518 try: 519 ODR(data, p, beta0=[0.1, 13], errfile=error_file_path, 520 rptfile=report_file_path).run() 521 ODR(data, p, beta0=[0.1, 13], errfile=error_file_path, 522 rptfile=report_file_path, overwrite=True).run() 523 finally: 524 # remove output files for clean up 525 shutil.rmtree(tmp_dir) 526 527