Merge pull request #58 from ClimateImpactLab/findpolymin-1

fix error in _findpolymin

HISTORY.rst

@@ -5,12 +5,21 @@ History
 0.1.2 (Current version)
 -----------------------
 
+API changes
+~~~~~~~~~~~
+
 * :py:func:`impax.csvv.get_gammas` has been deprecated. Use :py:func:`impax.read_csvv` instead (:issue:`37`)
 * :py:meth:`~impax.csvv.Gammas._prep_gammas` has been removed, and :py:meth:`~impax.csvv.Gammas.sample` now
   takes no arguments and returns a sample by default. Seeding the random number generator is now left up to
   the user (:issue:`36`)
+
+
+Bug fixes
+~~~~~~~~~
+
 * fix py3k compatability issues (:issue:`39`)
 * fix travis status icon in README
+* add tests for :py:func:`impax.mins._minimize_polynomial`, fix major math errors causing a failure to find minima in :py:mod:`impax.mins` module, and clarify documentation (:issue:`58`)
 
 
 0.1.0 (2017-10-12)

impax/__init__.py

@@ -5,13 +5,13 @@
 
 __author__ = """Justin Simcock"""
 __email__ = 'jsimcock@rhg.com'
-__version__ = '0.1.1'
+__version__ = '0.1.2'
 
 
 _module_imports = (
     minimize_polynomial,
     construct_covars,
-    construct_weather, 
+    construct_weather,
     read_csvv,
     MultivariateNormalEstimator
 )

impax/estimate.py

@@ -1,9 +1,9 @@
 from __future__ import absolute_import
-import csv
-import xarray as xr
+
 import pandas as pd
 import numpy as np
-from scipy.stats import multivariate_normal as mn
+import csv
+from scipy.stats import multivariate_normal
 
 import warnings
 
@@ -20,7 +20,8 @@ def read_csvv(csvv_path):
     Returns
     -------
     estimator : MultivariateNormalEstimator
-        :py:class:`Gamma` object with median and VCV matrix indexed by prednames, covarnames, and outcomes
+        :py:class:`Gamma` object with median and VCV matrix indexed by
+        prednames, covarnames, and outcomes
 
     '''
 
@@ -40,38 +41,46 @@ def read_csvv(csvv_path):
                 data['prednames'] = [i.strip() for i in next(reader)]
             if row[0] == 'covarnames':
                 data['covarnames'] = [i.strip() for i in next(reader)]
-            if row[0] == 'outcome': 
-                data['outcome'] =[cv.strip() for cv in next(reader)]
+            if row[0] == 'outcome':
+                data['outcome'] = [cv.strip() for cv in next(reader)]
 
     index = pd.MultiIndex.from_tuples(
-        list(zip(data['outcome'], data['prednames'], data['covarnames'])), 
+        list(zip(data['outcome'], data['prednames'], data['covarnames'])),
         names=['outcome', 'prednames', 'covarnames'])
 
     g = MultivariateNormalEstimator(data['gamma'], data['gammavcv'], index)
 
-    return g 
+    return g
 
 
 def get_gammas(*args, **kwargs):
-    warnings.warn('get_gammas has been deprecated, and has been replaced with read_csvv', DeprecationWarning)
+    warnings.warn(
+        'get_gammas has been deprecated, and has been replaced with read_csvv',
+        DeprecationWarning)
+
     return read_csvv(*args, **kwargs)
 
 
 class MultivariateNormalEstimator(object):
     '''
-    Stores a median and residual VCV matrix for multidimensional variables with named indices
-    and provides multivariate sampling and statistical analysis functions
+    Stores a median and residual VCV matrix for multidimensional variables with
+    named indices and provides multivariate sampling and statistical analysis
+    functions
 
     Parameters
     ----------
-    coefficients: array 
-        length :math:`(m_1*m_2*\cdots*m_n)` 1-d :py:class:`numpy.ndarray` with regression coefficients
+    coefficients: array
+        length :math:`(m_1*m_2*\cdots*m_n)` 1-d :py:class:`numpy.ndarray` with
+        regression coefficients
 
     vcv: array
-        :math:`(m_1*m_2*\cdots*m_n) x (m_1*m_2*\cdots*m_n)` :py:class:`numpy.ndarray` with variance-covariance matrix for multivariate distribution
+        :math:`(m_1*m_2*\cdots*m_n) x (m_1*m_2*\cdots*m_n)`
+        :py:class:`numpy.ndarray` with variance-covariance matrix for
+        multivariate distribution
 
     index: Index
-        :py:class:`~pandas.Index` or :math:`(m_1*m_2*\cdots*m_n)` 1-d :py:class:`~pandas.MultiIndex` describing the multivariate space
+        :py:class:`~pandas.Index` or :math:`(m_1*m_2*\cdots*m_n)` 1-d
+        :py:class:`~pandas.MultiIndex` describing the multivariate space
 
     '''
 
@@ -95,20 +104,24 @@ def median(self):
     def sample(self, seed=None):
         '''
         Sample from the multivariate normal distribution
-        
+
         Takes a draw from a multivariate distribution and returns
         an :py:class:`xarray.DataArray` of parameter estimates.
 
         Returns
         ----------
         draw : DataArray
-            :py:class:`~xarray.DataArray` of parameter estimates drawn from the multivariate normal
+            :py:class:`~xarray.DataArray` of parameter estimates drawn from the
+            multivariate normal
 
         '''
         if seed is not None:
-            warnings.warn(
-                'Sampling with a seed has been deprecated. In future releases, this will be up to the user.',
+            warnings.warn((
+                'Sampling with a seed has been deprecated. In future releases,'
+                ' this will be up to the user.'),
                 DeprecationWarning)
             np.random.seed(seed)
 
-        return pd.Series(mn.rvs(self.coefficients, self.vcv), index=self.index).to_xarray()
+        return pd.Series(
+            multivariate_normal.rvs(self.coefficients, self.vcv),
+            index=self.index).to_xarray()

impax/impax.py

@@ -3,10 +3,7 @@
 import xarray as xr
 import pandas as pd
 import numpy as np
-from toolz import memoize
 from impax.mins import minimize_polynomial
-import time
-
 
 
 def construct_weather(**weather):
@@ -27,7 +24,6 @@ def construct_weather(**weather):
             variables, with variables concatenated along the
             new `prednames` dimension
 
-
     '''
     prednames = []
     weather_data = []
@@ -43,19 +39,20 @@ def construct_weather(**weather):
 
     return xr.concat(weather_data, pd.Index(prednames, name='prednames'))
 
+
 def construct_covars(add_constant=True, **covars):
     '''
     Helper function to construct the covariates dataarray
 
     Parameters
     -----------
     add_constant : bool
-        flag indicating whether a constant term should be added. The constant term will have the
-        same shape as the other covariate DataArrays
-    
+        flag indicating whether a constant term should be added. The constant
+        term will have the same shape as the other covariate DataArrays
+
     covars: dict
-        dictionary of covariate name, covariate (``str`` path or :py:class:`xarray.DataArray`)
-        pairs
+        dictionary of covariate name, covariate (``str`` path or
+        :py:class:`xarray.DataArray`) pairs
 
     Returns
     -------
@@ -66,7 +63,7 @@ def construct_covars(add_constant=True, **covars):
     '''
     covarnames = []
     covar_data = []
-    for covar, path  in covars.items():
+    for covar, path in covars.items():
         if hasattr(path, 'dims'):
             covar_data.append(path)
 
@@ -77,24 +74,24 @@ def construct_covars(add_constant=True, **covars):
         covarnames.append(covar)
 
     if add_constant:
-        ones = xr.DataArray(np.ones(shape=covar_data[0].shape),
+        ones = xr.DataArray(
+            np.ones(shape=covar_data[0].shape),
             coords=covar_data[0].coords,
             dims=covar_data[0].dims)
         covarnames.append('1')
         covar_data.append(ones)
-        
+
     return xr.concat(covar_data, pd.Index(covarnames, name='covarnames'))
-        
+
 
 class Impact(object):
     '''
     Base class for computing an impact as specified by the Climate Impact Lab
-    
+
     '''
 
     min_function = NotImplementedError
 
-
     def impact_function(self, betas, weather):
         '''
         computes the dot product of betas and annual weather by outcome group
@@ -118,65 +115,69 @@ def impact_function(self, betas, weather):
             overrides `impact_function` method in Impact base class
 
         '''
-    
+
         return (betas*weather).sum(dim='prednames')
-
 
-    def compute(self,
+    def compute(
+            self,
             weather,
             betas,
             clip_flat_curve=True,
             t_star=None):
         '''
-        Computes an impact for a unique set of gdp, climate, weather and gamma coefficient inputs.
-        For each set of these, we take the analytic minimum value between two points, 
-        save t_star to disk and compute analytical min for function m_star for a givene covariate set
-        This operation is called for every adaptation scenario specified in the run script
+        Computes an impact for a unique set of gdp, climate, weather and gamma
+        coefficient inputs. For each set of these, we take the analytic minimum
+        value between two points, save t_star to disk and compute analytical
+        min for function m_star for a given covariate set.
+
+        This operation is called for every adaptation scenario specified in the
+        run script.
 
         Parameters
         ----------
-        
+
         weather: DataArray
           weather :py:class:`~xarray.DataArray`
 
         betas: DataArray
           covarname by outcome :py:class:`~xarray.DataArray`
-         
+
         clip_flat_curve: bool
           flag indicating that flat-curve clipping should be performed
           on the result
 
         t_star: DataArray
-          :py:class:`xarray.DataArray` with minimum temperatures used for clipping
+          :py:class:`xarray.DataArray` with minimum temperatures used for
+          clipping
 
 
         Returns
         -------
-          :py:class `~xarray.Dataset` of impacts by hierid by outcome group 
+          :py:class `~xarray.Dataset` of impacts by hierid by outcome group
 
         '''
 
-        #Compute Raw Impact
+        # Compute Raw Impact
         impact = self.impact_function(betas, weather)
 
         if clip_flat_curve:
 
-            #Compute the min for flat curve adaptation
+            # Compute the min for flat curve adaptation
             impact_flatcurve = self.impact_function(betas, t_star)
 
-            #Compare values and evaluate a max
+            # Compare values and evaluate a max
             impact = xr.ufuncs.maximum((impact - impact_flatcurve), 0)
 
         impact = self.postprocess_daily(impact)
 
-        #Sum to annual
+        # Sum to annual
         impact = impact.sum(dim='time')
 
-        impact_annual = self.postprocess_annual(impact) 
+        impact_annual = self.postprocess_annual(impact)
 
         return impact_annual
 
-    def get_t_star(self,betas, bounds, t_star_path=None):
+    def get_t_star(self, betas, bounds, t_star_path=None):
         '''
         Read precomputed t_star
 
@@ -190,10 +191,10 @@ def get_t_star(self,betas, bounds, t_star_path=None):
             values between which to evaluate function
 
         path: str
-          place to load t-star from 
+          place to load t-star from
 
         '''
-        
+
         try:
             with xr.open_dataarray(t_star_path) as t_star:
                 return t_star.load()
@@ -204,17 +205,17 @@ def get_t_star(self,betas, bounds, t_star_path=None):
         except (IOError, ValueError):
             try:
                 os.remove(t_star_path)
-            except:
+            except (IOError, OSError):
                 pass
 
-        #Compute t_star according to min function
+        # Compute t_star according to min function
         t_star = self.compute_t_star(betas, bounds=bounds)
 
-        #write to disk
+        # write to disk
         if t_star_path is not None:
             if not os.path.isdir(os.path.dirname(t_star_path)):
                 os.makedirs(os.path.dirname(t_star_path))
-            
+
             t_star.to_netcdf(t_star_path)
 
         return t_star
@@ -237,8 +238,9 @@ class PolynomialImpact(Impact):
     @staticmethod
     def min_function(*args, **kwargs):
         '''
-        helper function to call minimization function for given mortality polynomial spec
-        mortality_polynomial implements findpolymin through `np.apply_along_axis`
+        helper function to call minimization function for given mortality
+        polynomial spec mortality_polynomial implements findpolymin through
+        `np.apply_along_axis`
 
         Parameters
         ----------
@@ -260,4 +262,3 @@ def min_function(*args, **kwargs):
         '''
 
         return minimize_polynomial(*args, **kwargs)
-