Changeset 2666 for trunk/python

Timestamp:

10/15/12 15:52:38 (12 years ago)

Author:

Malte Marquarding

Message:

Ticket #173: added setting of constraints on the fitter.

File:

• trunk/python/asapfitter.py (modified) (15 diffs)

trunk/python/asapfitter.py

@@ -3 +3 @@
 from asap.logging import asaplog, asaplog_post_dec
 from asap.utils import _n_bools, mask_and
-
+from numpy import ndarray
 
 class fitter:
@@ -26 +26 @@
         self._selection = None
         self.uselinear = False
+        self._constraints = []
 
     def set_data(self, xdat, ydat, mask=None):
@@ -73 +74 @@
         Set the function to be fit.
         Parameters:
-            poly:     use a polynomial of the order given with nonlinear least squares fit
-            lpoly:    use polynomial of the order given with linear least squares fit
+            poly:     use a polynomial of the order given with nonlinear 
+                      least squares fit
+            lpoly:    use polynomial of the order given with linear least 
+                      squares fit
             gauss:    fit the number of gaussian specified
             lorentz:  fit the number of lorentzian specified
             sinusoid: fit the number of sinusoid specified
         Example:
-            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
-            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
+            fitter.set_function(poly=3)  # will fit a 3rd order polynomial 
+                                         # via nonlinear method
+            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial 
+                                          # via linear method
             fitter.set_function(gauss=2) # will fit two gaussians
             fitter.set_function(lorentz=2) # will fit two lorentzians
@@ -117 +122 @@
             self.components = [ 3 for i in range(n) ]
             self.uselinear = False
+        elif kwargs.has_key('expression'):
+            self.uselinear = False
+            raise RuntimeError("Not yet implemented")
         else:
             msg = "Invalid function type."
@@ -122 +130 @@
 
         self.fitter.setexpression(self.fitfunc,n)
+        self._constraints = []
         self.fitted = False
         return
@@ -147 +156 @@
             raise RuntimeError(msg)
 
-        else:
-            if self.data is not None:
-                self.x = self.data._getabcissa(row)
-                self.y = self.data._getspectrum(row)
-                #self.mask = mask_and(self.mask, self.data._getmask(row))
-                if len(self.x) == len(self.mask):
-                    self.mask = mask_and(self.mask, self.data._getmask(row))
-                else:
-                    asaplog.push('lengths of data and mask are not the same. preset mask will be ignored')
-                    asaplog.post('WARN','asapfit.fit')
-                    self.mask=self.data._getmask(row)
-                asaplog.push("Fitting:")
-                i = row
-                out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (self.data.getscan(i),
-                                                                      self.data.getbeam(i),
-                                                                      self.data.getif(i),
-                                                                      self.data.getpol(i),
-                                                                      self.data.getcycle(i))
-                asaplog.push(out,False)
+        if self.data is not None:
+            self.x = self.data._getabcissa(row)
+            self.y = self.data._getspectrum(row)
+            #self.mask = mask_and(self.mask, self.data._getmask(row))
+            if len(self.x) == len(self.mask):
+                self.mask = mask_and(self.mask, self.data._getmask(row))
+            else:
+                asaplog.push('lengths of data and mask are not the same. '
+                             'preset mask will be ignored')
+                asaplog.post('WARN','asapfit.fit')
+                self.mask=self.data._getmask(row)
+            asaplog.push("Fitting:")
+            i = row
+            out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (
+                self.data.getscan(i),
+                self.data.getbeam(i),
+                self.data.getif(i),
+                self.data.getpol(i),
+                self.data.getcycle(i))
+            
+            asaplog.push(out, False)
+
         self.fitter.setdata(self.x, self.y, self.mask)
         if self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
@@ -174 +186 @@
         if len(fxdpar) and fxdpar.count(0) == 0:
              raise RuntimeError,"No point fitting, if all parameters are fixed."
+        if self._constraints:
+            for c in self._constraints:
+                self.fitter.addconstraint(c[0]+[c[-1]])
         if self.uselinear:
             converged = self.fitter.lfit()
@@ -234 +249 @@
             msg = "Please specify a fitting function first."
             raise RuntimeError(msg)
-        if (self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid") and component is not None:
+        if (self.fitfunc == "gauss" or self.fitfunc == "lorentz" 
+            or self.fitfunc == "sinusoid") and component is not None:
             if not self.fitted and sum(self.fitter.getparameters()) == 0:
                 pars = _n_bools(len(self.components)*3, False)
@@ -338 +354 @@
             raise ValueError(msg)
 
+
+    def add_constraint(self, xpar, y):
+        """Add parameter constraints to the fit. This is done by setting up
+        linear equations for the related parameters.
+
+        For example a two component gaussian fit where the amplitudes are
+        constraint by amp1 = 2*amp2
+        needs a constraint   
+
+            add_constraint([1, 0, 0, -2, 0, 0, 0], 0)
+
+        a velocity difference of v2-v1=17
+
+            add_constraint([0.,-1.,0.,0.,1.,0.], 17.)
+
+        """
+        self._constraints.append((xpar, y))
+        
+
     def get_area(self, component=None):
         """
@@ -381 +416 @@
         cerrs = errs
         if component is not None:
-            if self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid":
+            if self.fitfunc == "gauss" or self.fitfunc == "lorentz" \
+                    or self.fitfunc == "sinusoid":
                 i = 3*component
                 if i < len(errs):
@@ -463 +499 @@
                 out += "%s%s = %3.3f %s, " % (pnam[1], fix1, pars[i+1], aunit)
                 out += "%s%s = %3.3f %s\n" % (pnam[2], fix2, pars[i+2], aunit)
-                if len(area): out += "      area = %3.3f %s %s\n" % (area[i], ounit, aunit)
+                if len(area): out += "      area = %3.3f %s %s\n" % (area[i], 
+                                                                     ounit, 
+                                                                     aunit)
                 c+=1
                 i+=3
@@ -571 +609 @@
             ylab = self.data._get_ordinate_label()
 
-        colours = ["#777777","#dddddd","red","orange","purple","green","magenta", "cyan"]
+        colours = ["#777777","#dddddd","red","orange","purple","green",
+                   "magenta", "cyan"]
         nomask=True
         for i in range(len(m)):
@@ -604 +643 @@
             if isinstance(components,int): cs = [components]
             if plotparms:
-                self._p.text(0.15,0.15,str(self.get_parameters()['formatted']),size=8)
+                self._p.text(0.15,0.15,
+                             str(self.get_parameters()['formatted']),size=8)
             n = len(self.components)
             self._p.palette(3)
@@ -656 +696 @@
         asaplog.push("Fitting:")
         for r in rows:
-            out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (scan.getscan(r),
-                                                                   scan.getbeam(r),
-                                                                   scan.getif(r),
-                                                                   scan.getpol(r),
-                                                                   scan.getcycle(r))
+            out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (
+                scan.getscan(r),
+                scan.getbeam(r),
+                scan.getif(r),
+                scan.getpol(r),
+                scan.getcycle(r)
+                )
             asaplog.push(out, False)
             self.x = scan._getabcissa(r)
@@ -668 +710 @@
                 self.mask = mask_and(self.mask, self.data._getmask(row))
             else:
-                asaplog.push('lengths of data and mask are not the same. preset mask will be ignored')
+                asaplog.push('lengths of data and mask are not the same. '
+                             'preset mask will be ignored')
                 asaplog.post('WARN','asapfit.fit')
                 self.mask=self.data._getmask(row)