#! /usr/bin/env python # ========================================================================== # Estimate error for world-wide carbon footprint by bootstraping # # Copyright (C) 2021 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import copy import random import math # ============ # # Set OIR data # # ============ # def set_oir_data(): """ Set OIR data """ oir = [{'obs': 'cfht', 'num': 1, 'footprint': 84989, 'yearly': 2073, 'lifetime': 41}, {'obs': 'eso3.6m', 'num': 1, 'footprint': 79608, 'yearly': 1851, 'lifetime': 43}, {'obs': 'aat', 'num': 1, 'footprint': 205610, 'yearly': 4470, 'lifetime': 46}, {'obs': 'vlt', 'num': 4, 'footprint': 539655, 'yearly': 25698, 'lifetime': 21}, {'obs': 'geminis', 'num': 1, 'footprint': 97280, 'yearly': 4864, 'lifetime': 20}, {'obs': 'gtc', 'num': 1, 'footprint': 29880, 'yearly': 2716, 'lifetime': 11}] # Return return oir # ================== # # Set small OIR data # # ================== # def set_small_oir_data(): """ Set small OIR data """ oir = [{'obs': 'tarot', 'num': 3, 'footprint': 711, 'yearly': 32, 'lifetime': 22}, {'obs': 'pdm', 'num': 1, 'footprint': 1665, 'yearly': 29, 'lifetime': 57}, {'obs': 'nowt', 'num': 1, 'footprint': 415, 'yearly': 49, 'lifetime': 7}, {'obs': 'c2pu', 'num': 2, 'footprint': 780, 'yearly': 98, 'lifetime': 6}, {'obs': 'ogle', 'num': 1, 'footprint': 1604, 'yearly': 70, 'lifetime': 23}, {'obs': 'tillinghast', 'num': 1, 'footprint': 4312, 'yearly': 85, 'lifetime': 51}, {'obs': 'spm', 'num': 1, 'footprint': 4241, 'yearly': 85, 'lifetime': 50}, {'obs': 'kmtnet', 'num': 3, 'footprint': 6377, 'yearly': 856, 'lifetime': 5}, {'obs': 'boao', 'num': 1, 'footprint': 3827, 'yearly': 159, 'lifetime': 24}, {'obs': 'ohp', 'num': 1, 'footprint': 9763, 'yearly': 157, 'lifetime': 62}, {'obs': 'tbl', 'num': 1, 'footprint': 11435, 'yearly': 286, 'lifetime': 40}, {'obs': 'hct', 'num': 1, 'footprint': 4331, 'yearly': 228, 'lifetime': 19}, {'obs': 'xinglong', 'num': 1, 'footprint': 10137, 'yearly': 220, 'lifetime': 46}, {'obs': 'lijiang', 'num': 1, 'footprint': 5168, 'yearly': 431, 'lifetime': 12}] # Return return oir # ============== # # Set Radio data # # ============== # def set_radio_data(): """ Set Radio data """ radio = [{'obs': 'nro', 'num': 1, 'footprint': 26609, 'yearly': 700, 'lifetime': 38}, {'obs': 'noto', 'num': 1, 'footprint': 12096, 'yearly': 378, 'lifetime': 32}, {'obs': 'apex', 'num': 1, 'footprint': 14925, 'yearly': 995, 'lifetime': 15}, {'obs': 'jcmt', 'num': 1, 'footprint': 54194, 'yearly': 1642, 'lifetime': 33}, {'obs': 'gbt', 'num': 1, 'footprint': 75088, 'yearly': 3952, 'lifetime': 19}, {'obs': 'lmt', 'num': 1, 'footprint': 23221, 'yearly': 2637, 'lifetime': 6}] # Return return radio # ==================== # # Set Radio array data # # ==================== # def set_radio_array_data(): """ Set Radio array data """ radio = [{'obs': 'sma', 'num': 1, 'footprint': 14354, 'yearly': 844, 'lifetime': 17}, {'obs': 'iram', 'num': 1, 'footprint': 124740, 'yearly': 4158, 'lifetime': 30}, {'obs': 'lofar', 'num': 1, 'footprint': 66326, 'yearly': 7091, 'lifetime': 8}, {'obs': 'alma', 'num': 1, 'footprint': 535340, 'yearly': 56154, 'lifetime': 9}, {'obs': 'atca', 'num': 1, 'footprint': 50787, 'yearly': 1587, 'lifetime': 32}, {'obs': 'meerkat', 'num': 1, 'footprint': 37004, 'yearly': 6252, 'lifetime': 2}, {'obs': 'vla', 'num': 1, 'footprint': 178826, 'yearly': 4471, 'lifetime': 40}, {'obs': 'vlba', 'num': 1, 'footprint': 136194, 'yearly': 5044, 'lifetime': 27}, {'obs': 'eht', 'num': 1, 'footprint': 12580, 'yearly': 1144, 'lifetime': 11}] # Return return radio # ============== # # Set other data # # ============== # def set_other_data(): """ Set other data """ other = [{'obs': 'themis', 'num': 1, 'footprint': 5775, 'yearly': 275, 'lifetime': 21}, {'obs': 'mlso', 'num': 1, 'footprint': 16817, 'yearly': 306, 'lifetime': 55}, {'obs': 'sofia', 'num': 1, 'footprint': 464919, 'yearly': 48729, 'lifetime': 9}, {'obs': 'hess', 'num': 1, 'footprint': 49126, 'yearly': 2890, 'lifetime': 17}] # Return return other # ========================= # # Set solar mass-based data # # ========================= # def set_solar_mass_data(): """ Set solar mass-based data """ data = [{'obs': 'psp', 'num': 1, 'footprint': 34250, 'yearly': 3425}, {'obs': 'sdo', 'num': 1, 'footprint': 155000, 'yearly': 15500}, {'obs': 'soho', 'num': 1, 'footprint': 92500, 'yearly': 3700}] # Return return data # ========================= # # Set solar cost-based data # # ========================= # def set_solar_cost_data(): """ Set solar cost-based data """ data = [{'obs': 'psp', 'num': 1, 'footprint': 183456, 'yearly': 18346}, {'obs': 'sdo', 'num': 1, 'footprint': 121164, 'yearly': 12116}, {'obs': 'soho', 'num': 1, 'footprint': 205617, 'yearly': 8225}] # Return return data # ========================== # # Set plasma mass-based data # # ========================== # def set_plasma_mass_data(): """ Set plasma mass-based data """ data = [{'obs': 'mms', 'num': 1, 'footprint': 68000, 'yearly': 6800}, {'obs': 'cluster', 'num': 1, 'footprint': 240000, 'yearly': 12000}, {'obs': 'stereo', 'num': 1, 'footprint': 61900, 'yearly': 4421}, {'obs': 'ace', 'num': 1, 'footprint': 37600, 'yearly': 1635}, {'obs': 'geotail', 'num': 1, 'footprint': 50450, 'yearly': 1802}, {'obs': 'wind', 'num': 1, 'footprint': 62500, 'yearly': 2404}, {'obs': 'timed', 'num': 1, 'footprint': 33000, 'yearly': 1833}] # Return return data # ========================== # # Set plasma cost-based data # # ========================== # def set_plasma_cost_data(): """ Set plasma cost-based data """ data = [{'obs': 'mms', 'num': 1, 'footprint': 147501, 'yearly': 14750}, {'obs': 'cluster', 'num': 1, 'footprint': 132207, 'yearly': 6610}, {'obs': 'stereo', 'num': 1, 'footprint': 86021, 'yearly': 6144}, {'obs': 'timed', 'num': 1, 'footprint': 36196, 'yearly': 2011}] # Return return data # ============================ # # Set plantary mass-based data # # ============================ # def set_planetary_mass_data(): """ Set planetary mass-based data """ data = [{'obs': 'maven', 'num': 1, 'footprint': 122700, 'yearly': 12270}, {'obs': 'mex', 'num': 1, 'footprint': 61150, 'yearly': 3597}, {'obs': 'insight', 'num': 1, 'footprint': 34700, 'yearly': 3470}, {'obs': 'msl', 'num': 1, 'footprint': 194650, 'yearly': 19465}, {'obs': 'mro', 'num': 1, 'footprint': 109000, 'yearly': 7786}, {'obs': 'juno', 'num': 1, 'footprint': 181250, 'yearly': 18125}] # Return return data # ============================ # # Set plantary cost-based data # # ============================ # def set_planetary_cost_data(): """ Set planetary cost-based data """ data = [{'obs': 'maven', 'num': 1, 'footprint': 89270, 'yearly': 8927}, {'obs': 'mex', 'num': 1, 'footprint': 52332, 'yearly': 3078}, {'obs': 'insight', 'num': 1, 'footprint': 99922, 'yearly': 9992}, {'obs': 'msl', 'num': 1, 'footprint': 362595, 'yearly': 36259}, {'obs': 'mro', 'num': 1, 'footprint': 129850, 'yearly': 9275}, {'obs': 'juno', 'num': 1, 'footprint': 151547, 'yearly': 15155}] # Return return data # ========================= # # Set astro mass-based data # # ========================= # def set_astro_mass_data(): """ Set astro mass-based data """ data = [{'obs': 'tess', 'num': 1, 'footprint': 16250, 'yearly': 1625}, {'obs': 'gaia', 'num': 1, 'footprint': 101700, 'yearly': 10170}, {'obs': 'wise', 'num': 1, 'footprint': 33050, 'yearly': 3005}, {'obs': 'hst', 'num': 1, 'footprint': 555500, 'yearly': 18517}, {'obs': 'xmm', 'num': 1, 'footprint': 190000, 'yearly': 9048}, {'obs': 'chandra', 'num': 1, 'footprint': 293000, 'yearly': 13952}, {'obs': 'nicer', 'num': 1, 'footprint': 18600, 'yearly': 1860}, {'obs': 'nustar', 'num': 1, 'footprint': 18000, 'yearly': 1800}, {'obs': 'astrosat', 'num': 1, 'footprint': 75750, 'yearly': 7575}, {'obs': 'swift', 'num': 1, 'footprint': 42150, 'yearly': 2634}, {'obs': 'integral', 'num': 1, 'footprint': 200000, 'yearly': 11111}, {'obs': 'fermi', 'num': 1, 'footprint': 215150, 'yearly': 17929}] # Return return data # ========================= # # Set astro cost-based data # # ========================= # def set_astro_cost_data(): """ Set astro cost-based data """ data = [{'obs': 'tess', 'num': 1, 'footprint': 38478, 'yearly': 3848}, {'obs': 'gaia', 'num': 1, 'footprint': 145114, 'yearly': 14511}, {'obs': 'wise', 'num': 1, 'footprint': 46855, 'yearly': 4260}, {'obs': 'hst', 'num': 1, 'footprint': 1125197, 'yearly': 37507}, {'obs': 'xmm', 'num': 1, 'footprint': 155845, 'yearly': 7421}, {'obs': 'chandra', 'num': 1, 'footprint': 575955, 'yearly': 27426}, {'obs': 'nicer', 'num': 1, 'footprint': 7374, 'yearly': 737}, {'obs': 'nustar', 'num': 1, 'footprint': 21799, 'yearly': 2180}, {'obs': 'astrosat', 'num': 1, 'footprint': 3751, 'yearly': 375}, {'obs': 'swift', 'num': 1, 'footprint': 39030, 'yearly': 2439}, {'obs': 'integral', 'num': 1, 'footprint': 58720, 'yearly': 3262}, {'obs': 'fermi', 'num': 1, 'footprint': 120881, 'yearly': 10073}] # Return return data # ============ # # Flatten data # # ============ # def flatten_data(data): """ Flatten data """ # Initialise flattened data flattened_data = [] # Append data for d in data: # Get attributes name = d['obs'] num = d['num'] footprint = d['footprint'] yearly = d['yearly'] # Append observatory for i in range(num): if num > 1: d['obs'] = '%s-%d' % (name, i+1) d['num'] = 1 d['footprint'] = float(footprint) / float(num) d['yearly'] = float(yearly) / float(num) flattened_data.append(copy.deepcopy(d)) # Return return flattened_data # ============== # # Bootstrap data # # ============== # def bootstrap_data(category, data, ntotal, nsamples=1000, fix_largest=False, uncfrac=0.8): """ Bootstrap data """ # Get number of observatories num = len(data) # Determine observatory with largest yearly footprint largest_footprint = 0.0 largest_yearly = 0.0 largest_obs = '' for d in data: footprint = float(d['footprint']) / float(d['num']) yearly = float(d['yearly']) / float(d['num']) if yearly > largest_yearly: largest_footprint = footprint largest_yearly = yearly largest_obs = d['obs'] # Initialise mean and rms mean_total = 0.0 mean_total_unc = 0.0 mean_yearly = 0.0 mean_yearly_unc = 0.0 var_total = 0.0 var_yearly = 0.0 # Loop over all samples for i in range(nsamples): # Initialise sums sum_total = 0.0 unc_total = 0.0 sum_yearly = 0.0 unc_yearly = 0.0 # If infrastructure with largest yearly write-off should be fixed then add it # now to sum if fix_largest: sum_total += largest_footprint unc_total += uncfrac*uncfrac*largest_footprint*largest_footprint sum_yearly += largest_yearly unc_yearly += uncfrac*uncfrac*largest_yearly*largest_yearly ninfra = ntotal - 1 else: ninfra = ntotal # Bootstrap ninfra infrastructures for k in range(ninfra): # Get index of infrastructure if fix_largest: inx = -1 while inx < 0: inx = random.randrange(0, num) if data[inx]['obs'] == largest_obs: inx = -1 continue else: inx = random.randrange(0, num) # Add total and yearly footprint sum_total += float(data[inx]['footprint']) unc_total += uncfrac*uncfrac*float(data[inx]['footprint'])*float(data[inx]['footprint']) sum_yearly += float(data[inx]['yearly']) unc_yearly += uncfrac*uncfrac*float(data[inx]['yearly'])*float(data[inx]['yearly']) # Compute means and variance mean_total += sum_total mean_total_unc += math.sqrt(unc_total) mean_yearly += sum_yearly mean_yearly_unc += math.sqrt(unc_yearly) var_total += sum_total * sum_total var_yearly += sum_yearly * sum_yearly # Compute mean and standard deviation var_total = (var_total - mean_total * mean_total / float(nsamples)) / float(nsamples) var_yearly = (var_yearly - mean_yearly * mean_yearly / float(nsamples)) / float(nsamples) mean_total /= float(nsamples) mean_yearly /= float(nsamples) mean_total_unc /= float(nsamples) mean_yearly_unc /= float(nsamples) # Compute total error (bootstrap and FE-related) var_total += mean_total_unc*mean_total_unc var_yearly += mean_yearly_unc*mean_yearly_unc # Package results results = {'mean_total': mean_total, 'mean_total_unc': mean_total_unc, 'var_total': var_total, 'mean_yearly': mean_yearly, 'mean_yearly_unc': mean_yearly_unc, 'var_yearly': var_yearly} # Print result print('%-12s [%2d->%4d]: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year' % (category, num, ntotal, mean_total/1000.0, math.sqrt(var_total)/1000.0, mean_yearly/1000.0, math.sqrt(var_yearly)/1000.0)) # Return return results # ========= # # Bootstrap # # ========= # def bootstrap(nsamples=1000, fix_largest=False): """ Bootstrap """ # Set mass-based space missions mass_solar = flatten_data(set_solar_mass_data()) mass_plasma = flatten_data(set_plasma_mass_data()) mass_planetary = flatten_data(set_planetary_mass_data()) mass_astro = flatten_data(set_astro_mass_data()) # Set cost-based space missions cost_solar = flatten_data(set_solar_cost_data()) cost_plasma = flatten_data(set_plasma_cost_data()) cost_planetary = flatten_data(set_planetary_cost_data()) cost_astro = flatten_data(set_astro_cost_data()) # Set ground-based data oir = flatten_data(set_oir_data()) oir_small = flatten_data(set_small_oir_data()) radio = flatten_data(set_radio_data()) radio_array = flatten_data(set_radio_array_data()) other = flatten_data(set_other_data()) # Bootstrap space mass-based estimates print(' ************* Space missions (mass-based) *************') mass_solar_results = bootstrap_data('Solar', mass_solar, 3, nsamples=nsamples, fix_largest=fix_largest) mass_plasma_results = bootstrap_data('Plasma', mass_plasma, 13, nsamples=nsamples, fix_largest=fix_largest) mass_planetary_results = bootstrap_data('Planetary', mass_planetary, 21, nsamples=nsamples, fix_largest=fix_largest) mass_astro_results = bootstrap_data('Astro', mass_astro, 18, nsamples=nsamples, fix_largest=fix_largest) # Compute mass-based footprint of space missions value_mass_total = (mass_solar_results['mean_total']+ mass_plasma_results['mean_total']+ mass_planetary_results['mean_total']+ mass_astro_results['mean_total']) var_mass_total = (mass_solar_results['var_total']+ mass_plasma_results['var_total']+ mass_planetary_results['var_total']+ mass_astro_results['var_total']) value_mass_yearly = (mass_solar_results['mean_yearly']+ mass_plasma_results['mean_yearly']+ mass_planetary_results['mean_yearly']+ mass_astro_results['mean_yearly']) var_mass_yearly = (mass_solar_results['var_yearly']+ mass_plasma_results['var_yearly']+ mass_planetary_results['var_yearly']+ mass_astro_results['var_yearly']) # Print total print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year' % ('Total', value_mass_total/1000.0, math.sqrt(var_mass_total)/1000.0, value_mass_yearly/1000.0, math.sqrt(var_mass_yearly)/1000.0)) # Bootstrap space cost-based estimates print('') print(' ************* Space missions (cost-based) *************') cost_solar_results = bootstrap_data('Solar', cost_solar, 3, nsamples=nsamples, fix_largest=fix_largest) cost_plasma_results = bootstrap_data('Plasma', cost_plasma, 13, nsamples=nsamples, fix_largest=fix_largest) cost_planetary_results = bootstrap_data('Planetary', cost_planetary, 21, nsamples=nsamples, fix_largest=fix_largest) cost_astro_results = bootstrap_data('Astro', cost_astro, 18, nsamples=nsamples, fix_largest=fix_largest) # Compute cost-based footprint of space missions value_cost_total = (cost_solar_results['mean_total']+ cost_plasma_results['mean_total']+ cost_planetary_results['mean_total']+ cost_astro_results['mean_total']) var_cost_total = (cost_solar_results['var_total']+ cost_plasma_results['var_total']+ cost_planetary_results['var_total']+ cost_astro_results['var_total']) value_cost_yearly = (cost_solar_results['mean_yearly']+ cost_plasma_results['mean_yearly']+ cost_planetary_results['mean_yearly']+ cost_astro_results['mean_yearly']) var_cost_yearly = (cost_solar_results['var_yearly']+ cost_plasma_results['var_yearly']+ cost_planetary_results['var_yearly']+ cost_astro_results['var_yearly']) # Print total print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year' % ('Total', value_cost_total/1000.0, math.sqrt(var_cost_total)/1000.0, value_cost_yearly/1000.0, math.sqrt(var_cost_yearly)/1000.0)) # Bootstrap ground-based estimates print('') print(' ************* Ground-based observatories *************') oir_results = bootstrap_data('OIR (>3m)', oir, 37, nsamples=nsamples, fix_largest=fix_largest) oir_small_results = bootstrap_data('OIR (<3m)', oir_small, 1000, nsamples=nsamples, fix_largest=fix_largest) radio_results = bootstrap_data('Radio', radio, 74, nsamples=nsamples, fix_largest=fix_largest) radio_array_results = bootstrap_data('Radio Arrays', radio_array, 27, nsamples=nsamples, fix_largest=fix_largest) other_results = bootstrap_data('Others', other, 4, nsamples=nsamples, fix_largest=fix_largest) # Compute footprint of ground-based observatories value_ground_total = (oir_results['mean_total'] + oir_small_results['mean_total'] + radio_results['mean_total'] + radio_array_results['mean_total'] + other_results['mean_total']) var_ground_total = (oir_results['var_total'] + oir_small_results['var_total'] + radio_results['var_total'] + radio_array_results['var_total'] + other_results['var_total']) value_ground_yearly = (oir_results['mean_yearly'] + oir_small_results['mean_yearly'] + radio_results['mean_yearly'] + radio_array_results['mean_yearly'] + other_results['mean_yearly']) var_ground_yearly = (oir_results['var_yearly'] + oir_small_results['var_yearly'] + radio_results['var_yearly'] + radio_array_results['var_yearly'] + other_results['var_yearly']) # Print total print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year' % ('Total', value_ground_total/1000.0, math.sqrt(var_ground_total)/1000.0, value_ground_yearly/1000.0, math.sqrt(var_ground_yearly)/1000.0)) # Compute total for mass-based space estimates total_value_mass = value_mass_total + value_ground_total total_var_mass = var_mass_total + var_ground_total yearly_value_mass = value_mass_yearly + value_ground_yearly yearly_var_mass = var_mass_yearly + var_ground_yearly # Compute total for cost-based space estimates total_value_cost = value_cost_total + value_ground_total total_var_cost = var_cost_total + var_ground_total yearly_value_cost = value_cost_yearly + value_ground_yearly yearly_var_cost = var_cost_yearly + var_ground_yearly # Compute average estimates total_value_average = 0.5*(total_value_mass + total_value_cost) total_var_average = 0.5*(total_var_mass + total_var_cost) yearly_value_average = 0.5*(yearly_value_mass + yearly_value_cost) yearly_var_average = 0.5*(yearly_var_mass + yearly_var_cost) # Print total results print('') print(' *************** Total carbon footprint ***************') #print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % # ('Total (mass-based)', # total_value_mass/1000.0, math.sqrt(total_var_mass)/1000.0, # yearly_value_mass/1000.0, math.sqrt(yearly_var_mass)/1000.0, # yearly_value_mass/30000.0, math.sqrt(yearly_var_mass)/30000.0)) #print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % # ('Total (cost-based)', # total_value_cost/1000.0, math.sqrt(total_var_cost)/1000.0, # yearly_value_cost/1000.0, math.sqrt(yearly_var_cost)/1000.0, # yearly_value_cost/30000.0, math.sqrt(yearly_var_cost)/30000.0)) #print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % # ('Total (average)', # total_value_average/1000.0, math.sqrt(total_var_average)/1000.0, # yearly_value_average/1000.0, math.sqrt(yearly_var_average)/1000.0, # yearly_value_average/30000.0, math.sqrt(yearly_var_average)/30000.0)) # Compute space-based estimates total_space = 0.5 * (value_mass_total + value_cost_total) yearly_space = 0.5 * (value_mass_yearly + value_cost_yearly) total_space_error = 0.5 * math.sqrt(var_mass_total+var_cost_total) + 0.5*abs(value_mass_total - value_cost_total) yearly_space_error = 0.5 * math.sqrt(yearly_var_mass+yearly_var_cost) + 0.5*abs(value_mass_yearly - value_cost_yearly) # Print space mission and ground-based observatories print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % ('Space missions', total_space/1000.0, total_space_error/1000.0, yearly_space/1000.0, yearly_space_error/1000.0, yearly_space/30000.0, yearly_space_error/30000.0)) print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % ('Ground-based', value_ground_total/1000.0, math.sqrt(var_ground_total)/1000.0, value_ground_yearly/1000.0, math.sqrt(var_ground_yearly)/1000.0, value_ground_yearly/30000.0, math.sqrt(var_ground_yearly)/30000.0)) # Add half the difference between mass-based and cost-based to variance total_error_global = math.sqrt(total_var_average) + 0.5*abs(total_value_mass - total_value_cost) yearly_error_global = math.sqrt(yearly_var_average) + 0.5*abs(yearly_value_mass - yearly_value_cost) # Print total print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % ('Total', total_value_average/1000.0, total_error_global/1000.0, yearly_value_average/1000.0, yearly_error_global/1000.0, yearly_value_average/30000.0, yearly_error_global/30000.0)) # Put results in dictionary results = {'total': total_value_average, 'err_total': total_error_global, 'annual': yearly_value_average, 'err_annual': yearly_error_global} # Return results return results # =========== # # Entry point # # =========== # if __name__ == '__main__': # Initialise random number generator seed random.seed(0) # Perform bootstraping full = bootstrap() print('') print('') limited = bootstrap(fix_largest=True) # Compute final result total = 0.5 * (full['total'] + limited['total']) err_total = 0.5 * math.sqrt(full['err_total']*full['err_total'] + limited['err_total']*limited['err_total']) + \ 0.5 * abs(full['total'] - limited['total']) annual = 0.5 * (full['annual'] + limited['annual']) err_annual = 0.5 * math.sqrt(full['err_annual']*full['err_annual'] + limited['err_annual']*limited['err_annual']) + \ 0.5 * abs(full['annual'] - limited['annual']) # Print global print('') print('') print(' ********************** Synhesis **********************') print('%-23s: %5.0f +/- %4.0f kt CO2e %4.0f +/- %3.0f kt CO2e / year %4.1f +/- %3.1f t CO2e / year / astronomer' % ('Total', total/1000.0, err_total/1000.0, annual/1000.0, err_annual/1000.0, annual/30000.0, err_annual/30000.0))