diff --git a/src/scripts/costing/cost_estimation.py b/src/scripts/costing/cost_estimation.py index 02d0971955..9dc5a120c4 100644 --- a/src/scripts/costing/cost_estimation.py +++ b/src/scripts/costing/cost_estimation.py @@ -156,7 +156,8 @@ def get_discount_factor(year): # Compute the cumulative discount factor as the product of (1 + discount_rate) for all previous years discount_factor = 1 for y in range(_initial_year + 1, - year + 1): # only starting from initial year + 1 as the discount factor for initial year should be 1 + year + 1): # only starting from initial year + 1 as the discount factor for initial year + # should be 1 discount_factor *= (1 + _discount_rate.get(y, 0)) # Default to 0 if year not in dictionary return discount_factor else: @@ -282,6 +283,8 @@ def clean_equipment_name(name: str, equipment_drop_list = None) -> str: def estimate_input_cost_of_scenarios(results_folder: Path, resourcefilepath: Path, + suspended_results_folder: Path = None, + scaling_factor: Optional[float] = None, _draws: Optional[list[int]] = None, _runs: Optional[list[int]] = None, summarize: bool = False, @@ -298,6 +301,13 @@ def estimate_input_cost_of_scenarios(results_folder: Path, Path to the directory containing simulation output files. resourcefilepath : Path, optional Path to the resource files + suspended_results_folder: Path, optional + Path to the directory containing suspended simulation output files (using the suspend and resume functionality), + This is used to extract the scaling_factor to scale result to actual population size. If None, then the + 'scaling_factor' is obtained from the results_folder. + scaling_factor: float, optional + Scaling factor used to scale results to actual population size. If provided, this is used in place of the + scaling_factor logged in the suspended_results_folder or results_folder. _draws : list, optional Specific draws to include in the cost estimation. Defaults to all available draws. _runs : list, optional @@ -316,8 +326,10 @@ def estimate_input_cost_of_scenarios(results_folder: Path, Returns: ------- pd.DataFrame - A dataframe containing discounted costs disaggregated by category, sub-category, category-specific subgroup, year, draw, and run. - Note that if a discount rate is used, the dataframe will provide cost as the NPV during the first year of the dataframe + A dataframe containing discounted costs disaggregated by category, sub-category, category-specific subgroup, + year, draw, and run. + Note that if a discount rate is used, the dataframe will provide cost as the NPV during the first year of the + dataframe """ # Useful common functions @@ -343,6 +355,8 @@ def melt_model_output_draws_and_runs(_df, id_vars): _draws = range(0, info['number_of_draws']) if _runs is None: _runs = range(0, info['runs_per_draw']) + if suspended_results_folder is None: + suspended_results_folder = results_folder # Load cost input files # ------------------------ @@ -354,7 +368,8 @@ def melt_model_output_draws_and_runs(_df, id_vars): facility_id_levels_dict = dict(zip(mfl['Facility_ID'], mfl['Facility_Level'])) fac_levels = set(mfl.Facility_Level) - # If variable discount rate is provided, use the average across the relevant years for the purpose of annuitization of HR and equipment costs + # If variable discount rate is provided, use the average across the relevant years for the purpose of annuitization + # of HR and equipment costs def calculate_annuitization_rate(_discount_rate, _years): if isinstance(_discount_rate, (int, float)): # Single discount rate, return as is @@ -477,25 +492,69 @@ def merge_cost_and_model_data(cost_df, model_df, varnames): return merged_df # Get available staff count for each year and draw - def get_staff_count_by_facid_and_officer_type(_df: pd.Series) -> pd.Series: - """Summarise the parsed logged-key results for one draw (as dataframe) into a pd.Series.""" - _df = _df.set_axis(_df['date'].dt.year).drop(columns=['date']) - _df.index.name = 'year' + def get_staff_count_by_facid_and_officer_type(_df: pd.DataFrame) -> pd.Series: + """ + Convert logged staff dictionary output into tidy format, + summing staff counts across all clinic columns. + + Returns pd.Series indexed by: + (year, FacilityID, Officer) + """ + + df = _df.copy() + df["year"] = df["date"].dt.year + df = df.drop(columns=["date"]) + + clinic_cols = df.columns.difference(["year"]) + + long_frames = [] + + for clinic in clinic_cols: + expanded = df[[clinic, "year"]].copy() + expanded = expanded[expanded[clinic].notna()] + + expanded_dict = expanded[clinic].apply(pd.Series) + expanded_dict["year"] = expanded["year"].values + + long_frames.append(expanded_dict) + + # Combine all clinics + combined = pd.concat(long_frames, ignore_index=True) - def change_to_standard_flattened_index_format(col): - parts = col.split("_", 3) # Split by "_" only up to 3 parts - if len(parts) > 2: - return parts[0] + "=" + parts[1] + "|" + parts[2] + "=" + parts[ - 3] # Rejoin with "I" at the second occurrence - return col # If there's no second underscore, return the string as it is + # Melt to long format + long_df = ( + combined + .melt(id_vars=["year"], + var_name="facility_officer", + value_name="count") + .dropna(subset=["count"]) + ) + + # Split FacilityID and Officer + parts = long_df["facility_officer"].str.split("_Officer_", expand=True) - _df.columns = [change_to_standard_flattened_index_format(col) for col in _df.columns] + long_df["FacilityID"] = ( + parts[0] + .str.replace("FacilityID_", "", regex=False) + .astype(int) + ) + long_df["Officer"] = parts[1] + + # SUM ACROSS CLINICS HERE + result = ( + long_df + .groupby(["year", "FacilityID", "Officer"])["count"] + .sum() + .sort_index() + ) - return unflatten_flattened_multi_index_in_logging(_df).stack(level=[0, 1]) # expanded flattened axis + return result # Staff count by Facility ID available_staff_count_by_facid_and_officertype = extract_results( Path(results_folder), + suspended_results_folder=suspended_results_folder, + scaling_factor=scaling_factor, module='tlo.methods.healthsystem.summary', key='number_of_hcw_staff', custom_generate_series=get_staff_count_by_facid_and_officer_type, @@ -510,8 +569,14 @@ def change_to_standard_flattened_index_format(col): available_staff_count_by_facid_and_officertype['Facility_Level'] = available_staff_count_by_facid_and_officertype[ 'Facility_ID'].map(facility_id_levels_dict) idx = pd.IndexSlice - available_staff_count_by_level_and_officer_type = available_staff_count_by_facid_and_officertype.drop( - columns=[idx['Facility_ID']]).groupby([idx['year'], idx['Facility_Level'], idx['OfficerType']]).sum() + # Sum staff counts across facilities within each year, Facility_Level and OfficerType. The grouping keys are + # passed as Series and the (draw, run) value columns selected explicitly, because dropping labels from this + # unsorted column MultiIndex (which mixes strings and integers, so cannot be sorted) raises a PerformanceWarning + _staff_df = available_staff_count_by_facid_and_officertype + _key_cols = ['year', 'Facility_Level', 'OfficerType'] + _value_cols = [col for col in _staff_df.columns if col[0] not in _key_cols + ['Facility_ID']] + available_staff_count_by_level_and_officer_type = _staff_df[_value_cols].groupby( + [_staff_df[(col, '')].rename(col) for col in _key_cols]).sum() available_staff_count_by_level_and_officer_type = melt_model_output_draws_and_runs( available_staff_count_by_level_and_officer_type.reset_index(), id_vars=['year', 'Facility_Level', 'OfficerType']) @@ -519,22 +584,72 @@ def change_to_standard_flattened_index_format(col): 'Facility_Level'].astype(str) # make sure facility level is stored as string available_staff_count_by_level_and_officer_type = available_staff_count_by_level_and_officer_type.drop( available_staff_count_by_level_and_officer_type[available_staff_count_by_level_and_officer_type[ - 'Facility_Level'] == '5'].index) # drop headquarters because we're only concerned with staff engaged in service delivery + 'Facility_Level'] == '5'].index) # drop headquarters + # because we're only concerned with staff engaged in service delivery available_staff_count_by_level_and_officer_type.rename(columns={'value': 'staff_count'}, inplace=True) # Get list of cadres which were utilised in each run to get the count of staff used in the simulation - # Note that we still cost the full staff count for any cadre-Facility_Level combination that was ever used in a run, and - # not the amount of time which was used - def get_capacity_used_by_officer_type_and_facility_level(_df: pd.Series) -> pd.Series: - """Summarise the parsed logged-key results for one draw (as dataframe) into a pd.Series.""" - _df = _df.set_axis(_df['date'].dt.year).drop(columns=['date']) - _df.index.name = 'year' - return unflatten_flattened_multi_index_in_logging(_df).stack(level=[0, 1]) # expanded flattened axis + # Note that we still cost the full staff count for any cadre-Facility_Level combination that was ever used in a run, + # and not the amount of time which was used + def get_capacity_used_by_officer_type_and_facility_level( + _df: pd.DataFrame + ) -> pd.Series: + """ + Parse logging output and return a Series indexed by: + (year, OfficerType, FacilityLevel) + + Collapses (sums) across clinics. + Uses facility_id_levels_dict to map FacilityID → FacilityLevel. + """ + + # ---- 1. Set year index ---- + _df = _df.set_axis(_df["date"].dt.year).drop(columns=["date"]) + _df.index.name = "year" + + # ---- 2. Unflatten logging columns ---- + _df = unflatten_flattened_multi_index_in_logging(_df) + + # Expect columns like: + # ('Clinic', 'facID_and_officer') + + col_df = _df.columns.to_frame(index=False) + + # ---- 3. Extract OfficerType ---- + col_df["OfficerType"] = ( + col_df["facID_and_officer"] + .str.split("_Officer_") + .str[-1] + ) + + # ---- 4. Extract FacilityID ---- + col_df["FacilityID"] = ( + col_df["facID_and_officer"] + .str.split("_Officer_") + .str[0] + .str.replace("FacilityID_", "", regex=False) + .astype(int) + ) + + # ---- 5. Map to FacilityLevel ---- + col_df["FacilityLevel"] = col_df["FacilityID"].map(facility_id_levels_dict) + + # ---- 6. Rebuild MultiIndex (drop clinic level) ---- + _df.columns = pd.MultiIndex.from_frame( + col_df[["OfficerType", "FacilityLevel"]] + ) + + # ---- 7. Collapse across clinics ---- + _df = _df.groupby(level=["OfficerType", "FacilityLevel"], axis=1).sum() + + # ---- 8. Return stacked format ---- + return _df.stack(["OfficerType", "FacilityLevel"]) annual_capacity_used_by_cadre_and_level = extract_results( Path(results_folder), + suspended_results_folder=suspended_results_folder, + scaling_factor=scaling_factor, module='tlo.methods.healthsystem.summary', - key='Capacity_By_OfficerType_And_FacilityLevel', + key='Capacity_By_FacID_and_Officer', custom_generate_series=get_capacity_used_by_officer_type_and_facility_level, do_scaling=False, ) @@ -552,7 +667,10 @@ def get_capacity_used_by_officer_type_and_facility_level(_df: pd.Series) -> pd.S average_capacity_used_by_cadre_and_level[average_capacity_used_by_cadre_and_level['capacity_used'] != 0][ ['OfficerType', 'FacilityLevel', 'draw', 'run']] print( - f"Out of {average_capacity_used_by_cadre_and_level.groupby(['OfficerType', 'FacilityLevel']).size().count()} cadre and level combinations available, {list_of_cadre_and_level_combinations_used.groupby(['OfficerType', 'FacilityLevel']).size().count()} are used across the simulations") + f"Out of {average_capacity_used_by_cadre_and_level.groupby(['OfficerType', 'FacilityLevel']).size().count()} " + f"cadre and level combinations available, " + f"{list_of_cadre_and_level_combinations_used.groupby(['OfficerType', 'FacilityLevel']).size().count()} " + f"are used across the simulations") list_of_cadre_and_level_combinations_used = list_of_cadre_and_level_combinations_used.rename( columns={'FacilityLevel': 'Facility_Level'}) @@ -564,11 +682,13 @@ def get_capacity_used_by_officer_type_and_facility_level(_df: pd.Series) -> pd.S if (cost_only_used_staff): print( - "The input for 'cost_only_used_staff' implies that only cadre-level combinations which have been used in the run are costed") + "The input for 'cost_only_used_staff' implies that only cadre-level combinations which have been used in " + "the run are costed") staff_size_chosen_for_costing = used_staff_count_by_level_and_officer_type else: print( - "The input for 'cost_only_used_staff' implies that all staff are costed regardless of the cadre-level combinations which have been used in the run are costed") + "The input for 'cost_only_used_staff' implies that all staff are costed regardless of the cadre-level " + "combinations which have been used in the run are costed") staff_size_chosen_for_costing = available_staff_count_by_level_and_officer_type # Calculate various components of HR cost @@ -607,7 +727,8 @@ def calculate_npv_past_training_expenses_by_row(row, r=_discount_rate): if partial_year > 0: npv += annual_cost * partial_year * (1 + r) ** (1 + r) - # Add recruitment cost assuming this happens during the partial year or the year after graduation if partial year == 0 + # Add recruitment cost assuming this happens during the partial year or the year after graduation if + # partial year == 0 npv += row['recruitment_cost_per_person_recruited_usd'] * (1 + r) return npv @@ -619,34 +740,33 @@ def calculate_npv_past_training_expenses_by_row(row, r=_discount_rate): npv_values.append(npv) preservice_training_cost['npv_of_training_and_recruitment_cost'] = npv_values - preservice_training_cost['npv_of_training_and_recruitment_cost_per_recruit'] = preservice_training_cost[ - 'npv_of_training_and_recruitment_cost'] * \ - (1 / (preservice_training_cost[ - 'absorption_rate_of_students_into_public_workforce'] + - preservice_training_cost[ - 'proportion_of_workforce_recruited_from_abroad'])) * \ - (1 / preservice_training_cost[ - 'graduation_rate']) * (1 / - preservice_training_cost[ - 'licensure_exam_passing_rate']) - if _discount_rate == 0: # if the discount rate is 0, then the pre-service + recruitment cost simply needs to be divided by the number of years in tenure + preservice_training_cost['npv_of_training_and_recruitment_cost_per_recruit'] \ + = (preservice_training_cost['npv_of_training_and_recruitment_cost'] * + (1 / (preservice_training_cost['absorption_rate_of_students_into_public_workforce'] + + preservice_training_cost['proportion_of_workforce_recruited_from_abroad'])) * + (1 / preservice_training_cost['graduation_rate']) * + (1 /preservice_training_cost['licensure_exam_passing_rate'])) + if _discount_rate == 0: # if the discount rate is 0, then the pre-service + recruitment cost simply + # needs to be divided by the number of years in tenure preservice_training_cost['annuitisation_rate'] = preservice_training_cost[ 'average_length_of_tenure_in_the_public_sector'] else: preservice_training_cost['annuitisation_rate'] = 1 + (1 - (1 + annuitization_rate) ** ( -preservice_training_cost[ 'average_length_of_tenure_in_the_public_sector'] + 1)) / annuitization_rate - preservice_training_cost['annuitised_training_and_recruitment_cost_per_recruit'] = preservice_training_cost[ - 'npv_of_training_and_recruitment_cost_per_recruit'] / \ - preservice_training_cost[ - 'annuitisation_rate'] - - # Cost per student trained * 1/Rate of absorption from the local and foreign graduates * 1/Graduation rate * attrition rate - # the inverse of attrition rate is the average expected tenure; and the preservice training cost needs to be divided by the average tenure + preservice_training_cost['annuitised_training_and_recruitment_cost_per_recruit'] = \ + (preservice_training_cost['npv_of_training_and_recruitment_cost_per_recruit'] / + preservice_training_cost['annuitisation_rate']) + + # Cost per student trained * 1/Rate of absorption from the local and foreign graduates + # * 1/Graduation rate * attrition rate + # the inverse of attrition rate is the average expected tenure; and the preservice training cost needs to + # be divided by the average tenure preservice_training_cost['cost'] = preservice_training_cost[ 'annuitised_training_and_recruitment_cost_per_recruit'] * \ preservice_training_cost['staff_count'] * preservice_training_cost[ - 'annual_attrition_rate'] # not multiplied with attrition rate again because this is already factored into 'Annual_cost_per_staff_recruited' + 'annual_attrition_rate'] # not multiplied with attrition rate again + # because this is already factored into 'Annual_cost_per_staff_recruited' preservice_training_cost = preservice_training_cost[ ['draw', 'run', 'year', 'OfficerType', 'Facility_Level', 'cost']] @@ -680,7 +800,8 @@ def label_rows_of_cost_dataframe(_df, label_var, label): # Initialize HR with the salary data if (cost_only_used_staff): human_resource_costs = retain_relevant_column_subset( - label_rows_of_cost_dataframe(salary_for_staff, 'cost_subcategory', 'salary_for_cadres_used'), 'OfficerType') + label_rows_of_cost_dataframe(salary_for_staff, 'cost_subcategory', 'salary_for_cadres_used'), + 'OfficerType') # Concatenate additional cost categories additional_costs = [ (preservice_training_cost, 'preservice_training_and_recruitment_cost_for_attrited_workers'), @@ -689,7 +810,8 @@ def label_rows_of_cost_dataframe(_df, label_var, label): ] else: human_resource_costs = retain_relevant_column_subset( - label_rows_of_cost_dataframe(salary_for_staff, 'cost_subcategory', 'salary_for_all_staff'), 'OfficerType') + label_rows_of_cost_dataframe(salary_for_staff, 'cost_subcategory', 'salary_for_all_staff'), + 'OfficerType') # Concatenate additional cost categories additional_costs = [ (preservice_training_cost, 'preservice_training_and_recruitment_cost_for_attrited_workers'), @@ -741,6 +863,8 @@ def get_counts_of_items_requested(_df): cons_req = extract_results( results_folder, + suspended_results_folder=suspended_results_folder, + scaling_factor=scaling_factor, module='tlo.methods.healthsystem.summary', key='Consumables', custom_generate_series=get_counts_of_items_requested, @@ -762,7 +886,8 @@ def get_counts_of_items_requested(_df): # 2.1 Cost of consumables dispensed # --------------------------------------------------------------------------------------------------------------- # Multiply number of items needed by cost of consumable - # consumables_dispensed.columns = consumables_dispensed.columns.get_level_values(0).str() + "_" + consumables_dispensed.columns.get_level_values(1) # Flatten multi-level columns for pandas merge + # consumables_dispensed.columns = consumables_dispensed.columns.get_level_values(0).str() + "_" + + # consumables_dispensed.columns.get_level_values(1) # Flatten multi-level columns for pandas merge unit_costs['consumables'].columns = pd.MultiIndex.from_arrays( [unit_costs['consumables'].columns, [''] * len(unit_costs['consumables'].columns)]) cost_of_consumables_dispensed = consumables_dispensed.merge(unit_costs['consumables'], on=idx['Item_Code'], @@ -796,8 +921,20 @@ def get_counts_of_items_requested(_df): left_on='Item_Code', right_on='item_code', validate='m:1', how='left') + # Identify rows where excess_stock_proportion_of_dispensed is NaN + missing_excess_stock = ( + cost_of_excess_consumables_stocked + ['excess_stock_proportion_of_dispensed'] + .isna() + ) + + # Fill missing values with the average inflow-to-outflow ratio minus 1 + fill_value = average_inflow_to_outflow_ratio_ratio - 1 + cost_of_excess_consumables_stocked.loc[ - cost_of_excess_consumables_stocked.excess_stock_proportion_of_dispensed.isna(), 'excess_stock_proportion_of_dispensed'] = average_inflow_to_outflow_ratio_ratio - 1 # TODO disaggregate the average by program + missing_excess_stock, + 'excess_stock_proportion_of_dispensed' + ] = fill_value # TODO: disaggregate the average by program cost_of_excess_consumables_stocked[quantity_columns] = cost_of_excess_consumables_stocked[ quantity_columns].multiply(cost_of_excess_consumables_stocked[idx[price_column]], axis=0) cost_of_excess_consumables_stocked[quantity_columns] = cost_of_excess_consumables_stocked[ @@ -815,7 +952,8 @@ def get_counts_of_items_requested(_df): def melt_and_label_consumables_cost(_df, label): multi_index = pd.MultiIndex.from_tuples(_df.columns) _df.columns = multi_index - # Select 'Item_Code', 'year', and all columns where both levels of the MultiIndex are numeric (these are the (draw,run) columns with cost values) + # Select 'Item_Code', 'year', and all columns where both levels of the MultiIndex are numeric + # (these are the (draw,run) columns with cost values) selected_columns = [col for col in _df.columns if (col[0] in ['Item_Code', 'year']) or (isinstance(col[0], int) and isinstance(col[1], int))] _df = _df[selected_columns] # Subset the dataframe with the selected columns @@ -827,20 +965,23 @@ def melt_and_label_consumables_cost(_df, label): melted_df['consumable'] = melted_df['Item_Code'].map(consumables_dict) melted_df['cost_subcategory'] = label melted_df[ - 'Facility_Level'] = 'all' # TODO this is temporary until 'tlo.methods.healthsystem.summary' only logs consumable at the aggregate level + 'Facility_Level'] = 'all' + # TODO this is temporary until 'tlo.methods.healthsystem.summary' only logs consumable at the aggregate level melted_df = melted_df.rename(columns={'value': 'cost'}) return melted_df def disaggregate_separately_managed_medical_supplies_from_consumable_costs(_df, _consumables_dict, - # This is a dictionary mapping codes to names + # This is a dictionary mapping codes to names list_of_unique_medical_products): reversed_consumables_dict = {value: key for key, value in _consumables_dict.items()} # reverse dictionary to map names to codes new_df = _df.copy() new_df['item_code'] = new_df['consumable'].map(reversed_consumables_dict) - cost_of_consumables = new_df[~new_df['item_code'].isin(list_of_unique_medical_products)] - cost_of_separately_managed_medical_supplies = new_df[new_df['item_code'].isin(list_of_unique_medical_products)] + # Take explicit copies so that relabelling below modifies these dataframes and not a view of new_df + cost_of_consumables = new_df[~new_df['item_code'].isin(list_of_unique_medical_products)].copy() + cost_of_separately_managed_medical_supplies = new_df[ + new_df['item_code'].isin(list_of_unique_medical_products)].copy() cost_of_separately_managed_medical_supplies['cost_subcategory'] = cost_of_separately_managed_medical_supplies[ 'cost_subcategory'].replace( {'consumables_dispensed': 'separately_managed_medical_supplies_dispensed', @@ -849,24 +990,26 @@ def disaggregate_separately_managed_medical_supplies_from_consumable_costs(_df, columns='item_code') separately_managed_medical_supplies = [127, 141, 161] # Oxygen, Blood, IRS - cost_of_consumables_dispensed, cost_of_separately_managed_medical_supplies_dispensed = disaggregate_separately_managed_medical_supplies_from_consumable_costs( + cost_of_consumables_dispensed, cost_of_separately_managed_medical_supplies_dispensed = ( + disaggregate_separately_managed_medical_supplies_from_consumable_costs( _df=retain_relevant_column_subset( melt_and_label_consumables_cost(cost_of_consumables_dispensed, 'cost_of_consumables_dispensed'), 'consumable'), _consumables_dict=consumables_dict, - list_of_unique_medical_products=separately_managed_medical_supplies) - cost_of_excess_consumables_stocked, cost_of_separately_managed_medical_supplies_excess_stock = disaggregate_separately_managed_medical_supplies_from_consumable_costs( + list_of_unique_medical_products=separately_managed_medical_supplies)) + cost_of_excess_consumables_stocked, cost_of_separately_managed_medical_supplies_excess_stock = ( + disaggregate_separately_managed_medical_supplies_from_consumable_costs( _df=retain_relevant_column_subset( melt_and_label_consumables_cost(cost_of_excess_consumables_stocked, 'cost_of_excess_consumables_stocked'), 'consumable'), _consumables_dict=consumables_dict, - list_of_unique_medical_products=separately_managed_medical_supplies) + list_of_unique_medical_products=separately_managed_medical_supplies)) consumable_costs = pd.concat([cost_of_consumables_dispensed, cost_of_excess_consumables_stocked]) # 2.4 Supply chain costs # --------------------------------------------------------------------------------------------------------------- - # Assume that the cost of procurement, warehousing and distribution is a fixed proportion of consumable purchase costs + # Assume that the cost of procurement,warehousing and distribution is a fixed proportion of consumable purchase # The fixed proportion is based on Resource Mapping Expenditure data from 2018 resource_mapping_data = unit_costs['actual_expenditure_data'] # Make sure values are numeric @@ -922,7 +1065,8 @@ def disaggregate_separately_managed_medical_supplies_from_consumable_costs(_df, # -------------------------------------------- print("Now estimating Medical equipment costs...") - # Total cost of equipment required as per SEL (HSSP-III) only at facility IDs where it has been used in the simulation + # Total cost of equipment required as per SEL (HSSP-III) only at facility IDs where it has been used in the + # simulation # Get list of equipment used in the simulation by district and level def get_equipment_used_by_district_and_facility(_df: pd.Series) -> pd.Series: """Summarise the parsed logged-key results for one draw (as dataframe) into a pd.Series.""" @@ -934,6 +1078,8 @@ def get_equipment_used_by_district_and_facility(_df: pd.Series) -> pd.Series: list_of_equipment_used_by_draw_and_run = extract_results( Path(results_folder), + suspended_results_folder=suspended_results_folder, + scaling_factor=scaling_factor, module='tlo.methods.healthsystem.summary', key='EquipmentEverUsed_ByFacilityID', custom_generate_series=get_equipment_used_by_district_and_facility, @@ -943,62 +1089,68 @@ def get_equipment_used_by_district_and_facility(_df: pd.Series) -> pd.Series: list_of_equipment_used_by_draw_and_run[col] = list_of_equipment_used_by_draw_and_run[col].apply( ast.literal_eval) - # Initialize an empty DataFrame - equipment_cost_across_sim = pd.DataFrame() + # Collect equipment cost for each draw and run, then concatenate once + equipment_cost_frames = [] # Extract equipment cost for each draw and run for d in _draws: for r in _runs: print(f"Processing draw {d} and run {r} of equipment costs") # Extract a list of equipment which was used at each facility level within each district - equipment_used = {district: {level: [] for level in fac_levels} for district in list( - district_dict.values())} # create a dictionary with a key for each district and facility level - list_of_equipment_used_by_current_draw_and_run = list_of_equipment_used_by_draw_and_run[ - (d, r)].reset_index() - for dist in list(district_dict.values()): - for level in fac_levels: - equipment_used_subset = list_of_equipment_used_by_current_draw_and_run[ - (list_of_equipment_used_by_current_draw_and_run['District'] == dist) & ( - list_of_equipment_used_by_current_draw_and_run['Facility_Level'] == level)] - equipment_used_subset.columns = ['District', 'Facility_Level', 'EquipmentEverUsed'] - equipment_used[dist][level] = set().union(*equipment_used_subset['EquipmentEverUsed']) - equipment_used = pd.concat({ - k: pd.DataFrame.from_dict(v, 'index') for k, v in equipment_used.items()}, - axis=0) - full_list_of_equipment_used = set(equipment_used.values.flatten()) - full_list_of_equipment_used = set(filter(pd.notnull, full_list_of_equipment_used)) - - equipment_df = pd.DataFrame() - equipment_df.index = equipment_used.index - for item in full_list_of_equipment_used: - equipment_df[str(item)] = 0 - for dist_fac_index in equipment_df.index: - equipment_df.loc[equipment_df.index == dist_fac_index, str(item)] = equipment_used[ - equipment_used.index == dist_fac_index].isin([item]).any(axis=1) - # equipment_df.to_csv('./outputs/equipment_use.csv') - - equipment_df = equipment_df.reset_index().rename( - columns={'level_0': 'District', 'level_1': 'Facility_Level'}) - equipment_df = pd.melt(equipment_df, id_vars=['District', 'Facility_Level']).rename( - columns={'variable': 'Item_code', 'value': 'whether_item_was_used'}) + equipment_series = list_of_equipment_used_by_draw_and_run[(d, r)] + districts = list(district_dict.values()) + levels = list(fac_levels) + district_set = set(districts) + used_records = [ + (district, level, item) + for (district, level), items in equipment_series.items() + if district in district_set and level in fac_levels + for item in items + if pd.notna(item) + ] + used_membership = pd.DataFrame( + used_records, columns=['District', 'Facility_Level', 'Item_code'] + ).drop_duplicates() + all_items = used_membership['Item_code'].drop_duplicates() + equipment_df = pd.MultiIndex.from_product( + [districts, levels, all_items], + names=['District', 'Facility_Level', 'Item_code'], + ).to_frame(index=False) + equipment_df = equipment_df.merge( + used_membership.assign(whether_item_was_used=True), + on=['District', 'Facility_Level', 'Item_code'], + how='left', validate='1:1', + ) + equipment_df['whether_item_was_used'] = ( + equipment_df['whether_item_was_used'].eq(True) + ) equipment_df['Item_code'] = pd.to_numeric(equipment_df['Item_code']) # Merge the count of facilities by district and level equipment_df = equipment_df.merge(mfl[['District', 'Facility_Level', 'Facility_Count']], on=['District', 'Facility_Level'], how='left') equipment_df.loc[equipment_df.Facility_Count.isna(), 'Facility_Count'] = 0 - # Because levels 1b and 2 are collapsed together, we assume that the same equipment is used by level 1b as that recorded for level 2 + # Because levels 1b and 2 are collapsed together, we assume that the same equipment is used by level 1b as + # that recorded for level 2 def update_itemuse_for_level1b_using_level2_data(_df): - # Create a list of District and Item_code combinations for which use == True - list_of_equipment_used_at_level2 = \ - _df[(_df.Facility_Level == '2') & (_df['whether_item_was_used'] is True)][['District', 'Item_code']] - # Now update the 'whether_item_was_used' for 'Facility_Level' == '1b' to match that of level '2' - _df.loc[ - (_df['Facility_Level'] == '1b') & - (_df[['District', 'Item_code']].apply(tuple, axis=1).isin( - list_of_equipment_used_at_level2.apply(tuple, axis=1))), - 'whether_item_was_used' - ] = True + level2_items = pd.MultiIndex.from_frame( + _df.loc[ + _df["Facility_Level"].eq("2") + & _df["whether_item_was_used"].eq(True), + ["District", "Item_code"], + ] + ) + + row_items = pd.MultiIndex.from_frame( + _df[["District", "Item_code"]] + ) + + copy_use = ( + _df["Facility_Level"].eq("1b") + & row_items.isin(level2_items) + ) + + _df.loc[copy_use, "whether_item_was_used"] = True return _df @@ -1020,19 +1172,12 @@ def update_itemuse_for_level1b_using_level2_data(_df): equipment_cost['Quantity'] * equipment_cost[ 'Facility_Count'] equipment_cost['year'] = max(years) - 1 - if equipment_cost_across_sim.empty: - equipment_cost_across_sim = equipment_cost.groupby(['year', 'Facility_Level', 'Equipment_tlo'])[ - [item + '_annual_total' for item in categories_of_equipment_cost]].sum() - equipment_cost_across_sim['draw'] = d - equipment_cost_across_sim['run'] = r - else: - equipment_cost_for_current_sim = equipment_cost.groupby(['year', 'Facility_Level', 'Equipment_tlo'])[ - [item + '_annual_total' for item in categories_of_equipment_cost]].sum() - equipment_cost_for_current_sim['draw'] = d - equipment_cost_for_current_sim['run'] = r - # Concatenate the results - equipment_cost_across_sim = pd.concat([equipment_cost_across_sim, equipment_cost_for_current_sim], - axis=0) + equipment_cost_for_current_sim = equipment_cost.groupby( + ['year', 'Facility_Level', 'Equipment_tlo'] + )[[item + '_annual_total' for item in categories_of_equipment_cost]].sum() + equipment_cost_frames.append(equipment_cost_for_current_sim.assign(draw=d, run=r)) + + equipment_cost_across_sim = pd.concat(equipment_cost_frames, axis=0) equipment_costs = pd.melt(equipment_cost_across_sim.reset_index(), id_vars=['draw', 'run', 'Facility_Level', 'Equipment_tlo'], # Columns to keep @@ -1043,7 +1188,8 @@ def update_itemuse_for_level1b_using_level2_data(_df): # Assume that the annual costs are constant each year of the simulation equipment_costs = pd.concat([equipment_costs.assign(year=year) for year in years]) - # TODO If the logger is updated to include year, we may wish to calculate equipment costs by year - currently we assume the same annuitised equipment cost each year + # TODO If the logger is updated to include year, we may wish to calculate equipment costs by year + # (currently we assume the same annuitised equipment cost each year) equipment_costs = equipment_costs.reset_index(drop=True) equipment_costs = equipment_costs.rename(columns={'Equipment_tlo': 'Equipment'}) equipment_costs = prepare_cost_dataframe(equipment_costs, _category_specific_group='Equipment', @@ -1137,8 +1283,8 @@ def update_itemuse_for_level1b_using_level2_data(_df): # Define a function to summarize cost data from -# Note that the dataframe needs to have draw as index and run as columns. if the dataframe is long with draw and run as index, then -# first unstack the dataframe and subsequently apply the summarize function +# Note that the dataframe needs to have draw as index and run as columns. if the dataframe is long with draw and run as +# index, then first unstack the dataframe and subsequently apply the summarize function def summarize_cost_data(_df, _metric: Literal['mean', 'median'] = 'mean') -> pd.DataFrame: """ @@ -1194,8 +1340,8 @@ def estimate_projected_health_spending(resourcefilepath: Path, """ Estimate total projected health spending for a simulation period. - Combines health spending per capita projections (Dieleman et al, 2019) with simulated population estimates to calculate - total health expenditure, optionally applying a discount rate and summarizing across runs. + Combines health spending per capita projections (Dieleman et al, 2019) with simulated population estimates to + calculate total health expenditure, optionally applying a discount rate and summarizing across runs. Parameters: ---------- @@ -1424,7 +1570,8 @@ def do_stacked_bar_plot_of_cost_by_category(_df: pd.DataFrame, if (_disaggregate_by_subgroup is True): for name, df in dfs.items(): dfs[name] = df.copy() # Choose the dataframe to modify - # If sub-groups are more than 10 in number, then disaggregate the top 10 and group the rest into an 'other' category + # If sub-groups are more than 10 in number, then disaggregate the top 10 and group the rest into an + # 'other' category if (len(dfs[name]['cost_subgroup'].unique()) > 10): # Calculate total cost per subgroup subgroup_totals = dfs[name].groupby('cost_subgroup')['cost'].sum() @@ -1870,7 +2017,8 @@ def wrap_text(text, width=15): if (len(_df['cost_subgroup'].unique()) > 10): # Step 2: Group all other consumables into "Other" other_cost = _df.iloc[10:]["cost"].sum() - top_10 = pd.concat([top_10, pd.DataFrame([{"cost_subgroup": "Other", "cost": other_cost}])], ignore_index=True) + top_10 = pd.concat([top_10, pd.DataFrame([{"cost_subgroup": "Other", "cost": other_cost}])], + ignore_index=True) # Prepare data for the treemap total_cost = top_10["cost"].sum() @@ -2015,7 +2163,8 @@ def generate_multiple_scenarios_roi_plot(_monetary_value_of_incremental_health: # Initialize an empty DataFrame to store values for each 'run' all_run_values = pd.DataFrame() - # Create an array of implementation costs ranging from 0 to the max value of max ability to pay for the current draw + # Create an array of implementation costs ranging from 0 to the max value of max ability to pay for the current + # draw implementation_costs = np.linspace(0, max_ability_to_pay_for_implementation.loc[draw_index].max(), 50) # Add fixed values for ROI ratio calculation additional_costs = np.array([1_000_000_000, 3_000_000_000]) @@ -2121,7 +2270,8 @@ def generate_multiple_scenarios_roi_plot(_monetary_value_of_incremental_health: # Replace specific x-ticks with % of health spending values if _projected_health_spending: xtick_labels[ - 1] = f'{xticks[1]:,.0f}\n({xticks[1] / (_projected_health_spending / 1e6) :.2%} of \n projected total \n health spend)' + 1] = (f'{xticks[1]:,.0f}\n({xticks[1] / (_projected_health_spending / 1e6) :.2%} of \n projected total ' + f'\n health spend)') for i, tick in enumerate(xticks): if (i != 0) & (i != 1): # Replace for 4000 xtick_labels[i] = f'{tick:,.0f}\n({tick / (_projected_health_spending / 1e6) :.2%})' diff --git a/src/tlo/analysis/utils.py b/src/tlo/analysis/utils.py index 1580a352b1..c6714ded1f 100644 --- a/src/tlo/analysis/utils.py +++ b/src/tlo/analysis/utils.py @@ -295,6 +295,8 @@ def extract_results(results_folder: Path, index: str = None, custom_generate_series=None, do_scaling: bool = False, + suspended_results_folder: Path = None, + scaling_factor: Optional[float] = None, ) -> pd.DataFrame: """Utility function to unpack results. @@ -307,16 +309,23 @@ def extract_results(results_folder: Path, `custom_generate_series`. Optionally, with `do_scaling=True`, each element is multiplied by the scaling_factor recorded in the simulation. + If the suspend-and-resume functionality is used, scaling factor may be available in the folder where the logs + of the suspended run are stored. If the 'scaling_factor' is provided, this value will be used for scaling. + Otherwise, if suspended_results_folder is provided, the scaling_factor will be extracted from the demography log + of this suspended results folder. If neither is provided, the function will attempt to extract this from the + demography log of the (resume) results folder - this contains the scaling_factor only if suspend-and-resume + functionality is NOT used. Note that if runs in the batch have failed (such that logs have not been generated), these are dropped silently. """ - def get_multiplier(_draw, _run): + + def get_multiplier(results_folder, _draw, _run): """Helper function to get the multiplier from the simulation. Note that if the scaling factor cannot be found a `KeyError` is thrown.""" return load_pickled_dataframes( - results_folder, _draw, _run, 'tlo.methods.population' - )['tlo.methods.population']['scaling_factor']['scaling_factor'].values[0] + results_folder, _draw, _run, 'tlo.methods.demography' + )['tlo.methods.demography']['scaling_factor']['scaling_factor'].values[0] if custom_generate_series is None: # If there is no `custom_generate_series` provided, it implies that function required selects the specified @@ -352,7 +361,12 @@ def generate_series(dataframe: pd.DataFrame) -> pd.Series: 'Custom command does not generate a pd.Series' ) if do_scaling: - res[draw_run] = output_from_eval * get_multiplier(draw, run) + if scaling_factor is not None: + res[draw_run] = output_from_eval * scaling_factor + elif suspended_results_folder is not None: + res[draw_run] = output_from_eval * get_multiplier(suspended_results_folder, 0, 0) + else: + res[draw_run] = output_from_eval * get_multiplier(results_folder, draw, run) else: res[draw_run] = output_from_eval @@ -386,7 +400,7 @@ def check_info_value_changes(df): prev_info = row["Info"] return problems - + def remove_events_for_individual_after_death(df): rows_to_drop = []