99from scipy .optimize import differential_evolution
1010from scipy .special import gamma
1111
12- def stretch_exp_fit (TRPL , t , Tc = (0 ,1e5 ), Beta = (0 ,1 ), A = (0 ,1e6 )):
12+ def stretch_exp_fit (TRPL , t , Tc = (0 ,1e5 ), Beta = (0 ,1 ), A = (0 ,1e6 ), noise = ( 0 , 1e6 ) ):
1313
14- def exp_stretch (t , tc , beta , a ):
15- return (a * np .exp (- ((1.0 / tc ) * t ) ** beta ))
14+ def exp_stretch (t , tc , beta , a , noise ):
15+ return (( a * np .exp (- ((1.0 / tc ) * t ) ** beta )) + noise )
1616
1717 def avg_tau_from_exp_stretch (tc , beta ):
1818 return (tc / beta ) * gamma (1.0 / beta )
@@ -25,13 +25,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
2525 tc = params [0 ]
2626 beta = params [1 ]
2727 a = params [2 ]
28+ noise = params [3 ]
2829
29- PL_sim = exp_stretch (t ,tc ,beta ,a )
30+ PL_sim = exp_stretch (t ,tc ,beta ,a , noise )
3031
31- Resid = np .sqrt ( np . sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (PL_sim )** 2 ) ))
32+ Resid = np .sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (TRPL )** 2 ))
3233 return Resid #returns the difference between the PL data and simulated data
3334
34- bounds = [Tc , Beta , A ]
35+ bounds = [Tc , Beta , A , noise ]
3536
3637 result = differential_evolution (residuals , bounds )
3738 return result .x
@@ -41,20 +42,21 @@ def residuals(params):#params are the parameters to be adjusted by differential
4142 tc = p [0 ]
4243 beta = p [1 ]
4344 a = p [2 ]
45+ noise = p [3 ]
4446
45- PL_fit = exp_stretch (t ,tc ,beta ,a )
47+ PL_fit = exp_stretch (t ,tc ,beta ,a , noise )
4648
4749 avg_tau = avg_tau_from_exp_stretch (tc ,beta )
4850
49- return tc , beta , a , avg_tau , PL_fit
51+ return tc , beta , a , avg_tau , PL_fit , noise
5052
51- def double_exp_fit (TRPL , t , tau1_bounds = (0 ,100 ), a1_bounds = (0 ,1e5 ), tau2_bounds = (0 ,100 ), a2_bounds = (0 ,1e5 )):
53+ def double_exp_fit (TRPL , t , tau1_bounds = (0 ,1000 ), a1_bounds = (0 ,1e6 ), tau2_bounds = (0 ,10000 ), a2_bounds = (0 ,1e5 ), noise = ( 0 , 1e6 )):
5254
5355 def single_exp (t , tau , a ):
54- return (a * np .exp (- ((1.0 / tau )* t ) ))
56+ return (a * np .exp (- ((1.0 / tau )* t )))
5557
56- def double_exp (t , tau1 , a1 , tau2 , a2 ):
57- return ((single_exp (t , tau1 , a1 )) + (single_exp (t , tau2 , a2 )))
58+ def double_exp (t , tau1 , a1 , tau2 , a2 , noise ):
59+ return ((single_exp (t , tau1 , a1 )) + (single_exp (t , tau2 , a2 )) + noise )
5860
5961 def avg_tau_from_double_exp (tau1 , a1 , tau2 , a2 ):
6062 return (((tau1 * a1 ) + (tau2 * a2 ))/ (a1 + a2 ))
@@ -68,13 +70,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
6870 a1 = params [1 ]
6971 tau2 = params [2 ]
7072 a2 = params [3 ]
73+ noise = params [4 ]
7174
72- PL_sim = double_exp (t ,tau1 , a1 , tau2 , a2 )
75+ PL_sim = double_exp (t ,tau1 , a1 , tau2 , a2 , noise )
7376
74- Resid = np .sqrt ( np . sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (PL_sim )** 2 ) ))
77+ Resid = np .sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (TRPL )** 2 ))
7578 return Resid #returns the difference between the PL data and simulated data
7679
77- bounds = [tau1_bounds , a1_bounds , tau2_bounds , a2_bounds ]
80+ bounds = [tau1_bounds , a1_bounds , tau2_bounds , a2_bounds , noise ]
7881
7982 result = differential_evolution (residuals , bounds )
8083 return result .x
@@ -85,17 +88,18 @@ def residuals(params):#params are the parameters to be adjusted by differential
8588 a1 = p [1 ]
8689 tau2 = p [2 ]
8790 a2 = p [3 ]
91+ noise = p [4 ]
8892
89- PL_fit = double_exp (t , tau1 , a1 , tau2 , a2 )
93+ PL_fit = double_exp (t , tau1 , a1 , tau2 , a2 , noise )
9094
9195 avg_tau = avg_tau_from_double_exp (tau1 , a1 , tau2 , a2 )
9296
93- return tau1 , a1 , tau2 , a2 , avg_tau , PL_fit
97+ return tau1 , a1 , tau2 , a2 , avg_tau , PL_fit , noise
9498
95- def single_exp_fit (TRPL , t , tau_bounds = (0 ,1000 ), a_bounds = (0 ,1e5 )):
99+ def single_exp_fit (TRPL , t , tau_bounds = (0 ,10000 ), a_bounds = (0 ,1e6 ), noise = ( 0 , 1e6 )):
96100
97- def single_exp (t , tau , a ):
98- return (a * np .exp (- ((1.0 / tau )* t ) ))
101+ def single_exp (t , tau , a , noise ):
102+ return (a * np .exp (- ((1.0 / tau )* t ) ) + noise )
99103
100104 def Diff_Ev_Fit_singleExp (TRPL ):
101105 TRPL = TRPL
@@ -104,13 +108,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
104108 #Variable Rates
105109 tau = params [0 ]
106110 a = params [1 ]
111+ noise = params [2 ]
107112
108- PL_sim = single_exp (t , tau , a )
113+ PL_sim = single_exp (t , tau , a , noise )
109114
110- Resid = np .sqrt ( np . sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (PL_sim )** 2 ) ))
115+ Resid = np .sum (((PL_sim - TRPL )** 2 )/ (np .sqrt (TRPL )** 2 ))
111116 return Resid #returns the difference between the PL data and simulated data
112117
113- bounds = [tau_bounds , a_bounds ]
118+ bounds = [tau_bounds , a_bounds , noise ]
114119
115120 result = differential_evolution (residuals , bounds )
116121 return result .x
@@ -119,8 +124,9 @@ def residuals(params):#params are the parameters to be adjusted by differential
119124
120125 tau = p [0 ]
121126 a = p [1 ]
127+ noise = p [2 ]
122128
123- PL_fit = single_exp (t , tau , a )
129+ PL_fit = single_exp (t , tau , a , noise )
124130
125- return tau , a , PL_fit
131+ return tau , a , PL_fit , noise
126132
0 commit comments