u

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

src/yapsut/montecarlo_fitting.py

+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