From 8da0eb122e93bf9defb6169a27baecc4000dd1a7 Mon Sep 17 00:00:00 2001 From: Matthew Gidden Date: Wed, 6 Sep 2017 14:49:12 +0200 Subject: [PATCH] Added second projection model which minimized population threshold differences and moves theil difference into constraint --- salamanca/models/project.py | 135 ++++++++++++++++++++++++++++++++++-- tests/test_project.py | 35 +++++++++- 2 files changed, 162 insertions(+), 8 deletions(-) diff --git a/salamanca/models/project.py b/salamanca/models/project.py index 45b290d..96b49d3 100644 --- a/salamanca/models/project.py +++ b/salamanca/models/project.py @@ -24,22 +24,34 @@ def gdp_sum_rule(m): return sum(m.data['n_frac'][idx] * m.i[idx] for idx in m.idxs) == m.data['I'] -def threshold_lo_rule(m, f=1.0, b=0.95, relative=True): +def theil_lo_rule(m, b=0.05): + lhs = model_T_b(m) + model_T_w(m, income_from_data=False) + rhs = m.data['T'] + return lhs >= (1 - b) * rhs + + +def theil_hi_rule(m, b=0.05): + lhs = model_T_b(m) + model_T_w(m, income_from_data=False) + rhs = m.data['T'] + return lhs <= (1 + b) * rhs + + +def threshold_lo_rule(m, f=1.0, b=0.05, relative=True): i = m.data['I'] x = f * i if relative else f rhs = below_threshold(x, i, m.data['T']) lhs = sum(m.data['n_frac'][idx] * below_threshold(x, m.i[idx], m.t[idx]) for idx in m.idxs) - return lhs >= b * rhs + return lhs >= (1 - b) * rhs -def threshold_hi_rule(m, f=1.0, b=1.05, relative=True): +def threshold_hi_rule(m, f=1.0, b=0.05, relative=True): i = m.data['I'] x = f * i if relative else f rhs = below_threshold(x, i, m.data['T']) lhs = sum(m.data['n_frac'][idx] * below_threshold(x, m.i[idx], m.t[idx]) for idx in m.idxs) - return lhs <= b * rhs + return lhs <= (1 + b) * rhs def theil_diff_hi_rule(m, idx, b): @@ -134,6 +146,15 @@ def theil_total_sum_obj(m): return (m.data['T'] - T_b - T_w) ** 2 +def population_below_obj(m, f=1.0, relative=True): + i = m.data['I'] + x = f * i if relative else f + p_nat = below_threshold(x, i, m.data['T']) + p_sum = sum(m.data['n_frac'][idx] * below_threshold(x, m.i[idx], m.t[idx]) + for idx in m.idxs) + return (p_nat - p_sum) ** 2 + + class Model(object): """Base class for Projection Models""" @@ -256,10 +277,12 @@ class Model1(Model): - GDP sum - constrained CDF - with optionally (assuming 10 year time steps): + with optionally (assuming equally spaced time steps): - - maximum theil diffusion of 1% per year - - maximum income share diffusion of 2% per year + - maximum theil diffusion % per year + - maximum income share diffusion % per year + - incomes track with national income + - theils track with national theil """ def construct(self, diffusion={}, direction={}): @@ -337,6 +360,104 @@ def construct(self, diffusion={}, direction={}): return self +class Model2(Model): + """ + Minimize L2 norm of population under income threshold subject to: + + - GDP sum + - constrained Theil + + + with optionally (assuming equally spaced time steps): + + - maximum theil diffusion % per year + - maximum income share diffusion % per year + - incomes track with national income + - theils track with national theil + """ + + def construct(self, diffusion={}, direction={}, + theil_within=0.05, threshold=1.0): + self.model = m = mo.ConcreteModel() + # Model Data + m.data = self.model_data + # Sets + m.idxs = mo.Set(initialize=m.data['idxs']) + # Variables + m.i = mo.Var(m.idxs, within=mo.PositiveReals, + bounds=(m.data['i_min'], m.data['i_max'])) + m.t = mo.Var(m.idxs, within=mo.PositiveReals, + bounds=(m.data['t_min'], m.data['t_max'])) + # Constraints + m.gdp_sum = mo.Constraint(rule=gdp_sum_rule, doc='gdp sum = gdp') + m.theil_lo = mo.Constraint( + rule=lambda m: theil_lo_rule(m, b=theil_within), + doc='Theil within 5%' + ) + m.theil_hi = mo.Constraint( + rule=lambda m: theil_hi_rule(m, b=theil_within), + doc='Theil within 5%' + ) + # optional constraints + b = diffusion.pop('income', False) + if b: + b = 0.2 if b is True else b + m.i_hi = mo.Constraint( + m.idxs, + rule=lambda m, idx: income_diff_hi_rule(m, idx, b), + doc='income within 20% from past', + ) + m.i_lo = mo.Constraint( + m.idxs, + rule=lambda m, idx: income_diff_lo_rule(m, idx, b), + doc='income within 20% from past', + ) + b = diffusion.pop('share', False) + if b: + b = 0.2 if b is True else b + m.s_hi = mo.Constraint( + m.idxs, + rule=lambda m, idx: share_diff_hi_rule(m, idx, b), + doc='income share within 20% from past', + ) + m.s_lo = mo.Constraint( + m.idxs, + rule=lambda m, idx: share_diff_lo_rule(m, idx, b), + doc='income share within 20% from past', + ) + b = diffusion.pop('theil', False) + if b: + b = 0.1 if b is True else b + m.t_hi = mo.Constraint( + m.idxs, + rule=lambda m, idx: theil_diff_hi_rule(m, idx, b), + doc='income share within 20% from past', + ) + m.t_lo = mo.Constraint( + m.idxs, + rule=lambda m, idx: theil_diff_lo_rule(m, idx, b), + doc='income share within 20% from past', + ) + if direction.pop('income', False): + m.i_dir = mo.Constraint( + m.idxs, + rule=income_direction_rule, + doc='income must track with national values', + ) + if direction.pop('theil', False): + m.t_dir = mo.Constraint( + m.idxs, + rule=theil_direction_rule, + doc='theil must track with national values', + ) + # Objective + m.obj = mo.Objective( + rule=lambda m: population_below_obj(m, f=threshold), + sense=mo.minimize + ) + return self + + class Runner(object): """Simple class that runs all projection values for a given full projection instance diff --git a/tests/test_project.py b/tests/test_project.py index 4bf8272..218676d 100644 --- a/tests/test_project.py +++ b/tests/test_project.py @@ -4,7 +4,7 @@ import pandas as pd from salamanca import ineq -from salamanca.models.project import Model, Model1 +from salamanca.models.project import Model, Model1, Model2 from salamanca.models.project import Runner from pandas.testing import assert_frame_equal @@ -130,6 +130,39 @@ def test_Model1_diffusion_result(): assert_frame_equal(obs, exp) +def test_Model2_result(): + natdata, subdata = data() + model = Model2(natdata, subdata) + model.construct() + model.solve() + df = model.result() + + # this is a regression test, results tested 09-06-17 + obs = df[['gini', 'i']] + exp = pd.DataFrame({ + 'i': [7.80113, 6.549360], + 'gini': [0.348809, 0.347297] + }, index=pd.Index(['foo', 'bar'], name='name')) + assert_frame_equal(obs, exp) + + +def test_Model2_diffusion_result(): + natdata, subdata = data() + model = Model2(natdata, subdata) + diffusion = {'theil': True} + model.construct(diffusion=diffusion, theil_within=0.1) + model.solve() + df = model.result() + + # this is a regression test, results tested 09-06-17 + obs = df[['gini', 'i']] + exp = pd.DataFrame({ + 'i': [6.39826598759, 7.33847538198], + 'gini': [0.47774761513, 0.285297024826] + }, index=pd.Index(['foo', 'bar'], name='name')) + assert_frame_equal(obs, exp) + + def test_runner_horizon(): # test on actual data natdata, subdata = data()