From 3f16e7c18538f662f18e3f9a890796b920cf763e Mon Sep 17 00:00:00 2001
From: "Michele.Maris" <michele.maris@inaf.it>
Date: Tue, 11 Jun 2024 20:20:58 +0200
Subject: [PATCH] u
---
src/yapsut/__init__.py | 1 +
src/yapsut/montecarlo_fitting.py | 187 +++++++++++++++++++++++++++++++
2 files changed, 188 insertions(+)
create mode 100644 src/yapsut/montecarlo_fitting.py
diff --git a/src/yapsut/__init__.py b/src/yapsut/__init__.py
index 4e2fca0..86a86b4 100644
--- a/src/yapsut/__init__.py
+++ b/src/yapsut/__init__.py
@@ -18,6 +18,7 @@ from .template_subs import discoverKeys, templateFiller
from .periodic_stats import periodic_stats, periodic_centroid
from .stats import CumulativeOfData
from .colored_noise import gaussian_colored_noise
+from .montecarlo_fitting import EnsembleFitting_Base
# planck legacy is under editing
#from .planck_legacy import cosmological_dipole
diff --git a/src/yapsut/montecarlo_fitting.py b/src/yapsut/montecarlo_fitting.py
new file mode 100644
index 0000000..f975a60
--- /dev/null
+++ b/src/yapsut/montecarlo_fitting.py
@@ -0,0 +1,187 @@
+import numpy as np
+import pandas as pd
+
+class EnsembleFitting_Base :
+ """This class implements Ensemble fitting with curve_fit for (a,b,c,Rdark) out of the input curve.
+ """
+ @property
+ def initial_values(self) :
+ """return the initial values"""
+ return self._initial_values
+ #
+ @initial_values.setter
+ def initial_values(self,this) :
+ """return the initial values"""
+ self._initial_values=this
+ #
+ @property
+ def param_names(self) :
+ return self._param_names
+ #
+ @property
+ def Model(self) :
+ return self._Model
+ #
+ @Model.setter
+ def Model(self,Model,param_names) :
+ self._param_names=param_names
+ self._Model=this
+ #
+ @property
+ def fit_success(self) :
+ return self._success
+ #
+ @property
+ def size(self) :
+ return self._size
+ #
+ @property
+ def fit(self) :
+ return self._fit
+ #
+ def copy(self) :
+ from copy import deepcopy
+ return copy(self)
+ #
+ def __len__(self) :
+ return self._size
+ #
+ def __init__(self,X,Xerr,Y,Yerr) :
+ """It takes in input X, Y and their errors"""
+ self._X=X
+ self._Xerr=Xerr if not Xerr is None else X*0
+ #
+ self._Y=Y
+ self._Yerr=Yerr if not Yerr is None else Y*0
+ #
+ self._size=len(X)
+ #
+ self._fitting_clean()
+ #
+ def _fitting_clean(self) :
+ self._mc=None
+ self._initial_values=None
+ self._Model=None
+ self._param_names=None
+ self._success=False
+ self._fit=None
+ self._fitting_kargs={}
+ #
+ @property
+ def fitting_kargs(self) :
+ """access to the fitting arguments dictionary"""
+ return self._fitting_kargs
+ #
+ @property
+ def fitting_result(self) :
+ """returns the whole last fitting result"""
+ return self._fitting_result
+ #
+ def fitting_bounds(self) :
+ """returns fitting bounds"""
+ return self._fitting_args['bounds']
+ #
+ @property
+ def residuals(self) :
+ if not self.fit_success :
+ return np.nan*self._X
+ return self._Y-self._Model(self._X,*self._fit)
+ #
+ @property
+ def chisq(self) :
+ if not self.fit_success :
+ return np.nan
+ r=(self.residuals/self._Yerr)**2
+ return r.sum()
+ #
+ def mc_shot(self) :
+ """creates a montecarlo realization of the I,P data.
+ Assumes uniform error for I and gaussian error for P
+ """
+ n=self._size
+ X=np.random.uniform(low=-0.5,high=0.5,size=n)*self._Xerr+self._X
+ Y=np.random.normal(size=n)*self._Yerr+self._Y
+ return X,Y
+ #
+ def montecarlo(self,nmc,force_positive=True,as_pandas=True,max_nmc=-3) :
+ """creates a montecarlo serie of parameters
+ output: mc realization of parameters
+ """
+ out=[]
+ n=len(self)
+ imc=nmc
+ mcount=abs(max_nmc)*nmc if max_nmc < 0 else max_nmc
+ if nmc == 0 : mcount = nmc
+ while imc > 0 and mcount > 0:
+ mcount-=1
+ #
+ X,Y=self.mc_shot()
+ fit,fit_success=self.fitting(X,Y,self._Yerr)
+ #
+ if fit_success :
+ if self.valid_fit(fit) :
+ imc+=-1
+ out.append(fit)
+ self._mc=np.array(out)
+ if as_pandas :
+ out=pd.DataFrame(self._mc,columns=self._param_names)
+ return out
+ #
+ @property
+ def montecarlo_ensemble_averages(self) :
+ """returns the ensemble averaged values from the montecarlo simulation"""
+ try :
+ return np.array([np.median(k) for k in self._mc.T])
+ except :
+ return np.ones(len(self._param_names))+np.nan
+ #
+ def fitting(self,X,Y,sgm) :
+ """
+ the fitter function, must be specialized according to the this template
+
+def fitting(self,X,Y,sgm) :
+ from scipy.optimize import curve_fit
+ #
+ # resets all the variables to be generated
+ self._success=False
+ self._fit=None
+ self._fitting_result=None
+ #
+ # perfom the fit
+ try :
+ self._fitting_result=curve_fit(self.Model,X,Y,sigma=sgm,p0=self._initial_values,full_output=True,**self._fitting_kargs)
+ #
+ # saves the fit result and the flags success
+ self._fit=self._fitting_result[0]
+ self._success=self._fitting_result[-1] in [1,2,3,4]
+ except :
+ self._success=False
+ #
+ # returns the result
+ return self._fit,self._success
+
+ """
+ from scipy.optimize import curve_fit
+ #
+ # resets all the variables to be generated
+ self._success=False
+ self._fit=None
+ self._fitting_result=None
+ #
+ # perfom the fit
+ try :
+ self._fitting_result=curve_fit(self.Model,X,Y,sigma=sgm,p0=self._initial_values,full_output=True,**self._fitting_kargs)
+ #
+ # saves the fit result and the flags success
+ self._fit=self._fitting_result[0]
+ self._success=self._fitting_result[-1] in [1,2,3,4]
+ except :
+ self._success=False
+ #
+ # returns the result
+ return self._fit,self._success
+ #
+ def valid_fit(self,fit) :
+ """test wether the fit is valid: a specialized function for the application"""
+ #return True or False
+ return True
--
GitLab