andreatramacere / jetset

JetSeT a framework for self-consistent modeling and fitting of astrophysical relativistic jets
BSD 3-Clause "New" or "Revised" License
30 stars 15 forks source link

Latest version (1.3.0) issue with "minuit" minimizer #84

Open Ayon2001 opened 3 weeks ago

Ayon2001 commented 3 weeks ago

BL_Lac_Two_Zone_EC.pdf [Uploading BL_Lac_Two_Zone_EC.pdf…]() @andreatramacere

### Latest version (1.3.0) issue with "minuit" minimizer

Why this issue is happening? Till yesterday code is working fine. From today, the same code is giving this error. "lsb" minimizer is working properly. Issue is happening with "minuit" minimizer.

**_25 model_minimizer_lsb=ModelMinimizer('minuit')
---> 26 best_fit=model_minimizer_lsb.fit(composite_jet,sed_data,1E11,1E29,fitname='SSC-best-fit-minuit',repeat=2)_**

File ~/mambaforge/envs/jetset_latest/lib/python3.10/site-packages/jetset/minimizer.py:465, in ModelMinimizer.fit(self, fit_model, sed_data, nu_fit_start, nu_fit_stop, fitname, fit_workplace, loglog, silent, get_conf_int, max_ev, use_fake_err, use_UL, skip_minimizer, repeat)
    463     self.asymm_errors = self.minimizer.asymm_errors
    464 self.covar=self.minimizer.covar
--> 465 self.corr=self.minimizer.corr
    467 self.reset_to_best_fit()
    468 self.minimizer._fit_stats()

File ~/mambaforge/envs/jetset_latest/lib/python3.10/site-packages/jetset/minimizer.py:699, in Minimizer.corr(self)
    696 @property
    697 def corr(self):
    698     if self.covar is not None:
--> 699         v = np.sqrt(np.diag(self.covar))
    700         outer_v = np.outer(v, v)
    701         correlation = self.covar / outer_v

File ~/mambaforge/envs/jetset_latest/lib/python3.10/site-packages/numpy/lib/twodim_base.py:303, in diag(v, k)
    301     return diagonal(v, k)
    302 else:
--> 303     raise ValueError("Input must be 1- or 2-d.")

ValueError: Input must be 1- or 2-d.
andreatramacere commented 3 weeks ago

It is possible that there is an issue with the covariance matrix due to the changes in the fit range, maybe you have included data with zero error, you were not using previously. I will wrap the correlation computation within a try/except clause, to just raise a warning in the next version. If you do not want to wait for the next release you can install from source when ready. I guess you added this data:

15 GHz MOJAVE DATA

1.5E10 0 1.6446E-12 1.6446E-13

4.3E10 0 6.755945E-9 0

which has zero error. Please, se a reasonable error and let me know if the issue is still there

Ayon2001 commented 3 weeks ago

@andreatramacere

No, I had previously used the same data with no changes.

Additionally, Now I attempted to apply reasonable error bars, I encountered the same error.

Is it possible to solve the issue now ?

andreatramacere commented 3 weeks ago

It is strange that the frequentist fit changes behaviour re-running it, unless you did not change the data or the workflow. I can do a quick. I can do a quick fix you can install from source, otherwise, you can do as follows 1) run the minuit minimizer without the repeat argument i.e.: model_minimizer_lsb.fit(composite_jet,sed_data,1E11,1E29,fitname='SSC-best-fit-minuit') 2) repeat that step the times you want e.g. 2 or 3

Since the creation of the correlation matrix does not impact on the fit results you can simply replace your old instruction with

n_repeat=2
for i in range(n_repeat):
    try:
        best_fit=model_minimizer_lsb.fit(composite_jet,sed_data,1E11,1E29,fitname='SSC-best-fit-minuit',repeat=2)
    except:
        model_minimizer_lsb.reset_to_best_fit()
        model_minimizer_lsb.minimizer._fit_stats()
        best_fit=model_minimizer_lsb.get_fit_results(composite_jet,1E11,1E29,fitname='SSC-best-fit-minuit',silent=False)
Ayon2001 commented 3 weeks ago

@andreatramacere Thank you! It works. The parameters and fit range are very sensitive.

I just want to ask if it's possible to change the legends(such as font size etc) in the plot. I'd like to change the line colors and styles. Also, is it possible to assign different colors to different data categories (other than the spectral indices as mentioned in the plot), such as Fermi, X-ray, and Optical data, so that it's easier to understand which data is from which telescope, easily understandable visually from the final fit results?

How can one do that?