fix wellbeing calculation

mikhailsirenko · Aug 23, 2023 · f3dd986 · f3dd986
1 parent b37ab8f
commit f3dd986
Show file tree

Hide file tree

Showing 5 changed files with 32 additions and 19 deletions.
diff --git a/unbreakable/analysis/analyzer.py b/unbreakable/analysis/analyzer.py
@@ -203,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',

diff --git 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
@@ -300,9 +300,8 @@ def calc_tot_wellbeing_loss(households: pd.DataFrame, consump_util: float) -> fl
     Returns:
         float: Total wellbeing loss.
     '''
-    x = households['aeexp'].multiply(
-        households['popwgt']) / households['popwgt'].sum()
+    x = (households['aeexp'].multiply(households['popwgt'])
+         ).sum() / households['popwgt'].sum()
     W = x**(-consump_util)
-    wellbeing = households['wellbeing']
-    wellbeing_loss = - wellbeing / W
-    return wellbeing_loss.sum()
+    wellbeing_loss = - households['wellbeing'].sum() / W
+    return wellbeing_loss
diff --git a/unbreakable/modules/integrator.py b/unbreakable/modules/integrator.py
@@ -199,10 +199,9 @@ def calculate_wellbeing(households: pd.DataFrame, consump_util: float, discount_
         # 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,8 +235,20 @@ def calculate_wellbeing(households: pd.DataFrame, consump_util: float, discount_
                                 < poverty_line_adjusted, 'weeks_pov'] += 1
 
         # Integrate well-being
-        affected_households['wellbeing'] = affected_households['c_t_unaffected']**(1 - consump_util)\
-            / (1-consump_util) * dt\
+        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)

diff --git 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)