Double source implementation (#87)

* Initial commit on double source

* Initial commit on double source**2

* Double source on all models

* Double source on all models

* Fix trajectories bug

* Fix import bug

* Updated pyLIMA_plots to include a plot_parameters_table and plot_distribution, also ML_fit now produces the bokeh distribution and bokeh parameters in the html output

* Updated pyLIMA_plots to include a plot_parameters_table and plot_distribution, also ML_fit now produces the bokeh distribution and bokeh parameters in the html output

* Likelihood fit correction

* Fix DoubleSource

* Fix some likelihood bugs

* Fix test

* Add likelihood test and fix tests

* Fix likelihood functions

* Fix likelihoods etc.

* Fix likelihood

* Extend space models graphs

* Small adjustement

* Add teff/tstar

* Rebranded fancy parameters

* Rebranded fancy parameters**2

* Correct Fancy

* Transform pyLIMA_parameters

* Correct pyLIMA_parameters

* FIxing grids and CROIN

* Fix rescale parameters + MCMC_chains

* Fix pyLIMASS MCMC:

* Replace xallarap with double source atttribute in models

* Fix unit tests

* Fix astrometry fits

* Implement ftotal

* Changing Isochrones to more Fe sampling

* Adding the old Isochrones

* Fix tests

* Fix LINT

* Fix LINT**2

---------

Co-authored-by: [email protected] <[email protected]>

examples/pyLIMA_example_1.py

@@ -96,7 +96,7 @@
 ### Let's try another fit with the differential evolution (DE) algorithm. 
 ### This will take longer... 
 
-my_fit2 = DE_fit.DEfit(pspl)
+my_fit2 = DE_fit.DEfit(pspl,loss_function='chi2')
 my_fit2.fit()
 
 ### Look at the results:
@@ -112,7 +112,7 @@
 
 ### You can still use the FITTING ALGORITHM that you imported previously. 
 ### Let's just use DE_fit for this:
-my_fit3 = DE_fit.DEfit(fspl)
+my_fit3 = DE_fit.DEfit(fspl, loss_function='chi2')
 my_fit3.fit()
 
 ### Let's see some plots. You can zoom close to the peak to see what is going on.
@@ -128,7 +128,7 @@
 your_event.telescopes[1].ld_gamma = 0.5
 
 ### Fit again:
-my_fit4 = DE_fit.DEfit(fspl)
+my_fit4 = DE_fit.DEfit(fspl, loss_function='chi2')
 my_fit4.fit()
 
 ### And plot it. Then zoom at the peak again.

examples/pyLIMA_example_6.py

@@ -66,6 +66,7 @@
 
 SLP.bounds[1] = [0.0,15]
 SLP.update_priors()
+
 ### Find maximum likelihood region
 
 results = SLP.de(popsize=10)
@@ -96,4 +97,3 @@
 
 plt.show()
 ### This concludes tutorial 6.
-

pyLIMA/astrometry/astrometric_positions.py

@@ -47,15 +47,15 @@ def astrometric_positions_of_the_source(telescope, pyLIMA_parameters, time_ref=N
     time = telescope.astrometry['time'].value
 
     if time_ref is None:
-        time_ref = pyLIMA_parameters.t0
+        time_ref = pyLIMA_parameters['t0']
 
-    ref_N = getattr(pyLIMA_parameters, 'position_source_N_' + telescope.name)
-    ref_E = getattr(pyLIMA_parameters, 'position_source_E_' + telescope.name)
-    mu_N = pyLIMA_parameters.mu_source_N
-    mu_E = pyLIMA_parameters.mu_source_E
+    ref_N = pyLIMA_parameters['position_source_N_' + telescope.name]
+    ref_E = pyLIMA_parameters['position_source_E_' + telescope.name]
+    mu_N = pyLIMA_parameters['mu_source_N']
+    mu_E = pyLIMA_parameters['mu_source_E']
 
     earth_vector = telescope.Earth_positions_projected['astrometry']
-    parallax_source = pyLIMA_parameters.pi_source
+    parallax_source = pyLIMA_parameters['pi_source']
     Earth_projected = earth_vector * parallax_source  # mas
 
     if telescope.astrometry['ra'].unit == 'deg':
@@ -135,21 +135,25 @@ def lens_astrometric_positions(model, telescope, pyLIMA_parameters, time_ref=Non
     source_EN = source_astrometric_positions(telescope, pyLIMA_parameters, shifts=None,
                                              time_ref=time_ref)
 
-    source_relative_to_lens = model.source_trajectory(telescope, pyLIMA_parameters,
-                                                      data_type='astrometry')
+    (source1_trajectory_x, source1_trajectory_y,
+    source2_trajectory_x, source2_trajectory_y,
+    dseparation, dalpha) = model.sources_trajectory(
+        telescope, pyLIMA_parameters,
+        data_type='astrometry')
 
     source_relative_to_lens_EN = xy_shifts_to_NE_shifts(
-        (source_relative_to_lens[0], source_relative_to_lens[1])
-        , pyLIMA_parameters.piEN, pyLIMA_parameters.piEE)
+        (source1_trajectory_x, source1_trajectory_y)
+        , pyLIMA_parameters['piEN'], pyLIMA_parameters['piEE'])
 
-    lens_EN = np.array(source_relative_to_lens_EN) * pyLIMA_parameters.theta_E
+    source_relative_to_lens_EN = (np.array(source_relative_to_lens_EN) *
+                                  pyLIMA_parameters['theta_E']) #mas
 
     if telescope.astrometry['ra'].unit == 'deg':
 
-        lens_EN = source_EN - lens_EN / 3600. / 1000.
+        lens_EN = source_EN - source_relative_to_lens_EN / 3600. / 1000.
 
     else:
 
-        lens_EN = source_EN - lens_EN / telescope.pixel_scale
+        lens_EN = source_EN - source_relative_to_lens_EN / telescope.pixel_scale
 
     return lens_EN

pyLIMA/fits/DE_fit.py

@@ -5,6 +5,7 @@
 import scipy
 from pyLIMA.fits.ML_fit import MLfit
 from tqdm import tqdm
+from pyLIMA.priors import parameters_priors
 
 
 class DEfit(MLfit):
@@ -49,6 +50,9 @@ def objective_function(self, fit_process_parameters):
     def fit(self, initial_population=[], computational_pool=None):
 
         start_time = python_time.time()
+        # Safety, recompute in case user changes boundaries after init
+        self.priors = parameters_priors.default_parameters_priors(
+            self.priors_parameters)
 
         if computational_pool:
 
@@ -60,7 +64,7 @@ def fit(self, initial_population=[], computational_pool=None):
 
         if initial_population == []:
 
-            init = 'latinhypercube'
+            init = 'sobol'
 
         else:
 
@@ -71,10 +75,10 @@ def fit(self, initial_population=[], computational_pool=None):
         differential_evolution_estimation = scipy.optimize.differential_evolution(
             self.objective_function,
             bounds=bounds,
-            mutation=(0.5, 1.5), popsize=int(self.DE_population_size),
+            mutation=(0.5, 1.0), popsize=int(self.DE_population_size),
             maxiter=self.max_iteration, tol=0.00,
             atol=1, strategy=self.strategy,
-            recombination=0.5, polish=False, init=init,
+            recombination=0.7, polish=False, init=init,
             disp=self.display_progress, workers=worker)
 
         self.trials = np.array(self.trials)
@@ -131,19 +135,14 @@ def fit_type(self):
 
     def objective_function(self, fit_process_parameters):
 
-        likelihood, pyLIMA_parameters = self.model_likelihood(fit_process_parameters)
-        likelihood *= -1
-
-        # Priors
-        priors = self.get_priors(fit_process_parameters)
-
-        likelihood += -priors
-
-        return likelihood
+        objective = self.standard_objective_function(fit_process_parameters)
+        return objective
 
     def fit(self, initial_population=[], computational_pool=None):
 
         start_time = python_time.time()
+        # Safety, recompute in case user changes boundaries after init
+        self.priors = parameters_priors.default_parameters_priors(self.fit_parameters)
 
         if computational_pool:
 
@@ -171,7 +170,7 @@ def fit(self, initial_population=[], computational_pool=None):
             strategy=self.strategy,
             recombination=0.5, polish=False,
             init=init, disp=self.display_progress,
-            workers=worker)
+            workers=worker,updating='deferred')
 
         if initial_population == []:
 

pyLIMA/fits/DREAM_fit.py

@@ -369,7 +369,7 @@ def fit(self, initial_population=[], computational_pool=None):
         self.acceptance = np.array(all_acceptance)
         DEMC_population = np.copy(self.population)
         self.Z = Z_prime
-
+        breakpoint()
         print(self.swap / self.max_iteration)
         computation_time = python_time.time() - start_time
         print(sys._getframe().f_code.co_name, ' : ' + self.fit_type() + ' fit SUCCESS')

pyLIMA/fits/GRIDS_fit.py

@@ -2,8 +2,8 @@
 import time as python_time
 
 import numpy as np
-import scipy.optimize as so
 from pyLIMA.fits.ML_fit import MLfit
+from pyLIMA.fits import DE_fit
 from tqdm import tqdm
 
 
@@ -21,7 +21,7 @@ class GRIDfit(MLfit):
     grid_resolution : int, the resolution of the grid for each grid parameters
     """
     def __init__(self, model, rescale_photometry=False, rescale_astrometry=False,
-                 telescopes_fluxes_method='polyfit', DE_population_size=2,
+                 telescopes_fluxes_method='polyfit', DE_population_size=5,
                  max_iteration=2000,
                  fix_parameters=[], grid_resolution=10):
         """The fit class has to be intialized with an event object."""
@@ -34,33 +34,11 @@ def __init__(self, model, rescale_photometry=False, rescale_astrometry=False,
         self.max_iteration = max_iteration
         self.fix_parameters = fix_parameters
         self.grid_resolution = grid_resolution
+        self.intervals = []
 
     def fit_type(self):
         return "Grids"
 
-    def reconstruct_fit_process_parameters(self, moving_parameters, fix_parameters):
-        """
-        """
-
-        fit_process_parameters = []
-
-        ind_fix = 0
-        ind_move = 0
-
-        for key in list(self.fit_parameters.keys()):
-
-            if key not in self.fix_parameters:
-
-                fit_process_parameters.append(moving_parameters[ind_move])
-                ind_move += 1
-
-            else:
-
-                fit_process_parameters.append(fix_parameters[ind_fix])
-                ind_fix += 1
-
-        return np.array(fit_process_parameters)
-
     def construct_the_hyper_grid(self):
         """Define the grid. ON CONSTRUCTION.
         """
@@ -69,11 +47,10 @@ def construct_the_hyper_grid(self):
 
         for parameter_name in self.fix_parameters:
             parameter_range = self.fit_parameters[parameter_name][1]
+            self.intervals.append((parameter_range[1]-parameter_range[0])/(self.grid_resolution))
 
-            parameters_on_the_grid.append(
-                np.linspace(parameter_range[0], parameter_range[1],
-                            self.grid_resolution))
-
+            parameters_on_the_grid.append( np.arange(parameter_range[0], parameter_range[1],
+                                                     self.intervals[-1]))
         parameters_on_the_grid = np.array(parameters_on_the_grid)
         params = map(np.asarray, parameters_on_the_grid)
         grid = np.broadcast_arrays(
@@ -85,76 +62,78 @@ def construct_the_hyper_grid(self):
 
     def objective_function(self, fit_process_parameters):
 
-        likelihood = -self.model_likelihood(fit_process_parameters)
-
-        # Priors
-        priors = self.get_priors(fit_process_parameters)
+        objective = self.standard_objective_function(fit_process_parameters)
 
-        likelihood += -priors
-
-        return likelihood
+        return objective
 
     def fit_on_grid_pixel(self, *fixed_parameters):
 
+        fixed_parameters, computational_pool = fixed_parameters[0]
+
         fixed_parameters = np.ravel(fixed_parameters)
-        differential_evolution_estimation = so.differential_evolution(
-            self.objective_function,
-            bounds=self.bounds,
-            mutation=(0.5, 1.5),
-            popsize=self.DE_population_size,
-            args=([fixed_parameters]),
-            maxiter=self.max_iteration, tol=0.00,
-            atol=1.0, strategy='rand1bin',
-            recombination=0.5, polish=True,
-            init='latinhypercube',
-            disp=False)
-
-        fitted_parameters = differential_evolution_estimation['x']
-        best_model = self.reconstruct_fit_process_parameters(fitted_parameters,
-                                                             fixed_parameters)
-        best_model = np.append(best_model, differential_evolution_estimation['fun'])
+        defit = DE_fit.DEfit(self.model,DE_population_size= self.DE_population_size,display_progress=False,
+                             strategy='best1bin', loss_function='chi2',max_iteration=self.max_iteration)
+
+        for key in self.fit_parameters:
+
+            defit.fit_parameters[key][1] = self.fit_parameters[key][1]
+
+        for ind, key in enumerate(self.fix_parameters):
+
+            defit.fit_parameters[key][1] = [fixed_parameters[ind]+self.intervals[ind]/2, fixed_parameters[ind]+self.intervals[ind]/2]
+
+        defit.fit(computational_pool=computational_pool)
+        fitted_parameters = defit.fit_results['best_model']
+        best_model = np.append( fitted_parameters,self.objective_function(fitted_parameters))
         return best_model
 
     def fit(self, computational_pool=None):
 
-        hyper_grid = self.construct_the_hyper_grid()
 
+        hyper_grid = self.construct_the_hyper_grid()
         start_time = python_time.time()
 
-        self.bounds = [self.fit_parameters[key][1] for key in self.fit_parameters.keys()
-                       if key not in self.fix_parameters]
+        self.bounds = [self.fit_parameters[key][1] for key in self.fit_parameters.keys()]
+
+        #if computational_pool is not None:
+
+            #with computational_pool as pool, tqdm(total=len(hyper_grid)) as pbar:
+
+            #    res = [pool.apply_async(self.fit_on_grid_pixel, args=(hyper_grid[i],),
+            #                            callback=lambda _: pbar.update(1)) for i in
+            #           range(len(hyper_grid))]
 
-        if computational_pool is not None:
+            #    population = np.array([r.get() for r in res])
 
-            with computational_pool as pool, tqdm(total=len(hyper_grid)) as pbar:
+        #else:
+        #    population = []
 
-                res = [pool.apply_async(self.fit_on_grid_pixel, args=(hyper_grid[i],),
-                                        callback=lambda _: pbar.update(1)) for i in
-                       range(len(hyper_grid))]
+        #    for j in tqdm(range(len(hyper_grid))):
+        #        new_step = self.fit_on_grid_pixel([hyper_grid[j]])
+        #        population.append(new_step)
 
-                population = np.array([r.get() for r in res])
+        ###This is not clean, can not have a parallelization of the grid because of the Manage List of the DE fit
 
-        else:
-            population = []
+        population = []
 
-            for j in tqdm(range(self.max_iteration)):
-                new_step = self.fit_on_grid_pixel([hyper_grid[j]])
-                population.append(new_step)
+        for j in tqdm(range(len(hyper_grid))):
+            new_step = self.fit_on_grid_pixel([hyper_grid[j], computational_pool])
+            population.append(new_step)
 
-            GRIDS_population = np.array(population)
+        GRIDS_population = np.array(population)
 
-            computation_time = python_time.time() - start_time
-            print(sys._getframe().f_code.co_name,
-                  ' : ' + self.fit_type() + ' fit SUCCESS')
+        computation_time = python_time.time() - start_time
+        print(sys._getframe().f_code.co_name,
+              ' : ' + self.fit_type() + ' fit SUCCESS')
 
-            best_model_index = GRIDS_population[:, -1].argmin()
+        best_model_index = GRIDS_population[:, -1].argmin()
 
-            print('best_model:', GRIDS_population[best_model_index, :-1],
-                  self.loss_function, GRIDS_population[best_model_index, -1])
+        print('best_model:', GRIDS_population[best_model_index, :-1],
+              self.loss_function, GRIDS_population[best_model_index, -1])
 
-            self.fit_results = {'best_model': GRIDS_population[best_model_index, :-1],
-                                self.loss_function: GRIDS_population[
-                                    best_model_index, -1],
-                                'fit_time': computation_time,
-                                'GRIDS_population': GRIDS_population}
+        self.fit_results = {'best_model': GRIDS_population[best_model_index, :-1],
+                            self.loss_function: GRIDS_population[
+                                best_model_index, -1],
+                            'fit_time': computation_time,
+                            'GRIDS_population': GRIDS_population}
 

pyLIMA/fits/LM_fit.py

@@ -27,6 +27,7 @@ def __init__(self, model, telescopes_fluxes_method='fit', loss_function='chi2'):
                          loss_function=loss_function)
 
         self.guess = []
+        #self.priors = None
 
     def fit_type(self):
         return "Levenberg-Marquardt"