Merge pull request #55 from mikhailsirenko/develop

Develop

config/SaintLucia.yaml

@@ -6,7 +6,7 @@ constants:
     pov_bias_rnd_distr: uniform
     pov_bias_rnd_high: 1.5
     pov_bias_rnd_low: 0.5
-  consumption_utility: 1.5
+  consump_util: 1.5
   country: Saint Lucia
   ident_affected_params:
     delta_pct: 0.05
@@ -86,7 +86,7 @@ levers:
 
 # ------------------------------- Uncertainties ------------------------------ #
 uncertainties:
-  consumption_utility:
+  consump_util:
   - 1.0
   - 1.5
   discount_rate:

unbreakable/analysis/analyzer.py

@@ -5,6 +5,7 @@
 import ast
 import geopandas as gpd
 from tqdm import tqdm
+import yaml
 
 
 def prepare_outcomes(results: tuple, add_policies: bool, add_uncertainties: bool) -> pd.DataFrame:
@@ -18,38 +19,16 @@ def prepare_outcomes(results: tuple, add_policies: bool, add_uncertainties: bool
     Returns:
         pd.DataFrame: Outcomes.
     '''
-    # * Note that we specify all outcomes in `get_outcomes` function in `write.py`
-    # * Here we just read them in the same sequence that they are written
-    outcome_names = [
-        'total_population',
-        'total_asset_loss',
-        'total_consumption_loss',
-        'tot_exposed_asset',
-        'tot_asset_surv',
-        'expected_loss_fraction',
-        'n_affected_people',
-        'annual_average_consumption',
-        'poverty_line_adjusted',
-        'district_pml',
-        'n_poor_initial',
-        'n_poor_affected',
-        'n_new_poor',
-        'initial_poverty_gap',
-        'new_poverty_gap_initial',
-        'new_poverty_gap_all',
-        'annual_average_consumption_loss',
-        'annual_average_consumption_loss_pct',
-        'r',
-        'mean_recovery_rate',
-        'weighted_vuln_quint',
-        'weighted_vuln_dec',
-        'years_in_poverty',
-    ]
-
-    uncertainty_names = ['consumption_utility',
-                         'discount_rate',
-                         'income_and_expenditure_growth',
-                         'poverty_bias']
+    # Read outcome names from a yaml file
+    with open("../unbreakable/analysis/outcome_names.yaml", "r") as f:
+        outcome_names = yaml.safe_load(f)
+
+    # TODO: Read uncertainty names from results
+    # uncertainty_names = ['consump_util',
+    #                      'discount_rate',
+    #                      'income_and_expenditure_growth',
+    #                      'poverty_bias']
+    uncertainty_names = []
 
     experiments, _ = results
     experiments['random_seed'] = experiments['random_seed'].astype(int)
@@ -172,16 +151,6 @@ def prepare_outcomes(results: tuple, add_policies: bool, add_uncertainties: bool
             i += 1  # increase row index of the outcomes dataframe
     outcomes = pd.DataFrame(outcomes, columns=columns)
 
-    # Convert numeric columns to numeric
-    if add_policies:
-        numeric_columns = outcomes.columns[5:-4].tolist()
-        outcomes[numeric_columns] = outcomes[numeric_columns].apply(
-            pd.to_numeric)
-    else:
-        numeric_columns = outcomes.columns[4:-4].tolist()
-        outcomes[numeric_columns] = outcomes[numeric_columns].apply(
-            pd.to_numeric)
-
     # Rename a district
     outcomes['district'].replace(
         {'AnseLaRayeCanaries': 'Anse-La-Raye & Canaries'}, inplace=True)
@@ -196,13 +165,6 @@ def prepare_outcomes(results: tuple, add_policies: bool, add_uncertainties: bool
     outcomes = outcomes.assign(n_new_poor_increase_pp=outcomes['n_new_poor'].div(
         outcomes['total_population']).multiply(100))
 
-    outcomes['pct_poor_before'] = outcomes['n_poor_initial'].div(
-        outcomes['total_population'])
-    outcomes['pct_poor_after'] = outcomes['n_new_poor'].add(
-        outcomes['n_poor_initial']).div(outcomes['total_population'])
-    outcomes['pct_poor_increase'] = outcomes['pct_poor_after'].sub(
-        outcomes['pct_poor_before'])
-
     # Move years_in_poverty column to the end of the data frame
     outcomes = outcomes[[c for c in outcomes if c not in [
         'years_in_poverty']] + ['years_in_poverty']]
@@ -241,6 +203,7 @@ def get_spatial_outcomes(outcomes: pd.DataFrame, outcomes_of_interest: list = []
 
     if len(outcomes_of_interest) == 0:
         outcomes_of_interest = ['total_asset_loss',
+                                'tot_exposed_asset',
                                 'total_consumption_loss',
                                 'n_affected_people',
                                 'n_new_poor',
@@ -314,7 +277,7 @@ def get_policy_effectiveness_tab(outcomes: pd.DataFrame) -> pd.DataFrame:
     return df
 
 
-def get_weeks_in_poverty_tab(outcomes: pd.DataFrame, max_years: int = 10) -> pd.DataFrame:
+def get_weeks_in_poverty_tab(outcomes: pd.DataFrame) -> pd.DataFrame:
     '''Get the average across scenarios number of weeks in poverty for each district.
 
     Args:

unbreakable/analysis/outcome_names.yaml

@@ -0,0 +1,24 @@
+- total_population
+- total_asset_loss
+- total_consumption_loss
+- tot_exposed_asset
+- tot_asset_surv
+- expected_loss_frac
+- n_affected_people
+- annual_average_consumption
+- poverty_line_adjusted
+- district_pml
+- n_poor_initial
+- n_poor_affected
+- n_new_poor
+- initial_poverty_gap
+- new_poverty_gap_initial
+- new_poverty_gap_all
+- annual_average_consumption_loss
+- annual_average_consumption_loss_pct
+- r
+- mean_recovery_rate
+- tot_wellbeing_loss
+- weighted_vuln_quint
+- weighted_vuln_dec
+- years_in_poverty

unbreakable/analysis/visualizer.py

@@ -137,13 +137,14 @@ def raincloud_plot(outcomes: pd.DataFrame, savefig: bool,  x_columns: list = [],
                                                                                    f'M={df[x_column].median():.2f}',
                                                                                    horizontalalignment='left', size='small', color='black')
 
-            initial_poverty_gap = df['initial_poverty_gap'].iloc[0]
+            # initial_poverty_gap = df['initial_poverty_gap'].iloc[0]
             # Add initial poverty gap as in the legend to the plot
-            if x_column == 'new_poverty_gap_all' or x_column == 'new_poverty_gap_initial':
-                ax[districts.index(district) // 3, districts.index(district) % 3].text(0.025, 0.9,
-                                                                                       f'Poverty gap before disaster={initial_poverty_gap:.2f}',
-                                                                                       horizontalalignment='left', size='small', color='black',
-                                                                                       transform=ax[districts.index(district) // 3, districts.index(district) % 3].transAxes)
+            # if x_column == 'new_poverty_gap_all' or x_column == 'new_poverty_gap_initial':
+            #     ax[districts.index(district) // 3, districts.index(district) % 3].text(0.025, 0.9,
+            #                                                                            f'Poverty gap before disaster={initial_poverty_gap:.2f}',
+            #                                                                            horizontalalignment='left', size='small', color='black',
+            #                                                                            transform=ax[districts.index(district) // 3, districts.index(district) % 3].transAxes)
+
         fig.tight_layout()
         if savefig:
             plt.savefig(

unbreakable/data/writer.py

@@ -3,9 +3,10 @@
 
 import pandas as pd
 import numpy as np
+import yaml
 
 
-def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_recover) -> dict:
+def get_outcomes(households: pd.DataFrame, tot_exposed_asset: float, expected_loss_frac: float, years_to_recover: int, consump_util: float) -> dict:
     '''Calculate outcomes of interest from the simulation model.
 
     Args:
@@ -50,13 +51,14 @@ def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_rec
     annual_average_consumption_loss, annual_average_consumption_loss_pct = calculate_average_annual_consumption_loss(
         affected_households, years_to_recover)
 
-    r = calculate_resilience(affected_households)
+    tot_wellbeing_loss = calc_tot_wellbeing_loss(households, consump_util)
+    r = calculate_resilience(affected_households, tot_wellbeing_loss)
 
     # Get the weighted average vulnerability by consumption quintile and decile
     weighted_vuln_quint = get_weighted_vuln(affected_households, quintile=True)
     weighted_vuln_dec = get_weighted_vuln(affected_households, quintile=False)
 
-    return {
+    outcomes = {
         'total_population': total_population,
         'total_asset_loss': total_asset_loss,
         'total_consumption_loss': total_consumption_loss,
@@ -77,11 +79,18 @@ def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_rec
         'annual_average_consumption_loss_pct': annual_average_consumption_loss_pct,
         'r': r,
         'mean_recovery_rate': mean_recovery_rate,
+        'tot_wellbeing_loss': tot_wellbeing_loss,
         'weighted_vuln_quint': weighted_vuln_quint,
         'weighted_vuln_dec': weighted_vuln_dec,
         'years_in_poverty': years_in_poverty
     }
 
+    # Save outcome names in a yaml file to pick up in preprocessing
+    with open('analysis/outcome_names.yaml', 'w') as f:
+        yaml.dump(list(outcomes.keys()), f)
+
+    return outcomes
+
 
 def find_poor(households: pd.DataFrame, poverty_line: float, years_to_recover: int) -> tuple:
     '''Get the poor at the beginning of the simulation and the poor at the end of the simulation
@@ -226,13 +235,14 @@ def calculate_average_annual_consumption_loss(affected_households: pd.DataFrame,
     return annual_average_consumption_loss, annual_average_consumption_loss_pct
 
 
-def calculate_resilience(affected_households: pd.DataFrame) -> float:
+def calculate_resilience(affected_households: pd.DataFrame, tot_wellbeing_loss: float) -> float:
     '''Calculate socio-economic resilience of affected households.
 
     Socio-economic resilience is a ratio of asset loss to consumption loss.
 
     Args:
         affected_households (pd.DataFrame): Affected households.
+        tot_wellbeing_loss (float): Total wellbeing loss.
 
     Returns:
         float: Socio-economic resilience
@@ -247,7 +257,8 @@ def calculate_resilience(affected_households: pd.DataFrame) -> float:
         r = np.nan
 
     else:
-        r = total_asset_damage / total_consumption_loss
+        # r = total_asset_damage / total_consumption_loss
+        r = total_asset_damage / tot_wellbeing_loss
 
     return r
 
@@ -277,3 +288,20 @@ def get_weighted_vuln(affected_households: pd.DataFrame, quintile: bool) -> dict
         pop_by_d = df.groupby('decile').sum(numeric_only=True)[['popwgt']]
         average_v_by_d = v_weighted_by_d['v_weighted'].div(pop_by_d['popwgt'])
         return average_v_by_d.to_dict()
+
+
+def calc_tot_wellbeing_loss(households: pd.DataFrame, consump_util: float) -> float:
+    '''Calculate wellbeing loss.
+
+    Args:
+        households (pd.DataFrame): Households.
+        consump_util (float): Consumption utility.
+
+    Returns:
+        float: Total wellbeing loss.
+    '''
+    x = (households['aeexp'].multiply(households['popwgt'])
+         ).sum() / households['popwgt'].sum()
+    W = x**(-consump_util)
+    wellbeing_loss = - households['wellbeing'].sum() / W
+    return wellbeing_loss

unbreakable/model.py

@@ -50,7 +50,7 @@ def model(**params) -> dict:
 
     # Uncertainties
     poverty_bias = params['poverty_bias']
-    consumption_utility = params['consumption_utility']
+    consump_util = params['consump_util']
     discount_rate = params['discount_rate']
     lambda_increment = params['lambda_increment']
     income_and_expenditure_growth = params['income_and_expenditure_growth']
@@ -99,16 +99,16 @@ def model(**params) -> dict:
                                 .pipe(calculate_exposure, poverty_bias, calc_exposure_params)
                                 .pipe(identify_affected, ident_affected_params)
                                 .pipe(apply_policy, my_policy)
-                                .pipe(calculate_recovery_rate, consumption_utility, discount_rate, lambda_increment, years_to_recover)
-                                .pipe(calculate_wellbeing, consumption_utility, discount_rate, income_and_expenditure_growth, years_to_recover, add_income_loss, cash_transfer))
+                                .pipe(calculate_recovery_rate, consump_util, discount_rate, lambda_increment, years_to_recover)
+                                .pipe(calculate_wellbeing, consump_util, discount_rate, income_and_expenditure_growth, years_to_recover, add_income_loss, cash_transfer))
 
         if save_households:
             Path(f'../experiments/households/').mkdir(parents=True, exist_ok=True)
             households.to_csv(
                 f'../experiments/households/{district}_{random_seed}.csv')
 
         array_outcomes = np.array(list(get_outcomes(
-            households, tot_exposed_asset, expected_loss_frac, years_to_recover).values()))
+            households, tot_exposed_asset, expected_loss_frac, years_to_recover, consump_util).values()))
 
         # * To check whether we have different households affected in different runs
         # if district == 'Castries':