Merge pull request #58 from esheldon/wpixels

Wpixels

Author: esheldon

ngmix/__init__.py

@@ -1,9 +1,15 @@
 __version__="v1.2"
 
 from . import gmix
-from .gmix import GMix
-from .gmix import GMixModel
-from .gmix import GMixCoellip
+from .gmix import (
+    GMix,
+    GMixModel,
+    GMixBDF,
+    GMixCoellip,
+
+    GMixList,
+    MultiBandGMixList,
+)
 
 from . import gmix_ndim
 from .gmix_ndim import GMixND
@@ -25,7 +31,7 @@
 from .gexceptions import GMixRangeError, GMixFatalError, GMixMaxIterEM
 
 from . import fitting
-from .fitting import print_pars
+from .fitting import print_pars, format_pars
 from . import simplex
 
 from . import galsimfit

ngmix/bootstrap.py

@@ -2693,10 +2693,57 @@ def _get_lm_fitter_class(self):
         from .fitting import LMComposite
         return LMComposite
 
+class BDRunner(MaxRunner):
+    """
+    wrapper to generate guesses and run the BD fitter a few times
+    """
+    def __init__(self,
+                 obs,
+                 max_pars,
+                 guesser,
+                 prior=None):
+        self.obs=obs
+
+        self.max_pars=max_pars
+
+        self.method=max_pars.get('method','lm')
+        if self.method == 'lm':
+            self.send_pars=max_pars['lm_pars']
+        else:
+            self.send_pars=max_pars
+
+        self.prior=prior
+
+        self.guesser=guesser
+
+    def _go_lm(self, ntry=1):
+        fitclass=self._get_lm_fitter_class()
+
+        for i in xrange(ntry):
+            guess=self.guesser()
+            fitter=fitclass(
+                self.obs,
+                lm_pars=self.send_pars,
+                prior=self.prior,
+            )
+
+            fitter.go(guess)
+
+            res=fitter.get_result()
+            if res['flags']==0:
+                break
+
+        res['ntry'] = i+1
+        self.fitter=fitter
+
+    def _get_lm_fitter_class(self):
+        from .fitting import LMBD
+        return LMBD
+
 
 class BDFRunner(MaxRunner):
     """
-    wrapper to generate guesses and run the psf fitter a few times
+    wrapper to generate guesses and run the BDF fitter a few times
     """
     def __init__(self,
                  obs,

ngmix/fitting.py

@@ -180,8 +180,7 @@ def _set_totpix(self):
         totpix=0
         for obs_list in self.obs:
             for obs in obs_list:
-                shape=obs.image.shape
-                totpix += shape[0]*shape[1]
+                totpix += obs.pixels.size
 
         self.totpix=totpix
 
@@ -849,8 +848,7 @@ def _set_totpix(self):
 
         totpix=0
         for obs in self.obs:
-            shape=obs.image.shape
-            totpix += shape[0]*shape[1]
+            totpix += obs.pixels.size
 
         self.totpix=totpix
 
@@ -1674,13 +1672,17 @@ def go(self, guess):
 
         self._make_lists()
 
+        bounds = self._get_bounds()
+
         result = run_leastsq(
             self._calc_fdiff,
             guess,
             self.n_prior_pars,
+            bounds=bounds,
             **self.lm_pars
         )
 
+
         result['model'] = self.model_name
         if result['flags']==0:
             result['g'] = result['pars'][2:2+2].copy()
@@ -1692,6 +1694,15 @@ def go(self, guess):
 
         self._result=result
 
+    def _get_bounds(self):
+        """
+        get bounds on parameters
+        """
+        bounds=None
+        if self.prior is not None:
+            if hasattr(self.prior,'bounds'):
+                bounds=self.prior.bounds
+        return bounds
 
     run_max=go
     run_lm=go
@@ -1728,7 +1739,7 @@ def _setup_data(self, guess):
             self._init_gmix_all(guess)
         except ZeroDivisionError:
             raise GMixRangeError("got zero division")
-
+        
     def get_band_pars(self, pars_in, band):
         """
         Get linear pars for the specified band
@@ -1896,6 +1907,77 @@ def _make_model(self, band_pars):
         gm0=gmix.GMixCM(self.fracdev, self.TdByTe, band_pars)
         return gm0
 
+class LMBD(LMSimple):
+    """
+    exp+dev model with Td/Te and fracdev free
+    """
+    def __init__(self, obs, **keys):
+        super(LMBD,self).__init__(obs, 'bd', **keys)
+
+        self._band_pars=zeros(8)
+
+    def _set_n_prior_pars(self):
+        # center1 + center2 + shape + T + TdByTe + fracdev + fluxes
+        if self.prior is None:
+            self.n_prior_pars=0
+        else:
+            self.n_prior_pars=1 + 1 + 1 + 1 + 1 + 1 + self.nband
+
+    def get_gmix(self, band=0):
+        """
+        Get a gaussian mixture at the "best" parameter set, which
+        definition depends on the sub-class
+        """
+        res=self.get_result()
+        pars=self.get_band_pars(res['pars'], band)
+        return self._make_model(pars)
+
+    def get_band_pars(self, pars_in, band):
+        """
+        Get linear pars for the specified band
+        """
+
+        pars=self._band_pars
+
+        pars[0:7] = pars_in[0:7]
+        pars[7] = pars_in[7+band]
+
+        return pars
+
+
+    def _make_model(self, band_pars):
+        """
+        incoming parameters are
+            [c1,c2,g1,g2,T,TdByTe,fracdev,F]
+        """
+        return gmix.GMixModel(band_pars,'bd')
+
+    def _set_flux(self, res):
+        if self.nband==1:
+            res['flux'] = res['pars'][7]
+            res['flux_err'] = sqrt(res['pars_cov'][7,7])
+        else:
+            res['flux'] = res['pars'][7:]
+            res['flux_cov'] = res['pars_cov'][7:, 7:]
+            res['flux_err'] = sqrt(diag(res['flux_cov']))
+
+    '''
+    def _setup_data(self, guess):
+        super(LMBD,self)._setup_data(guess)
+
+        import images
+        gm = self._gmix_all[0][0]
+        print('gm:',gm)
+        im = gm.make_image(
+            self.obs[0][0].image.shape,
+            jacobian=self.obs[0][0].jacobian,
+        )
+        print(im.std())
+        images.multiview(im)
+        stop
+    '''
+
+
 class LMBDF(LMSimple):
     """
     exp+dev model with Tdev/Texp=1 but fracdev free
@@ -1979,7 +2061,8 @@ def run_leastsq(func, guess, n_prior_pars, **keys):
     xtol:
         Relative error desired in solution. 1.0e-6
     """
-    from scipy.optimize import leastsq
+    #from scipy.optimize import leastsq
+    from .leastsqbound import leastsqbound as leastsq
 
     npars=guess.size
     k_space=keys.pop('k_space',False)
@@ -3363,6 +3446,15 @@ def print_pars(pars, stream=stdout, fmt='%8.3g',front=None):
     if pars is None:
         stream.write('%s\n' % None)
     else:
-        fmt = ' '.join( [fmt+' ']*len(pars) )
-        stream.write(fmt % tuple(pars))
+        #fmt = ' '.join( [fmt+' ']*len(pars) )
+        #stream.write(fmt % tuple(pars))
+        s = format_pars(pars, fmt=fmt)
+        stream.write(s)
         stream.write('\n')
+
+def format_pars(pars, fmt='%8.3g'):
+    """
+    make a nice string of the pars with no line breaks
+    """
+    fmt = ' '.join( [fmt+' ']*len(pars) )
+    return fmt % tuple(pars)

ngmix/gmix.py

@@ -35,13 +35,16 @@
 # this is for backward compatibility
 from .gmix_ndim import GMixND
 
-def make_gmix_model(pars, model):
+def make_gmix_model(pars, model, **kw):
     """
     get a gaussian mixture model for the given model
     """
 
+    model = get_model_num(model)
     if model==GMIX_COELLIP:
         return GMixCoellip(pars)
+    elif model==GMIX_BDF:
+        return GMixBDF(pars=pars, **kw)
     elif model==GMIX_FULL:
         return GMix(pars=pars)
     else:
@@ -831,6 +834,10 @@ class GMixModel(GMix):
     """
     def __init__(self, pars, model):
 
+        assert model != 'bdf','use GMixBDF for bdf model'
+        self._do_init(pars, model)
+
+    def _do_init(self, pars, model):
         self._model      = _gmix_model_dict[model]
         self._model_name = _gmix_string_dict[self._model]
 
@@ -928,17 +935,34 @@ def __repr__(self):
 class GMixBDF(GMixModel):
     """
     Gaussian mixture representing a bulge+disk with
-    fixed size ratio Td/Te=1
+    fixed size ratio
+
+    Parameters
+    ----------
+    pars: sequence
+        [c1,c2,g1,g2,T,flux]
+    TdByTe: number, optional
+        Optionally set TdByTe.  Defaults to 1.0
+
+        Notes: a value of 1.0 is not the most common value for real world
+        galaxies.  It is more typically ~0.3 when fitting to COSMOS galaxies.
+        But 1.0 provides much more stable fits generally and does not reduce
+        accuracy much.
     """
-    def __init__(self, pars):
-        super(GMixBDF,self).__init__(pars,'bdf')
+    def __init__(self, pars=None, TdByTe=1.0):
+        assert pars is not None,'send pars='
+        assert TdByTe is not None,'send TdByTe='
+
+        self._TdByTe = float(TdByTe)
+        super(GMixBDF,self)._do_init(pars,'bdf')
 
     def copy(self):
         """
         Get a new GMix with the same parameters
         """
         return GMixBDF(
-            self._pars,
+            pars=self._pars.copy(),
+            TdByTe=self._TdByTe,
         )
 
     def _fill(self, pars):
@@ -958,12 +982,14 @@ def _fill(self, pars):
         self._fill_func(
             gm,
             self._pars,
+            self._TdByTe,
         )
 
 
     def __repr__(self):
         rep=super(GMixBDF,self).__repr__()
         rep = [
+            'TdByTe:  %g' % self._TdByTe,
             rep,
         ]
         return '\n'.join(rep)
@@ -1060,6 +1086,8 @@ def copy(self):
 GMIX_SERSIC=8
 GMIX_CM=9
 
+GMIX_BD=10
+
 _gmix_model_dict={
     'full':       GMIX_FULL,
     GMIX_FULL:    GMIX_FULL,
@@ -1073,6 +1101,10 @@ def copy(self):
     GMIX_DEV:     GMIX_DEV,
     'bdc':        GMIX_BDC,
     GMIX_BDC:     GMIX_BDC,
+
+    'bd':        GMIX_BD,
+    GMIX_BD:     GMIX_BD,
+
     'bdf':        GMIX_BDF,
     GMIX_BDF:     GMIX_BDF,
 
@@ -1099,6 +1131,10 @@ def copy(self):
     'dev':'dev',
     GMIX_BDC:'bdc',
     'bdc':'bdc',
+
+    GMIX_BD:'bd',
+    'bd':'bd',
+
     GMIX_BDF:'bdf',
     'bdf':'bdf',
 
@@ -1120,10 +1156,12 @@ def copy(self):
     GMIX_DEV:6,
 
     GMIX_CM:6,
+
+    GMIX_BD:8,
+
     GMIX_BDF:7,
 
     GMIX_BDC:8,
-    GMIX_BDF:7,
     GMIX_SERSIC:7,
 }
 
@@ -1138,10 +1176,10 @@ def copy(self):
     'dev':10,
 
     GMIX_CM:16,
+    GMIX_BD:16,
     GMIX_BDF:16,
 
     GMIX_BDC:16,
-    GMIX_BDF:16,
     GMIX_SERSIC:4,
 
     'em1':1,