Skip to content

Commit

Permalink
Bug fixing
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebachelet
ebachelet committed Oct 3, 2024
1 parent 71b55a6 commit 3a1fb70
Show file tree
Hide file tree
Showing 10 changed files with 200 additions and 140 deletions.
82 changes: 41 additions & 41 deletions examples/pyLIMA_example_3.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,8 @@
"outputs": [],
"source": [
"### Import the required libraries.\n",
"%matplotlib widget\n",
"\n",
"import matplotlib.pyplot as plt\n",
"from pyLIMA.fits import DE_fit\n",
"from pyLIMA.models import DSPL_model\n",
"from pyLIMA.models import PSPL_model\n",
"from pyLIMA.outputs import pyLIMA_plots\n",
"### Import the simulator to be used for generating the simulated light curve\n",
Expand Down Expand Up @@ -83,8 +80,8 @@
"metadata": {},
"outputs": [],
"source": [
"CTIO_I = simulator.simulate_a_telescope(name='CTIO_I', time_start=2457365.5, \n",
" time_end=2457965.5, sampling=4, \n",
"CTIO_I = simulator.simulate_a_telescope(name='CTIO_I', time_start=2457365.5,\n",
" time_end=2457965.5, sampling=4,\n",
" location='Earth', camera_filter='I',\n",
" uniform_sampling=True, astrometry=False)"
]
Expand Down Expand Up @@ -163,7 +160,7 @@
"outputs": [],
"source": [
"pspl_parameters = simulator.simulate_microlensing_model_parameters(pspl)\n",
"print (pspl_parameters)"
"print(pspl_parameters)"
]
},
{
Expand Down Expand Up @@ -255,10 +252,10 @@
"outputs": [],
"source": [
"plt.close('all')\n",
"plt.errorbar(CTIO_I.lightcurve_magnitude['time'].value-2450000,\n",
"plt.errorbar(CTIO_I.lightcurve_magnitude['time'].value - 2450000,\n",
" CTIO_I.lightcurve_magnitude['mag'].value,\n",
" yerr = CTIO_I.lightcurve_magnitude['err_mag'].value,\n",
" fmt = '.', label=CTIO_I.name)\n",
" yerr=CTIO_I.lightcurve_magnitude['err_mag'].value,\n",
" fmt='.', label=CTIO_I.name)\n",
"\n",
"plt.gca().invert_yaxis()\n",
"plt.legend()\n",
Expand Down Expand Up @@ -315,43 +312,44 @@
"metadata": {},
"outputs": [],
"source": [
"SAAO_I = simulator.simulate_a_telescope(name='SAAO_I', time_start=2457575.5, \n",
" time_end=2457625.5, sampling=2.5, \n",
"SAAO_I = simulator.simulate_a_telescope(name='SAAO_I', time_start=2457575.5,\n",
" time_end=2457625.5, sampling=2.5,\n",
" location='Earth', camera_filter='I',\n",
" uniform_sampling=False, altitude=400, \n",
" longitude = 20.659279, \n",
" latitude = -32.3959, \n",
" bad_weather_percentage=20.0 / 100, \n",
" moon_windows_avoidance=20, minimum_alt=15, \n",
" uniform_sampling=False, altitude=400,\n",
" longitude=20.659279,\n",
" latitude=-32.3959,\n",
" bad_weather_percentage=20.0 / 100,\n",
" moon_windows_avoidance=20, minimum_alt=15,\n",
" astrometry=False)\n",
"\n",
"SSO_I = simulator.simulate_a_telescope(name='SSO_I', time_start=2457535.5, \n",
" time_end=2457645.5, sampling=2.5, \n",
"SSO_I = simulator.simulate_a_telescope('SSO_I', time_start=2457535.5,\n",
" time_end=2457645.5, sampling=2.5,\n",
" location='Earth', camera_filter='I',\n",
" uniform_sampling=False, altitude=1165, \n",
" longitude = 149.0685, \n",
" latitude = -31.2749, \n",
" bad_weather_percentage=35.0 / 100, \n",
" moon_windows_avoidance=20, minimum_alt=15, \n",
" uniform_sampling=False, altitude=1165,\n",
" longitude=149.0685,\n",
" latitude=-31.2749,\n",
" bad_weather_percentage=35.0 / 100,\n",
" moon_windows_avoidance=20, minimum_alt=15,\n",
" astrometry=False)\n",
"\n",
"CTIO_I = simulator.simulate_a_telescope(name='CTIO_I', time_start=2457365.5, \n",
" time_end=2457965.5, sampling=4.5, \n",
"CTIO_I = simulator.simulate_a_telescope('CTIO_I', time_start=2457365.5,\n",
" time_end=2457965.5, sampling=4.5,\n",
" location='Earth', camera_filter='I',\n",
" uniform_sampling=False, altitude=1000, \n",
" longitude = -109.285399, \n",
" latitude = -27.130, \n",
" bad_weather_percentage=10.0 / 100, \n",
" moon_windows_avoidance=30, minimum_alt=30, \n",
" uniform_sampling=False, altitude=1000,\n",
" longitude=-109.285399,\n",
" latitude=-27.130,\n",
" bad_weather_percentage=10.0 / 100,\n",
" moon_windows_avoidance=30, minimum_alt=30,\n",
" astrometry=False)\n",
"\n",
"CTIO_V = simulator.simulate_a_telescope(name='CTIO_V', time_start=2457365.5, \n",
" time_end=2457965.5, sampling=24.5, \n",
"CTIO_V = simulator.simulate_a_telescope('CTIO_V', time_start=2457365.5,\n",
" time_end=2457965.5, sampling=24.5,\n",
" location='Earth', camera_filter='V',\n",
" uniform_sampling=False, altitude=1000, \n",
" longitude = -109.285399, \n",
" latitude = -27.130, bad_weather_percentage=10.0 / 100, \n",
" moon_windows_avoidance=30, minimum_alt=30, \n",
" uniform_sampling=False, altitude=1000,\n",
" longitude=-109.285399,\n",
" latitude=-27.130,\n",
" bad_weather_percentage=10.0 / 100,\n",
" moon_windows_avoidance=30, minimum_alt=30,\n",
" astrometry=False)"
]
},
Expand Down Expand Up @@ -442,7 +440,7 @@
"metadata": {},
"outputs": [],
"source": [
"dspl = DSPL_model.DSPLmodel(your_event2)"
"dspl = PSPL_model.PSPLmodel(your_event2, double_source=['Static',2457500])"
]
},
{
Expand All @@ -452,7 +450,7 @@
"source": [
"Now that the MODEL is there, we need to set the relevant parameters.\n",
"\n",
"The parameters are drawn uniformly from the bounds defined but you can also set them manually. Please consult the documentation for more details on the parameters of the MODEL you want to use. For the DSPL example, dspl_parameters = [to, uo, delta_to, delta_uo, tE, q_fluxr_1, q_fluxr2, ...]\n",
"The parameters are drawn uniformly from the bounds defined but you can also set them manually. Please consult the documentation for more details on the parameters of the MODEL you want to use. For the DSPL example, dspl_parameters = [to, uo, tE,delta_to, delta_uo, q_fluxr_1, q_fluxr2, ...]\n",
"where q_fluxr_* is the flux ratio in each observing band."
]
},
Expand All @@ -464,7 +462,7 @@
"outputs": [],
"source": [
"dspl_parameters = simulator.simulate_microlensing_model_parameters(dspl)\n",
"print (dspl_parameters)"
"print(dspl_parameters)"
]
},
{
Expand Down Expand Up @@ -590,7 +588,9 @@
"metadata": {},
"outputs": [],
"source": [
"dspl_parameters[0:7] = [2457760.216627234, 0.8605811108889658, 143.4484970433387, -0.6046788112617074, 116.43231096591524, 0.15157064165919296, 0.18958495421162946]"
"dspl_parameters[0:7] = [2457760.216627234, 0.8605811108889658, 116.43231096591524, 143.4484970433387,\n",
" -0.6046788112617074, 0.15157064165919296,\n",
" 0.18958495421162946]"
]
},
{
Expand Down Expand Up @@ -872,7 +872,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.15"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
14 changes: 6 additions & 8 deletions examples/pyLIMA_example_3.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
### Import the required libraries.
import matplotlib.pyplot as plt
from pyLIMA.fits import DE_fit
from pyLIMA.models import DSPL_model
from pyLIMA.models import PSPL_model
from pyLIMA.outputs import pyLIMA_plots
### Import the simulator to be used for generating the simulated light curve
Expand Down Expand Up @@ -166,13 +165,13 @@
### link it to the EVENT you prepared.
### We will use the double-source point-lens (DSPL) model for this example.

dspl = DSPL_model.DSPLmodel(your_event2)
dspl = PSPL_model.PSPLmodel(your_event2, double_source=['Static',2457500])

### Now that the MODEL is there, we need to set the relevant parameters.
### The parameters are drawn uniformly from the bounds defined but you can
### also set them manually. Please consult the documentation for more
### details on the parameters of the MODEL you want to use. For the DSPL example,
### dspl_parameters = [to, uo, delta_to, delta_uo, tE, q_fluxr_1, q_fluxr2, ...]
### dspl_parameters = [to, uo, tE, elta_to, delta_uo, q_fluxr_1, q_fluxr2, ...]
### where q_fluxr_* is the flux ratio in each observing band.
dspl_parameters = simulator.simulate_microlensing_model_parameters(dspl)
print(dspl_parameters)
Expand Down Expand Up @@ -219,8 +218,8 @@
### Here's how you can do that. Let's fix the DSPL parameters to some values where
### the binary source model produces two clear peaks, and then just adjust the flux
### parameters.
dspl_parameters[0:7] = [2457760.216627234, 0.8605811108889658, 143.4484970433387,
-0.6046788112617074, 116.43231096591524, 0.15157064165919296,
dspl_parameters[0:7] = [2457760.216627234, 0.8605811108889658, 116.43231096591524, 143.4484970433387,
-0.6046788112617074, 0.15157064165919296,
0.18958495421162946]

### The order of the parameters is:
Expand Down Expand Up @@ -265,11 +264,11 @@
print(dspl_parameters)
parameter_commentary = ['Time of minimum impact parameter for source 1',
'minimum impact parameter for source 1',
'angular Einstein radius crossing time',
'difference of time of minimum impact parameter between the '
'two sources',
'difference of minimum impact parameters between the two '
'sources',
'angular Einstein radius crossing time',
'flux ratio in I between source 1 and source 2',
'flux ratio in V between source 1 and source 2',
'source flux of source 1 for telescope CTIO_I (survey '
Expand All @@ -291,7 +290,7 @@
### Let's try to fit this now! (This can take a while!)
### You can check the first tutorial again for a detailed explanation if needed.

my_fit = DE_fit.DEfit(dspl, display_progress=True, strategy='best1bin')
my_fit = DE_fit.DEfit(dspl, display_progress=True, loss_function='likelihood',strategy='best1bin')
my_fit.fit()

my_fit.fit_results['best_model']
Expand All @@ -303,5 +302,4 @@
pyLIMA_plots.plot_lightcurves(dspl, my_fit.fit_results['best_model'])
pyLIMA_plots.plot_lightcurves(dspl, dspl_parameters)
plt.show()

### This concludes tutorial 3.
8 changes: 5 additions & 3 deletions pyLIMA/fits/DE_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,11 +90,13 @@ def fit(self, initial_population=[], computational_pool=None):
print('DE converge to parameters : = ',
differential_evolution_estimation['x'].astype(str))

results = differential_evolution_estimation['x']
results_log_likelihood = differential_evolution_estimation['fun']

DE_population = np.c_[self.trials_parameters,self.trials_objective,self.trials_priors]

best = np.where(np.all(differential_evolution_estimation['x']==DE_population[:,:len(self.fit_parameters)]
,axis=1))[0][0]
results = DE_population[best,:-2]
results_log_likelihood = differential_evolution_estimation['fun']

computation_time = python_time.time() - start_time
print(sys._getframe().f_code.co_name, ' : ' + self.fit_type() + ' fit SUCCESS')

Expand Down
18 changes: 10 additions & 8 deletions pyLIMA/fits/MCMC_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,9 +40,11 @@ def objective_function(self, fit_process_parameters):

if limits_check is not None:

#self.trials_parameters.append(fit_process_parameters.tolist())
#self.trials_priors.append(np.inf)
#self.trials_objective.append(np.inf)
bad_parameters = np.zeros(len(self.priors_parameters))
bad_parameters[:len(self.fit_parameters)] = fit_process_parameters
self.trials_parameters.append(bad_parameters.tolist())
self.trials_priors.append(np.inf)
self.trials_objective.append(np.inf)

return -limits_check #i.e. -np.inf

Expand Down Expand Up @@ -177,12 +179,12 @@ def reconstruct_chains(self, mcmc_samples, mcmc_prob):

for unique_values in unique_sample[0]:

index = np.where(self.trials_parameters[:, :len(self.fit_parameters)][:, -1]
== unique_values[-1])[0][0]
index = np.where(np.all(self.trials_parameters[:, :len(self.fit_parameters)]
== unique_values,axis=1))[0][0]

unique_trials.append(self.trials_parameters[index].tolist())
unique_objective.append(self.trials_objective[index].tolist())
unique_priors.append(self.trials_priors[index].tolist())
unique_trials.append(self.trials_parameters[index].tolist())
unique_objective.append(self.trials_objective[index].tolist())
unique_priors.append(self.trials_priors[index].tolist())

MCMC_chains_with_fluxes[:,j][:,:-2] = np.array(unique_trials)[unique_sample[1].ravel()]
MCMC_chains_with_fluxes[:,j][:,-2] = np.array(unique_objective)[unique_sample[1].ravel()]
Expand Down
1 change: 0 additions & 1 deletion pyLIMA/fits/ML_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,6 @@ def define_parameters(self, include_telescopes_fluxes=True):
fit_parameters_indexes.append(ind)
fit_parameters_boundaries.append(theboundaries)


if self.rescale_photometry:

for telescope in self.model.event.telescopes:
Expand Down
1 change: 1 addition & 0 deletions pyLIMA/fits/TRF_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ def fit(self):

for telescope in self.model.event.telescopes:
n_data = n_data + telescope.n_data('flux')
n_data = n_data + telescope.n_data('astrometry')

if self.model.Jacobian_flag != 'Numerical':

Expand Down
Loading

0 comments on commit 3a1fb70

Please sign in to comment.