diff --git a/config/SaintLucia.yaml b/config/SaintLucia.yaml index 10c7d26..1777f52 100644 --- a/config/SaintLucia.yaml +++ b/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: diff --git a/unbreakable/analysis/analyzer.py b/unbreakable/analysis/analyzer.py index d892de1..f1937dd 100644 --- a/unbreakable/analysis/analyzer.py +++ b/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: diff --git a/unbreakable/analysis/outcome_names.yaml b/unbreakable/analysis/outcome_names.yaml new file mode 100644 index 0000000..7455465 --- /dev/null +++ b/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 diff --git a/unbreakable/analysis/visualizer.py b/unbreakable/analysis/visualizer.py index b504cd0..ba9af96 100644 --- a/unbreakable/analysis/visualizer.py +++ b/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( diff --git a/unbreakable/data/writer.py b/unbreakable/data/writer.py index 1b219f1..ee7c68c 100644 --- a/unbreakable/data/writer.py +++ b/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 diff --git a/unbreakable/model.py b/unbreakable/model.py index 95918fb..8426e52 100644 --- a/unbreakable/model.py +++ b/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,8 +99,8 @@ 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) @@ -108,7 +108,7 @@ def model(**params) -> dict: 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': diff --git a/unbreakable/modules/integrator.py b/unbreakable/modules/integrator.py index f81e61b..3057d6c 100644 --- a/unbreakable/modules/integrator.py +++ b/unbreakable/modules/integrator.py @@ -6,7 +6,7 @@ import pickle -def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float, discount_rate: float, lambda_increment: float, years_to_recover: int) -> pd.DataFrame: +def calculate_recovery_rate(households: pd.DataFrame, consump_util: float, discount_rate: float, lambda_increment: float, years_to_recover: int) -> pd.DataFrame: '''Calculates the recovery rate for each affected household. Recovery rate is a function of the household vulnerability (v), @@ -14,7 +14,7 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float Args: households (pd.DataFrame): Households data frame. - consumption_utility (float): Consumption utility. + consump_util (float): Consumption utility. discount_rate (float): Discount rate. lambda_increment (float): Lambda increment for the integration. years_to_recover (int): Number of years to recover. @@ -42,7 +42,7 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float # Calculate the recovery rate for each affected household affected_households['recovery_rate'] = affected_households['v'].apply(lambda x: integrate_and_find_recovery_rate( - x, results, consumption_utility, discount_rate, average_productivity, lambda_increment, years_to_recover)) + x, results, consump_util, discount_rate, average_productivity, lambda_increment, years_to_recover)) # Set recovery rate to zero for unaffected households households['recovery_rate'] = 0 @@ -54,13 +54,13 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float return households -def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumption_utility: float, discount_rate: float, average_productivity: float, lambda_increment: float, years_to_recover: int) -> float: +def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consump_util: float, discount_rate: float, average_productivity: float, lambda_increment: float, years_to_recover: int) -> float: '''Find recovery rate (lambda) given the value of `v` (household vulnerability). Args: v (float): Household vulnerability. results (pd.DataFrame): Data frame to store integration results. - consumption_utility (float): Consumption utility. + consump_util (float): Consumption utility. discount_rate (float): Discount rate. average_productivity (float): Average productivity. lambda_increment (float): Lambda increment for the integration. @@ -92,7 +92,7 @@ def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumptio factor = average_productivity + _lambda part1 = average_productivity - \ factor * v * np.e**(-_lambda * _t) - part1 = part1**(-consumption_utility) + part1 = part1**(-consump_util) part2 = _t * factor - 1 part3 = np.e**(-_t * (discount_rate + _lambda)) dwdlambda += part1 * part2 * part3 * dt @@ -107,12 +107,12 @@ def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumptio _lambda += lambda_increment -def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, discount_rate: float, income_and_expenditure_growth: float, years_to_recover: int, add_income_loss: bool, cash_transfer: dict = {}) -> pd.DataFrame: +def calculate_wellbeing(households: pd.DataFrame, consump_util: float, discount_rate: float, income_and_expenditure_growth: float, years_to_recover: int, add_income_loss: bool, cash_transfer: dict = {}) -> pd.DataFrame: '''Calculates consumption loss and well-being for each affected household. Args: households (pd.DataFrame): Households data frame. - consumption_utility (float): Consumption utility. + consump_util (float): Consumption utility. discount_rate (float): Discount rate. income_and_expenditure_growth (float): Income and expenditure growth. years_to_recover (int): Number of years to recover. @@ -134,9 +134,8 @@ def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, di 'net_consumption_loss_NPV', 'c_t', 'c_t_unaffected', - 'w_final', 'weeks_pov', - 'w_final2'] + 'wellbeing'] # Set values to zero households[columns] = 0 @@ -200,10 +199,9 @@ def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, di # Check if any of the households has a negative consumption # It may happen if a household have very low savings or expenditure if (affected_households['c_t'] < 0).any(): - # If so, then set the consumption to 0 - affected_households.loc[affected_households['c_t'] < 0, 'c_t'] = 0 - # * Previously it was set to 1 - # affected_households.loc[affected_households['c_t'] < 1, 'c_t'] = 1 + # If so, then set the consumption to 1 + # We must have 1 to to avoid -inf in the wellbeing integration + affected_households.loc[affected_households['c_t'] < 0, 'c_t'] = 1 # Consumption after the disaster should be lower than or equal to consumption before the disaster if (affected_households['c_t'] > affected_households['c_t_unaffected']).any(): @@ -236,13 +234,24 @@ def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, di affected_households.loc[affected_households['c_t'] < poverty_line_adjusted, 'weeks_pov'] += 1 - # Integrate wellbeing - # !!!: Fix wellbeing integration, currently it is nan - affected_households['w_final'] += dt * (affected_households['c_t'])**(1 - consumption_utility) * \ - np.e**(-discount_rate * _t) / (1 - consumption_utility) - - affected_households['w_final2'] += affected_households['c_t_unaffected']**(1-consumption_utility)/(1-consumption_utility)*dt*( - (1-((affected_households['c_t_unaffected'] - affected_households['c_t'])/affected_households['c_t_unaffected'])*np.e**(-affected_households['recovery_rate']*_t))**(1-consumption_utility)-1)*np.e**(-discount_rate*_t) + # Integrate well-being + a = affected_households['c_t_unaffected']**(1 - consump_util) + b = (1 - consump_util) * dt + c = ((1 - ((affected_households['c_t_unaffected'] - affected_households['c_t']) / affected_households['c_t_unaffected']) + * np.e**(-affected_households['recovery_rate'] * _t))**(1 - consump_util) - 1) + d = np.e**(-discount_rate * _t) + e = a.div(b).mul(c).mul(d) + f = affected_households['c_t_unaffected']**(1 - consump_util)\ + / (1 - consump_util) * dt\ + * ((1 - ((affected_households['c_t_unaffected'] - affected_households['c_t']) / affected_households['c_t_unaffected']) + * np.e**(-affected_households['recovery_rate'] * _t))**(1 - consump_util) - 1)\ + * np.e**(-discount_rate * _t) + + affected_households['wellbeing'] += affected_households['c_t_unaffected']**(1 - consump_util)\ + / (1 - consump_util) * dt\ + * ((1 - ((affected_households['c_t_unaffected'] - affected_households['c_t']) / affected_households['c_t_unaffected']) + * np.e**(-affected_households['recovery_rate'] * _t))**(1 - consump_util) - 1)\ + * np.e**(-discount_rate * _t) # Use to examine individual consumption recovery # Save consumption recovery value at the time _t diff --git a/unbreakable/run.py b/unbreakable/run.py index caebcb8..7b2ce4c 100644 --- a/unbreakable/run.py +++ b/unbreakable/run.py @@ -40,13 +40,14 @@ for district in constants['districts']] # Specify the number of scenarios and policies - n_scenarios = 5 + n_scenarios = 48 n_policies = 0 + # Perform the experiments + # results = perform_experiments( # models=my_model, scenarios=n_scenarios, policies=n_policies) - # Perform the experiments with MultiprocessingEvaluator(my_model, n_processes=12) as evaluator: results = evaluator.perform_experiments( scenarios=n_scenarios, policies=n_policies)