test_odr.py 22 KB

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