OpenMPCD
test_Fit.py
1 def test_fit():
2  from MPCDAnalysis.Fit import Fit
3  import pytest
4 
5 
6  def oneParam(x, p):
7  return p * x
8 
9  def twoParam(x, p, c):
10  return p * x + c
11 
12 
13  fit = Fit()
14 
15 
16  with pytest.raises(TypeError):
17  fit.fit(1, [1], [1])
18 
19  with pytest.raises(TypeError):
20  fit.fit(oneParam, 1, [1])
21 
22  with pytest.raises(TypeError):
23  fit.fit(oneParam, [1], 1)
24 
25  with pytest.raises(TypeError):
26  fit.fit(oneParam, [1], [1], 1)
27 
28 
29  with pytest.raises(ValueError):
30  fit.fit(oneParam, [1, 1], [1])
31 
32  with pytest.raises(ValueError):
33  fit.fit(oneParam, [1], [1, 1])
34 
35  with pytest.raises(ValueError):
36  fit.fit(oneParam, [1], [1], [1, 1])
37 
38 
39  fit.fit(oneParam, [1, 2], [2, 4])
40  fit.fit(oneParam, [1, 2], [2, 4], [0.1, 0.2])
41 
42  import numpy
43 
44  fit.fit(oneParam, numpy.array([1, 2]), [2, 4])
45  fit.fit(oneParam, [1, 2], numpy.array([2, 4]))
46  fit.fit(oneParam, numpy.array([1, 2]), numpy.array([2, 4]))
47  fit.fit(oneParam, numpy.array([1, 2]), [2, 4], [0.1, 0.2])
48  fit.fit(oneParam, [1, 2], numpy.array([2, 4]), [0.1, 0.2])
49  fit.fit(oneParam, [1, 2], [2, 4], numpy.array([0.1, 0.2]))
50  fit.fit(oneParam, numpy.array([1, 2]), numpy.array([2, 4]), [0.1, 0.2])
51  fit.fit(oneParam, numpy.array([1, 2]), [2, 4], numpy.array([0.1, 0.2]))
52  fit.fit(oneParam, [1, 2], numpy.array([2, 4]), numpy.array([0.1, 0.2]))
53  fit.fit(
54  oneParam,
55  numpy.array([1, 2]), numpy.array([2, 4]), numpy.array([0.1, 0.2]))
56 
57  fit.fit(twoParam, [1, 2], [2, 4])
58 
59 
60 
61 def test_fit_bounds():
62  from MPCDAnalysis.Fit import Fit
63  from MPCDAnalysis.Exceptions import OutOfRangeException
64  import pytest
65 
66 
67  def oneParam(x, p):
68  if p < 0 or p > 7:
69  raise OutOfRangeException()
70  return p * x
71 
72  def twoParam(x, p, c):
73  if p < 0 or p > 5:
74  raise OutOfRangeException()
75  if c < 3 or c > 10:
76  raise OutOfRangeException()
77  return p * x + c
78 
79 
80  fit = Fit()
81 
82 
83  with pytest.raises(TypeError):
84  fit.fit(oneParam, [1], [1], lowerBounds = 1)
85  with pytest.raises(TypeError):
86  fit.fit(oneParam, [1], [1], upperBounds = 1)
87  with pytest.raises(TypeError):
88  fit.fit(oneParam, [1], [1], lowerBounds = 1, upperBounds = 1)
89 
90  with pytest.raises(TypeError):
91  fit.fit(oneParam, [1], [1], lowerBounds = [1, 2])
92  with pytest.raises(TypeError):
93  fit.fit(oneParam, [1], [1], upperBounds = [1, 2])
94 
95  with pytest.raises(TypeError):
96  fit.fit(twoParam, [1, 2], [1, 2], lowerBounds = [1])
97  with pytest.raises(TypeError):
98  fit.fit(twoParam, [1, 2], [1, 2], upperBounds = [1])
99 
100 
101  with pytest.raises(OutOfRangeException):
102  fit.fit(oneParam, [1, 2], [-2, -4], lowerBounds = None)
103  fit.fit(oneParam, [1, 2], [-2, -4], lowerBounds = [0])
104 
105 
106  with pytest.raises(OutOfRangeException):
107  fit.fit(twoParam, [1, 2], [-3, -5], lowerBounds = None)
108  with pytest.raises(OutOfRangeException):
109  fit.fit(twoParam, [1, 2], [-3, -5], lowerBounds = [None, 3])
110  with pytest.raises(OutOfRangeException):
111  fit.fit(twoParam, [1, 2], [-3, -5], lowerBounds = [0, None])
112  fit.fit(twoParam, [1, 2], [-3, -5], lowerBounds = [0, 3])
113 
114  with pytest.raises(OutOfRangeException):
115  fit.fit(
116  twoParam, [1, 2], [30, 50],
117  lowerBounds = [0, 3], upperBounds = [None, 10])
118 
119  with pytest.raises(OutOfRangeException):
120  fit.fit(
121  twoParam, [1, 2], [30, 50],
122  lowerBounds = [0, 3], upperBounds = [5, None])
123 
124  fit.fit(
125  twoParam, [1, 2], [30, 50],
126  lowerBounds = [0, 3], upperBounds = [5, 10])
127 
128 
129 def test_getFitParameterCount():
130  from MPCDAnalysis.Fit import Fit
131  import pytest
132 
133 
134  def oneParam(x, p):
135  return p * x
136 
137  def twoParam(x, p, c):
138  return p * x + c
139 
140 
141  fit = Fit()
142 
143 
144  with pytest.raises(RuntimeError):
145  fit.getFitParameterCount()
146 
147 
148  fit.fit(oneParam, [1, 2], [2, 4])
149  assert fit.getFitParameterCount() == 1
150  assert fit.getFitParameterCount() == 1
151 
152  fit.fit(oneParam, [1, 2], [2, 4], [1.0, 0.1])
153  assert fit.getFitParameterCount() == 1
154  assert fit.getFitParameterCount() == 1
155 
156  fit.fit(twoParam, [1, 2], [2, 4])
157  assert fit.getFitParameterCount() == 2
158  assert fit.getFitParameterCount() == 2
159 
160  fit.fit(oneParam, [1, 2], [2, 4])
161  assert fit.getFitParameterCount() == 1
162  assert fit.getFitParameterCount() == 1
163 
164 
165 
166 def test_getFitParameter():
167  from MPCDAnalysis.Fit import Fit
168  import pytest
169 
170 
171  def oneParam(x, p):
172  return p * x
173 
174  def twoParam(x, p, c):
175  return p * x + c
176 
177 
178  fit = Fit()
179 
180 
181  with pytest.raises(RuntimeError):
182  fit.getFitParameter(0)
183 
184 
185  fit.fit(oneParam, [1, 2], [2, 4])
186 
187  with pytest.raises(TypeError):
188  fit.getFitParameter(0.0)
189 
190  with pytest.raises(ValueError):
191  fit.getFitParameter(-1)
192 
193  with pytest.raises(ValueError):
194  fit.getFitParameter(1)
195 
196  assert fit.getFitParameter(0) == 2
197 
198 
199  fit.fit(twoParam, [1, 2], [2, 4])
200  assert fit.getFitParameter(0) == pytest.approx(2)
201  assert fit.getFitParameter(1) == pytest.approx(0, abs = 3e-12)
202 
203  fit.fit(twoParam, [1, 2], [3, 5])
204  assert fit.getFitParameter(0) == pytest.approx(2)
205  assert fit.getFitParameter(1) == pytest.approx(1)
206 
207 
208  import scipy.optimize
209 
210  testcases = \
211  [
212  [
213  twoParam,
214  [1, 2],
215  [3.1, 5.5],
216  None
217  ],
218  [
219  twoParam,
220  [1, 2],
221  [3.1, 5.5],
222  [0.1, 0.5]
223  ],
224  ]
225 
226  for testcase in testcases:
227  f, x, y, err = testcase
228 
229  expected = \
230  scipy.optimize.curve_fit(
231  f, x, y, sigma = err, absolute_sigma = True)
232 
233  fit.fit(f, x, y, err)
234 
235  for i in range(0, fit.getFitParameterCount()):
236  assert fit.getFitParameter(i) == pytest.approx(expected[0][i])
237 
238 
239 
240 def test_getFitParameterVariance():
241  from MPCDAnalysis.Fit import Fit
243 
244  import numpy
245  import pytest
246 
247 
248  def oneParam(x, p):
249  return p * x
250 
251  def twoParam(x, p, c):
252  return p * x + c
253 
254 
255  fit = Fit()
256 
257 
258  with pytest.raises(RuntimeError):
259  fit.getFitParameterVariance(0)
260 
261 
262  fit.fit(oneParam, [1, 2], [2, 4])
263 
264  with pytest.raises(TypeError):
265  fit.getFitParameterVariance(0.0)
266 
267  with pytest.raises(ValueError):
268  fit.getFitParameterVariance(-1)
269 
270  with pytest.raises(ValueError):
271  fit.getFitParameterVariance(1)
272 
273 
274  import scipy.optimize
275 
276  testcases = \
277  [
278  [
279  twoParam,
280  [1, 2],
281  [2, 4],
282  None
283  ],
284  [
285  twoParam,
286  [1, 2],
287  [3, 5],
288  None
289  ],
290  [
291  twoParam,
292  [1, 2],
293  [3.1, 5.5],
294  None
295  ],
296  [
297  twoParam,
298  [1, 2],
299  [3.1, 5.5],
300  [0.1, 0.5]
301  ],
302 
303  [
304  twoParam,
305  [1, 2, 3],
306  [2, 4, 6],
307  None
308  ],
309  [
310  twoParam,
311  [1, 2, 3],
312  [3, 5, 7],
313  None
314  ],
315  [
316  twoParam,
317  [1, 2, 5],
318  [3.1, 5.5, 9.9],
319  None
320  ],
321  [
322  twoParam,
323  [1, 2, 3.14],
324  [3.1, 5.5, 7.77],
325  [0.1, 0.5, 0.123]
326  ],
327  ]
328 
329  for testcase in testcases:
330  f, x, y, err = testcase
331 
332  paramCount = \
334 
335  estimateAvailable = paramCount < len(x) or err is not None
336 
337  import warnings
338  with warnings.catch_warnings():
339  if not estimateAvailable:
340  warnings.simplefilter("ignore", scipy.optimize.OptimizeWarning)
341 
342  expected = \
343  scipy.optimize.curve_fit(
344  f, x, y, sigma = err,
345  absolute_sigma = err is not None)
346 
347  fit.fit(f, x, y, err)
348 
349  for i in range(0, fit.getFitParameterCount()):
350  val = expected[1][i][i]
351  assert fit.getFitParameterVariance(i) == pytest.approx(val)
352 
353  if not estimateAvailable:
354  assert fit.getFitParameterVariance(i) == numpy.inf
355 
356 
357  fit.fit(oneParam, [1, 2, 3], [3, 6, 9])
358  assert fit.getFitParameter(0) == 3
359  assert fit.getFitParameterVariance(0) == 0
360 
361  fit.fit(oneParam, [2, 3], [6, 9])
362  assert fit.getFitParameter(0) == 3
363  assert fit.getFitParameterVariance(0) == 0
364 
365  fit.fit(oneParam, [3], [9])
366  assert fit.getFitParameter(0) == 3
367  assert fit.getFitParameterVariance(0) == numpy.inf
368 
369 
370 
371 def test_getFitParameterStandardDeviation():
372  from MPCDAnalysis.Fit import Fit
373  import pytest
374 
375 
376  def oneParam(x, p):
377  return p * x
378 
379  def twoParam(x, p, c):
380  return p * x + c
381 
382 
383  fit = Fit()
384 
385 
386  with pytest.raises(RuntimeError):
387  fit.getFitParameterStandardDeviation(0)
388 
389 
390  fit.fit(oneParam, [1, 2], [2, 4])
391 
392  with pytest.raises(TypeError):
393  fit.getFitParameterStandardDeviation(0.0)
394 
395  with pytest.raises(ValueError):
396  fit.getFitParameterStandardDeviation(-1)
397 
398  with pytest.raises(ValueError):
399  fit.getFitParameterStandardDeviation(1)
400 
401 
402  import numpy
403 
404  testcases = \
405  [
406  [
407  twoParam,
408  [1, 2],
409  [2, 4],
410  None
411  ],
412  [
413  twoParam,
414  [1, 2],
415  [3, 5],
416  None
417  ],
418  [
419  twoParam,
420  [1, 2],
421  [3.1, 5.5],
422  None
423  ],
424  [
425  twoParam,
426  [1, 2],
427  [3.1, 5.5],
428  [0.1, 0.5]
429  ],
430  ]
431 
432  for testcase in testcases:
433  f, x, y, err = testcase
434 
435  fit.fit(f, x, y, err)
436 
437  for i in range(0, fit.getFitParameterCount()):
438  val = numpy.sqrt(fit.getFitParameterVariance(i))
439  assert fit.getFitParameterStandardDeviation(i) == pytest.approx(val)
440 
441 
442  fit.fit(oneParam, [1, 2, 3], [3, 6, 9])
443  assert fit.getFitParameter(0) == 3
444  assert fit.getFitParameterStandardDeviation(0) == 0
445 
446  fit.fit(oneParam, [2, 3], [6, 9])
447  assert fit.getFitParameter(0) == 3
448  assert fit.getFitParameterStandardDeviation(0) == 0
449 
450  fit.fit(oneParam, [3], [9])
451  assert fit.getFitParameter(0) == 3
452  assert fit.getFitParameterStandardDeviation(0) == numpy.inf
453 
454 
455 
456 def test_getFitParameterCovariance():
457  from MPCDAnalysis.Fit import Fit
459 
460  import numpy
461  import pytest
462 
463 
464  def oneParam(x, p):
465  return p * x
466 
467  def twoParam(x, p, c):
468  return p * x + c
469 
470 
471  fit = Fit()
472 
473 
474  with pytest.raises(RuntimeError):
475  fit.getFitParameterCovariance(0, 0)
476 
477 
478  fit.fit(oneParam, [1, 2], [2, 4])
479 
480  with pytest.raises(TypeError):
481  fit.getFitParameterCovariance(0.0, 0)
482 
483  with pytest.raises(TypeError):
484  fit.getFitParameterCovariance(0, 0.0)
485 
486  with pytest.raises(TypeError):
487  fit.getFitParameterCovariance(0.0, 0.0)
488 
489  with pytest.raises(ValueError):
490  fit.getFitParameterCovariance(-1, 0)
491 
492  with pytest.raises(ValueError):
493  fit.getFitParameterCovariance(0, -1)
494 
495  with pytest.raises(ValueError):
496  fit.getFitParameterCovariance(-1, -1)
497 
498  with pytest.raises(ValueError):
499  fit.getFitParameterCovariance(0, 1)
500 
501  with pytest.raises(ValueError):
502  fit.getFitParameterCovariance(1, 0)
503 
504  with pytest.raises(ValueError):
505  fit.getFitParameterCovariance(1, 1)
506 
507 
508  import scipy.optimize
509 
510  testcases = \
511  [
512  [
513  twoParam,
514  [1, 2],
515  [2, 4],
516  None
517  ],
518  [
519  twoParam,
520  [1, 2],
521  [3, 5],
522  None
523  ],
524  [
525  twoParam,
526  [1, 2],
527  [3.1, 5.5],
528  None
529  ],
530  [
531  twoParam,
532  [1, 2],
533  [3.1, 5.5],
534  [0.1, 0.5]
535  ],
536 
537  [
538  twoParam,
539  [1, 2, 3],
540  [2, 4, 6],
541  None
542  ],
543  [
544  twoParam,
545  [1, 2, 3],
546  [3, 5, 7],
547  None
548  ],
549  [
550  twoParam,
551  [1, 2, 5],
552  [3.1, 5.5, 9.9],
553  None
554  ],
555  [
556  twoParam,
557  [1, 2, 3.14],
558  [3.1, 5.5, 7.77],
559  [0.1, 0.5, 0.123]
560  ],
561  ]
562 
563  for testcase in testcases:
564  f, x, y, err = testcase
565 
566  paramCount = \
568 
569  estimateAvailable = paramCount < len(x) or err is not None
570 
571 
572  fit.fit(f, x, y, err)
573 
574  import warnings
575  with warnings.catch_warnings():
576  if not estimateAvailable:
577  warnings.simplefilter("ignore", scipy.optimize.OptimizeWarning)
578 
579  expected = \
580  scipy.optimize.curve_fit(
581  f, x, y, sigma = err,
582  absolute_sigma = err is not None)
583 
584  for i in range(0, fit.getFitParameterCount()):
585  for j in range(0, fit.getFitParameterCount()):
586  cov_ij = fit.getFitParameterCovariance(i, j)
587  cov_ji = fit.getFitParameterCovariance(j, i)
588 
589  assert cov_ij == pytest.approx(cov_ji)
590  assert cov_ij == pytest.approx(expected[1][i][j])
591 
592  if not estimateAvailable:
593  assert cov_ij == numpy.inf
594  assert cov_ji == numpy.inf
595 
596 
597 
598 def test_multipleIndependentDataFits():
599  from MPCDAnalysis.Fit import Fit
600  import pytest
601 
602 
603  def twoParam(x, p, c):
604  return p * x + c
605 
606 
607  testcases = \
608  [
609  [
610  [1, 2],
611  [2, 4],
612  None
613  ],
614  [
615  [1, 2],
616  [3, 5],
617  None
618  ],
619  [
620  [1, 2],
621  [3.1, 5.5],
622  None
623  ],
624  [
625  [1, 2],
626  [3.1, 5.5],
627  [0.1, 0.5]
628  ],
629  ]
630 
631  results = Fit.multipleIndependentDataFits(twoParam, testcases)
632 
633  assert isinstance(results, list)
634 
635  for index, result in enumerate(results):
636  assert isinstance(result, dict)
637  assert result['data'] == testcases[index]
638 
639  fit = Fit()
640  fit.fit(
641  twoParam,
642  testcases[index][0], testcases[index][1], testcases[index][2])
643  assert result['fit'] == fit
644 
645 
646 
647 def test_getBestOutOfMultipleIndependentDataFits():
648  from MPCDAnalysis.Fit import Fit
649  import pytest
650 
651 
652  def twoParam(x, p, c):
653  return p * x + c
654 
655 
656  testcases = \
657  [
658  [
659  [1, 2],
660  [2, 4],
661  None
662  ],
663  [
664  [1, 2],
665  [3, 5],
666  None
667  ],
668  [
669  [1, 2],
670  [3.1, 5.5],
671  None
672  ],
673  [
674  [1, 2],
675  [3.1, 5.5],
676  [0.1, 0.5]
677  ],
678  ]
679 
680 
681  def comp(old, new):
682  oldVariance = old["fit"].getFitParameterVariance(0)
683  newVariance = new["fit"].getFitParameterVariance(0)
684  return newVariance < oldVariance
685 
686  assert \
687  Fit.getBestOutOfMultipleIndependentDataFits(twoParam, [], comp) is None
688 
689  fits = Fit.multipleIndependentDataFits(twoParam, testcases)
690 
691  best = fits[0]
692  for i in range(1, len(fits)):
693  if comp(best, fits[i]):
694  best = fits[i]
695 
696  result = \
697  Fit.getBestOutOfMultipleIndependentDataFits(twoParam, testcases, comp)
698  assert best == result
699 
700  minVariance = 1e300
701  for fit in fits:
702  minVariance = min(minVariance, fit["fit"].getFitParameterVariance(0))
703  assert minVariance == result["fit"].getFitParameterVariance(0)
704 
705 
706 
707  xData = []
708  yData = []
709  for x in range(-3, 3 + 1):
710  xData.append(x)
711  yData.append(x * x - 1)
712 
713 
714  def oneParam(x, k):
715  return x * k
716 
717  def discardFrontGenerator(xData, yData, minNumElements):
718  for i in range(0, len(xData) - minNumElements + 1):
719  yield [xData[i:], yData[i:], None]
720 
721 
722  def comp(old, new):
723  oldVariance = old["fit"].getFitParameterVariance(0)
724  newVariance = new["fit"].getFitParameterVariance(0)
725  return newVariance < oldVariance
726 
727 
728  result = \
729  Fit.getBestOutOfMultipleIndependentDataFits(
730  oneParam,
731  discardFrontGenerator(xData, yData, 1),
732  comp)
733 
734  #expected results are taken reference data:
735  assert \
736  result["data"] == [xData[len(xData) - 4:], yData[len(yData) - 4:], None]
737  assert result["fit"].getFitParameter(0) == 2.1428571428596355
738  assert result["fit"].getFitParameterVariance(0) == 0.2312925188464976
739 
740 
741  result = \
742  Fit.getBestOutOfMultipleIndependentDataFits(
743  twoParam,
744  discardFrontGenerator(xData, yData, 2),
745  comp)
746 
747  #expected results are taken reference data:
748  assert \
749  result["data"] == [xData[len(xData) - 3:], yData[len(yData) - 3:], None]
750  assert result["fit"].getFitParameter(0) == 4.0000000000000009
751  assert result["fit"].getFitParameterVariance(0) == 0.33333333333587728
752  assert result["fit"].getFitParameter(1) == -4.3333333333333339
753  assert result["fit"].getFitParameterVariance(1) == 1.5555555698350492
754 
755 
756 
757 def test___eq__():
758  from MPCDAnalysis.Fit import Fit
759  import pytest
760 
761  empty1 = Fit()
762  empty2 = Fit()
763 
764  with pytest.raises(TypeError):
765  empty1 == 1
766 
767  assert empty1 == empty1
768  assert empty2 == empty2
769 
770  assert empty1 == empty2
771  assert empty2 == empty1
772 
773  f1 = Fit()
774  f1.fit(lambda x, p: p * x, [1], [2])
775 
776  f2 = Fit()
777  f2.fit(lambda x, p: p * x, [1], [2])
778 
779  f3 = Fit()
780  f3.fit(lambda x, p: p * x, [1], [2], [3])
781 
782  f4 = Fit()
783  f4.fit(lambda x, p: p * x, [1], [3])
784 
785  f5 = Fit()
786  f5.fit(lambda x, p: p * x, [2], [2])
787 
788 
789  assert f1 == f1
790  assert f1 == f2
791 
792  assert not f1 == f3
793  assert not f1 == f4
794  assert not f1 == f5
795 
796  assert not f3 == f4
797  assert not f3 == f5
798 
799  assert not f4 == f5
800 
801 
802 
803 def test___neq__():
804  from MPCDAnalysis.Fit import Fit
805  import pytest
806 
807  empty1 = Fit()
808  empty2 = Fit()
809 
810  with pytest.raises(TypeError):
811  empty1 != 1
812 
813  assert not empty1 != empty1
814  assert not empty2 != empty2
815 
816  assert not empty1 != empty2
817  assert not empty2 != empty1
818 
819  f1 = Fit()
820  f1.fit(lambda x, p: p * x, [1], [2])
821 
822  f2 = Fit()
823  f2.fit(lambda x, p: p * x, [1], [2])
824 
825  f3 = Fit()
826  f3.fit(lambda x, p: p * x, [1], [2], [3])
827 
828  f4 = Fit()
829  f4.fit(lambda x, p: p * x, [1], [3])
830 
831  f5 = Fit()
832  f5.fit(lambda x, p: p * x, [2], [2])
833 
834 
835  assert not f1 != f1
836  assert not f1 != f2
837 
838  assert f1 != f3
839  assert f1 != f4
840  assert f1 != f5
841 
842  assert f3 != f4
843  assert f3 != f5
844 
845  assert f4 != f5
MPCDAnalysis.Exceptions
Definition: Exceptions.py:1
MPCDAnalysis.Utilities
Definition: Utilities.py:1
MPCDAnalysis.Fit
Definition: Fit.py:1
MPCDAnalysis.Utilities.getNumberOfArgumentsFromCallable
def getNumberOfArgumentsFromCallable(f)
Definition: Utilities.py:214