correct fitting with some models

ptiming_ana/_version.py

@@ -12,5 +12,5 @@
 __version_tuple__: VERSION_TUPLE
 version_tuple: VERSION_TUPLE
 
-__version__ = version = '0.3.dev120+g19d6b46.d20240115'
-__version_tuple__ = version_tuple = (0, 3, 'dev120', 'g19d6b46.d20240115')
+__version__ = version = '0.3.dev122+g8b32d86.d20240124'
+__version_tuple__ = version_tuple = (0, 3, 'dev122', 'g8b32d86.d20240124')

ptiming_ana/phaseogram/lightcurve.py

@@ -109,18 +109,21 @@ def draw_fitting(self,pulsar_phases,color,density=False,label=None):
         if pulsar_phases.fitting.shift==0:
             x=np.linspace(0,1,100000)
         else:
-            x=np.linspace(2*pulsar_phases.fitting.shift,1+2*pulsar_phases.fitting.shift,100000)
+            x=np.linspace(pulsar_phases.fitting.shift,1+pulsar_phases.fitting.shift,100000)
 
         #Calculate distribution y-points using a model and its fitted params
         hoff=np.mean(self.lc[0])
         if pulsar_phases.fitting.model=='dgaussian':
             y=double_gaussian(x, *pulsar_phases.fitting.params[0:7])
+
+        elif pulsar_phases.fitting.model=='tgaussian':
+            y=triple_gaussian(x, *pulsar_phases.fitting.params)
 
         elif pulsar_phases.fitting.model=='asym_dgaussian':
-            assymetric_gaussian_pdf_vec=np.vectorize(assymetric_double_gaussian)
-            y=assymetric_gaussian_pdf_vec(x, *pulsar_phases.fitting.params)
+            assymetric_double_gaussian_vec=np.vectorize(assymetric_double_gaussian,excluded = [1,2,3,4,5,6,7,8,9,10])
+            y=assymetric_double_gaussian_vec(x, *pulsar_phases.fitting.params)
 
-        elif pulsar_phases.fitting.model=='lorentzian':
+        elif pulsar_phases.fitting.model=='double_lorentz':
             y=double_lorentz(x, *pulsar_phases.fitting.params)
 
         elif pulsar_phases.fitting.model=='gaussian':
@@ -137,8 +140,8 @@ def draw_fitting(self,pulsar_phases,color,density=False,label=None):
         else:
             plt.plot(x,y,color=color,label='Fit')
 
-        plt.plot(x+1,y,color=color)
-        plt.plot(x-1,y,color=color)
+        plt.plot(x+1,y,color=color,linewidth=2)
+        plt.plot(x-1,y,color=color,linewidth=2)
 
         #Add labels
         plt.xlabel('Pulsar phase',fontsize=10)

ptiming_ana/phaseogram/models.py

@@ -26,13 +26,12 @@ def triple_gaussian(x, Bkg, mu, sigma,mu_2,sigma_2,mu_3,sigma_3, A,B,C):
 @nb.njit(**kwd)
 def assymetric_gaussian_pdf(x,mu,sigma1,sigma2):
     if x<=mu:
-            return(2 / np.sqrt(2 * np.pi) / (sigma1+sigma2) * np.exp(-(x - mu) ** 2 / 2. / sigma1 ** 2))
+            return(2 / np.sqrt(2 * np.pi) / (abs(sigma1)+abs(sigma2)) * np.exp(-(x - mu) ** 2 / 2. / sigma1 ** 2))
     else:
-            return(2 / np.sqrt(2 * np.pi) / (sigma1+sigma2) * np.exp(-(x - mu) ** 2 / 2. / sigma2 ** 2))
+            return(2 / np.sqrt(2 * np.pi) / (abs(sigma1)+abs(sigma2)) * np.exp(-(x - mu) ** 2 / 2. / sigma2 ** 2))
 
 @nb.njit(**kwd)
 def assymetric_double_gaussian(x, mu, sigma1,sigma2,mu_2,sigma1_2,sigma2_2,A,B,C):
-    #assymetric_gaussian_pdf_vec=np.vectorize(assymetric_gaussian_pdf)
     return (A+B*assymetric_gaussian_pdf(x,mu,sigma1,sigma2)+C*assymetric_gaussian_pdf(x,mu_2,sigma1_2,sigma2_2))
 
 

ptiming_ana/phaseogram/pfitting.py

@@ -63,6 +63,8 @@ def est_initial_values(self,pulsar_phases):
                 #Set different initial values for different models
                 if self.model=='tgaussian':
                     regions_names=['P1','P2','P3']
+                elif self.model =='gaussian':
+                    regions_names=['P2']
                 else:
                     regions_names=['P1','P2']
 
@@ -200,9 +202,10 @@ def fit_Binned(self,pulsar_phases):
 
 
             elif self.model=='asym_dgaussian':
-
-                c = cost.LeastSquares(bin_height,np.sqrt(bin_height),bin_centres,assymetric_double_gaussian)
-                minuit = Minuit(c, *self.init)
+                assymetric_double_gaussian_vec=np.vectorize(assymetric_double_gaussian,excluded = [1,2,3,4,5,6,7,8,9,10])
+                custom_adgaussian = lambda x, mu,sigma1,sigma2,mu_2,sigma1_2,sigma2_2,B,C: assymetric_double_gaussian_vec(x, mu,sigma1,sigma2,mu_2,sigma1_2,sigma2_2,self.init[-1],B,C)
+
+                params,pcov_l=curve_fit(custom_adgaussian,bin_centres,bin_height,sigma = np.sqrt(bin_height), p0=self.init[:-1])
                 self.parnames=['mu', 'sigma1','sigma2','mu_2','sigma1_2','sigma2_2','A','B','C']
 
 
@@ -217,15 +220,18 @@ def fit_Binned(self,pulsar_phases):
                 self.parnames=['A','mu', 'sigma','mu_2','sigma_2','mu_3','sigma_3','B','C','D']
 
 
-            elif self.model=='double_lorentzian':
-                c = cost.LeastSquares(bin_height,np.sqrt(bin_height),bin_centres,double_lorentzian)
-                minuit = Minuit(c, *self.init)
+            elif self.model=='double_lorentz':
+                custom_lorentzian = lambda x, mu_1, gamma_1,mu_2,gamma_2,B,C: double_lorentz(x, mu_1, gamma_1,mu_2,gamma_2,self.init[-1],B,C)
+
+                params,pcov_l=curve_fit(custom_lorentzian,bin_centres,bin_height,sigma = np.sqrt(bin_height), p0=self.init[:-1])
+
                 self.parnames=['mu_1', 'gamma_1','mu_2','gamma_2','A','B','C']
 
 
             elif self.model=='gaussian':
-                c = cost.LeastSquares(bin_height,np.sqrt(bin_height),bin_centres,gaussian)
-                minuit = Minuit(c, *self.init)
+                custom_gaussian = lambda x,mu, sigma,B: gaussian(x, mu,sigma,self.init[-1],B)
+
+                params,pcov_l=curve_fit(custom_gaussian,bin_centres,bin_height,sigma = np.sqrt(bin_height), p0=self.init[:-1])
                 self.parnames=['mu', 'sigma','A','B']