Add 6th calibration model that minimizes difference with national value

salamanca/models/calibrate_ineq.py

@@ -65,9 +65,18 @@ def threshold_hi_rule(m, idx, f=1.0, b=1.1, relative=True):
     return lhs <= b * rhs
 
 
-# def theil_sum_rule(m, b=0.0):
-#     t_w = sum(m.t[idx] * m.data['g'][idx] for idx in m.idxs) / m.data['G']
-#     return abs(t_w - m.data['T_w']) <= b
+def theil_sum_lo_rule(m, b=0.05):
+    lhs = sum(m.t[idx] * m.data['g'][idx] for idx in m.idxs)
+    rhs = m.data['T_w'] * m.data['G']
+    return lhs >= (1 - b) * rhs
+
+
+def theil_sum_hi_rule(m, b=0.05):
+    lhs = sum(m.t[idx] * m.data['g'][idx] for idx in m.idxs)
+    rhs = m.data['T_w'] * m.data['G']
+    return lhs <= (1 + b) * rhs
+
+
 def theil_sum_rule(m):
     return sum(m.t[idx] * m.data['g'][idx] for idx in m.idxs) == \
         m.data['T_w'] * m.data['G']
@@ -115,6 +124,17 @@ def threshold_obj(m, factors=[1.0], weights=[1.0], relative=True):
         for (f, w), idx in itertools.product(zip(factors, weights), m.idxs))
 
 
+def national_threshold_obj(m, f=1.0, relative=True):
+    I = m.data['I']
+    x = f * I if relative else f
+    lhs = sum(
+        below_threshold(x, m.data['i'][idx], m.t[idx])
+        for idx in m.idxs
+    )
+    rhs = below_threshold(x, I, m.data['T'])
+    return (lhs - rhs) ** 2
+
+
 class Model(object):
     """Base class for Inequality Calibration Models"""
 
@@ -165,6 +185,7 @@ def _setup_model_data(self, natdata, subdata):
             'I': ndf[i],
             'G': ndf[n] * ndf[i],
             'T_w': T_w,
+            'T': T,
         }
 
     def _check_model_data(self):
@@ -374,3 +395,36 @@ def construct(self):
         # Objective
         m.obj = mo.Objective(rule=l2_norm_obj, sense=mo.minimize)
         return self
+
+
+class Model6(Model):
+    """
+    Minimize low threshold L2-norm under position, std, and approx theil sum.
+    """
+
+    def construct(self):
+        self.model = m = mo.ConcreteModel()
+        # Model Data
+        m.data = self.model_data
+        # Sets
+        m.idxs = mo.Set(initialize=m.data['idxs'])
+        # Variables
+        m.t = mo.Var(m.idxs, within=mo.PositiveReals,
+                     bounds=(1e-5, ineq.MAX_THEIL))
+        # Constraints
+        m.position = mo.Constraint(m.idxs, rule=position_rule,
+                                   doc='ordering between provinces must be maintained')
+        m.theil_sum_hi = mo.Constraint(rule=theil_sum_hi_rule,
+                                       doc='')
+        m.theil_sum_lo = mo.Constraint(rule=theil_sum_lo_rule,
+                                       doc='')
+        m.spacing = mo.Constraint(m.idxs, rule=spacing_rule,
+                                  doc='')
+        m.std_lo = mo.Constraint(m.idxs, rule=std_diff_lo_rule,
+                                 doc='')
+        m.std_hi = mo.Constraint(m.idxs, rule=std_diff_hi_rule,
+                                 doc='')
+        # Objective
+        rule = lambda m: national_threshold_obj(m, f=0.25)
+        m.obj = mo.Objective(rule=rule, sense=mo.minimize)
+        return self

tests/test_calibrate_ineq.py

@@ -4,7 +4,8 @@
 import pandas as pd
 
 from salamanca import ineq
-from salamanca.models.calibrate_ineq import Model, Model1, Model2, Model3, Model4, Model4b
+from salamanca.models.calibrate_ineq import Model, Model1, Model2, Model3, \
+    Model4, Model4b, Model5, Model6
 
 from utils import assert_almost_equal, assert_array_almost_equal
 
@@ -182,3 +183,27 @@ def test_Model4b_result():
     obs = model.result()['gini'].values
     exp = [0.48849954, 0.0277764]
     assert_array_almost_equal(obs, exp, eps=1e-3)
+
+
+def test_Model5_result():
+    natdata, subdata = data()
+    model = Model5(natdata, subdata)
+    model.construct()
+    model.solve()
+
+    # ginis in original order
+    obs = model.result()['gini'].values
+    exp = [0.45663393, 0.19798727]
+    assert_array_almost_equal(obs, exp, eps=1e-3)
+
+
+def test_Model6_result():
+    natdata, subdata = data()
+    model = Model6(natdata, subdata)
+    model.construct()
+    model.solve()
+
+    # ginis in original order
+    obs = model.result()['gini'].values
+    exp = [0.4624727, 0.21394398]
+    assert_array_almost_equal(obs, exp, eps=1e-3)