From 6de068317dbe908415836fbe21987f936d811cd0 Mon Sep 17 00:00:00 2001
From: Andreas Larsen <44347092+LarsenAndreas@users.noreply.github.com>
Date: Fri, 10 Jan 2025 14:45:32 +0100
Subject: [PATCH 1/3] Rewrite and Restructure
---
.gitignore | 4 +
EOS.py | 1211 -
EOS/EOS.py | 335 +
EOS/LPP.py | 111 +
__init__.py => EOS/__init__.py | 0
EOS/constraints/__init__.py | 0
EOS/constraints/capacity.py | 9 +
EOS/constraints/misc.py | 55 +
EOS/constraints/stereo.py | 53 +
EOS/constraints/strips.py | 23 +
EOS/data/city_locations.xlsx | Bin 0 -> 1443149 bytes
EOS/dataGenerator.py | 337 +
EOS/docs/README.md | 316 +
cloud_dist.png => EOS/docs/cloud_dist.png | Bin
cloud_gen.png => EOS/docs/cloud_gen.png | Bin
EOS/docs/generated_requests_start.png | Bin 0 -> 58011 bytes
.../docs/single_scenario_map.PNG | Bin
EOS/scoring.py | 88 +
EOS/solvers/__init__.py | 0
EOS/solvers/solver_cplex.py | 11 +
EOS/solvers/solver_dag.py | 196 +
EOS/solvers/solver_glpk.py | 28 +
EOS/solvers/solver_gurobi.py | 59 +
EOS/solvers/solver_pulp.py | 63 +
EOS/solvers/solver_random.py | 37 +
EOS/solvers/solver_vns.py | 63 +
EOS/tests/default.py | 0
EOS/utility/__init__.py | 0
EOS/utility/api.py | 160 +
EOS/utility/plotting.py | 152 +
EOS/utility/utils.py | 164 +
EOS_example.html | 135429 ---------------
ID_create.py | 11 -
LPP_data_setup.py | 784 -
README.md | 298 -
all_requests.html | 135411 --------------
construct_performance_df.py | 267 -
distance_matrix.py | 41 -
easy_funcs.py | 180 -
electre.py | 274 -
electre_parallel.py | 79 -
gen_locoations.py | 75 -
generated_requests_start.png | Bin 136563 -> 0 bytes
plot_requests.py | 111 -
sat_path.html | 126304 --------------
schedule_rel_criteria.py | 314 -
scipy_sparce_to_spmatrix.py | 5 -
worldloc.xlsx | 100001 -----------
48 files changed, 2264 insertions(+), 500795 deletions(-)
delete mode 100644 EOS.py
create mode 100644 EOS/EOS.py
create mode 100644 EOS/LPP.py
rename __init__.py => EOS/__init__.py (100%)
create mode 100644 EOS/constraints/__init__.py
create mode 100644 EOS/constraints/capacity.py
create mode 100644 EOS/constraints/misc.py
create mode 100644 EOS/constraints/stereo.py
create mode 100644 EOS/constraints/strips.py
create mode 100644 EOS/data/city_locations.xlsx
create mode 100644 EOS/dataGenerator.py
create mode 100644 EOS/docs/README.md
rename cloud_dist.png => EOS/docs/cloud_dist.png (100%)
rename cloud_gen.png => EOS/docs/cloud_gen.png (100%)
create mode 100644 EOS/docs/generated_requests_start.png
rename single_scenario_map.PNG => EOS/docs/single_scenario_map.PNG (100%)
create mode 100644 EOS/scoring.py
create mode 100644 EOS/solvers/__init__.py
create mode 100644 EOS/solvers/solver_cplex.py
create mode 100644 EOS/solvers/solver_dag.py
create mode 100644 EOS/solvers/solver_glpk.py
create mode 100644 EOS/solvers/solver_gurobi.py
create mode 100644 EOS/solvers/solver_pulp.py
create mode 100644 EOS/solvers/solver_random.py
create mode 100644 EOS/solvers/solver_vns.py
create mode 100644 EOS/tests/default.py
create mode 100644 EOS/utility/__init__.py
create mode 100644 EOS/utility/api.py
create mode 100644 EOS/utility/plotting.py
create mode 100644 EOS/utility/utils.py
delete mode 100644 EOS_example.html
delete mode 100644 ID_create.py
delete mode 100644 LPP_data_setup.py
delete mode 100644 README.md
delete mode 100644 all_requests.html
delete mode 100644 construct_performance_df.py
delete mode 100644 distance_matrix.py
delete mode 100644 easy_funcs.py
delete mode 100644 electre.py
delete mode 100644 electre_parallel.py
delete mode 100644 gen_locoations.py
delete mode 100644 generated_requests_start.png
delete mode 100644 plot_requests.py
delete mode 100644 sat_path.html
delete mode 100644 schedule_rel_criteria.py
delete mode 100644 scipy_sparce_to_spmatrix.py
delete mode 100644 worldloc.xlsx
diff --git a/.gitignore b/.gitignore
index b6e4761..7aa385e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -127,3 +127,7 @@ dmypy.json
# Pyre type checker
.pyre/
+
+
+# Custom
+*.html
\ No newline at end of file
diff --git a/EOS.py b/EOS.py
deleted file mode 100644
index 7953a29..0000000
--- a/EOS.py
+++ /dev/null
@@ -1,1211 +0,0 @@
-#overall functions
-
-def scenario(seconds_gran = 20, number_of_requests_0 = 1000, NORAD_ids = [38755, 40053], weather_real = False, simplify = True,
- schedule_start = [2021,7,21,9,40],
- hours_horizon = 8,
- max_off_nadir_angle = 30, #degrees
- height_satellite = 694, #km
- rotation_speed = 30/12, #degrees per second #per https://directory.eoportal.org/web/eoportal/satellite-missions/s/spot-6-7
- cam_resolution = 1, #m^2 per pixel
- capacity_limit = 1000000, #in mega byte
- satellite_swath = 3600,
- map_generation=True): #compute maps?
- # NORAD_ids = [38012, 39019]
- # seconds_gran =10
- # number_of_requests_0 =30
- # weather_real = False
- ## entire setup ##
- import datetime
- #from satellite_tle import fetch_tle_from_celestrak
- import requests
-
-
-
- TLEs = list()
- for i in range(0, len(NORAD_ids)):
- #tle = requests.get('https://www.celestrak.com/satcat/tle.php?CATNR={}'.format(NORAD_ids[i]), #old version
- tle = requests.get('https://www.celestrak.com/NORAD/elements/gp.php?CATNR={}&FORMAT=TLE'.format(NORAD_ids[i]), #new version
- verify=True,
- timeout=20)
- TLE = tle.text.split('\r\n')
- TLEs.append(TLE)
- #NORAD_ids = [38012, 39019, 38755, 40053]
-
- ################### SCHEDULE INPUT ##################
- #starting time for horizon
- import random
- random.seed(42)
-
- #time_stochastic = random.uniform(0,24)
- start_schedule = datetime.datetime(schedule_start[0],schedule_start[1],schedule_start[2],schedule_start[3],schedule_start[4]) #- datetime.timedelta(hours=int(time_stochastic))
- #start_schedule = datetime.datetime(2020,5,8,12,50)
- #datetime.datetime.utcnow() - datetime.timedelta(days=total_days+day) - datetime.timedelta(hours = 18)
- #hours ahead to schedule
- hours_ahead = hours_horizon
- print(start_schedule, 'to', start_schedule + datetime.timedelta(hours=int(hours_ahead)))
- #Note, Dates always use Universal Time, NOT local time
- #granularity of scheduling (discretization),
- #i.e. time segmentation
- #number of satellite
- number_of_satellites = len(NORAD_ids)
- #compute maps?
- map_generation=True
- ####################################################
-
- ################## SATELLITE INPUT #################
- #TLE of satellite 1
- # =============================================================================
- # name_1 = "pleiades 1A";
- # line1_1 = "1 38012U 11076F 19303.06642174 -.00000056 00000-0 -23405-5 0 9997";
- # line2_1 = "2 38012 98.1883 15.9149 0001494 90.8918 269.2452 14.58560404418937"
- #
- # #TLE of satellite 2
- # name_2 = "pleiades 1B";
- # line1_2 = "1 39019U 12068A 19303.10054518 -.00000046 00000-0 -91078-7 0 9993";
- # line2_2 = "2 39019 98.1870 15.9140 0001220 88.0086 272.1258 14.58558050367792"
- #
- # #TLE of satellite 3
- # name_3 = "SPOT 6";
- # line1_3 = "1 38755U 12047A 19303.08272605 -.00000060 00000-0 -30989-5 0 9990";
- # line2_3 = "2 38755 98.1980 8.4889 0000909 81.1977 278.9329 14.58563161380015"
- #
- # #TLE of satellite 4
- # name_4 = "SPOT 7";
- # line1_4 = "1 40053U 14034A 19303.11654248 .00001178 00000-0 26330-3 0 9995";
- # line2_4 = "2 40053 98.2135 8.3331 0001054 68.4868 291.6410 14.58559748283956"
- # =============================================================================
-
- #SATELLITE SPECS assumed to be the same
- #no energy constraint - ease computation
- #satellite cameras reachability
- #max_off_nadir_angle = 30 #degrees
- #height_satellite = 694 #km
- #rotation_speed = 30/12 #degrees per second #per https://directory.eoportal.org/web/eoportal/satellite-missions/s/spot-6-7
- #cam_resolution = 1 #m^2 per pixel
- #capacity_limit = 1000000 #in mega byte
- #satellite_swath = 3600
- ########################################################
- #allowed_weather = 50
- #succesful_weather = 10
-
- ############# DATA GENERATION INPUT ####################
- #number of requests from customers through # days!
- number_of_requests = 0
- #request already in database
- #number_of_requests_0 = 0
- ########################################################
-
- #packages
- import ephem
- from math import degrees, floor
- import folium
- import pandas as pd
- import numpy as np
- np.random.seed(42)
-
- #import timeit
- #data_time_start = timeit.default_timer()
-
- ################### DATA GENERATION ####################
- #### schedule relative data ###
- from EOSpython.schedule_rel_criteria import schedule_criteria
- day=1
- total_days = 1
- df = schedule_criteria(number_of_requests, total_days, number_of_requests_0, satellite_swath)
- #df.info()
- avg_pri_df = [[list(np.mean(df[df["priority"] == j], axis = 0))[i] for i in [3,4,5,6,7,8,11]] for j in list(range(1,8))]
- #np_avg_pri_df = np.array(avg_pri_df)
- col_names_avg = list(df.columns[[4,5,6,7,8,9,12]])
- df_avg_pri = pd.DataFrame(
- {"measures": col_names_avg,
- "pri 1": avg_pri_df[0],
- "pri 2": avg_pri_df[1],
- "pri 3": avg_pri_df[2],
- "pri 4": avg_pri_df[3],
- "pri 5": avg_pri_df[4],
- "pri 6": avg_pri_df[5],
- "pri 7": avg_pri_df[6],
- }
- )
- df_avg_pri.iloc[:,:4]
- df_avg_pri.iloc[:,4:]
-
-
- #### prelimenary data analysis begins - COLLECT VALID DATA
- #changing waiting time accordingly
- DF_i = df[(df["day"] <= day) & (df["acquired"] == 0)]
- pd.options.mode.chained_assignment = None
-
- added_wait = list()
- for i in range(0,DF_i.shape[0]):
- added_wait.append(random.randint(1,14)) #age distribution randomly uniform 1 to 14!!!!!!
- DF_i["waiting time"] = day - DF_i["day"] + np.array(added_wait)
- #stereo requests into multiple requests
- #DF_i["stereo"]>1
-
-
- #data_time_end_generation = data_time_start - timeit.default_timer()
-
-
- ##Start map generation
- if map_generation == True:
- m = folium.Map(location=[20, 0], zoom_start=2) #tiles= 'Cartodb Positron'
- from EOSpython.plot_requests import plot_requests, plot_requests2
- plot_requests(m, df = DF_i, name = "request location", radius = 5)
- #file is called all_requests.html
-
-
-
-
-
- ### satellite path computation
-
- sat_names = list()
- tle_rec = list()
- for i in range(0, len(TLEs)):
- sat_names.append(TLEs[i][0])
- tle_rec.append(ephem.readtle(TLEs[i][0], TLEs[i][1], TLEs[i][2]))
-
- increment = datetime.timedelta(seconds=seconds_gran)
- number_of_acq_points = floor((hours_ahead*60*60)/seconds_gran) #convert hours to seconds
- time_slots = list()
- location_slots = list()
- for i in range(0, len(NORAD_ids)):
- location_slots.append(list())
- for i in range(0,number_of_acq_points):
- time = start_schedule + (i)*increment
- time_slots.append(time)
- for k in range(0, number_of_satellites):
- tle_rec[k].compute(time)
- location_slots[k].append([degrees(tle_rec[k].sublat), degrees(tle_rec[k].sublong)])
-
- ### DATA COMPUTATION (satellite relative)
- from EOSpython.distance_matrix import distance_matrix
- distance = distance_matrix(location_slots, DF_i, max_off_nadir_angle, height_satellite, number_of_satellites)
- print("number of reachable attempts:", np.sum(~np.isnan(distance)))
-
- #data_time_end_path = data_time_start - timeit.default_timer()
-
- if weather_real == True:
- weather = True
- generate_weather = False
- else:
- weather = False
- generate_weather = True
-
- from EOSpython.construct_performance_df import construct_performance_df
- performance_df = construct_performance_df(DF_i, seconds_gran, location_slots, time_slots,
- distance, height_satellite, hours_ahead,
- weather = weather, generate_weather = generate_weather)
- print("number of attempts within thresholds:", performance_df.shape[0])
- print("number of requests", len(np.unique(performance_df['ID'])))
- #plot satellite path
- if map_generation == True:
- sat = folium.FeatureGroup(name = "Satellite path")
- colors = list(["royalblue", "black", "darkolivegreen", "green"])
- for k in range(0,len(location_slots)):
- for i in range(0,len(location_slots[k])):
- popup_text = "Satellite: {}
time: {}
possible acquisitions:
{}"
- popup_text = popup_text.format(sat_names[k],
- time_slots[i].strftime('%H:%M:%S-%m/%d/%Y'),
- list(np.where(performance_df["time"] == time_slots[i])[0]))
- sat.add_child(folium.CircleMarker(location=location_slots[k][i], radius = 4, opacity = 1-(i/(1.3*number_of_acq_points)),popup=folium.Popup(popup_text), color = colors[k]))
- m.add_child(sat)
- #m.add_child(folium.LayerControl())
- m.save("sat_path.html")
-
- if map_generation == True:
- plot_requests2(m, df = performance_df, name = "request location")
- #file is called all_requests.html
-
- #data_time_end_preprocessing = data_time_start - timeit.default_timer()
-
-
- #performance_df.info()
- ###################################################################
- ########################### LPP SETUP #############################
- ###################################################################
-
- from EOSpython.LPP_data_setup import LPP_data_multi
- LPP = LPP_data_multi(performance_df, number_of_acq_points, time_slots, location_slots,
- height_satellite, rotation_speed, seconds_gran,
- capacity_limit, cam_resolution, simplify)
- #LPP.LHS
- #LPP.RHS
- #LPP.eLHS
- #LPP.eRHS
-
- #data_time_end_LPP = data_time_start - timeit.default_timer()
- #start = timeit.default_timer()
- ##############################################################
- ##### Stochastic multicriteria acceptability analysis ########
- ##############################################################
-
-
- #1# create different weights
-
- #2# get electre iii score for all sets of w
- #2# get topsis score for all sets of w
- # take average for all sets of scores
-
- #3# insert score in LPP an find solution
- #3# analyse relationship between weight input and solution
-
- # how often is each alternative chosen
- # when alternative is chosen what is the average weight distribution?
- # based on certain weight dis input - what is closest schedule(s)
-
- scenario.LPP = LPP
- scenario.df = DF_i
- scenario.pf_df = LPP.performance_df
- scenario.m = m
- return(EOSscenario)
-
-
-
-
-
-
-
-
-def solve(x_data, scoring_method=2, solution_method="DAG",
- criteria_weights = np.array([0,0,0,0,0,0,1,0]),
- threshold_parameters= np.array([[0,0,1000],
- [0,0,40],
- [0,0,40],
- [0,0,15],
- [0,0,4],
- [0,0,20000],
- [0,0,1],
- [0,0,1]]),
- alpha = 1):
- #import ephem
- #from math import degrees, floor
- #import folium
- import pandas as pd
- import numpy as np
- import time
- #performance_df = x_data.pf_df
- #performance_df.info()
- #performance_df.columns[[7,11,12,13,15,16,17,18,19]]
-
- #### Call generated data frames
- #performance_df.to_csv(r'C:\Users\allex\Desktop\SPECIALE\data\10reach.csv', index = ";", header=True)
- #read diff data
- #reach10_data = pd.read_csv(r'C:\Users\allex\Desktop\SPECIALE\data\10reach.csv', sep = ",")
- #performance_df = reach10_data.iloc[:,1:]
-
- #naming convention
- LPP = x_data.LPP
- DF_i = x_data.df
- performance_df = x_data.pf_df
-
- #IDENTIFY CRITERIA TO INCLUDE IN SCORING PROCEDURE
- dat = np.array(performance_df.iloc[:,[7,11,12,13,15,17,18,19]].transpose()) #which criteria is important in performance df
- #NOTE FOR PAPER 3 - 16th col is not included as criteria
-
- #from fuzzy_topsis import fuzzytopsis
- from EOSpython.easy_funcs import topsis
- from EOSpython.electre_parallel import parallelectre
- from cvxopt import matrix
- from cvxopt import glpk
- #glpk.options["show_progress"] = True
- #glpk.options["maxiters"] = 1000
-
- binVars = range(dat.shape[1])
-
- ###### setup test environment #######
- #scoring_method = 0 #0 = airbus, 1 = TOPSIS, 2 = ELECTRE, 3 = naive scoring method
- SMAA_version = 0 #0 = no usage, 1 = version 1 , 1 = version 2
- MC_runs = 1
- #(area, distance, angle, sun elevation, cloud cover, pri, type, price, age, uncertainty)
-# =============================================================================
-# q = np.array([1000, 10, 2, 1, 0, 0, 0, 0, 1, 0]) #indifference
-# p = np.array([1500, 25, 4, 5, 5, 0.8, 0.5, 3000, 5, 0.2]) #preferred
-# #v = np.array([20000, 100, 10, 10, 15, 0.3, 2, 40000, 6, 0.3]) #veto
-# v = np.array([100000, 250, 40, 70, 60, 1.2, 0.8, 100000, 10, 0.8]) #veto
-# =============================================================================
- #####################################
-
-
- if SMAA_version == 2:
- #create weights
- MonteCarlo_runss = MC_runs
- weights = np.zeros((dat.shape[0], MonteCarlo_runss))
- i=0
- for k in list(np.array(list(range(1,1000, int(1000/(MonteCarlo_runss)))))/1000):
- weights[:,i] = np.random.dirichlet(np.ones(dat.shape[0])*float(k),size=1)
- i=i+1
- #np.std(weights, axis = 0)
-
- if SMAA_version == 1:
- MonteCarlo_runss = 1
-
- if SMAA_version == 0:
- MonteCarlo_runss = 1
- weights = np.zeros((dat.shape[0], MonteCarlo_runss))
- weights[:,0] = np.array([1]*dat.shape[0])/dat.shape[0] #equal
- #w_other = (1-0.5)/(dat.shape[0]-1)
- #customer type 0.5
- #weights[:,0] = np.array([w_other,w_other,w_other,w_other, 0.5, w_other,w_other,w_other])
- #uncertainty 0.5
- #weights[:,0] = np.array([w_other,w_other,w_other,w_other, w_other, w_other,w_other,0.5])
-
- #scenarios
- #(area, angle, sun elevation, cloud cover, pri, price, age, uncertainty)
- #weights[:,0] = np.array([0,0,0,0,0,0,0,1,0,0]) #profit
- #weights[:,0] = np.array([0,0,0,0,0,0,0,0,1,0]) #lead time
- #weights[:,0] = np.array([0,0,0,0,1,0,0,0,0,0]) #cloud cover
- #weights[:,0] = np.array([0.25, 0, 0, 0, 0.25, 0.25, 0.25, 0, 0, 0])
- if type(threshold_parameters) is not np.ndarray:
- q = np.array([0, 2, 2, 0, 0,0,0,0])
- p = np.array([50,5,5,5,1, 1000, 0, 2])
- v = np.array([1000, 40, 40, 15, 2, 10000, 13, 5])
-# =============================================================================
-# The airbus mimiced parameters
-# q = np.array([0,0, 2, 2.424, 4.021, 0,0,0,1,0])
-# p = np.array([1500,0,4,5,5,0.8, 0.5, 3000, 5, 2])
-# v = np.array([100000, 0, 40, 70, 15.319, 0.812, 0.704, 100000, 10, 5])
-# =============================================================================
- else:
- q = threshold_parameters[:,0]
- p = threshold_parameters[:,1]
- v = threshold_parameters[:,2]
-
-
- if type(criteria_weights) is not np.ndarray:
- weights[:,0] = np.array([1]*dat.shape[0])/dat.shape[0] #equal
- else:
- weights[:,0] = criteria_weights
-
- #weights[:,0] = np.array([0.0005, 0, 0.0004, 0.0002, 0.2885727, 0.13339262, 0.87803468, 0, 0.0005, 0])
- objective_np = np.array([1,0,1,0,0,0,1,1,1]) #note for paper 3 - remove customer critiera 6th
-
- schedules = np.zeros((performance_df.shape[0], MonteCarlo_runss))
- for i in range(0, MonteCarlo_runss):
-
- if SMAA_version == 1:
- MonteCarlo_runs_v1 = MC_runs
- score_v1 = np.zeros((performance_df.shape[0], MonteCarlo_runs_v1))
- #len_x = len(dat[4,:])
- for h in range(0, MonteCarlo_runs_v1):
- dat[3,:] = dat[3,:] + np.multiply(np.random.uniform(-1,1,1)*20, dat[8,:]) #len_x instead of 1
- dat[3,:][dat[3,:] > 100] = 100
- dat[3,:][dat[3,:] < 0] = 0
- #global score with topsis or electre
- if (scoring_method == 0):
- #priority
- priority_airbus = [100000000,1000000,10000,6600,3300,120,1]
- weather_airbus = (100-dat[3,:])/100
- s_1cell = 3600
- score_airbus = np.zeros((dat.shape[1]))
- for s in range(0,dat.shape[1]):
- score_airbus[s] = priority_airbus[int(dat[4,s])-1] * (1+ 4*weather_airbus[s] + 2*(dat[0,s]/s_1cell))
- score = score_airbus
-
- if (scoring_method == 1):
- FT = topsis(dat, objective_np, weights[:,i])
- score = FT.score
- print('scored!')
- #note topsis can score a request with 0, if it is the worst global alternative..
- score = score + 0.000001
-
- if (scoring_method == 2):
- FT = parallelectre(dat, q, p, v, objective_np, weights[:,i])
- score = np.mean(FT.score, axis = 1)
- print('scored!')
-
- if (scoring_method == 3):
- objective = objective_np
- w_n = np.zeros((dat.shape))
- score = np.zeros((dat.shape[1]))
- for j in range(0, dat.shape[0]):
- w_n[j,:] = dat[j,:]/max(dat[j,:])
- if objective[j] == 0:
- w_n[j,:] = 1 - w_n[j,:]
- for j in range(0,dat.shape[1]):
- score[j] = w_n[:,j] @ weights[:,0]
- print('scored!')
- #save v1 scoring
- score_v1[:,h] = score
- score = np.mean(score_v1, axis = 1)
-
- else:
- #global score with topsis, electre, naive
-
- #priority modification for airbus score
-
-
- if (scoring_method == 1):
- FT = topsis(dat, objective_np, weights[:,i])
- score = FT.score
- #note topsis can score a request with 0, if it is the worst global alternative..
- score = score + 0.000001
-
- if (scoring_method == 2):
- FT = parallelectre(dat, p, q, v, objective_np, weights[:,i])
- score = np.mean(FT.score, axis = 1)
-
- if (scoring_method == 3):
- objective = objective_np
- w_n = np.zeros((dat.shape))
- score = np.zeros((dat.shape[1]))
- for j in range(0, dat.shape[0]):
- w_n[j,:] = dat[j,:]/max(dat[j,:])
- if objective[j] == 0:
- w_n[j,:] = 1 - w_n[j,:]
- for j in range(0,dat.shape[1]):
- score[j] = w_n[:,j] @ weights[:,i]
-
- #political valuation
- score = score**alpha
- print('scoring complete')
-
-
-
-
- if solution_method == "gurobi":
- import gurobipy as grb
- opt_model = grb.Model(name="BLP_Model")
- #LPP solution
- f = -score
-
- set_I = range(1, len(f)+1)
- # if x is Binary
- x_vars = {(i): opt_model.addVar(vtype=grb.GRB.BINARY,
- name="x_{0}".format(i)) for i in set_I}
-
- # <= constraints
- set_J = range(1, len(LPP.RHS)+1)
- a = {(j,i) : np.array(matrix(LPP.LHS))[j-1,i-1] for j in set_J for i in set_I}
- b = {(j): LPP.RHS[j-1] for j in set_J}
- c = {(j,i) : LPP.eLHS[j-1,i-1] for j in range(1,LPP.eLHS.shape[0]) for i in set_I}
- d = {(j): LPP.eRHS[j-1] for j in range(1,LPP.eLHS.shape[0])}
-
- #LPP.eLHS
- #LPP.eRHS
- # <= constraints
- constraints1 = {j :
- opt_model.addConstr(
- lhs=grb.quicksum(a[j,i] * x_vars[i] for i in set_I),
- sense=grb.GRB.LESS_EQUAL,
- rhs=b[j],
- name="constraint1_{0}".format(j))
- for j in set_J}
-
- # == constraints
- constraints2 = {j :
- opt_model.addConstr(
- lhs=grb.quicksum(c[j,i] * x_vars[i] for i in set_I),
- sense=grb.GRB.EQUAL,
- rhs=b[j],
- name="constraint2_{0}".format(j))
- for j in range(1,LPP.eLHS.shape[0])}
-
-
- #OBJECTIVE
- f = {(i): -score[i-1] for i in set_I}
- objective = grb.quicksum(x_vars[i] * f[i] for i in set_I)
-
- # for minimization
- opt_model.ModelSense = grb.GRB.MINIMIZE
- opt_model.setObjective(objective)
-
- #solve
- opt_model.optimize()
-
- #assign solution
- opt_df = pd.DataFrame.from_dict(x_vars, orient="index", columns = ["variable_object"])
- opt_df.reset_index(inplace=True)
-
- opt_df["solution_value"] = opt_df["variable_object"].apply(lambda item: item.X)
-
- #save x
- schedules[:,i] = np.array(opt_df["solution_value"])
-
- opt_model
- print(opt_model.solve_details.time, opt_model.solve_details.status)
-
-
-
-
-
-
- if solution_method == "PuLP":
- import pulp as plp
- opt_model = plp.LpProblem(name="BLP_Model")
-
- #LPP solution
- f = -score
-
- set_I = range(1, len(f)+1)
- x_vars = {(i): plp.LpVariable(cat=plp.LpBinary, name="x_{0}_1".format(i)) for i in set_I}
-
-
- # <= constraints
- set_J = range(1, len(LPP.RHS)+1)
- a = {(j,i) : np.array(matrix(LPP.LHS))[j-1,i-1] for j in set_J for i in set_I}
- b = {(j): LPP.RHS[j-1] for j in set_J}
- c = {(j,i) : LPP.eLHS[j-1,i-1] for j in range(1,LPP.eLHS.shape[0]) for i in set_I}
- d = {(j): LPP.eRHS[j-1] for j in range(1,LPP.eLHS.shape[0])}
-
- #LPP.eLHS
- #LPP.eRHS
- constraints1 = {j : opt_model.addConstraint(
- plp.LpConstraint(
- e=plp.lpSum(a[j,i] * x_vars[i] for i in set_I),
- sense=plp.LpConstraintLE,
- rhs=b[j],
- name="constraint1_{0}".format(j)))
- for j in set_J}
-
- # == constraints
- constraints2 = {j : opt_model.addConstraint(
- plp.LpConstraint(
- e=plp.lpSum(c[j,i] * x_vars[i] for i in set_I),
- sense=plp.LpConstraintEQ,
- rhs=d[j],
- name="constraint2_{0}".format(j)))
- for j in range(1,LPP.eLHS.shape[0])}
-
- #OBJECTIVE
- f = {(i): -score[i-1] for i in set_I}
- objective = plp.lpSum(x_vars[i] * f[i] for i in set_I)
-
- # for minimization
- opt_model.sense = plp.LpMinimize
- opt_model.setObjective(objective)
-
- # solving with CBC
- opt_model.solve()
-
- #assign solution
- import pandas as pd
- opt_df = pd.DataFrame.from_dict(x_vars, orient="index", columns = ["variable_object"])
- opt_df.reset_index(inplace=True)
-
- opt_df["solution_value"] = opt_df["variable_object"].apply(lambda item: item.varValue)
-
-
- #save x
- schedules[:,i] = np.array(opt_df["solution_value"])
-
- print(plp.LpStatus[opt_model.status])
-
-
- if solution_method == "cplex":
- import docplex.mp.model as cpx
- opt_model = cpx.Model(name="BLP_Model")
- #LPP solution
- f = -score
-
- set_I = range(1, len(f)+1)
- x_vars = {(i,0): opt_model.binary_var(name="x_{0}_0".format(i)) for i in set_I}
-
- # <= constraints
- set_J = range(1, len(LPP.RHS)+1)
- a = {(j,i) : np.array(matrix(LPP.LHS))[j-1,i-1] for j in set_J for i in set_I}
- b = {(j): LPP.RHS[j-1] for j in set_J}
- c = {(j,i) : LPP.eLHS[j-1,i-1] for j in range(1,LPP.eLHS.shape[0]) for i in set_I}
- d = {(j): LPP.eRHS[j-1] for j in range(1,LPP.eLHS.shape[0])}
-
- #LPP.eLHS
- #LPP.eRHS
- constraints1 = {j :
- opt_model.add_constraint(
- ct=opt_model.sum(a[j,i] * x_vars[i,0] for i in set_I) <= b[j],
- ctname="constraint1_{0}".format(j)) for j in set_J}
-
- # == constraints
- constraints2 = {j :
- opt_model.add_constraint(
- ct=opt_model.sum(c[j,i] * x_vars[i,0] for i in set_I) == d[j],
- ctname="constraint2_{0}".format(j)) for j in range(1,LPP.eLHS.shape[0])}
-
- #OBJECTIVE
- f = {(i): -score[i-1] for i in set_I}
- objective = opt_model.sum(x_vars[i,0] * f[i]
- for i in set_I)
-
- #solve with cplex cloud
- opt_model.time_limit = 60 #1e75
-
- API_cplexcloud = 'api_a60ae489-9091-4c39-b7e2-90135566c662'
- cplex_cloud_url = 'https://api-oaas.docloud.ibmcloud.com/job_manager/rest/v1/'
- opt_model.solve(url=cplex_cloud_url, key=API_cplexcloud)
-
- #assign solution
- opt_df = pd.DataFrame.from_dict(x_vars, orient="index", columns = ["variable_object"])
- opt_df.index = pd.MultiIndex.from_tuples(opt_df.index, names=["column_i", "column_j"])
- opt_df.reset_index(inplace=True)
-
- opt_df["solution_value"] = opt_df["variable_object"].apply(lambda item: item.solution_value)
-
- #save x
- schedules[:,i] = np.array(opt_df["solution_value"])
-
- print(opt_model.solve_details.time, opt_model.solve_details.status)
-
-
- if solution_method == "GLPK":
- #LPP solution
- f = -score
- start_time = time.time()
-
- output_cvxopt = glpk.ilp(c = matrix(f), G = LPP.LHS, h = matrix(LPP.RHS), # A = matrix(LPP.eLHS), b = matrix(LPP.eRHS),
- I = set(binVars), B = set(binVars))
- x = np.array(output_cvxopt[1])
-
- end_time = time.time() - start_time
-
- schedules[:,i] = np.squeeze(x)
-
- print(f@x, end_time)
- #global airbus score
- #GS = np.dot(np.squeeze(x), score_airbus)
- #f@x
- #np.sum(x)/len(x)
-
-
- if solution_method[0:3] == "DAG":
- #note, currently all attempts are denoted by a digit - we do another abstraction
- #to represent each choice as the edge between two attempts, that is each attempt
- #is a node, while all the feasible rotations between attempts are the edges of
- #the graph. All edges are represented by the inverted infeasibility matrix.
- #NOW SORT the NODES after satellite and time:
- #performance_df.info()
-
-
- ##TOPOLOOGICAL SORTING - relative to satellite and time
- new_index_list = list(range(0,len(performance_df)))
- for i in range(0,len(performance_df)-1):
- if not (performance_df.iloc[i]["satellite"] == performance_df.iloc[i+1]["satellite"] and performance_df.iloc[i]["time"] <= performance_df.iloc[i+1]["time"]):
- #something is maybe wrong with topological sorting
- if (performance_df.iloc[i]["time"] >= performance_df.iloc[i+1]["time"] and performance_df.iloc[i]["satellite"] == performance_df.iloc[i+1]["satellite"]):
- #interchange index i and i+1
- new_index_list[i+1] = i
- new_index_list[i] = i+1
-
- #check for changes needed
- numberofsatellites = 2
- check_1 = any(np.diff(new_index_list) == 0)
- check_2 = np.sum(np.diff(np.array(performance_df['satellite']))!=0) != numberofsatellites-1 #number of satellites
- #check_3 = any(np.array(np.diff(np.array(performance_df['time'], dtype = 'timedelta64[ns]')), dtype = 'int64')<0)
- if check_1 or check_2: # or check_3:
- #if changes are required modify performance_df, score, F, etc. - unlikely as this sorting should have been done prior.
- print('topological sorting error with', numberofsatellites, 'satellites')
-
-
- #create edges matrix (no loops)
- edges0 = LPP.F==0
- edges = np.triu(edges0, 1) #remove loops and directed (cannot go backwards)
-
- #Algorithm should not include interdependent sets, that is, attempts representing the same request
- #inter_attempts = LPP.B
-
-
-
- #stereo and strips interdependencies
-
- #stereo
- stereo = LPP.stereo
-
- #strips
-
- all_strip_reqs_id = performance_df['ID'].iloc[np.where(performance_df['strips']>=1)]
- unique_strip_reqs = np.unique(all_strip_reqs_id)
- strips = np.zeros((len(unique_strip_reqs), len(performance_df)))
- strips_num_acq = list()
- for ii in range(0,len(unique_strip_reqs)):
- index_i = all_strip_reqs_id.iloc[np.where(all_strip_reqs_id == unique_strip_reqs[ii])[0]].index
- strips[ii, index_i] = 1
- strips_num_acq.append(max(performance_df[['stereo','strips']].iloc[index_i[0]]+np.array([1,0])))
-
-
- depth_edges_Search = solution_method[3:]
- if depth_edges_Search == '':
- depth_edges_Search = 25
- else:
- depth_edges_Search = int(depth_edges_Search)
-
- start_time = time.time()
-
- #remove function
- def REMOVE(path, I, edges, stereo, strips, strips_num_acq,score):
- path[I] = 0
- if any(path) == False:
- return(path)
- else:
- number_before_I = int(np.sum(path[0:I]))
- if number_before_I > 0:
- from_a = np.where(path)[0][number_before_I-1] #the number before variable starts from 1
- else:
- from_a = 0
- number_after_I = int(np.sum(path[I:]))
- if number_after_I > 0:
- from_b = np.where(path)[0][number_before_I]
- else:
- return(path)
-
- search_range = [from_a, from_b]
-
- into=np.where(edges[search_range[0],:])[0]
- feasible_nodes = into[np.where(edges[into, search_range[1]])[0]]
- feasible_nodes = np.setdiff1d(feasible_nodes, I)
- if len(feasible_nodes) == 0:
- return(path)
- else:
- longest_path_tempp = [[] for ii in range(0,len(feasible_nodes))]
- weight_of_path_tempp = np.zeros((len(feasible_nodes)))
- for iii in range(0,len(feasible_nodes)):
- path_temp = np.copy(path)
- path_temp[feasible_nodes[iii]] = 1
- if all(strips @ path_temp <= strips_num_acq) and all(stereo @ path_temp == 0) :
- longest_path_tempp[iii] = path_temp
- weight_of_path_tempp[iii] = score[feasible_nodes[iii]]
- if np.max(weight_of_path_tempp) == 0:
- return(path)
- else:
- return(longest_path_tempp[np.argmax(weight_of_path_tempp)])
-
- #insert node in path function
- def INSERT(x, I, keep, edges, stereo, strips, strips_num_acq,score):
- x[I] = 1
- if any(strips @ x > strips_num_acq):
- #remove least contributing interdep node
- interdep = np.where(strips[np.where(strips[:,I])[0],:])[1]
- interdep_i_acq = np.intersect1d(interdep,np.where(x)[0])
- if keep != []:
- interdep_i_acq = np.setdiff1d(interdep_i_acq,keep)
- interdep = np.delete(interdep_i_acq,np.where(interdep_i_acq==I)[0]) #single
- x = REMOVE(x, interdep[np.argmin(score[interdep])], edges, stereo, strips, strips_num_acq, score)
-
- if any(stereo @ x != 0):
- stereo_idx = np.where(stereo[np.where(stereo @ x != 0)[0], :]==1)[1] #single
- x = REMOVE(x, stereo_idx[np.argmin(score[stereo_idx])], edges, stereo, strips, strips_num_acq, score)
-
- r1 = np.squeeze(edges[np.where(x[0:I])[0],I])
- r2 = np.squeeze(edges[I,np.where(x[I:])[0] + I])
- feasible_maneuver = np.concatenate((r1.reshape((r1.size)), r2.reshape((r2.size))))
- if np.sum(feasible_maneuver == False) > 1: #because I is infeasible with it self
- infeas = np.where(x)[0][np.where((np.where(x)[0] != I)*1 + ~feasible_maneuver.astype(bool) == 2)[0]]
-
- if len(infeas) < 2:
- x = REMOVE(x, infeas[0], edges, stereo, strips, strips_num_acq, score)
- else:
- infeas_lims = infeas[[0,len(infeas)-1]]
- x[infeas] = 0
- for infs in infeas_lims:
- x = REMOVE(x, infs, edges, stereo, strips, strips_num_acq, score)
-
- return(x)
-
- #EXAMPLE from paper:
-# edges = np.array([[0,1,1,1,1,1,1,1,1,1],
-# [0,0,1,1,1,1,1,1,1,1],
-# [0,0,0,0,1,1,1,1,1,1],
-# [0,0,0,0,1,1,1,1,1,1],
-# [0,0,0,0,0,1,1,1,1,1],
-# [0,0,0,0,0,0,1,1,1,1],
-# [0,0,0,0,0,0,0,1,1,1],
-# [0,0,0,0,0,0,0,0,0,1],
-# [0,0,0,0,0,0,0,0,0,1],
-# [0,0,0,0,0,0,0,0,0,0]])
-# strips =np.array([[1,0,0,0,0,1,0,0,0,0],
-# [0,1,1,0,0,0,1,1,0,0],
-# [0,0,0,1,1,0,0,0,1,1]])
-# strips_num_acq = np.array([1,2,2])
-# stereo = np.array([[0,0,0,1,0,0,0,0,-1,0],
-# [0,0,0,0,1,0,0,0,0,-1]])
-# score = np.array([1,1,2,3,2,2,2,2,2,3])
-
- #x = [0, 0, 1, 0, 0, 1, 1, 0, 1, 0]
- #I=3
- #INSERT(path, I, edges, stereo, strips, strips_num_acq)
-
- ##EXTENDED LONGEST PATH ALGORITHM
- longest_path_to_node = [[] for i in range(0,len(edges))]
- weight_of_path = np.zeros((len(edges)))
-
- for i in range(0, len(edges)):
- incomming_neighbours = list(np.where(edges[0:i,i])[0])
- #if zero longest path is just it se lf.
- if len(incomming_neighbours) == 0:
- if all(stereo[:,i] == 0): #non-stereo
- #longest_path_to_node[i].append(-1)
- longest_path_to_node[i].append(i)
- weight_of_path[i] = score[i]
- else:
- weight_of_path[i] = 0
- else:
- ##initiate loop to find largest path not including an interdependent node
- sort_neighbour = np.argsort(weight_of_path[incomming_neighbours])[::-1]
- vertice_which = np.array(incomming_neighbours)[sort_neighbour]
- #naive shorting of incoming neighbours - parameter to how deep it should investigate
- depth = min(depth_edges_Search,len(incomming_neighbours))
- longest_path_temp = [[] for jj in range(0,depth)]
- weight_of_path_temp = np.zeros((depth))
- for j in range(0,depth):
- path = np.concatenate((longest_path_to_node[vertice_which[j]], [i]))
- x = np.zeros((len(edges)))
- x[path.astype(int)] = 1
- if all(strips @ x<= strips_num_acq):
- if all(stereo @ x == 0):
- if j == 0:
- longest_path_to_node[i] = np.where(x)[0]
- weight_of_path[i] = score @ x
- break
- else:
- longest_path_temp[j] = x
- weight_of_path_temp[j] = score @ x
- else:
- #ADD stereo or path is not possible
- stereo_set = np.where(stereo[np.where(stereo[:,i] == -1)[0],:] == 1)[1]
- if len(stereo_set) == 0:
- weight_of_path_temp[j] = 0
- continue
- else:
- #print(i,j,x,path,stereo_set)
- x_stereo = INSERT(x, stereo_set[0], i, edges, stereo, strips, strips_num_acq,score)
- longest_path_temp[j] = x_stereo
- weight_of_path_temp[j] = score @ x_stereo
- else:
- #remove least contributing interdep node
- if all(stereo[:,i] == 0): #non-stereo
- np.intersect1d(np.where(strips[np.where(strips @ x > strips_num_acq)[0],:])[1], np.where(x[:i])[0])
- interdep_i = np.where(strips[np.where(strips[:,i])[0],:])[1]
- interdep_i_acq = np.intersect1d(interdep_i,np.where(x)[0])
- interdep_i_acq = np.setdiff1d(interdep_i_acq, i)
-
- min_interdep_i = np.argmin(score[interdep_i_acq])
- min_ilegal_node = interdep_i_acq[min_interdep_i]
- longest_path_temp[j] = REMOVE(x, min_ilegal_node, edges, stereo, strips, strips_num_acq,score)
- else:
- ind_stereo = np.where(stereo[np.where(stereo[:,i]==-1)[0],:i] == 1)[1]
- if len(ind_stereo) == 0:
- #legal path to this one is not feasible
- weight_of_path_temp[j] = 0
- continue
- else:
- #we now know that another set of stereo attempts are performed of the same stereo request as node i is trying to acquire
- #so that pair has to be terminated and the not-included should be included:
- #np.where(performance_df['index'][i] == performance_df['index'])[0]
- stereo_pairs = np.intersect1d(ind_stereo, path)
- if len(stereo_pairs) == 0: #the other pair is not included
- int_stereo_set = np.intersect1d(np.where(strips[np.where(strips @ x > strips_num_acq)[0],:])[1], np.where(x[:i])[0])
- x[int_stereo_set] = 0 #investigated after adding the others
- ##
- if len(ind_stereo) > 1: #mistake if this can happen
- print(ind_stereo)
- ##
-
- new_stereo = int(ind_stereo)
- #x_temp2 = np.zeros(len(x))
- #length_temp2 = np.zeros(len(ind_stereo))
- #for i_ss in range(0,len(ind_stereo)):
- x = INSERT(x, new_stereo, i, edges, stereo, strips, strips_num_acq, score)
- #length_temp2 = score @ x_temp2
- #x = x_temp2[:,length_temp2]
- for i_s in int_stereo_set:
- x = REMOVE(x, i_s, edges, stereo, strips, strips_num_acq, score) #investigates the first of the prior removed stereo attempt
- longest_path_temp[j] = x
- else:
- min_stereo = stereo_pairs[np.argmin(score[stereo_pairs])] #modified
- longest_path_temp[j] = REMOVE(x, min_stereo, edges, stereo, strips, strips_num_acq,score)
- weight_of_path_temp[j] = longest_path_temp[j] @ score
-
- if len(longest_path_to_node[i]) == 0:
- longest_path_to_node[i] = np.where(longest_path_temp[np.argmax(weight_of_path_temp)])[0]
- weight_of_path[i] = np.max(weight_of_path_temp)
-
-
- x = np.zeros((len(edges)))
- x[longest_path_to_node[np.argmax(weight_of_path)]] = 1
-
- #x = longest_path_to_node
-
- end_time = time.time() - start_time
- #schedules[:,i] = np.squeeze(x)
- print('acq, objvalue, runtime:', np.sum(x), -np.max(weight_of_path), end_time)
-
-# ##EXTENDED LONGEST PATH ALGORITHM
-# longest_path_to_node = [[] for i in range(0,len(performance_df))]
-# weight_of_path = np.zeros((len(performance_df)))
-#
-#
-# for i in range(0, len(performance_df)):
-# #i+=1
-# incomming_neighbours = list(np.where(edges[max(0,i-depth_edges_Search):i,i])[0])
-# #if zero longest path is just it se lf.
-# if len(incomming_neighbours) == 0:
-# #longest_path_to_node[i].append(-1)
-# longest_path_to_node[i].append(i)
-# weight_of_path[i] = score[i]
-# else:
-# #find all interdependent attempts relative to current node
-# which_id = np.where(inter_attempts[:,i]==1)[0]
-# inter_i = np.where(inter_attempts[which_id,:])[1]
-# #relative to strips and stereo
-# allowed_stereo = set(np.where(stereo[np.where(stereo[:,i] == -1)[0],:i])[0])
-# strips_interdependent_i = np.where(strips[:,i]==1)[0]
-#
-#
-# ##initiate loop to find largest path not including an interdependent node
-# max_path_weight = 0
-# max_path = list()
-# for j in range(0,min(depth_independent_Search,len(incomming_neighbours))):
-# max_neighbour = np.argmax(weight_of_path[incomming_neighbours])
-# vertice_which = incomming_neighbours[max_neighbour]
-# #check if vertice is already included (interdependent other node)
-# s1 = set(longest_path_to_node[vertice_which])
-# s2 = set(inter_i)
-# intersection = s1.intersection(s2)
-# if j == 0 and len(intersection) == 0:
-# max_path = list(s1)
-# break
-# #check if any interdependent can be omitted due to the stereo and strip allowing constraints?
-# # Note, intersection is the similarity between interdependent attempts and already included relative to the current investigated attempt.
-# # If we remove attempts from intersection, they are removed from interdependent list, and thereby added to the final list checked for max path.
-#
-# #stereo - just remove from intersection if they are allowed by stereo constraint.
-# intersection = intersection - allowed_stereo
-#
-# #strips
-# if len(strips_interdependent_i) > 0:
-# strips_interdependent = set(np.where(strips[strips_interdependent_i,:])[1])
-# strips_allowed = list(intersection.intersection(strips_interdependent))
-# #all strips are allowed, if number of strips does not exceed the allowed constrained number.
-# if len(strips_allowed) < strips_num_acq[int(strips_interdependent_i)]:
-# intersection = intersection - set(strips_allowed)
-# #locate least benefitting attempt in strips (interdepedent) - remove that! (note, neglects possible profit from alternative where old attempt where removed)
-# else:
-# min_strip_i = np.argmin(weight_of_path[strips_allowed])
-# min_strip = strips_allowed[min_strip_i]
-# #remove the least contributing strip from intersec
-# intersection = intersection - (set(strips_allowed) - set([min_strip]))
-# s1 = s1 - set([min_strip])
-# #As a rule of thumb the intersection will at most include one illegal node.
-# #we therefore search if there is any legal nodes that connects the same pair
-# #of nodes as the intersection, if so we include that
-# #if len(s1) > 1 and len(intersection)==1:
-# # s1.add(i)
-# # s1_check = np.array(list(s1))
-# # argsort = np.argsort(np.abs(s1_check-np.array(list(intersection))))
-# # idx_min = np.where(np.isin(argsort,[1,2]))
-# # s1_id = s1_check[idx_min]
-# # s1.remove(i)
-# # #find other connecting edge:
-# # connecting0 = np.where(1*edges[s1_id[0],s1_id[0]:(s1_id[1]+1)] + 1*edges[s1_id[0]:(s1_id[1]+1), s1_id[1]] == 2)[0]
-# # connecting = np.array(range(s1_id[0], s1_id[1]+1))[connecting0]
-# #
-# # #check if not already included
-# # not_included = np.where(np.isin(connecting, list(s1)) == False)[0]
-# # if len(not_included) > 0:
-# # #continue check and add node to that s1
-# # connecting = connecting[not_included]
-# # #legal?
-# # legal = connecting[np.where(np.isin(connecting, strips_interdependent_i) == False)[0]]
-# # #highest
-# # highest_connection = int(legal[np.argmax(score[legal])])
-# # s1.add(highest_connection)
-# # #Note, should not be added for future nodes as only relevant for this one
-#
-# #if no edges connects - path is as it is
-#
-# #path weight without previous interdependent node
-# max_path0 = list(s1-intersection)
-# max_path_weight0 = np.sum(score[max_path0])
-# if max_path_weight0 > max_path_weight:
-# max_path = max_path0
-# max_path_weight = max_path_weight0
-# #delete the already checked vertice
-# del incomming_neighbours[max_neighbour]
-#
-# longest_path_to_node[i] = max_path
-# longest_path_to_node[i].append(i)
-# weight_of_path[i] = np.sum(score[longest_path_to_node[i]])
-#
-#
-# x = np.zeros((len(performance_df)))
-# x[longest_path_to_node[np.argmax(weight_of_path)]] = 1
-#
-# #x = longest_path_to_node
-#
-# end_time = time.time() - start_time
-# #schedules[:,i] = np.squeeze(x)
-# print(depth_independent_Search, depth_edges_Search, np.max(weight_of_path), end_time)
-
-
- #SMAA
-
-
-
- if solution_method == "VNS":
- #object fct
- c = score
- #<=
- G = np.array(matrix(LPP.LHS))
- h = LPP.RHS
- #=
- A = LPP.eLHS
- b = LPP.eRHS
- #number of runs
- max_time = 120 #sec
- max_iter_N = 1000 #per neighbourhood
-
-
- #define neighbourhoods
- string_list = list()
- BBB = list(performance_df["request location"])
- string_list = list(set(list(map(str, BBB))))
- number_of_reach_areas = len(string_list)
- import ast
- string_list_np = np.array([ast.literal_eval(n) for n in string_list])
- pfloc_np = np.array([xi for xi in performance_df["request location"]])
- N_i = np.zeros((number_of_reach_areas, len(pfloc_np)))
- N_i_rhs = np.zeros((number_of_reach_areas))
- for i in range(0,number_of_reach_areas):
- N_i[i,:] = np.sum((pfloc_np == string_list_np[i,:]), axis = 1) == 2
- one_index_reach = np.where(np.sum(pfloc_np == string_list_np[i], axis = 1) == 2)[0][0]
- N_i_rhs[i] = max(performance_df["stereo"][one_index_reach], performance_df["strips"][one_index_reach])
-
-
-
- #begin requestbased neighbourhood search
- from random import randint
- import datetime
- VNS_time_max = datetime.datetime.now() + datetime.timedelta(seconds = max_time)
- fx_min = 0
- x_i = np.zeros((c.shape[0]))
- x_opt = np.zeros((c.shape[0]))
- N = 0
- while N < N_i.shape[0] and VNS_time_max > datetime.datetime.now():
- N_loc = np.where(N_i[N,:])[0]
- N_size = len(N_loc)
- x_i = list(x_opt)
- for i in range(0,max_iter_N):
- no_ones = randint(1,min(N_i_rhs[N],N_size))
- random_x = np.array([1] * no_ones + [0] * int(N_size-no_ones))
- np.random.shuffle(random_x)
- x_i = np.array(x_i)
- x_i[N_loc] = random_x
- fx_i = c@x_i
- if ((G@x_i <= h).all() and (A@x_i == b).all() and fx_i < fx_min):
- fx_min = fx_i
- x_opt = list(x_i)
- N = -1
- print(fx_min)
- break
-
- x = np.array(x_opt)
-
-
- if solution_method == "random":
- start_time = time.time()
- #generate 2000 random solutions and test objective function for each
- #choose interval that illustrates distribution of ones and zeros
- N = 1000
- binary_int = [0.03,0.18]
- binary_int_dist = binary_int[1]-binary_int[0]
- len_p = len(performance_df)
- X_test = np.zeros((len_p, N))
- for i in range(0, N):
- dist_1 = binary_int[0] + i/N * (binary_int_dist)
- X_test[:,i] = np.random.choice([0, 1], size=(len_p,), p=[1-dist_1, dist_1])
-
- #test if the solutions are valid
- valid_x = np.zeros((N))
- for i in range(0,N):
- if (all(np.array(matrix(LPP.LHS))@X_test[:,i] <= LPP.RHS) and all(np.array(matrix(LPP.eLHS))@X_test[:,i] == LPP.eRHS)):
- valid_x[i] = 1
-
- #find maximal score
- sum_valid = int(np.sum(valid_x))
- if sum_valid > 0:
- score_test = np.transpose(np.ones((sum_valid, len(score))) * score)
- score_sum = np.sum(X_test[:,np.where(valid_x)[0]]*score_test, axis = 0)
- max_i = np.argmax(score_sum)
- x = X_test[:,np.where(valid_x)[0]][:,max_i]
- print(score_sum[max_i])
- else:
- x = np.zeros(len_p)
- print('problem too complex')
-
- end_time = time.time() - start_time
- print(end_time)
-
-
-
- solve.x = np.squeeze(x)
- solve.score = score
- solve.time = end_time
- return(solve)
-
-
-
-def visualize(x_data, x_res, name_of_html = 'EOSpython', color = 'black'):
- #plot multi sat scenario and solution
- import numpy as np
- import folium
- #import copy
- performance_df = x_data.pf_df
- schedules = x_res.x
- a = x_data.m
- idx_acqs = np.where(schedules)[0]
- points_entire = list()
-
- for i in range(0,int(np.sum(schedules))):
- index_line = idx_acqs[i]
- points = list()
- points.append(tuple(performance_df.iloc[index_line, 4]))
- points.append(tuple(performance_df.iloc[index_line, 5]))
- #print(index_line, points)
-
- points_entire.append(points)
-
- folium.PolyLine(points_entire, color=color, weight=2.5, opacity=1).add_to(a)
-
- a.save(name_of_html + ".html")
-
-
-def evaluate(x_data, x_res):
- import numpy as np
- import pandas as pd
-
- #scenario generation specific evalation metrics
- reqss_m = int(len(np.unique(x_data.pf_df['ID'])))
- atts_m = int(len(x_data.pf_df))
- constraintss_m = int(x_data.LPP.eRHS.shape[0] + x_data.LPP.RHS.shape[0])
- angles_m = np.mean(x_data.pf_df['angle'])
- areas_m = np.mean(x_data.pf_df['area'])
- prices_m = np.mean(x_data.pf_df['price'])
- sun_elevations_m = np.mean(x_data.pf_df['sun elevation'])
- ccs_m = np.mean(x_data.pf_df['cloud cover real'])
- prio_m = np.mean(x_data.pf_df['priority'])
-
-
- #solution specific metrics
- acq = int(np.sum(x_res.x))
- profit = np.sum(x_data.pf_df['price'].iloc[np.where(x_res.x)])
- avg_cloud = np.mean(x_data.pf_df['cloud cover real'].iloc[np.where(x_res.x)])
- cloud_good = int(np.sum(x_data.pf_df['cloud cover real'].iloc[np.where(x_res.x)]<10))
- cloud_bad = int(np.sum(x_data.pf_df['cloud cover real'].iloc[np.where(x_res.x)]>30))
- avg_angle = np.mean(x_data.pf_df['angle'].iloc[np.where(x_res.x)])
- angle_good = int(np.sum(x_data.pf_df['angle'].iloc[np.where(x_res.x)]<=10))
- angle_bad = int(np.sum(x_data.pf_df['angle'].iloc[np.where(x_res.x)]>=30))
- avg_priority = np.mean(x_data.pf_df['priority'].iloc[np.where(x_res.x)])
- priority1 = int(np.sum(x_data.pf_df['priority'].iloc[np.where(x_res.x)]==1))
- priority2 = int(np.sum(x_data.pf_df['priority'].iloc[np.where(x_res.x)]==2))
- priority3 = int(np.sum(x_data.pf_df['priority'].iloc[np.where(x_res.x)]==3))
- priority4 = int(np.sum(x_data.pf_df['priority'].iloc[np.where(x_res.x)]==4))
- sunelevation = np.mean(x_data.pf_df['sun elevation'].iloc[np.where(x_res.x)])
- totalarea = np.sum(x_data.pf_df['area'].iloc[np.where(x_res.x)])
-
- EOSevaluate.scenario = pd.DataFrame({'metric':['requests', 'attempts','constraints', 'avg angle', 'avg area', 'avg price', 'avg sun elevation', 'avg cloud cover', 'avg priority'],
- 'value':[reqss_m, atts_m, constraintss_m, angles_m, areas_m, prices_m, sun_elevations_m, ccs_m, prio_m]})
-
-
- EOSevaluate.solution = pd.DataFrame({'metric':['acquisitions', 'total profit', 'avg cloud cover', 'cloud cover < 10', 'cloud cover > 30', 'avg angle', 'angle < 10', 'angle > 30', 'avg priority', 'priority 1', 'priority 2', 'priority 3', 'priority 4', 'avg sun elevation', 'total area'],
- 'value':[acq, profit, avg_cloud, cloud_good, cloud_bad, avg_angle, angle_good, angle_bad, avg_priority, priority1, priority2, priority3, priority4, sunelevation, totalarea]})
-
- return(EOSevaluate)
diff --git a/EOS/EOS.py b/EOS/EOS.py
new file mode 100644
index 0000000..b515461
--- /dev/null
+++ b/EOS/EOS.py
@@ -0,0 +1,335 @@
+import time
+from datetime import datetime, timedelta, timezone
+from math import floor
+
+import numpy as np
+import pandas as pd
+from dataGenerator import genCustomerRequests, genFeasibleRequests
+from scoring import electre_III, topsis
+from LPP import createLPPFormulation
+from skyfield.api import EarthSatellite, load
+from tqdm import tqdm
+from utility.api import getTLE
+from utility.plotting import plotRequests, plotSatPaths, plotSolutions
+from utility.utils import getFeasibleAngles
+
+from solvers import solver_cplex, solver_dag, solver_glpk, solver_gurobi, solver_pulp, solver_vns
+
+rng = np.random.default_rng(seed=42)
+
+
+def scenario(
+ db: pd.DataFrame,
+ NORAD_ids: list,
+ horizon_t0: str, # "YYYY-MM-DD HH:MM:SS"
+ n_requests: int = 1000,
+ dt_acq: int = 20, # Seconds
+ horizon_dt: int = 8, # Hours
+ ona_max: float = 30, # degrees
+ angular_velocity: float = 2.5, # degrees per second
+ cam_resolution: float = 1, # m^2 per pixel
+ **kwargs,
+) -> dict:
+ """Constructs a scenario for the Multi-Satellite Image Acquisition Scheduling Problem.
+
+ Args:
+ db (pd.DataFrame): Customer database. Will automatically generate on if none is provided to allow for testing.
+ NORAD_ids (list): Satellite Catalog Numbers specifying which satellites to utilize. (https://en.wikipedia.org/wiki/Satellite_Catalog_Number)
+ horizon_t0 (str): When to start calculating the satellite paths. Must be "YYYY-MM-DD HH:MM:SS".
+ n_requests (int): The number of requests to generate. Only utilized when no customer database is provided.
+ dt_acq (int, optional): Time between image attempts in seconds, i.e. discretization of the satellite path. Defaults to 20.
+ horizon_dt (int, optional): Time to schedule into the future in hourse. Defaults to 8.
+ ona_max (float, optional): Maximum Off Nadir Angle in degrees, to ensure quality of image in degrees (https://www.euspaceimaging.com/what-is-ona-off-nadir-angle-in-satellite-imagery/). Defaults to 30.
+ angular_velocity (float, optional): The agility of the satellite in degrees per second. Defaults to 2.5.
+ cam_resolution (float, optional): Image resolution in m2 per pixel. Defaults to 1.
+
+ Kwargs:
+ timezone (datetime.timezone, optional). Timezone for positional calculations. Defaults to "UTC".
+ sat_data (list[skyfield.sgp4lib.EarthSatellite,], optional). Allows the custom satellites to be specified using the SkyField library. Overwrites the automatic construction from the NORADS IDs.
+ filepath_requests_map (str, optional). Where to save the map of the requests. Defaults to "requests_map".
+ filepath_satellite_map (str, optional). Where to save the map of the satellite locations. Defaults to "satellite_map".
+ See
+ - genCustomerRequest,
+ - getFeasbileAngles,
+ - genFeasibleRequests,
+ - plotRequests,
+ - plotSatPaths
+
+ for further parameters.
+
+
+ Returns:
+ Dictionary of
+ "LPP": The LPP formulation. See createLPPFormulation for keys.
+ "df": The customer database.
+ "m": The folium map.
+ "sats": The NORAD IDs.
+ """
+ t_start = time.perf_counter()
+ print("## Starting Scenario Generation Process... ##\n")
+
+ if db is None:
+ db = genCustomerRequests(n_requests, **kwargs)
+
+ # Horizon start time
+ horizon_t0 = datetime.strptime(horizon_t0, "%Y-%m-%d %H:%M:%S").replace(tzinfo=kwargs.get("timezone", timezone.utc))
+ print(f"{horizon_t0} to {horizon_t0 + timedelta(hours=horizon_dt)}")
+
+ # Fetch TLE data
+ ts = load.timescale()
+
+ # Create the satellite objects
+ sat_data = kwargs.get("sat_data", [EarthSatellite(**getTLE(id), ts=ts) for id in tqdm(NORAD_ids, desc="Fetching TLEs")])
+
+ n_acq_points = floor((horizon_dt * 60 * 60) / dt_acq) # Points in the horizon
+ t_acq = [horizon_t0 + timedelta(seconds=i * dt_acq) for i in range(n_acq_points)] # Discretize horizon
+
+ ona_feasible = getFeasibleAngles(sat_data, t_acq, db, ona_max, **kwargs)
+
+ req_feasible = genFeasibleRequests(sat_data, t_acq, db, angle_matrix=ona_feasible, **kwargs)
+
+ m = None
+ if kwargs.get("plot_request_map"):
+ m = plotRequests(db, fname=kwargs.get("filepath_requests_map", "requests_map"), **kwargs)
+
+ # plot satellite path
+ if kwargs.get("plot_satellite_map"):
+ if m is None:
+ m = plotRequests(db, fname=kwargs.get("filepath_requests_map", "requests_map"), **kwargs)
+ # m = folium.Map(location=[20, 0], zoom_start=2)
+
+ plotSatPaths(
+ m,
+ req_feasible,
+ sat_data,
+ t_acq,
+ fname=kwargs.get("filepath_satellite_map", "satellite_map"),
+ **kwargs,
+ )
+
+ LPP = createLPPFormulation(req_feasible, t_acq, cam_resolution=cam_resolution, angular_velocity=angular_velocity)
+
+ print(f"Scenario generation completed in {time.perf_counter() - t_start:.2f}s")
+
+ return {
+ "LPP": LPP,
+ "df": db,
+ "m": m,
+ "sats": NORAD_ids,
+ }
+
+
+def solve(
+ LPP,
+ scoring_method: str = "TOPSIS", # "TOPSIS", "ELECTRE", or "WSA"
+ solution_method: str = "DAG", # "gurobi", "PuLP", "cplex", "VNS", "random", or "DAG"
+ **kwargs,
+):
+
+ t_start = time.perf_counter()
+
+ t_start_scoring = time.perf_counter()
+ print(f'Scoring using "{scoring_method}"...')
+ # IDENTIFY CRITERIA TO INCLUDE IN SCORING PROCEDURE
+ criteria_cols = kwargs.get("criteria_cols", ["area", "angle", "sun_elevation", "cc_est", "priority", "price", "age", "uncertainty"])
+ criteria_weights = kwargs.get("criteria_weights", {col: 1 for col in criteria_cols})
+ criteria_max = kwargs.get("criteria_max", ["area", "sun_elevation", "price", "age"])
+
+ data = LPP["req_feasible"][criteria_cols] # We only want to base the score on the specifed criteria
+
+ if scoring_method == "TOPSIS":
+ score = topsis(data, criteria_max, criteria_weights)
+
+ elif scoring_method == "ELECTRE":
+ criteria_thresholds = kwargs.get("criteria_thresholds", None)
+ if criteria_thresholds is None:
+ raise Exception("Criteria thresholds must be provided for ELECTRE method!")
+
+ score = electre_III(data, criteria_thresholds, criteria_max, criteria_weights)
+
+ elif scoring_method == "WSA": # Weight Space Analysis https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=10436690&tag=1
+ raise Exception("WSA not implemented yet!")
+
+ LPP["req_feasible"]["score"] = score # Add the score to the dataframe
+ print(f"Scoring Completed in {timedelta(seconds=time.perf_counter() - t_start_scoring)}s")
+
+ t_start_solving = time.perf_counter()
+ print(f'Solving using "{solution_method}"...')
+ match solution_method.lower():
+ case "gurobi":
+ x = solver_gurobi.solve(score=score, **LPP)
+ case "pulp":
+ x = solver_pulp.solve(score=score, **LPP)
+ case "cplex":
+ x = solver_cplex.solve(score=score, **LPP)
+ case "glpk":
+ x = solver_glpk.solve(score=score, bin_vars=range(len(LPP["req_feasible"])), **LPP)
+ case "dag":
+ x = solver_dag.solve(
+ score=score,
+ edge_search_depth=kwargs.get("edge_search_depth", 25), # Only used for solution_method="DAG"
+ **LPP,
+ )
+ case "vna":
+ x = solver_vns.solve(score, LPP, LPP["req_feasible"], max_time=120, max_iter_N=1000)
+ case "random":
+ raise Exception('"random" solver not yet implemented!')
+ case _:
+ raise Exception(f'Solution method "{solution_method}" not recognized! Please use "gurobi", "PuLP", "cplex", "GLPK", "DAG", or "VNS"')
+
+ print(f"Solution found in {timedelta(seconds=time.perf_counter() - t_start)}s")
+
+ print(f"Total Time: {timedelta(seconds=time.perf_counter() - t_start)}s")
+
+ return {"solution": np.squeeze(x), "score": score}
+
+
+def evaluate(x_data, x_res, **kwargs):
+ # scenario generation specific evalation metrics
+ reqss_m = int(len(np.unique(x_data.pf_df["ID"])))
+ atts_m = int(len(x_data.pf_df))
+ constraintss_m = int(x_data.LPP.eRHS.shape[0] + x_data.LPP.RHS.shape[0])
+ angles_m = np.mean(x_data.pf_df["angle"])
+ areas_m = np.mean(x_data.pf_df["area"])
+ prices_m = np.mean(x_data.pf_df["price"])
+ sun_elevations_m = np.mean(x_data.pf_df["sun elevation"])
+ ccs_m = np.mean(x_data.pf_df["cloud cover real"])
+ prio_m = np.mean(x_data.pf_df["priority"])
+
+ # solution specific metrics
+ acq = int(np.sum(x_res.x))
+ profit = np.sum(x_data.pf_df["price"].iloc[np.where(x_res.x)])
+ avg_cloud = np.mean(x_data.pf_df["cloud cover real"].iloc[np.where(x_res.x)])
+ cloud_good = int(np.sum(x_data.pf_df["cloud cover real"].iloc[np.where(x_res.x)] < 10))
+ cloud_bad = int(np.sum(x_data.pf_df["cloud cover real"].iloc[np.where(x_res.x)] > 30))
+ avg_angle = np.mean(x_data.pf_df["angle"].iloc[np.where(x_res.x)])
+ angle_good = int(np.sum(x_data.pf_df["angle"].iloc[np.where(x_res.x)] <= 10))
+ angle_bad = int(np.sum(x_data.pf_df["angle"].iloc[np.where(x_res.x)] >= 30))
+ avg_priority = np.mean(x_data.pf_df["priority"].iloc[np.where(x_res.x)])
+ priority1 = int(np.sum(x_data.pf_df["priority"].iloc[np.where(x_res.x)] == 1))
+ priority2 = int(np.sum(x_data.pf_df["priority"].iloc[np.where(x_res.x)] == 2))
+ priority3 = int(np.sum(x_data.pf_df["priority"].iloc[np.where(x_res.x)] == 3))
+ priority4 = int(np.sum(x_data.pf_df["priority"].iloc[np.where(x_res.x)] == 4))
+ sunelevation = np.mean(x_data.pf_df["sun elevation"].iloc[np.where(x_res.x)])
+ totalarea = np.sum(x_data.pf_df["area"].iloc[np.where(x_res.x)])
+
+ evaluate.scenario = pd.DataFrame(
+ {
+ "metric": [
+ "requests",
+ "attempts",
+ "constraints",
+ "avg angle",
+ "avg area",
+ "avg price",
+ "avg sun elevation",
+ "avg cloud cover",
+ "avg priority",
+ ],
+ "value": [
+ reqss_m,
+ atts_m,
+ constraintss_m,
+ angles_m,
+ areas_m,
+ prices_m,
+ sun_elevations_m,
+ ccs_m,
+ prio_m,
+ ],
+ }
+ )
+
+ evaluate.solution = pd.DataFrame(
+ {
+ "metric": [
+ "acquisitions",
+ "total profit",
+ "avg cloud cover",
+ "cloud cover < 10",
+ "cloud cover > 30",
+ "avg angle",
+ "angle < 10",
+ "angle > 30",
+ "avg priority",
+ "priority 1",
+ "priority 2",
+ "priority 3",
+ "priority 4",
+ "avg sun elevation",
+ "total area",
+ ],
+ "value": [
+ acq,
+ profit,
+ avg_cloud,
+ cloud_good,
+ cloud_bad,
+ avg_angle,
+ angle_good,
+ angle_bad,
+ avg_priority,
+ priority1,
+ priority2,
+ priority3,
+ priority4,
+ sunelevation,
+ totalarea,
+ ],
+ }
+ )
+
+ return evaluate
+
+
+if __name__ == "__main__":
+ from pprint import pprint
+
+ params = dict(
+ db=None,
+ NORAD_ids=[38755, 40053], # Spot 6 and 7 - Assumed to be heterogenous an capable of acquiring the customer requests in the database,
+ horizon_t0="2024-08-29 20:40:00",
+ n_requests=1000,
+ simplify=True,
+ dt_acq=20,
+ horizon_dt=8, # Planning horizon in hours
+ max_off_nadir_angle=30,
+ sat_height=694,
+ rotation_speed=30 / 12,
+ cam_resolution=1,
+ capacity_limit=1000000,
+ use_cloud_cover_data=False,
+ plot_request_map=True,
+ plot_satellite_map=True,
+ )
+ print("## Scenario Params ##")
+ pprint(params)
+ print("#" * len("## Scenario Params ##") + "\n")
+ s = scenario(**params)
+
+ print("## Solution Params ##")
+ params = dict(
+ scoring_method="TOPSIS",
+ solution_method="PULP",
+ alpha=1,
+ criteria_thresholds=pd.DataFrame(
+ {
+ "area": {"q": 0, "p": 50, "v": 1000},
+ "angle": {"q": 2, "p": 5, "v": 40},
+ "sun_elevation": {"q": 2, "p": 5, "v": 40},
+ "cc_est": {"q": 0, "p": 5, "v": 15},
+ "price": {"q": 0, "p": 1000, "v": 10000},
+ "priority": {"q": 0, "p": 1, "v": 2},
+ "age": {"q": 0, "p": 1, "v": 5},
+ "uncertainty": {"q": 0, "p": 2, "v": 5},
+ }
+ ),
+ edge_search_depth=25,
+ )
+ pprint(params)
+ print("#" * len("## Solution Params ##") + "\n")
+ solution = solve(s["LPP"], **params)
+
+ plotSolutions(req_feasible=s["LPP"]["req_feasible"], schedules=solution["solution"], m=s["m"])
+
+ print("Finished!")
diff --git a/EOS/LPP.py b/EOS/LPP.py
new file mode 100644
index 0000000..be11191
--- /dev/null
+++ b/EOS/LPP.py
@@ -0,0 +1,111 @@
+from datetime import datetime, timedelta
+
+import numpy as np
+import pandas as pd
+from cvxopt import matrix, spmatrix
+from scipy import sparse
+from skyfield.sgp4lib import EarthSatellite
+from tqdm import tqdm
+
+from constraints import capacity as ConCapacity
+from constraints import stereo as ConStereo
+from constraints import strips as ConStrips
+from constraints import misc as ConMisc
+
+
+def scipy_sparse_to_spmatrix(A: np.ndarray) -> spmatrix:
+ # Makes one of the solvers faster (GLPK?) if we convert the matrix to a spmatrix
+ coo = A.tocoo()
+ SP = spmatrix(coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=A.shape)
+ return SP
+
+
+# LPP data computation - multi satellite (no battery, higher dim location_slots, etc)
+# TODO: Make this a class
+def createLPPFormulation(
+ req_feasible: pd.DataFrame, # Dataframe containing all feasible request-satellite-time combinations
+ t_acq: list[datetime,],
+ cam_resolution: float,
+ angular_velocity=2.5, # degrees per second
+ stereo_angle: float = 15,
+ stereo_error: float = 5,
+ **kwargs,
+):
+ pbar = tqdm(total=9, desc="Formulating the LPP Problem")
+
+ #### Stereo Requests ####
+ feasible_pais = ConStereo.getFeasiblePairs(req_feasible, stereo_angle, stereo_error, **kwargs)
+ pbar.update()
+
+ constraint_matrix_stereo, req_feasible = ConStereo.getFeasibilityMatrix(req_feasible, feasible_pais, **kwargs)
+ pbar.update()
+
+ #### Maximum one attempt per timestep per satellite ####
+ constraint_matrix_A = ConMisc.getFeasibilityMatrixAcquisitions(req_feasible, t_acq, **kwargs)
+ pbar.update()
+
+ #### Maximum h attempt per request ####
+ B, B_rhs = ConMisc.getFeasbilityMatrixCombined(req_feasible, **kwargs)
+ pbar.update()
+
+ #### Rotational Constraints ####
+ F = ConMisc.getFeasibilityMatrixAngles(req_feasible, angular_velocity, **kwargs)
+ pbar.update()
+
+ # TODO: COMMENT AND REWRITE THIS
+ constraint_matrix_A = constraint_matrix_A[np.sum(constraint_matrix_A, axis=1) != 0, :] # Time
+
+ # Possibility of simplifying FF here!
+
+ pbar.update()
+
+ ## CAPACITY constraint
+ K = ConCapacity.getFeasibilityMatrix(req_feasible, cam_resolution, **kwargs)
+ pbar.update()
+
+ # STRIPS constraint
+ strips_constraint = ConStrips.getFeasbilityMatrix(req_feasible)
+ pbar.update()
+
+ #### OPTIMIZATION
+ # less than constraint
+ # LHS
+ LESS_THAN_Matrix = np.concatenate((B, constraint_matrix_A, F), axis=0) # (A_constraint)
+
+ # RHS
+ # one acq per request except for stereo -> it is there
+ F_rhs = np.ones((F.shape[0] + constraint_matrix_A.shape[0], 1)) # added A
+ rhs_ABF = np.concatenate((B_rhs, np.squeeze(F_rhs)), axis=0)
+
+ # equal to constraint
+ eLHS = constraint_matrix_stereo # np.concatenate((S_constraint, Strips_constraint), axis = 0)
+ eRHS = np.zeros(eLHS.shape[0])
+
+ # drop empty rows
+ non_empty_rows = ~(np.sum(LESS_THAN_Matrix, axis=1) == 0)
+ LHS_leq = LESS_THAN_Matrix[non_empty_rows, :]
+ RHS_leq = rhs_ABF[non_empty_rows]
+
+ # convert LHS matrix to sparce matrix
+ sparcity = np.sum(LHS_leq == 0) / (LHS_leq.shape[0] * LHS_leq.shape[1]) # level of sparsity
+ if sparcity >= 0.75:
+ b = sparse.csr_matrix(LHS_leq) # Compressed sparse row matrix
+ LHS_leq = scipy_sparse_to_spmatrix(b)
+ if sparcity < 0.75:
+ LHS_leq = matrix(LHS_leq)
+ #################### LPP DONE ################################
+
+ pbar.update()
+ pbar.close()
+
+ return {
+ "B": B,
+ "stereo": constraint_matrix_stereo,
+ "strips": strips_constraint,
+ "req_feasible": req_feasible,
+ "F": F,
+ "LHS": LHS_leq,
+ "RHS": RHS_leq,
+ "eLHS": eLHS,
+ "eRHS": eRHS,
+ }
diff --git a/__init__.py b/EOS/__init__.py
similarity index 100%
rename from __init__.py
rename to EOS/__init__.py
diff --git a/EOS/constraints/__init__.py b/EOS/constraints/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/EOS/constraints/capacity.py b/EOS/constraints/capacity.py
new file mode 100644
index 0000000..c557e6a
--- /dev/null
+++ b/EOS/constraints/capacity.py
@@ -0,0 +1,9 @@
+import numpy as np
+import pandas as pd
+
+
+def getFeasibilityMatrix(req_feasible: pd.DataFrame, cam_resolution: int = 1, compression_factor: int = 2, **kwargs) -> np.ndarray:
+ K = (req_feasible["area"] * (1 / cam_resolution)) / compression_factor
+ K = np.array(K).reshape((1, len(K)))
+
+ return K
diff --git a/EOS/constraints/misc.py b/EOS/constraints/misc.py
new file mode 100644
index 0000000..fc3e797
--- /dev/null
+++ b/EOS/constraints/misc.py
@@ -0,0 +1,55 @@
+import numpy as np
+import pandas as pd
+from tqdm import tqdm
+from datetime import timedelta
+
+
+# Maximum one attempt per timestep per satellite
+def getFeasibilityMatrixAcquisitions(req_feasible: pd.DataFrame, t_acq: list[timedelta,], **kwargs) -> np.ndarray:
+ A = []
+ for i, (sat_name, requests) in enumerate(tqdm(req_feasible.groupby(["sat_name"]), desc="Creating Single Attempt Constraint Matrix")):
+ mat = np.zeros((len(t_acq), len(req_feasible)))
+ for j, t in enumerate(requests["acq_time"]):
+ t_idx = t_acq.index(t) # Find the first index of the acquisition time
+ mat[t_idx, requests.index[j]] = 1
+ A.append(mat)
+ constraint_matrix_A = np.concatenate(A, axis=0)
+ constraint_matrix_A = constraint_matrix_A[np.sum(constraint_matrix_A, axis=1) != 0, :] # time
+
+ return constraint_matrix_A
+
+
+# Maximum h attempt per request
+def getFeasbilityMatrixCombined(req_feasible: pd.DataFrame, **kwargs) -> np.ndarray:
+ req_unique_latlon = req_feasible["req_latlon"].unique()
+ B = np.zeros(shape=(len(req_unique_latlon), len(req_feasible)))
+ B_rhs = np.ones(len(req_unique_latlon))
+ for i, pos in enumerate(tqdm(req_unique_latlon, desc="Creating h Attempt Constraint Matrix")):
+ mask = (req_feasible["req_latlon"] == pos).to_numpy()
+ B[i, :] = mask
+ idx_first = np.where(mask)[0][0]
+ B_rhs[i] = max(req_feasible["stereo"][idx_first] + 1, req_feasible["strips"][idx_first])
+ return B, B_rhs
+
+
+def getFeasibilityMatrixAngles(req_feasible: pd.DataFrame, angular_velocity: float, **kwargs) -> np.ndarray:
+ FF = np.identity(len(req_feasible))
+ for i, (sat_name, requests) in enumerate(tqdm(req_feasible.groupby(["sat_name"]), desc="Creating Rotational Contraint Matrix")):
+ for idx_i, row_i in tqdm(list(requests.iterrows()), postfix=f"Satellite {sat_name[0]}"):
+ vec_i = row_i["sat_xyz"] - row_i["req_xyz"]
+ vec_i = vec_i / np.linalg.norm(vec_i, ord=2)
+
+ req_finished = row_i["acq_time"] + timedelta(seconds=row_i["duration"]) # Timestamp for when the request has been processed
+ feasible = requests[requests["acq_time"] > req_finished] # Only consider acquisitions after the request has been processed
+ feasible_r = list(feasible.sort_values(by="acq_time", ascending=False).iterrows()) # Reverse order to get the latest acquisitions first
+ for j, (idx_j, row_j) in enumerate(feasible_r):
+ vec_j = row_j["sat_xyz"] - row_j["req_xyz"]
+ vec_j = vec_j / np.linalg.norm(vec_j, ord=2)
+
+ angle = np.degrees(np.acos(vec_i @ vec_j))
+ t_turn = timedelta(seconds=angle / angular_velocity)
+ if t_turn + req_finished > row_j["acq_time"]:
+ idx_infeasible = [idx for idx, _ in feasible_r[j:]]
+ FF[idx_i, idx_infeasible] = 1
+ break
+ return FF
diff --git a/EOS/constraints/stereo.py b/EOS/constraints/stereo.py
new file mode 100644
index 0000000..6f0b45a
--- /dev/null
+++ b/EOS/constraints/stereo.py
@@ -0,0 +1,53 @@
+import numpy as np
+from tqdm import tqdm
+import pandas as pd
+
+
+def getFeasiblePairs(req_feasible: pd.DataFrame, stereo_angle: float = 15, stereo_error: float = 5, **kwargs) -> np.ndarray:
+ # First we check the angle requirements for the stereo images
+ stereo_proj_min = np.cos(np.radians(stereo_angle - stereo_error)) # Ensures we don't have to calculate acos for every pair
+ stereo_proj_max = np.cos(np.radians(stereo_angle + stereo_error)) # because we can't use build in angle in req_feasible
+ # stereo_angle = 15, stereo_error=2.5 => 12.5-17.5 degrees
+
+ feasible_pais = []
+ req_stereo = req_feasible.loc[req_feasible["stereo"] == True].reset_index(drop=True)
+ for req_id, group in tqdm(req_stereo.groupby("ID"), desc="Computing Feasible Stereo Pairs"):
+ if len(group) == 1:
+ continue
+
+ group = group.sort_values(by="acq_time")
+ for i, (idx_i, row_i) in enumerate(group.iterrows()):
+ vec_i = row_i["sat_xyz"] - row_i["req_xyz"]
+ vec_i = vec_i / np.linalg.norm(vec_i, ord=2)
+ for j, (idx_j, row_j) in enumerate(group.iloc[i + 1 :].iterrows()):
+ vec_j = row_j["sat_xyz"] - row_j["req_xyz"] # row_j["req_xyz"] is the same as row_i["req_xyz"]
+ vec_j = vec_j / np.linalg.norm(vec_j, ord=2)
+
+ proj = vec_i @ vec_j
+ if stereo_proj_max < proj < stereo_proj_min: # Large angle yields small dot product, and vice versa
+ feasible_pais.append((idx_i, idx_j))
+
+ return feasible_pais
+
+
+def getFeasibilityMatrix(req_feasible: pd.DataFrame, feasible_pais: list[tuple[int, int]], **kwargs) -> pd.DataFrame:
+
+ # We need to ensure that every stereo request is acquired by at least two satellites, however, these pairs essentially correspond to a single decision. That means that if the feasible pair [1, 7] is taken, then we need to move the pair [1,13] in the matrix. If we don't do this, we force the solver to choose both [1, 7] and [1, 13], as the matrix-vector product need to sum to 0 everywhere for the Ax = b formulation. Consequently, we need to ensure that no two stereo pairs share an index in the stereo constraint matrix
+ idx_used = set() # Track repeated indices
+ rows_new = [] # Store df rows that should be repeated => Allows single call to pd.concat
+
+ constraint_matrix_stereo = np.zeros(shape=(len(feasible_pais), len(req_feasible)))
+ for i, pair in enumerate(tqdm(feasible_pais, desc="Creating Stereo Constraint Matrix")):
+ for idx in pair:
+ if idx in idx_used:
+ constraint_matrix_stereo = np.insert(constraint_matrix_stereo, -1, 0, axis=1) # Append column
+ constraint_matrix_stereo[i, -1] = 1
+ rows_new.append(req_feasible.iloc[idx])
+ else:
+ constraint_matrix_stereo[i, idx] = 1
+ idx_used.add(idx)
+
+ req_feasible = pd.concat([req_feasible, pd.DataFrame(rows_new)], ignore_index=True).reset_index(drop=True)
+ req_feasible = req_feasible.sort_values(by=["sat_name", "acq_time"]).reset_index(drop=True)
+
+ return constraint_matrix_stereo, req_feasible
diff --git a/EOS/constraints/strips.py b/EOS/constraints/strips.py
new file mode 100644
index 0000000..a529a93
--- /dev/null
+++ b/EOS/constraints/strips.py
@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def getFeasbilityMatrix(req_feasible):
+ strips_perf_df = np.where((req_feasible["stereo"] < 2)[req_feasible["strips"] >= 2])[0]
+ strips_IDs = list(np.unique(req_feasible["ID"][strips_perf_df]))
+ number_of_strips_IDs = list()
+ Strips = list()
+ for i in range(0, len(strips_IDs)):
+ range_strip = list(np.where(req_feasible["ID"] == strips_IDs[i])[0])
+ strips_for_i = req_feasible["strips"][range_strip[0]]
+ if len(range_strip) < strips_for_i:
+ continue
+ for a in range(0, int(len(range_strip) - strips_for_i)):
+ Strips.append(range_strip[a:])
+ number_of_strips_IDs.append(strips_for_i)
+
+ Strips_constraint = np.zeros((len(Strips), len(req_feasible)))
+ for i in range(len(Strips)):
+ Strips_constraint[i, Strips[i][0]] = number_of_strips_IDs[i]
+ Strips_constraint[i, Strips[i][1:]] = -1
+
+ return Strips_constraint
diff --git a/EOS/data/city_locations.xlsx b/EOS/data/city_locations.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..08d7a3057bd54587bc73e0fa9ee8788910a638e0
GIT binary patch
literal 1443149
zcmZ5`by!qg7cV*F%m9kizzn4*NVoLRihu$FN=tW2Hw+R3h$zwnN~feq3eqXvAV_yN
z_uzZK_kQ2~gWutqXFcnzz4mX#Ui+v$#>Kmhg@r{3{4rt~=_!e-Vq;-N<6&Xl0)Bak
zw6}9HwR17h@Nh77*5h=ywJChxW%q-dOz!wHynXTkkBhn{Wq=&K>rUFulfxE-~Z>gj2s&$^enLHBmv
zKVN4I4N}M%no>V1U+h~)5W&JXj2_zPd{N^O&8KwKktYKeDnxBBKTW?BzWqv
z1-l=vpiglvUk7b_>w|#^D{WSl{)O92$|F=Ly`vUVO;drAw%Yj*I1^H!Ck*2h6gp$3
z673_5P1p02;#S{PSvXT6n-X1(G+0kIR2bCApCUV7kF6aC)iLXf^$J?c|`tW
z<5zQ~^#I$t>f9US_tQ>SLy|!o7t!<{g%k-w%HH5t7mE#E_Rj5F5;VTpYX9T^y)>nX
zRlx6$aj>vRfTw|zsf{xi=M5IVYkJzuO)9nREw^}{Xs1Ys`SKm{dHOEDL?eqGT}4!8
zOVHN#+ZZX*C`$sg+qR0
zY02#J3A1u;M_sfIrD8kCGxj9iwNhLqfgOFjQ?<)d=6>JMYM5Ii`XH}p+vTjoem{&wcv4OSi!lHmOIn+!y&+rt10|EHNaI`Hc>d$p!DPdqVx*SpVx0S=X04
zD}X1y0Un_Rp6)hW|2XB1y_2=Ey}k7fC;#Iht7dL8+3ic>(qmmW3g!E;&i6xJO&;I|
zmsdWqQL?DAygVi+QlTfr5Tvz4j(sBdzTVTQhh=X7T+j(IfW%7U9Dm%Uva$4&VJzZ
zpnZ}%Lc~b8?>D?%DYfT2X_p91Sepdb0gmIy@!E^Bkw+>0z7CzEf?pCQ8o2k?ey=hA
zuYXkA7(xVdDY0I>w8NtQ$3Hioakelub#dnU`{228QDoiA!G|5TS#P!tnpV7go+R(3
zr}JYfD(-GPLbVfF}O-=
zxcc+vdUuek;p%j`;j$UJqv?Axt?GBa(C2ryKYMiu_dRL2>bpLBbBUV0Je!R6J0Fbp
zJ3Y0R1Z1zbRU2-U?E0NH`(0ghU7sAz`kt*;&8|w$**9DP0!Nef*H_04*O$e97gcFj
zXA8v!lGo>|euI*yAJhC!*P|O&uTQ!fuC7M?uK#4iugfn{S7#q*P3|8a`(5t#`JK<7
zOU`fmonG(votzq6{ps_&SWWZ0TDZPg2W;Ftmt6dqcD1))f41*;Sp{e=mb_SiUmciV
z@(u2PH}F0C=yx)oHn%++j5~*#yTlg959A`8fMjzSNZ_G|#J-Lp{nv=XL
zy_i3~#$3nyT@;+$dEsxc+vkXXM!8)bh^A~9oqcnmSbH9O9h`Q3($nWP{|DYV2>6)>
zyK?o-!W)4-&2LKiuRLwLP0luFU-0-*8s*p?i01IjT>I8sUb>uK4jRlvF!)M}U;9q@
zNuom|DRtv#vNi1|{QT_k6-NCan4gza*L|I3(xre3q-K!GNyhSI4uu^o47avyV}W$5oTF
zvz!eES+lbYzJ6X;mz(F?%f$na3-ubFnqA%oenuZRWD+V47FUtaUSIj1pP&CZRxMu6
zm%O;Tx;Q_-I))x8zX`}|C^}j_KEAx9i6>KRhVF9~YJ|MW&fSjNdZ&!3-#uj&^C%`p$n
zE49>2Ua=5U?G{sE96pPdLR?p+Q5EwK(&ryB4Go$;v^8HV7P~^v`$^3LWrR}12t6~{wW2E>QIXW6Miv#OKn+8;U6v)#ox6q
zh0}+2Lb~SpN~(zUbTPf3BG0A3ql#eXoX~7-aD*Y)ynN(MNUo*^ozRae$6)KjvNP=5
zGJJu%Li3N5Y7)*e*N{(K~ZVuel2v~9X!w7=D-1c;6UEI8#a%V
zC#8a@5d7hLu9L-kus@L9FEx8#h3&?puwowA5(2%$VnW`sOZfYYvq>#K+i%!)hIBT~
z^g+XlMDbD#cVS$Arv4@JmYcAyfEh*e4Au
zcC6{Sk(HF5dJ?K3P-#l0lp^UU7AS0ie&JFjqtf}76^ZT
zuFQY=TAiZ9=cR$g&Yc#VrTWHy>zEee31KuqM~Y?-4*_0m-fdpeD8<#`>nEYIl{Km
z{C>48Y1QLj8L6-D>=MnN4CC{QyjCH{X5D@60{^i}`@<13ufd?v&lLQ8T+>B`sphH%
zlWQPM@|ih!Qaak^ldxd#>nW5&O{haMQn)EDQjRd2IXH-QPI!zr%Yj8kf$+s-&WA7j
zlc2BLNn1!^J*03CWActtD^u{TE!CjIc-jnQHMo*HZFPeN(~(cyiqL!5hME(vFL>^D
zZ|ka3%bLrD7{lnpa`5M;>Am)Mq`wn_XMGyk!hSucImS*d!cJbodU^pkym3{SwF)V-
z#92uaA`}j#q~AzTGsv-aT9KPV**mo@I-lJ2peQArby=M2=V#qSWUhBnXu!<>e4$bc
znnBfTjb~Oz&$vX*@GIEosXn+fkc3`Q>46Y+XIi%;q->gh^0~9-FZUez>rEc(y(JfT
z-70OJC3Dy
z&^yxtG@51EBKLb~#*@$s&Q=wZx(QA-+ckuzGcSRuvzcVbEIzK#4RCb3AXL#H0!)J
zYrexSuJ^~Yn~ryrig&%9V(WW@2TVSchRxQ(rF<=MJ~!c<_a}3`THSR-_YT|fby5$p
zI#+tb8%}8(JkYa|fVsO#yD
zi)mew$I>#eRfJgPAJL7A$+JIYSgyA&8V2u~rp&Ath^CXfo
z)inl}!Ij9kcY0L&KBVSv=Ed6_iU2ctlGizr&+Cg`-oMS!-~7n+_1u#8@D~zA+fSI^
z=hVHU)VyVP|V8Md0Qk97eJLJB0NL-!rR&({M
z=IU;vS1ammOlT5%QQfwpNxSvi#hia~GzQ-gf5&o}PG`(NsDUmzfft>zTOOYPg73N<
z870DX$)o7+vL?@Tl6U0S7)JRde$;3ueapgXJ*p`2LaJ7l!L?v*bST?V#xt
zpeoM9$X_LXd-VDiYORj0-Ig3#2nGujAKH+C4CJ{v&_TD-8CQYtX7(Yee{0+2!+VGoo0Z^FO
zMc%WS@iyX&y=dTNN9?D_kS*As^!$
zP0xA$zJ+<@+(N&2ptg6QR)68pTjm*6Q*G+#?==_o?v1t4M@jj8ydSuDyO_9#KKiqm
zryoirj9PQ2^@ST8=2jMR8hePV4ky)&S;>vNMT>L^8;h*8sf(unxj#;%R2Ti{khW_b
z*2OL<9wo+~7Ya1;YFAFN%*`N@ngz
z*T%0?>DR1TkTMe3{JG%&OujZby5_$)XmKoEs~2?gs_Vpx^In@y8&h~y-QA-Ix_~TX
zV*I3;!1@>zuU^@bWBNh#A5W)E*C~Zvx2WH?=-*(tE5{G_N5!p5&bIwRU2a1n1Q#67
z>2oG|&@;mY(C!+ut!DS{p}QV%4*}1i(y)an^k-b^@4hRAqYjMcI$wCDqDp*ev^Jf4
z-0LFGxo3)uywNJftG0^cEDz$p}$Y~EpZM>zV^`K<4
zQV+GHdG=ICXLxQyqHMZt&n@cDK4v%n5XdBkGrhD`m``{z1bG9?l?xw|v1{r*rg8j^
zgttsQob}9{IqVnITbUB796K8^<9fs1zdh|QVjW>6N!Vp&YKO6b9m@mTK^_gDeG#cT
zPmf;%R0Fl}dmD)@dn1?qg~3k02BLX|F1?59wTC%BZg?D=Tg(&>cwgmvx>EHZ7k)D(
z9@ovTz?R_a*N<}
z?|6qx19bV4ws$#V@WnE)e1{p2QKJ_oRxRwgQBe;|6?F9vpqjx}@v-!)N9PVR#q-{r
zyY8+xJGHFmM&Ha^``|N#1Uc*@xwpqS5?5^8-Z?`-e>`7n(Cm@c>`@#odS||Nbm5%(
zhgJJS;?9e``twf-e7qG|PF}ihv)j!c4WauIsQnoRfOZ~|MO{XXh{*_DduP|V?^XLA
zSkQI;l{J-j4l{9PyVkH}FFbky&v&dH8bH9_A1s9Fy?GGz2Uoz)e%#y_
z-g_>P+^6>T1sovVd~Gx`w0?JF^`-JAG8`JIzj9T1p8nYIoFOHT&z}t33)fK#ENvI%IcneOoc2S)fOrAaAKU;ke
z@KGrV1CqMd!QC_Mmffc1N16vA4+>iILh>6OTO27L;T`^z$tO{LJ=M1B7`5Q6-t>=7
z>KmQ9&-i1M%g(}qvRpeeykIpl*;B>xU@P;aT7LJN-H|}*)6q7ju;_gkALcZF-*epT2~_&T%JhR_KDUjEH^agjSE^F1-qRFWAskT
zAd??8E6V;X5OU~DRY(6`p#8ag?~3+PUD1GW6-CzbbNGy0S(A(Q+YtOAyXBp@f0WMJ{ok(fg3;dPopv>glkTi8%zz
z--{>nd@7@gqSWu>E4;^pJcrSn6BXCt~y=OXY{w0PwemOLwWjWjAQ<2@B
zkioB5Y(yMc!n)nwG)EIx7fR6u2kp5ixoZ-B*HIPr+Xq`Dw&sIx-X)Tn#F!G0t5H}z
zaU4sasOI@1v7XB=e%L=4Y}@=S&G=!GwCePsw(W<|r~TWg=lV)YJ3K#kO_tDWqesiC
zo_FLM?`(K#gZ&d}&tq?k?jD6?5y(;I$BkUDfBzi%&f#@9Zbx@p1z>R60mp^^!?8m`
z4Vo6SCgJpxqmZNecr&HXzGi0)zv@J>J8vtz|CH3DKbe;1xm53Y-B#u4X3$?RU$(Yx
zXRmJ4J4yBIbV3<_SdVHNREA5{zTyCZX8S3XF%Seqc$&$->c{=l|6HISIubQzeC=&m
z%EfZ>K#L78^rES}0sQP~a!~jCIWP>eCEwRchn3vdme@>|Kc4W7b*GN2e0dO@r0H69
zSBcz|hNqW#fr+ozao+A=P8pi&sU12jK}e&`ywarCm#(bkhbnx$2H|y?YC)51HDdp^
zJLp%2rfk2o1}C}dIHW;u&7R8)XBtGCS(|E$47C4NZH?I{`=%J{{JHwm^a0o`=G9Qy
zdC-|@L(w=5bwqIHS0Bt5JGZH+2dh~#kKgS|f#UZ0zlFs)w^!4Nv-@E=t>jLVJ=84B
zdhhV0jH|^?F3n=IXid+$v$ZS7iGiK!n-iM;O{V!DXgY)E7^lkc$JoP41>e`=2k(mv
zAQ@R-k~OnsLOrZRICXzm#8pERs^K_f2^z_mOZ1aFL2C)h+nY?#
zBljygvtxQi=M;A(J9Q=^Mw7W^J;Eo&Hpq=@WH@dDY+wfEok=F>mHVBX**V>=Gma}9
z4c0bDH|xO(-WByT$SF>mrxr*nQC=wbA*a}bxgWeT>hRWLRnAmO3!0
zTHi@u28!C)t@2Em3Y~DiT@ao$L8>$r!evg%S;$lC2C9D)-8f!S}#_jJM#sa6H$$_r0fS7MV-b#hb(gm;XwB-RZ(1cUID
z2+&Fd*vbSDF~Cu#kR5&kw<-==uadO}y6z@(j0F8-gNQBIE8U*0K7`R=K+vFc`XIJM
zQV@+A)!`0?01SV;Hd9b`UyW^-n*h1t`|S?L8UNrkf`^)y8mM0h@RM&vN!hq}HoN>|z-
z=_#juNL@<(maX|%DzFf@tuP3y2-jPXfbTo!nn*FZ{6JcEy)F_O{5#SIvpJ$tmm(ok
zz|6`kg;Rp7P!i-{@>ljJAZvfih2###=+nL6j88ro3Ha)Vw%Cf?7s|&~E%d@4^ma73
z92#pxgnfcOnNuBZY&Ak={U2`HXtZ8n>W@Pz;VaK$v)=N)Vl=&lFUCUvv7++rC4qh8
z?EKW^CJl*4qP6Rlrn=?Xrx@HiHYI@-sdc0bC;o#}0)OeBfL0FjM
z$D#MD$Is35$Gi-3qO*0PifP2lI|HAn8be`A5S(QQtm`HU8Tu`ud748xd}P)Sh?f!}
zUg@gpq{;T;bEoGYjFwZx3ELkzXnzMN>m}u~!WeO1@KALY>#TFy`%zjrj`eb`l~bN8Tv>h~056099YPguq(A5XrfmNT
z#0TZ7xBf1W&XrVXtkP-{RE;MNx8|G?w!`^Nl{_
z^t$SOj>U#1?A&4Ao(!WN{I9a2F@^7Iz@6qZ7~1#Qm^vDn-|`<-Fnhrk4iDwq6o74Q~Y=(Kv`k
zWD}V|VdhXA3n=V~y~OD$h|ij9ZTvD|t2te|K@HW*dhl+Wq{@q9dLlB(Hu(DBTNLiDY7XoSJJwDMsVOR>2hnX4Gt!9?5V=EP9#9-lC`^AF
zATZoPGBsM0Hp!!G#jmX4P)~LIuHul{hQgDEv=+E99vlw@7w(+_qC#@cab=vs*7&mD
zI!b}Sdx+lu(A1iH*dl%=j)g_SflY#tL&Bl9j?KJv7heS4pIT>E;@S+YWK=9m=;NK7
zlD8irv&r1})uBlg1BJyxao#~;uG=?hk1&~Q@<_mC-;{#u&muQL;xnm_g%Q=J;9D{N
zr;V|vi1gFOfn&hI-2ro_U&l3@40%30$$`vL+0}gSc}h{TKhqs~*I&XRSb{K2!r^v3
zuxri;mlg7)`H@P_HirJ3bD7D4&B4+D3aiLpZE^0n1ZzQ3Yq$gp-($=6Kyc#O$`|&}
zPllFb{*X%d&Vvh9sy>=ka3@p*OZ`~KCD;g(+Q23F1F(D##I(F_hG1gQYuG1W_xQ5|
z>Y*6q9CpXiIi_NnlWZ7N)jR5}1X23xtgl=XBK!
zqe^@7dk$D;8p$sZh>J$`Ghk)40m%XcJy`Nu?=3K_Ees-75Gm)40shG%kr!}^7`X5|
z0}E{LfXqW6f;@0VtRaa3{ZYlmxl#p0v1N|~BEQ|0H(y^BP{KLDRX7Op2XOPRy-p)6
zRud+XxwdnTSlZhpLgJwCjY)AI&&oLe%f{GCMEYgpz(>6T#Wb15dSoYw#AOrn8dR%5vbym~ZC{Al)e4W8P)XR*Pq21E+4PX)NBJV6=`vv9
zGGw8XN{r&#FytaA?$rDifxx7Oibkg6cd6^jX81Js5cEqni0BXsJA&dILm8U&j
z@GrCP$4}qR8J?HQ2qZZ#{s5X3Zaolq{3HN*Ep2!muyEe<;&|4u_4=pE9Ea?XVV{Y{
zqpY{CJ87;TP}a3P4A^6AWe4J2S(BowWYqq1!hS+dTB+&!%3znYyw-U3&}UFj^f{jO
zb{QRh)(7ahyd=~y9Jo60qYc5tP8|A;opw^+UFXXLu2lHfm=
z?vK5_L+tR>qsVx-3q6j=OgfQM55LT=?
z=Jdt_-9|4a?cUc_Q^H1xYzW;*%0st1+=7BB2qCwqiTU*b6C6wStlnL|C1x2IRk@xS
z$O;LG&Q57GV6KUp8Donc7fflyX&S&`oWp4%$8Hv^`;9O2h!orPzMP5As1G|g(~;}!
zfCT$P{t6ls$UFiE#0Rf#S=Y0=b#9$V>Aj$q+t=ji9=8!!JPDY?3M3;8oI?aS
zH;oA;&;%O6EzSjGs<21bh98;>9y#^z=SnSdsb~(M#Yv+QKlro2q#Y2_ge=kyC5=GT
zzUW^Z!ljkAeIIBC3W#jrl1-wwdgxM-DH3#G``9J`mskct92f+eS$QbbEkI7T(PN#->Tf$nE*zr=yzDZqmL
z0i=xdG+Qz3AMEqO`!3%Ptvdwci>$AP1;3nI|RufCE3J**!>fJ=~lo)hPL?0d*3mrrrn
z@D$jB{Q+dG1=1+NB?n*0ckD?NiCWrv-wbfjFah&4_dm*9^ySB<3}mL1k)jODB){2^
zAm0x({%iBkN`@Z2LP~a#Mjk0m7d>XBQ^h7KA%m5Y;d~{7#qZn{F>i<9ECs*UJkhqy
z=%ur=Jm`}?^nBqgf?0Yjg|4=gYyu
zfpp+T67nOd`>qT@wDRKFl4L3shar4|$0}Vvg)gPoYsi_#V%9fYAE$@T<%E)zgmR~D
zHa!uHbb9KQ&ZS!Q!!hm_S=9NDLbA_xrnQr=T5Y~RT?w7r2_-xED;X&q>DeoTyXBub
z{-tQ%H86FM!9|6tg>dj);31+37PtVD!J-cgqUlA^&l;N1%ZuNePPZ%vt!9!((TS%m
z1m#(#Q63zOJkBgvV4SB9#AT2{Fa&1O-h@w(BZjWTSIW($>-O}a`Lvt1LDi8*=E!Hc
z1)>u&*eMy#85wM5N4~d{K63w=v?T1(onO5~|05!M%PG&pVL-Kz7Mmw7E&*3TD%f8N
zPeG(Mwv3m6ia0~>&mGrZCV5v<&xs;blB=;>e&o;cOQCs0ATFZ}f-x|N9(enwAQIA2
z>^nU^Exz_oeBEd>>g?P~RJ3@b$GZs?9m3|^!^r%@<`Pi*eU~717ADw}mZA>JpeRxs
zXPM^t-luXJ;-81M${osW)vtVi;@d&l#?Zp0Eoy<;O;;^`O%Cp@mSkdyATsfrOKJ%Yu;L*
z!LG~t_zw3xQ^o?5@f=eIGaRtG8kghVc@Fry<-)In7oo3>3@o-^mIf+on5$Wv&2Td^
zXD~AxBbYNx{;w?1HU5N*o>OyBLy<=zBdxZBfnOqq+5XJ-)jg7GlW%F1Z`m7a2;-x5
ztbfqtbi}o`+;8*LVz`82_|+K|jlg-Xao%9@@(f)%GjnKb4%LV3E4k#RF!i2+28$O?
z@`70Ykd
zVBVMI`qd*gM81bGH$jbjK%JZiKqKfvnNFl|qtg@k;wH}f?^h1>9)h-hf;KO--I%P3
z9$wmQ`CyCQWdhx05?wL!7L1T4&_>@6b|VCXt!RYQVti6m<{y^0+F27nrirAw#l$4a
z!~_Cx@kZ7IkTnW+TaX@kEussHs3=|l3)^k8ywg-SwW45P5@lclF);mjj=}B~Zb=6L
zRJU%EA4${wVY^yyb8cfZWMHx`60|6YeT$k&l-gK$(T`DVm}ah9CepZ96qkZ6@H1T$
zp2nYgJGyRb^$|4<=0}M?f`h}*!C?WxVLKXT{(OWy;UephVQHudjfG;Q?o&l4p#mQv
z8_fy35F!vcJ%}6^L_VL}iDgb~eA`)zhIE>9^L?Y|^p87O=JC?ZK|1daxlMbt5~Rom
zq{)ye`y}%Ngx#bqIZOOFXEV1KMy$2yB@>$oVk%DV=$X>#-m#{@VG_k*0^u8QJ}aiO1#BiT((Uhd4JZv
z`UeA3>g!DPq}cL}jjTSBS$AXTAq}rxJBWZR>3~JP(FltBnLCKUgXLm)J1
zpzS9ySy&lJ?nxtRr4dC9V0$P~Xzt+>oRs)#=CNqbGA8wtKNWS`F)I5icv)+;otZW0
z%YzNiz=a!Q!84fYfS{eCL!=&q^g>JjltdGf&y3o_6nbpJIq!nS#%De_KnEP5JOBoO
z?H4(Uz&mpcc6Ef+n~hpO8F2b~#WxHDZ9nV04Ys9QJ^GM<(1(j9vB!t6tV2kbAtWh)+RbzV5!Sf1YS4&7t?>hvDCzm9
zE(R9xXuY@3wqN_d)lfAifgO_J9+J_I0B-m_fyFh1v=*RPHye1!IYhlla=05N&}@~e
zoBnEG1n9zPgx45BYmA}=^HHfkE3ki$M;LM?-L-dHd4-H>sqa4+Qwy;ZKxqlauPBwE
zR!UIRB`C$dCBqPLS<&qKtH*}@R&+}2`w^F|x4a1N$Pb8i5|{~K5)7e07T__2#-{c(
z5qpKXB;u;b9#a1LsI*`9E&3$t-jTx?;)~t!VI3N<7XmvX!#yHf9fuW^pxm1nXj4Zz
zztEY;R&A091=kUu>Sqw>fT{0}Azj9r)gJ82cnT7*Z!rF{gd~ct8@*ApD#eWgGUqIG
zf-#rge}&-Vrz$+JFo;(ewpZ8zFoD`K8gfp=?cwdnk}Sbk1tmAbaxB-m9~@+jTe`LX
zzYU86oByKeooi-WHi3lMJPWaF67*%J*f@4IoRJz7v-qiu{Iw8`=N$&|4#V~iOSPf?
zQ&8Ikk1Z&Jp
z=9`RsN9%#MH0oa1nuomBTsYA&8Bl3rx(rsDLe={^VH;D(VUs(dQZ_kWPZlNC4{Yf`
z-2<=bbn8Y~pE^<4&;8wnm2aWyy`6L$W6dh)-Mg-_Jz&C_&9eS_1|E%OvUT|5TZ8|4
z;nt|bNZAjD-)M0SapObbp3s+-6cgUZA%@z@B_wS|u2(uf?0?;hBppL0SBFXibKs3J
zq?bhN#Bhm$V7qgF*$V;hQeV>+tytu0?aL?sVS=zDI4CtFmaw`;6aeImjBk=d-kM1^SoSos}Lzz2UGcv1%Ma$CJ0}QA*V`FBj){r
zglaEL#Ennp7Z#<;-hz*3h9iwXod0XY!K>yR^FI)`4Gd*bUdPZyK}q8BDc&kD=9W--
z_PmZ>n0V=28W_bLahJ>0_dA0#Q{%;jOHez*8;fp{f@)Qp)o|FYTamH@0_k5TYQ?Tt
zj)43G(0t$@QNRy*3PHIN3_Zmz+r%@8(V4i15W#?NM?%`
zEqtTxub@iup*q;vb#Iq?G`wNanerppBMj{k7T^)aU`2gc5ZTluRz7~eMl;ES?n$Yz
z0;E2ckh)2~%;pDw?BfK~N&<=+K+H`#^iSU&mm~HRW55FrOba<#L
znbt!6uc1?UFIuH~ALEZ~H$jG^s)G6zYq-kQl{JQ?Gp0@Syu|N{zvXf0Mv`_T57YoA
zn368mYV>p5I*7Lp6$}3i
z+eNYTR$mExsOw9ho)MU`NQ-GWZVdOII-)8V3ET~&r+y^1AP71$z<3CP`UmR!DV{SQ
zz4VeZeG1K^xGB?&pgl&;XaxufkA>E3k8f~3*+H%B
zpp1;?NW^D>^`8cLp5e~u9t~~g6r8wA95u@81x&M->?Zxfxx(`cgZPDE`-P?2-w<26
z4sBG~d0~bcq0l`)3i;5kROS1PGh|SV69do4gd3;BGur3?qP%gyz@6k4-YagtS{6hf
z*E_H~3|wPMqWsusQYda+Fp6Fu=~sKy|^muVwkfYR>GU*~_2`++im_F+6qq
z*SuOnKXhOK$^+oS-#l9l{3j}%k1@-I=$945cGz${oG@p4G0M7`58~fk99wfqd
z$=z<5EyR|1c3W`o)4-~^$Yyd`yIZ-Q&5l**t+sz#q`zw_ke)0I@UD395<`F!YBdiW
z1be5(LdzhJ<~Lb~fw6U)WJ{xD%aIM(49j@PITcni{pYqy#raAp!)N-D#!AWQU&He>
z&5VC1+)^!4Qh!D9{)(wCL#C~}}Vd(9
zwVIc$J+m(A&Zh|@t2V86HFc#sg{BNMoUa7gJV901D{d8
zJ_Rg0m|Y(I#q+K
z;^!JlC`^Psl(>|lYck7@viV&9a5%0SM?s=RN+ye!<}(Kcd0X&=2!@HTQQd}4Nswr&8hx++5eO3MM3er%EKm*vmniRwIU|#$TZtwm
z_Efo{&z*@Yv)V0Ni^#h~tNTjH|7AcvUy!dBlVFeg2Dy%o?0Y`+X5mB&#;(j(_G{n0
zf1BicKyiUOjTmp#}EQyxg4T#npa;gb+H#@4+_&^XsJ@1nuLC56arILiz`APPu8(TXtT)fh6%!&r7euS^_I6q5Ves-!WvCE?X3+)(G5n>1|wP_tLce2sb!z=Cy*3~
zSSZcInO5{^RseO>U`~;Mq?c+AY9$9nor4l)_ypoR;LZ~nmv>9)@^PrJ+fQG5O&6;|
zSE`=RaayKi6JWv+`nHG*C{>?4!q}G>(SvueBEU!V!7%xhy5o?S=7MY28U=m6pYdz|
zw^)q_rg$h?_`9&Ig~iJ>b$A(Apll17Xzcm7_Ss*IrHwzi$|~_KYJW){9y1TJ{bes*6{H;c2x7G|mbad6`yqR2ktY2E
z#T=u^93wt}pXX~e2Z`Mk@9j8D<_Uz;?ky~hCofP6SZfaUe8$|4i$wjGEHE%2=_gNV
zL<$NS;wu(Mx+AZt6)sU1n&j@-LxG>k}cHY)vxwclabI$FI^XVE^G_@5bTC6N{@`YCJQ
zs_+-nQMc1ki={U%zN03hO5##%^#zh-H;J5WukTNV_Ltz&U-~^TJ9h~bXF-XRE|}>r
zHUWKlw>YWiTsc1Sc}mfFA~tUG%9zi*)`#-gyj;t7DVPd2M)*5sI0)-JLTnN!;^R2N
z(IosjTW7l#6G-vy{VA%9X+qoj?6}X8|vIA)^pIDv%K+=qu|j_
ze>BX!XwY|iedtaJ{gJ;WY-8{T*n>El6haPe*KQxYNP*tM;s%uUpQ$>mux=2Td`T_DZ`_
z@5#tI$jR;l$af=Yi5-ya#NRyRpzFzJrXAyao7z8;SkR*MV}Zg43OpteJSKcRCXaR?
zF;^=ga+bS71C1{;8jZ9E0mgGB+5P$Hf9tF&X;+sm0P!ve@ys7VAvR3aPaLR7MaYzk
zW_+olQB5^IYP+zgJ}dY5rND9r#V@=-MM$6`cA(;phJh42AyW#O>{oS&+;fIUp%7#A
z9?qeSMTqSP*++Z4w+|i2WQoaSof4YAonO%P+-^C^OSuf6W;1Hei!cf-aQj#1w>}IJJ2NmrTDMxm`K=2eZh4X
zljXo|#MY4e6qClf^2+p)RmU9Pdk3-(N3#0>20(x-0q-Wji0y1koAPkCNr+6MgA*Ca
z*mPa^Y9I|=*njCY8p|^$H6EMV;bp|c4eV1R$U0`uMIV>u8zD8EU-Tuiejx(I!a%V&
z04{Ue4D;dhX+iAY>9DWWpz8njd_K({nTU
z660o@j1Lb;i6a0P@eB*`jv3Z_|4K*~n~bRM5gs_@3>m>gtSgo{fgN_j
z`g-<-Nz&a=IJTM{Yq@4Rw-pFZd(^vzFJ#cS-?#a1~?G_jdD}!V4&@~1ro3PQy
zJ{h9(7}z;Cthlzl32dG?D5P5owKRv@txRi?zHP3d3_kEjCnBeb4kA{KohQnAsg%c1
zi({zU0M3O1mz($SU1UfjHH?|b_?I!KX$Lb@cG>CiU`3&WTV1XB&u}DJ2qZsaNwO#i
z0i&%Nxv_ba*7|rmXJNQrOel8j`H#(^=`Avwm}uR#g9W9?fFSrYfB3u1oADY=5WR1;
zID9(bkwrMS%=3-<-q&S}ck=}5Njo7<6eo-#Cye+2*33sG$8!OCUoyY46fE{y$-m82
zXp!u>ks_N(WPl0OSN$&=Na+r{B}l#K;G2X;pmAyg{#Nz&>61pffPl#Y_cnQ1GoyRb
z@OMb~Gb#8x6S~k8&%1DW0fHID^!NHU6uC+FtMf&6_CNW^+%m8nsxBG74OefQLNNlTxt-YSO@WCnXFqP-LYyc8ieH#d8vCoZ71caQlD45uj_
zze^B-tm+@e(7wKB&5()btqAc}#P(K1-$tW`yBXx9^$BV>bvoNv=^UN1tRS6M
zki-BwGy@;NFyo4>$`2XkbIeVcggzu_l9IDf#u!BNxO}zx5bUFf_E8M*QQQMMxecV~
z7^>iV$yui{x577(^monlHe+fk`T`wA<0qW|P#s(Ag2}^;fc=P3Zw!?Km49w?Ujy+z
z?B%qSdIgud5~lEehZ#bh}!m43waPM=k_D=)?!F_eI#|8Qv{0Jey;9hkFz2B0)i3
z!u3OUta7T%gv$t|`_zrpTnJzzw)i7`Kph0W^7_BBy&f$>x{lCX{re!d*+&Hp#XT>s
zw|y9?ZRF{($Hy!0I(>3}{5Q~U1K4;o(;ori%+7Yq0#%jRCje<~KMFByoj^RIeKSsb
zi&y-ufWu#wIT-@HVYgQ63n1dPsvvIxfukuk(n9hGT1uXeaLjtH72Mil6xm|L2e9O(
zzzP9l+!2plOleg{tJHCIj>6Hg%B(^QHjW_OqshWsKu-8iW&ngma+ZZz)vtN^h1OH#
z)5|wVG^0hWeMSrVE2KQX#Kb^5Vxji|#M}gXw{#u<-Ic9;CG&x3F(t
zFtUWK*@f&OG?t;nn2=?3!BD+w^ntdxGQopM{pZDkU`96M+$7`N{ow=|3
zy7%03p6A?iFUJ!7Yh|lJ5RbStFpTKg-aOV&e_6s6*in&lE+XfyHXN;>CvP<+p9e@;
zz10jzX*;T63!JYx^+7H=IQy;SDmr)-4OvC|_I*EE83i%kA6jc=ex*4LU-Untv8wl>
zR)841&0ZSCadqmC0~!V!-YGFgYxVSw58}!w&Xv^5y}IP=9iIH)35~KIZUtL!d{j<(
zK%hDxOgz@I{p;8^@?t;d7UaHX)@+I)Ew}pT7Hxg7ifC|@`akavKZn1L@vo5d)toV&
z@gm`hSZ;b|S(x$<9;!d0|7$_gSJA**F12ab51DxziC;`y(icJDNTa=)G-dP}{|u_|
zKyp1wWblUeusLP5&!frg>xGMGjeNJ_%a0>poHd3qmikhq`*6wLmX@FZ`BMMaI!-9`%TtchLo)rJ&y|LHcJgKju#1(o!XBkBJ
z?IPz|L^5Wy>A~Wn4f$A3l_Ic9)
z`e3dfUJlLsq})8YXHnwcy0Y9i>f$M4fGk&L6MW|SO{5dyTD8#JFPtvaY$;zAVea-1`6$@}Nx?nW
zenqoO{kg}+8m-yujr#byS`wX^UZ{knzke!Z1jKy%ILZyY!K
z$Umq|BIq9L2;s6>%L;2EWUvfXF-QnCEmm+3KbkxFOS
zzv+Ut?X^XW(Z8e4Z3dkZ?w^E{n`jMHE~V|`Jx_TqH?^=nm}M5caEHH+uL#;?pkIo$*V2NTmZ>QA+7VOg{{fognKpi*!pu%@KGlE
za~=4NkHxXB5d{2m96IlRp@}B{VKU@DOlF!ff?lLbIB)W$#jEWhh6913IrP;LzkBIZ
zjrfGZ?T=Tf!spP?Idr%;u*FQ~EgC<-hH9YjbNl+M?R(K?UrFl2PDa!&rPa&VLH1DN
zE1do-{OVT_A#Tz*EcT;Sgl;^LH&di!UB
z7=PgO+k+&TGM3_#i856&q$QvBcScGHr5OcJgHSr|hB|Y)I_~B=`HZfbvp{E%+hfCA
zD|t`%6o2yqv%Nfy^`FCaR?3VsGDz*A!}riofY^3HZKMamQ0By;{yQ!47%8&K{3$AIxz=FeBip>D3mWnM;G9sUyy
z1^BN%eEO06At23rr+WD|e7k*>>1oQ`z2~~hDheKV81pw-$G2FK04bphQNJ&1eq&7#
zR-23~xV8vYWq!?}FbGZt&?WefF0;jC*?klHpT>cAq~N_mJc#Q|pnKVScLVIWrzcl;
z5_2-AR%y10vObM(r-OT5p{FF20nK
z87n6QvKB?@11o1#8!vT3L)GA5XhnlXMeES8cQ;^t6t9sxu4?Jb{bM!mwdpO*S2tEv
z0xa1R7*(GH@HgagQiy3#h|vRVo6aY1^?(8pvq6g|c5}9oy3fX$ph^lUz2xJzx@qT|
zcqQJoOtZF&Ch+Rautwr^=9qrf^6`yB_SxGk5Jh`v)4EHvA8M}l8SDkc>$oTD%wcur
z41dBOZYDF|5um>Rj8d7%mL9SlI$OsC)wJ_$^i>z)aSJY}#Jz7LlLOwDO~yV0cmz=Eo~
zOL%KLh+Tws3Df|`%6pYexH`JQb
z)p9r2V)o_1ivOx*e(IH=GCoQ^c8jt8;fHV~^yT}f*9?Msh9#U@&G;ELf>!zji<`0h%AC
zUuy{n9YZDcij9yk)5LdNxR{5sHBGkKaCK%f=Q3Pl8BPzd8n6#bz&;A18eazFe!FDj
za(#ADvoil-UKCP%P34w@`3wW>|MX#E#q)M%i8G=1C98(AkQYt1#uz-mj4t1m71=6W
z*Df5~B3w5s@NTX8b3~`Gfz|;(rN7W^0VCtzVTotA>+^Thvr;ZfMT
zoQho>w!Er**=NU>a=ad#Rzr8{Z@7oOw3B6(0%muv9ZQIPBqd^?bzhNc2F=Q#+CcTf
z^TBPsz957C=acF7P}ht7Pw;aA&E!KgEQ+thCSgZI(n8+
z5K|GA`tf5EUqoGEMmylGk3LvDlJCqmD>-7rQ^~L?u$GRZo$oHmzz-dM5E!NRs;I`)
zJqHI{PJiC-bZM-*W!f}L6HA81k|_aR=G*DnZYL(c{Slqg6rTo7QvM_w&66w4b~tEM
zpx?b__GiU2AFiDGzIU`k&dLYXB0kI&x^!LnRPx5c8Zkdct7pJZk(pP_CD?UML|lmBfW;hA|`ppI;wOu^jQmpiEM
zwQ1*xY%qB=cY8Hh_DmgMV;XyK`aSqJ%V+Su)pkTupP{&lbN@cO!8R#6OEBZfi|d{YZ4kWp*x1pGO52@Y
zdyYce-B2ap!t8Sdok*wu)C^=h^O3N@eFaDY*rpa#QmM5?88?oBZC4
z8*1l_ZtOKdeDLvBznuVrLyb*1J;1uZsoK`$58o(d7j%_R1;r|XEC?6f+ha}5c|MM2
zD9}^wIbCh{WUxW#`PH@maVAhQ?|$}>QO{rPv7%j-ag*C(V_`oVQRPxvAQSzEWn_0<
zO$4QP;i||D_Uei=es9|p++03^gR*`px6K9oJp-e@7Ghb)<5-bcRy*Mhc2F5*6_P{M
z%jP>|DSt*l^|!64zizVj7bv!R5<
zNHRW>90%;5uiDdvYqNJF`F7(D>xnUyI-@c6F_eL#R4;aMbg?Je?dQCe68u3%%z28K
z?E1ayrqIYy-hy|vA3b|2Na{ZG<#2O0hdfv@M^aTycEdzv*nWd#bBi7;QuoMG|4DB{
zwP<87jQI5&zpK6+u^)tDJ@dE{7OvW(JWCPuPnQjqkIKhsHMRqO1}}cN`h+`j~Em)K&OL
zH1s3-`^R=WP_V)L{DPJ*)3y+maFM)^o~f3nq;28*6m4U-Z=TOD{NsRf*r&AOrbZkliIQV
zA$OKRF~odvZe(Sv#qRt$=||!3yE5bX!}vO2QnbvYX9l*XS7N!L%qx-SFl4CVaNwo=>~-%O;(;Y-XbHN;
zJne>N0<<{y@tgi!>D3V38CwtAzPOYQ&zr+3?G1|{1g3{A5U3XjE0%6zeK|nGsuf_J
zMNgj#U4IcDQ=la{jIc;$zE`FAG#djQw}2us&}Z7p$G$QnQ+@F===usD%V^SQM3QA+
z=j(tTo=Z)L43=xti9uw15IHW0yi+NE%j-2zag}5C{sEW1lx}@_juhn7_+Tox?1Udb
zMqq=EvqO;p;mlsn`)N{b)6r{(eJ?nnX;z-kYxZZ$;7V-0-*N7{lAU*}4Bp?iGMk}+
zO%bT42tutO0Ir@}eNQ)9SRVVC98rKiakS}_S$fkzqCxr-O{K@%l^pzN6|+bXC9n?wI|Lq<3jDej(wdjf-rOfNCT-gxegz?XNl4
ze*$ALA(c%51wBUbDrEcXkLom}6K#wb4+3G{D-0m6i)f`BBiDi1=B4
zj?epbe0bq;n#n+|nCV@Mts?AP`0RODv#TWv|AHkRD@Wlkg275qmWy`
zUXtmorUwJOHCPu`X9zRUiYWz>V>yy4w2&WV$F?7XNSU2aT|V7nXoWx%Ooj%NDTB%Q
zi@^IP?NDdso-GILZ4c$5CY&YDPvQl%Cyb&MPbDSf0nJ#vl$IrnQ7y8lK9Mrbd4
z#d*^dpAA&t!FoodP
zh-#@EPyyQ;6r#KT%_uEbh`J{bhnV)TzAjmGLG4(7V()eP(ZkWlK?C_%mhm_iBtQ}q@Nz&l_N+^}
zTKFYd_Ypt#Y-PM#ydJIp@I|=HQh|583Sba(SQYm_g+^WAhLf}X>imVVp_N-XQ{aLD
z+oyiT0ZV$}GLnoRFt85<>JOMVmhy2u-tCc>GJ|cYCW+7KuGW2;4b*@6@Br`FpC5Im
z8YwPlf>Pl!RAKI(KJ!P~1z7HS8kH?_AVeS*HT}Lhp#8SD*?C*bK-mX9!yhkD1v#J@
z9ne8jZO>}v^Rhbl<6sW$Th>kca`(3pVF3?8B+KxsL40+lU1rNfm)5?~M$
zsG-L*%(VyzfZa|u&({Vl?E$;Es&e-22d&PZpR8)@b5xVms!22e0|2S&V2^{@wr_dR
z=r%~3$)~QrwK9ED>pn$GUJx~+Qs$yp&Z1D}GCc9Z;&C8tur7EyzqR_)v43FiD{t>;bI
zi9QIt4+pSc8Ys=-c?c<2trDs=YXNXB0#H5*u
zqkh!rhF6f5;h&hzPV8U3?DRC#7GtGm`=G*f_w2R=1|5V!GXflLtDxt(Ll3^M%jW8Z
zYO{Nv>o>a@m8(@Q-1RH#H*$XZzwZx<6xA4FTZ23ui)=OUVk)yYt8l%f5>_1ETX&}0
z=_!pb0_uyP^hMz9fq9;1o(iQ)u&}x=;4#VjWtIv_N^+dcr4@-;sCg*^P7?^IRk-rq
z-Q@sde00=bV;ED^7q0GBN@Qto-MU~_L1GEz6Ntt}O$41_9LB(gF;oDb95FkARZRO@
zj392O;C^QS%ZFB|ay_HX#gGi84@1B6jtpdifC*M0rJV-^Mt2if9|V4=>2BVz@SlrP
zQ~$_oFPrh@y?0N%Sk`S5qj+VmWaTWZGFQ4XD6ji6)h_$7Fn!%$z?gm_daRv-
zC&l;K_${16#LC6W>m-S2=L(Iufw=19x2ghg-9dN>7j%$IyPVs)q_7{~_S74`DIexk
zpHb_~v7g1YkHKkCVbdDf@y&vOg6O+P71@m>nnsdZ^ilkvsH?W^t?s4>Z8EqmkhTb_
z?w5GH`PW{oP2AR4I!fvbnzRU}fJ$2}Subq4;#8A`$M`$7HXCf$;_I3vvqxVddt($P
z9!W%_gQDBz+}9=XE=Rk4%mgTMG{Id@#9}r!S~K)TfSQT{ehp^#(U>B6h9a2*AX`$%
zE?M<1LvqOa&%U}$Cgt%}?jG#2tWw+fbaPEBZ}!HCxy}>aCaDs7uoN`RL>HEl!1sxK
z>o#nB-urC57*EW=YMB=KmLt^v>E(`fLrJa6XZVt5_>wsQB0zP)vWeN`ye<0SuhRRj
zc>0k81#36ldT*y}l`@%5&!kO_t4l4bMU9JI{)NSmjqk)(jrDMOv0dL{#J8F!Xy`=F
zB|*k^ajiaMHqh5(4uE{1(R~v$9TB_`udOndE5+0GmDQojQHQF7`ZBGhNV9k~u4J_=
zENEFS|CGwB?nJd0w-mM6dpn;keSQ`-7o=D`QwL2FP?5b|F)Z$up}#<^nKLusE6y%T%Z}jj8Fi2)tHfR&^Pyhwsj)4sYj`B*F#8GysCdk
z&*>6}6!M02U<8^m0v+_TO%Pe1fnNCBH0BT{yDAe@++3ZIhEL^g^r$y5UQC)gN13%yzsd`QF_^mi0Ho5``^~SkSk}&n#9dT6oq2
z;>;sZ6(&Cstw^?t9>KrUZwTOh&RAzEu90IxPMRSnagaf~1veDH`Ig+r{cP&0b?n~O
zUxc@Td32|Njs2%a8LNV*2#BTFW8vXYBX1q>%;QhD!gu*EaeRPmMno%D;-BlIO@69*7
zY3J>IsD9dSBBZf2Vq>NGao+=>%xbFRii3>rFtB$Rs&|;ZM-?t*{08!d7I9@OPXsT%
zdUChr+owBrTX(D`#Y5reCMaYX-(p~IF;s6c!gfh@7x~|GE0p-&cHTE;?v+E=KC5L;
zTlG8M4nMJSkC5nrzX+k?R0a7POo0+;MhWz^A+KBlb8*hw%>;{)g4=q^&!ahrs^0ylmb&Y40|
zOL3)0Yfcn#Y+;rMaqCIq`AXV`b32Ws0+p6yogX4BzV$7A&mW`clFMq5M_sgaTEl??DtH!co#}0Mxs(G-Qp1vH0HUUU#z%A72|=P86XTtiEi0h2
zhvTbv?cRAgDw|a}9#B>Ebux>JGrq#WUdh2CU%h6BfI=%MEe#Y|CokS{&v#!g`Z!V!a<2#>ka+ng9JXiIr2g%D9nvcy=O3(-5g2{v!{$F9H|d`i;*VuO4Kq{3?q@F+cjF0yufvt
zQ1GT`?g+o>IRU2kW%4MbRcZB}7A)k=c=Ywjvf%qyrn0loNR6R_#!!r7s7mwRO_ONp|z)!w(d9sRiy+qTaC%a|Qr*NTLZK4TDj&3u2mY
zJF(jFu6wQIMQ1+o>yOE?buU%5zz4+qVF{a-IHq>bZyEdEv^G;HWGn9V#k<-|pCi)W
z7w%l7C1$3vjh|#vRKCF-IE7-ILRF^QxX_-=VruMnQl@qm=sq=8AN7c;*y%6<4Rn3)bxz&xm3?aAPPJHdcC2rzzj&TDobf3J_7p?)6ti5+
zZ43H0$CDw{X@5S;+nsryY9}A>kT;PFZftnW7?u&p1X!_TAXBNfeH`~wj9CDDNxL-f
zb|=N5>TGtnE#-aLL9)B^d++hb#(&-(huO*a?oB)K22$S+@jgi$=pWM;Fo}QUGyiqm
zJ3z)E%&}KHFYy8be*qD90kJj!slUzbUCxaWbP3xN^hWymtkqJDdib<=-S;T6TfU27
z>abvv#4wd=*~c~1QLUlmK2dm#Ui*kz_yt0*#x)>n?I>Mxmcx^z=1HRQB-I5GL+ft$*ES>4hat)TGv
zd*f!jaU1{x6G3I;NP-9h=?%ZZg{Pj&_bNTE7Zruc7#_dCnQ`uFw`ZRbl`I#%EVF;e
zsURibdG9Z`vL;~Wy}3v=K6M2KU9pvYO%t!)7~Fk2Xm+?zniqPJ9(G@@w#a$Ln
zOn+bFvl;FGr10hI%5Jk+k=gsXPKaMsTJf8Pcp7U2)EYr)jliE)d%k*Ep%{{o&l*oy
zy5yH&ioI#-!*J??P8Ge6EblqX_mA9N#Edbx84QjCpzTrfytE=Tq&Jn!248B+(mC!Y
z?OSM72r3O-VOZNyh|S1Uqe_;^!pd-UiGiyXnq>Ca7_9Mh5n)3%79SQ3-T~QbpT_IV
zM-s-nO;IvjhB8^YGF-|3lBHm(1{bhAX}-&4;#0b9G5)wMb_fnvts2
zaKXV43stAEB1avel!NVVX?;`S*9c^1q3M-)P-JOFqA??B217hLq_ddH``004bqJyd
z6ZsZrSrNq8G-*`%dUZy8nk!|N%Vdh{jd(f6CLgtQrLl2ulj}
zuPYdq$ZiROQ1k6brst#Vf*|v6q*s?SBKbnMJO$NQ8Qk44=P^=QU+5UiF|cyE#D+Wn
zdU2+ukxEk`X_5ESJ$NG)Vt3`G=lr>FSW$(sE^bB_$Dxb+${>@)U4c1CY@67A_`FJd
z*iR~2l@|jSG*9wxVf2vG{QpRn1vQVem@dU-p5Voh)~Ax>Pg$w8*PWq~;i8wxqL9fl
z#E`0Wal!GBXV91%B~@NN-@*nim|;kss-uVFektmn&vpN=EG4(V&nUh#xMT!w5>evZWZJq^w7F69L$%vO3WGvX_L9t@X>-++6#NFH
zyre@hG~9ta?oN+zOoJeV*oLcB2){ROw8@&yr?xXKiuGz4?$UDULb}1D2kFKz@CK`l
zQ>srV$Via{Ahzr~h4{>C9eQT;iTo6!D)+Qu5K)}FC~giQB9RH1w@Oj4s;60?UNNI3
z)4}fj={yX=n67L(sk9=+WRnF68eKEU+F}Z~oaQLW8l*Qw}DN)#YEI#FY
zrVL@jxwGNs0P-F6LYCU(v9qoYoz|?&UP}AYi-&URRRM%KZ|dz=n{Vwjs1y=W!HFmZ
zhBjPoZyIIkvyqVYd3rzFkh^)ZsraqyQH?|IGZ2P8$`A^iI|XhIAVVUPavsxBvlPZ(
zQgV&lHP@w=z@UYq7|H!zIqlmNw9zPqXjCvjinfCjJvvR_7+>ap`*@06TkNuC`(XXO
z!GX^4Am@|t5ZFn|EAlecFCu_subz^Q*v-D)9Mo74(}t
z80GwBJVh~T;J#B*$UjF{?;`Mf!3uP1ZF!`0PENe#rT=-<#X;
z5&<1cf5KCPN4kD2e0cdEuEERqL|ae
z`zcI}=u8&N*qX^M^c#}?zY!bEtvJCf9mk4@MJg85>dV~Kn2xls?^~5@+8G=~KtGD~7d5NVH2xbV*213prGb|A@;}j_@O5WxOU&Z&4zZ
zJ*uYh$deMwrR}HHT{tpaIWpWh&P6sU5bK^P5Wl`u=tqAwYNB$hzHx1&`8-`}qLu8v
z$(Mc;dLxvA5h~aSC1ZY`&MG{d7`=v7Z!M;zVxxUpeY@y;{9ykLg3nL<&1aR2t5`Tk
zhB!xtB*&-7TlL}`gT?|hUvF6Ol(HKA9>wJwOAMbXTN<>HlNX<_K%_h!gDUgI$mf8fJ114Q(DDnl1S;v6B-{Z^vWyOXB_@6l;^
z-i^<{-T>9@%+J3eL1TQ)ip?YnPrSk<%GeDaIAgaMSUW%$EfQKvE#4~N8RJ%~ZZ6%F
z7!ILSH^>yOxpwtZKrq!A(C$g@EwP&t+ri=4Q}oIe<=H_!xg%>`^9=9Q$n0F7i5d?_
zd15puGYHHX1Q@`OO{axl>ETAGCQL$(oW-!o3KNERw!sTl>~K1baMI6buHrhDN&?O=
z0Y?J-GKG5gj$_z|5z}mcPQzDh-RH$_>a5Oi()~^Iv^@bSnp$f^YaClKP3ujsBsQxP_;nW982@NyB;9k2oOKSC!XG7-4W@$=L
z3MHuX9{XYmVty*PSwTnom#pX
zau7;8H75rWSKl}5Ia`1Xq`MW;J@qfiTLxkYp(7|%yN=+mSfnGb4%JfTy|OD1Tku|7
zK%#~1yjUroq7&V!DPh|t=w}B&V-~*4vi3Q@sC5D_4i(qLr>jS|RMKk9o$C@x!
zjsuJ#2ytYbxLhqNxE6KK%TQOtFwlY_j2|JidP&(&V)x2aOovdczmoCd_YJM+5h{B)
zzdamj|L;bT`J%T@dM@7dW;**FC`+F3^{BPe`ZJZsxu(E(L%OFSHwuL)*I<;Uh7_N>
znJ=E-nL4|;2JAhyM_iZMc~{caB~s3ZL(~M&_hCO1TJS4
zS9NeQfNm+T`30({yFcW@(4^u?DpoAc9gCa8;$o|9If+it4Q9$~HyJSwW2G039t*w;?vRh-ezaJ*GFqx>cqtTnTq5Ck=;^Xk
z3YeqS1aGuWG}JJ9pT<6PL2Gk2#S$5Hu6MsMEv4-KVh*PqQ#W|znYuBIKt+?VUU4J>
z^P!++>lpvjgu);hFCq_;@w*^6$Z(P3gaq{o3HlR&_<%*iPgU4KLrW8F5>+$$YUJ*c
zqS>p0#N~>rFJ)njvgfJ%;QW4Yq#u0Ct-Td`!N}Q<4^$;
zab=si)Fw6E1tTExb1bk$L#M@AhbW@P=I$N&(yah+M0xZnOX|V1Gfo5wdm}n~BdEPm
zl1W7Qm&0PAdP2k~@5RkXu^2gzJ{Yq=?jHX00S?h4Z}_|DA8}WHn*Al!!$sGFQs_Yi13Vfo$9Sw)B&o1fKv^q|May-c
z)yp9`lzZW3q;U8$&HqRih)>x$kwqhSA$?z7F7VSS(>nLKd`vPB50lpA$A3t=0zzUyS~DdRuWvIN8YT;$quc)dE=JXj
zBgLMC{~^FQH%K<7|F>=nq_5Gry$=Yy$e^=bq27NzS%z{!^HeoFsju>LCyQ~-w2wNG
zJA=fXQOrFUv8uwSaZ0Y4VYd3$SplB$Wfg(BiU6x32wULl(cFW)*0u)k!*^Go2Qv2*p?zg%-p#Ag1^t|+acXISOxE7w?RIl=p_eTYv4z`y
z&UuwRU0RU4@nV*_iZ>g}9=`ndYd{w>y6Ri-)<(`Yt}DhmxeEII^3$(->UD4DQo7h^
z$hG~N#abIBSsMY7uVnbxBgoQQnPF3j-8hqvU%r_uAZhGmJQ-FRtX3HA)VDD4S&{Pm
zcElp**j6#P?bF;Lw|BVH)J>V=
zIj2}nLY2AW?|@#qO`T`(n{_JQ%H;a@RCd=g96la*xoI5xMaa!wMP@_d1kd6|i`)9W
zg=Y|NAAOlQtQ9%Q$_Q^oraK*RO(mL|TU|N7IV7!rn66NV47;Ly{$o%|XTGH$5{56D
zqX1%7Mo{n
z7TfOg%;cQf?%*j#v(B@Mm*R5mE?`;02?E}2H#E*YZlLm$!vz#B-(*Sd*U2+IhaEZ#
z3FY#?_rtdD(z{Z}pf%k^Narw*vhRQgoIVSsl7$Kf$k%qYW=&)Rfs3|)@bmXjBed`z
zbfvuTr|DLnSM3fId*SVw)LOXeTDj`lx$33?qw`B)_|e97?;?}+M`ElSA|X@p0^2rO
z*O*n0aqQvLi%!Zi0!tY|2w*Y@9zu*;c+1#MP#43#tGV{1XMTHmw*Tpgd8*;dISW%*
zOLCGW8E;9Z(+B+TndFp`lwjqg{#(^%NVwl^S?G^tgX*hh1E)@e?1hfTVp%Hsck-LZ
z2uF#Yr!>)VkR|i}lcgU_SCpUU9VIx3$*eyTj-R&QTdQWF6vU7&DP4U8XaDqV_;&=o
z%p8-It=z9GE;N!=%j^DW8?w-~By3mE$i`V#je)m3SjadQGWlYUSC9k5h}7I)X=RSn
zpJvItdxrDRvCm#4(E3^^<0@<24Uho)bN-!&EWX
z#E@~)|J`UT_7uhy>ANazyj0bdHr70|8A_)aCoCeDsCbSuhw6_?%o&=&H`H)H+oi28
zQIn+v+9zY+zXKsmZ*)12B^WCS6R6GLf@W|HGq_GLB(h6-&e;rJ);sHFBrMS3Z00KK
z3sIo+dBWwPp*96(maA@-s}79T#^+yKELIT++3R=Vb>>@r}pOF=c!TDiu
zHB6fwN3x3H!JRvyCTE6+IZA3olI|!r2P0g?#K?Rn92h{DjKP8B>vu$1>oUdP4-KfE
z<*h2l_3u9JZKUDS3Sgnk5fG0Wr?P?zTET&R733k^2{GhWM@q%ZxrNUnqC7jjC$3g*
zTuquDlBS)tu(;2W!pIy)$z4ay4ZiN=>$!rQ=o`1nj@h%iH_<2su-k6lLg4LI9;5=BW2toi;1L5UsbM-4rL>ZxK
zR5a(6padclt&N
z>UULvhSu`S9W>??@#gj>GW}~94Wh3Y|4jFle}n1~#V$dGm!MP-NBm7(jWLln!r`nF
zKhE;&lSIG>l8@~a(x?2`6?^-J(Fo^fgj)rO0hyHoE095?X$$63U!*57J&9$5o`~tt
zgI^H?oA*d?gmzVi%hsN%Gyag1a_gZA?
zA4Q}C%5@Y&#>t61J{Be7-O&Y$k}ChHDFoZtPAn33Y5-y3pYMBL9*nAao|mwHLP~-F
zPrp;Kf1CI&xGmLZ>9+iOt-d5X?2j=RkApnZ9^nk0lw@31H6}7$-Y~{Lzwo=NH@>av
zer}B)OXTN@fZPA9vG?Wipj8RILuoPZ+XwRHkwGlyUeV%cy`D#<(m*yUsLsL#&yuKz
zUp8uUa|KWG$L}llP~?-U&*bI>?lia6CJeHr5glZ3cl#*Lk?GEnq31v$NW{ZShh0~U
zNY1onaB{H_&S@Vq$#m+I;_)4wv(^5kq*}GAuQON+a*_oZZ$YL*9mVq0I~Pt9KbumU
zHz+*Mp5-qOreSnwZE0t}rWz+Amx>BcMX9hI&6h~U<5W@Zk0q6|MEgm(1R%&adD>fl
z@He|w!GATmIRB9!XZEgJxMz)Z?wSjS+}NJ*ZedW@P%Dd8HML>d1DG6S96a)nagh2B
z4v!`FqT&;!vr5xb9G+|7oPuN)81z*1V$Ior!u+cP0qHmz|MN!u&AqgX*ZM;Q8N`2$
zN3SsCguT*Bze@5BprgNzQn`)_zm5t>52Nv4pqmJ#Reu`
zu1^aze+21Gu03*P`?Rw2drQc|wjjSx#d$@2x)VrHtCj+3dOtSYD!`*OeSIeS
z4$OdC%JNKNWUm?2G1wlRaI%I^@^y6(7U7P4L3{29!D
zOPqlL=f{9s1<3Z-5Pow+vp~WPv`#(8Z}PvFe^5!$QPTjoAEnZU3u?nP08)^6Ll!vI
z<)yA{+Kz7D?0>~>#eHY3u`RSn$S5j>f{inU@gtVRSy#+i*UI`YGssTyVn9aT{rK?Q
z8J5oL%S<;c!d^Yq2Gb-c7m5lGMXBt8WLjOO<(J3EvMeRNjKf$#V6T!HFM5U|iaWG2
zIQJ;gQlY}j9dtU41j{SI%TFqn?22}*snwdBcoO*OyB@H5g%aUyozl?7M8Czv8oHH^
z!*RBcOTYvu(Zt&&vy=;o8WrTLoPTmIm`aBudDQzY;~&Y@#l(+AEUg7*i}TwXZ<7jD
zLNSl`#74DQ^iDB;Sd|x+n+T>-gA1y`HPk>IxK3a)XOwd2%S#o4>p5pBnqT}5T|B5s
zD-OO>afZfjmcoooXGVsakq^gNKTd;pZ&Newy|$2h?+D|ckqW2~XR{LIsy9C!OoevY
zOswJG{L{&;Kb=H8;;x@0WNECcdixGtLZwU9AYdV!PC<_QI1q300v74Xn{7l~nMeUXigwhD5?K?UvSQM|mqurim
z*zR^GnBa=DdHxP77%acdqhzL`?hrc#BAq=>H0B_~JQpvj8gISgJ#ReyoJS2L5#0@?
z-2ZH|AmcOD8IU<0;M!akA6=>P*dn}ozEkTjkpXZ3r51rq|4u;Zj
zcd%t(zr&`>IQtpL!85=(zIl|y#AD0E7x^M47Re7Lub%m=xV|`HSZN|070NUj+9KuN
zF6G`O%Wp#wwt*3gj@#0wf8=4!$w}s9yg7N>7F|2b
z`Sp(Z*hpguP5;w~du`lTl$J#2rTK8bG%~)F9h0N^E19Glac%my{2FsBqwwniU7c(F
zCovAIQ$mXC-opNiPu@hWr`%a;VdK*!p;cduu^m|`{&-H!>)WQMPth2RNF{%a?rZ#I
z$oJxW$VP&K;jd#9Ty>OOb<}_hu|*-T=#nx+;^|oxC$ZN97RSGpu4-n|II3D`X&F@x
zwwb-5z>w*{Asq}kN%ybwiKxA|bm7i@8$$2S39@3%{kz(bZn
zUzfH?d_$s6W-4;Xe5!UHwzPF2O~CitB}U2~7zJb;Jd}{}Kiw_FdSz1#{Ys%Qwb3NF
zcd`dsG#xC78F{8onI7>=V!!FDm?L{~lxj5=e*KhsRIE-dr+&q`VTqBhkjYcDvXAit
zD>Vu(h=OYXr2CVP%wN=t*4PiD);RXMos!oJcIXi!NzeGE!Z*U7?SBor4^^R0{4MDO
zWY*RV4V{>dX^){W2lUS81_RD2B$lfVBxT3lPN{GKpW>K{UDzrbAQs=xgSj2chZP<`
z?b#cbExZzOILcxIOEE#HnDFCR%cpn-eNxgkx~i`pEpW{88@!b
z6t;*$z~zndIs*(jSB$BC4KlR+e4k!+p1CmC_gET*qiD-M{1Ts&;9uhpyo7MNQv9R0
z9I1CDDQ@N)_)6bdKAhMaWTm83)~Ad_t#TZY}Y*QFSq}M?n`ouIa||uS(JH8P}2iWSAu$K@oMju6&WzX
z-vr1?oKV|~a2@ZwS0$*aJ3F(>KEoUvM-KfD0~AibL+7?Y-?=+T(GfPx`R2w)<^33d)F@T4ORdZ*TLSHFvX_Xp}BE%0@)kkvT&E{^PPr7b7N
zfJXO}>$!btB$0LZuKEAu>_EHk1p@2O-J{hW3)G7}fK!y2D=ar8;(Ne(SVjW`Mgw4@
z>yiaHH1%P?5yz}6pXn<2Ny98sLdldIU>tro?_%nJ1$(&%5$)7jE8Lu
zSmV2W6z5e+0cl@osOF%DbZ7qg(8mwyT{TSHLVNLS$Oeo_kM=Ec@P
z=p7?+6)v|&=x1-Ai{V&fgj3O2xuJ6v$#Fj*b
zq>t&{qv3j_9?lOTFc|=wKtWFr_UJB>>qYQZRw+^xa?Sc8tahBBas5f!PfE#M!XrZF
zbmXp)cVmKLdME%wBoVf`Y~9}jR~`IJsOMl!jXexmSQ#Em!8rbU+zQ{Dm(}=Y82jJrw3R~^W#35>H@6DnKU4m&=w`Pr1
zd5vqp1Z1NQFb}wV0H8)kehp7jH!7(c156Hs
zom~y_-8NW#^nc+^0q}|j^o^<>ISp~0Fb$s2l=wU~H#rn8r%w6xn6m&NnIDp33+ri2$VxWQZ7BsEZa2r%NSHY0`M?x-lB(6L=K;8=dhX}qd%KvVbkUK;yS
zO7mOvRqiJq6Ev}(Ty)jKh05VcM2ihT{I|rwn!3|o-NobF2Px$vuq4fg#Rqo=CIcNaX|hA$jUU7vDCNa;jB6qU|7abNF|J(
z5ozKZ3p5qEWAqj!gh{$$dK|5upqrNus*xkl`Z~p{HmyiuF+EEpRv*XQ=~=DTiyb99>P`>2qJ%P+T3M;tPbJ>7CoJ4y)!IAg&Ib*bF>fWBZ+Eb(VD#n~
z7{C@7Q0TNkx-{7{6^_opOrSd%~ewm1)B*`&aecAtF9t3mVx3S0roBUpF&&qbq9|y>dzJa3A5Gy!;I0;1qDHb3HAJJkbAWwcB57
zS?qjq$&o_izIF!e3K+dH(p@r`4EIW%7$V_oMD&_spFGAVI}a0UpsjB4
zGV5sOk|wixy_A9z!0nzkP}i?cd%qsd`TH#)Hr`*9n2de6?jC(&YWs3ZU%rjM5t!fi{GBc=~@*rGTn
zQ4qVl`#`6|#j)5ZSX!J%DnNs#su93InyMH-7s3s&K>+v`=?o04FF=N-m=v^@|
zbhlmU-=4_EiV_`70%^qdK;oqN8%Q5}p~I(^yQap}Fw$5IEx2C*LU3>Mc&<4{N4H~;
zlnFW90J$L^>t2LfE+C9&Irs77MUMv^h2r>Ogg3h}=$ImfmM9q%dZ&3#2^U^<^WM@_
zRWHM6ehCY8H3do?A^&E;m~gk-?n_M1J5r4&0fd18$Z6l^C_Kcw3lsgG+&`4nSWu8=c^ogySmV
z_Jskv4Mq<~zI^k8@qC7PZY#!mWDQ(u%ZsMd9P#TK{e={CKzHD#3X|xmD(vs(
z>AW&R(0djtuXo^4De5xQPv93Z$Ue>OE(r5r1~LS83x`{-Bg4_9ffEK
zzhp8FC-@J;Dn;ksk@ja<1>MUs;D4j!|9$=wUrz${rz*V);MYrF3)ozL-STI!ZRI0uHFJoh)=)mx9D$uMgF+tcny!5Bw+ocWP
z>EgxT;SnkeW!kc!uDuKz2*=E&|sU2>hajrB=9q
zqL#lt#j=8eP%tQ<3m#XXJ!E~FTz@7axDBb^ddXCa?I!dDX81)-(g~2Wp*bah2H9#i
zW+0}gfOy{t>#{G;@~+P3<{`RM8Q40Ft?wrnU$}fk%;f}VbLs{En?}_yBa{WCu}XB{
zSV0Kl{WN1o6G_0mGhf3N)~D6vyUw{y_~
z@BFtT^(sk@iVKhsYgrLxHX=wf_3i^5idGmj2m^E!PdtJ)h=bsf3Dglww8-YhqZ68w
z7twrepgBUL1o_->Km`AhRzQAu$6*QrPAg{^xFnRy^9(t3=nFdPJTs?je^np
zVPKnpJw{*5?8;(4k;8t|WMnV=Wxi_(Bm>e}u5|f5QLL)6rH&;CA+Ur(nqc&z$S2}2
zw_>9oNc+j_^+dhy9FOayq@T4fr-;40pt*b2IP2adPpI<`WYo=8PST5HriKIS)PoN1
zF`|%56(O%)G_hsxZ;5Y^oLv1Nz0?blB&8RD6g?=rI&
zK+o!@0Gh)IC$)ipJ#Wwm|uPsg9I~c4@ht;D0vv7KjaD+$EO-54-6J
z6#W$Ful`%`0Ol92IhkRV{S#Cg)1yZ^f0Q;YY}y?Inf-{=2D5U~cd}dPy{Ku5{H|y!
zBzp6qh$zQ334VtH);%y_TS45h8{=Ed58oKX4W57A079;gbHiNHsjFP%z0`0EDCj*5
zn3N#w2(-Z7Wr$Tj!}8tqRp)?-6ZW?(A362VdxI!+T#v{)@_-{K7zVb~ZHETvl+lGz
z5g%_A+5xt7(JnheN-}rHd3^as%OFT!STWT<7VCxw?SrtBX4My`W*c>sYB>O
zAcZZtVqhOdvE|3v$W0=~r$8(Rd;y^9_klB@Iy~;Qw^5$h;v!3KxnZQe3?;8pKvsX6
zT*4N8R0Sc}%`kc&3~couGPK#osg|mZ1N*FBhZLOBmCr%TqFHhOXfLn6Zf1`-q<`pv@z?6x}WvtYxXkyH9F8$7K
zmat4#mc30{e6sK|_x}ZzIcHX-w=>pQU*`it>9e9lMx>Q2doc5hNGxmUzg*7MF&8fs
z869+|lXh^_KrL|bs;c&L8^%FD+@oeu!v*}(K{o~nJ!nA60dJO~sg}1kPtY6Rr}bm7!t4<5^x!`J~K(Ac3kS}C?d;8qr}dPk7x`CH)tb6kGTsn
z-aA?eV+ki_lTw9)dVwbbbTcbCw}0`Ho4bIb##|e>$V@nJux311Qm7mkKj=;w?=p;)
zOc7zPh*(xcr0=MIUfI9QxS>ZRU65)Ujg3J4E_t9U-I+DoIdR`gtWc@&KX@%HU>f!V
zd5+NB{D3_8fxd07uwV1@1^#MVK;+*e$mW(Omg+(V!i}RCll4LbZc0MLV=1yfSGf0
zXOxR>i#Rj_$#CZlT3ZwtouxB!gbz>om`(iLZ`9kVz$ABjbP|W*LyEMfK4~Q$R7C`;
z;tf^#=uiaY`uA+Gp5GSM#;vQ*&Gr6R#vCf#Me&+yHW$lR05Y9X)&OE()zxTH0yb*gwG9oc~moRn?jX~q&&2#J3
zgz+1^M85rRWMEOuTf|8b7>0xu1}$DDIj52mOYN_tSc$jXtcWZxyJ%vg!7!@M&p#Qt
z{csA9^aV(00WuRc`XCU#T})`wc7OWex*>%0D52*N8*Dlf)X%I{UyrY1i!)&X`4I7?r|n
zCn~)RI|@~+FIEyB$f2HB;6pT8Y$g_juqX6D*Mnb3JdbuOm0nCKntyH}A|sGr+3Bzf8a_I1S5-hro*m?8S3jMy|B|2Kc-u$9`5(X#?dWnIA8!BWLGS
z6GCay%!vYbO#^MeAeO%%3|^f9Z)XQSkAG7f33+V3ey*l=OcNON-skB1q@h93C-4`9
z{r^C&s?xGNuxa*}dojv3oev~g^1XYuBimEs-FQZdlGG|!v``!5?_Pggl**7;*&y2f|RNlu=
z&BW??m&I8=ucpzDkHE845xE*qqf(|DMr0rbd|X%~AB&a+V$W2_obA#CQnyHD{wQ0{
zDbf-qZbH~MA(opE230o?;9k!abMt}A1v`sr{Yo}8{)>k`y=#TE8t)-vYIoKc%^UwM
z3V)BIq`WIo(wyEqO&MHgcTV4b^Dr-UZLW6BzF(CYG(>zmB7!L)ib?BPGRw?X+Ro`F
zrhrCXaM|_U+mFYg!tI4eN}wg{*WAFYn=YmVVCaqZ2;3`NrDg+_x7jV0lOI=apIsRf3QN)Ju;F^RW6Lnp;vN0=s$%`;9Sw6LJm|3FH?w{>*e
z@^yO+v!;%0I%EoGXx{kTU2O19iq__|f2G|iMcpGs-7l3nCFwH@wE9_ewcD^m*%^GL
za`Pp#-dm>H7iDHIud%%)Ks@3}dBjuo=;$af|FK+tAt{r(8ErxC`NsIhhL|2?@trAm
zybcU|Ne{Gy#`drtb1jMufV&D#^D@6rEcNvEc)<@IYRFH#Q{*ey@63$_Nx+6A;6k*%
zeuL&=o;MQ^B;r|h$MIbEp@jxDI^Y~t1
zCyH8+`(9|kiPk>9po6kwOxPb3OS
zZ;vx?X0RV4-4YDQhjsofq6JBKz+QoFU?<#|o{X^Bc>g9X(3paw!6k^29TK_)8~jie
zk*oPhDqCSiE*oGSA9i`xPK&Q1R7e?vbD9#k=-xLUV-ZD9H8^1Vf4BAU3h>O2kl}@!
zFo^07T^>u$EAprnReE0)JN!q$1fr^()4>k
zq@R({&q%+S&6$}f`EN>uG_4y+WlASaKw>P%gZfafR4rGTM|3yV{%-8@-B{%y&ge5*zYnT#oDf^r*k#w)A9rOknZN0^17j%aV=r2ESa}~y{6aDZ7~@)_
z@16`%ww4jkMIf0IikJsEzv%_ktDKI3Ei{SqwC0#vnmK-`$f$1IT{oEQiV1nBC_}9x
zL#-i0UAwTf5X&1lhkP=9tM@IWs<0qO9Ci8R2q;1yX#m4gX@$I1N7KUtcn|?RDFHka
zQt*IiJ{X#shz%i+6wFz6|pji8$Sg`z2y2b)>Hz@W;kIY2!J2FEW
zqxW_lYq5cGC@$fzN705w4Qxxy2^i*t6@Z)jdI7B}rx@mKDuFY)Wr_{br6)nL%x3O7
zJT+h|i%uEp9vSL>naqDMt}#oNDQ4f4dv-=ODzx`jG#+IX#(*6V{{bi}{%*%}aA~3E
zcB?kSfs_+^5Idh;bGQT~C|qbRBS&8!sfH^i!yy+TApAka`-2Lm5l|2OKld59`i9>0
zzgX=DTYKL%CBL`c?IbcyWwsq2%!3H#NeSke17y<}?lKr2a?3EzSuf?!;R{~et*xbP
z3B&FfZ6#Y$(6JumST7r}$i^i+#^o^Fuy)(oDwfZd1Zf*FkyscKzMSd^%Gc&N2T1w@
zB=i6&{OD-c1N(^(HPtf*MGtmevi4PEdsCXZ!i1_t9%eNX_dd*m#k9usu&hvd0BrB?DUmLUl6jfOw(07RP*cRKA9(0=vyUQ5)@n&C=
zY*^M?YuRwtgb3DzC;*U^@unjd|Amwbe$@&{f306CT79o|Gqm
zPj5=EJ0x&`M?JoW*_rST-faDfuNV!X8Ab_iP2a${aAC0}V6!FQvL&!Xs*iSCyRsLs
zR;X{TLV-oqF?{SSnR9ozXJHmqUaI0UHIW
zdV^F&o~i^UwXjp1ce&iV?Ze>8q<}tN6k#t)l2RyeLHJ4zkQRWscZfka%H23?G;}
z$Mf^cOk6k~s1!&0Rxm;sQ3P*XCyMa^8bGQhPM1_4ndD>~*be!ynW+oEYB477$uaHo
z>l`~OGdZVOYRyToDpK?Uj}nF2Am9fW0X$gD+0rgSQkTcM-aHeUd^+=dXXk}80aSVv
zy}+YSp|*Vmm>blSb9v-hjk1k+{eX7FWu1&`iyizMyjnkl%KfxEC4AC8Y0z^*
z<;qKPsd{s%QZ4{q68FUdR_7~ZVnj7*>#xVRzTM4A#%+stt`#PyDxl~E9wdbt;Qa(2
z!9>CDScK7v2FGoO52aq@B@+7<;g8#ktc#`J$i|IQ_9LpvGNoQt{w(OpRyT$fWi>dO
zym_mseJwWn%Mo7x-Mqf0ja(e+<NADDVTOe8S(ZKL7MtC#p;-+1hHMVo&NGQ{oOEaS<4
znzA!sJ2~~M;Xd?9lQ?E^FcHpor|GT2?Yokd
zF;|(nW?;=hd=OmcHhT!>6RBno>D>kbXn80vWL<;xFtvUw{Gf<
zc@Z~fOU);nHk)qZs5{8GOn`_(0KRD3i&dU-q+F67
zlz94lw4uDOM&QF~PilEf{4Z76irFu5oXc=yoO{It>~pT8p&EBz9d8m8{`VLt;KJj
zeu$}T)w~0t6ziGnT2}2f@V)6HM)7LOd!vLDDacQmVYdHDy;Uu1u~1A+
zd)$Frr}+_l6)x&Dwwht8r+3LxE7cmOy5G?#TrnXY*qLEWqoQ0bt1A1nm^d!7yma`Z
zk0M7mJ}1`ge1$)yKNh-0X-&!{x2_D6{w`{?PKW{So}r%aBc;ZS%mY^h8jHgh*B%%47
zRcC|HbH(=D0>KOw)`vQqlx+E-1mlf_dXovM-Unp
zW1q)t`3c-bxdyEY)E?7NC_G`WmV-!XivR
z!unV>*hRmxu1Sh^9_KlKZ@lTmv&h#Dp$?NpM2ige`wVf)nH-fZe(N|nmOT7Ew7qvo
zV;%me$_~xmNffP-Kt};@@#$%DN7aBOwYwnX
z9Uy_*?l0ST7I*QPzuB=peyj5MQ|Ij`9QgQ9IGeE)ef~Qfj}^dRZa)&b-RynS{bf6k
zSXkO<#yIR1XINSX6xWCrk;qlBLI0py*(|M#XUZ1(c;!k{C7>127cxxZ^#G2Jx^yH9YYDY6~9%0d(68tN!oF2MjC*Iuc#1cQH0Qm!z%
z!r>Qv70Z-3Ly{JKS1_Jmte){TSjvxU83AClW`#IqCiir=XV@*Bmr96?*SN7hzzuz>
zzNdXGopOXOT$ac{#Tm+BUn3r
zrPDBqbukkjeo=5nT{s#rQW`s
zVqs(L=vlQI0U_8HS26p$IQ8lk(}BfnK9(-UU_KcAam-~XtquRWl+xXl&e12(2T`r9
z{0LJWxN>GyAAc#c|Hh6aOy>`r{w<69ol3Ov_Bl_&A4$xDliIB!^9)Sq(^qjWW)>Ox
zB_Cb0`#u!ELHxiKa{s|BtW-n&iI0a}h;dvG*K9Zk)F1mt74VWLGYdFw?VO;?KvhVG
z1eYA4`HPfy(vM~2aW8v_3Zf(JGP8@-?L?&tQ~}?!`9rWLlEVb`!`P(Vk8=IKrhs4O
zZ@AUzSd3OzH94Qd2Uv8}*rnpCSk5s+3f*aZ(3DXlO6QgSB|1Am+!+K?`-vJX<}b>?
z8t=k$Na-HoFHpb6s@@aHX@c4|HfhT@VQ3PO?~{@qJPjeZh_w*lD!n0?d~A{AuclLE
zR)WDw2exPBSy5qyYC1`%OfOw)I_YZz4u6#su1Hn0u*zWcpX;OnDfFrm
z)FcX8WotZ;5sjg_If|?v@DO*GgqU98Oni^}1-IV*l>ncQ=I`s3Lo)!cUQG!SPAy60
zuVoxIC)3CRgrCT0p}#_fj|E1!z74&z|8q`E3rW1m{c}wWB7Sxh$=h
z8~J5#5`qQ+d$w9|GXhFzvMtK28jcPPr!7(h4
zu|MxNI0x}qhT4=&qQ!E|K}rH-Fe7hkPFurk?K1AM=Q9C-i*zfncl3N2nAOrrtwOiD
z*NE@c6u_W5j0BZ`_vUG8d9^psqbmds=Ay}p64#ksNU|FM|V_lYoH;a1lVyFBtQ&j6=}YMp=kYcNJ|c9twB
zbSHtAl4eHl2i6cnkZeuSuP~6DU@GhdImPp9_xJTJ_0Y4*s;>nvW;%H(Z9juFcCt=^
z{d0ayL|x31YM`0PQa=%x)nW9W1A1d*QRnT$=Y65Ooj=hVsBS({!S{fMWQ~a|(>RMo
zkP&BzXO1&2
zG*mob$6)`KU5Crklm0!jU$)(q4DqG1fT*sig=Jk?tVNxYSPk;C=B0r|Tv%80%ZRhq
z2~}t5ZN6ifq>(}FqkpKa7h|fIU=klDaAm@6Vwb}c8jNwxN!DG5m)1I%nwIj5GC2+-
zdTbtZ4tMlWY4
zKE&}Th~T3CNx8z7O%EznKXvFkj_^7VRt5JZq5FQG^L-!NCQjbkeK1q6c!#pSZ2hb_
zSsYT^E*pN{sTH$t7d`-l%Wa>hpRR3TaVa=85*K)r4q7-Y@O$gvgn!bz-std@pBH*2
zt5RRi$ohu>$FVbeMwb5=_xzVtm27RR*R3&-p;hEP9s;%Pm)Wd-YLBC;7GZP$g>bf-
zgu)_Ftad?6$H@+AsWktTrLrvcs}I)5b*TP^SaShF^)~|DmeJ!P;O*=Z=jrE-p#;m`
zc8iSgy>9v@*S?#w#;=EF$zn&Zm!BoONLhEwuHYN+tBtA~?!4@Kg#zif=pRMt{y1pm
zXz7~rg66JYdl#>R+@<(8H}vX`SfB3swVZ$^Hq;oJy~d^+krVtAG&5w40X=pS*QOc}
z<>IBH+fKiM12eN)#rmumX3ZPe7y2n7jhy_eiKEmPVEyTTt*>(wt<9bEY`+?7RSx$!
zm`2(-s6lUgYsIH0y!KlEz9zq`V{&kp>XUlkl`CQgy1$Y**9Z7v?JDmZnHJ!SZv`zo
z3A=r71A4wo{$1YTUQ2|9kA!*=^Nozw>W5ObsIc%%LR8rD-c)!CfOB$f;N`WVwurQ6
z$a{Y+NL0s{oBx_D@{Isq@q#J^c~i%Z@YnElPUD%B_m$q_y|)N7bv#SYc|L1qW?g{5!hefCP^h025g%8tN;Pz)cM8~xMESb`L&E`Y(4(0imAx&ZUmY=Wvu$QU$
z2H>l_EMcFI0+*;Xlvxo?ry@Y`a5~m?+C<*k3ZNJ^MPDj9N+BAUQqWfcuAza6?zLEi
zRmfW<_3mfQ)4f^M#xjgN{=uXDG|1f}1&WzuyKe|i|7+{|)|bt;kJr5>qH}|e&l(^~
zbbHYI-or9OP&$#nT7ZLT3nn}7ABNpRGwu8=z3Je!$H5OL*goX===y{{&pbW*DfL^K
zfA;r1kJPHRfl2h^LN+YP*U6g9K0Bg3kWrFGl`N)
z;jM`o!G7>~;Q^JaqZ2`w(l@TurfgBgU?61qOJ{FfYoOP|r?U@Od;SiJ&bfo#
zrU(Nq*5(wJeQWH|e{j-nU*n{k5<4eLqUzgCH~d`~Yv7FO6R|0xvvjjhCI$qiQu!6;
zF>-CFUz-l{w{%GmS&rwO!ZV+C+1pe$+ZXdVQz)JH
zN~l5H_%cH0d-ttnYs{MTmCCS6D1vvt922}AxfYkF>6a^Ud0IA&T*Oi_nFz>$Kw>for=~qEy{b^GL5i)Wo
z!B!FIo52I#6_rz5j+Se2L^%SK;vpQaI!-0-89G;AHuDCTz
z%Dr|hQ1NwL_l?FMrX|ghqZKbfjW^UWovw!2OLK{bEr(4jd=X$WVe@;}UV#gw;oZ?q
z9gsAox|>nXx&*^z-~7Jqs>`Rn>JJ%+>>5gHZYQ<+R8dm334Ls!BlA9O+~YBgod0GBQ9#@;rZ4;l3Ktbko;=JI-02I*>P}BbQleoNiW~tp
zzhZ!DpBxk3?3}4CJ*MsSxtSdj)or~+1?3sp;v2_HxxKD#wlR(uinpsskGTNDx*q^$7`bv>Md#;M}OBrnS^ET6<
zN)8vqz`Ub;C`S8lch}zUE(`JtGtDR~&wK1%0Z)DQSmj!rek&sRdStkETH=vGc&*wE
z*2B_{6K^Ri%A4XATBR>Jf1j!T+4PC4&3*61Inqs*wCA>M?>|E?PYA=m$hyRHYEwCZ
zCM~dF_*m_f>w2T_^sVgx$?cWe-E$qEWOI_()KE)U%siaT{87?EYqCl|`W9ke<%d_B7LJfIsSz}iPPMz{5-qr
z1*ffRW2aufS6J9pelf3}AH#(5*g1I?(zm`#ZWFvNIfjmj$=O?8L3YS$d_sG?K!Obq
z&#po?Dp_`3z2r_hi|)iv6wC4nC|Z!raIxmy&6wK@5IdE+UIfZK;aYZ6o~vfJ%)5{i
zxK8R3r}omh&B-e%R8Hz-{V&jJYO83EkH~7@6M)nSqeicyPbjja;@+cI1i9}NLY_VH
z!e!XC=KInn;vlB6iS_{Ml7@qScIAX3le<^w(jKv~KLoDnEYe6-AHNSVU5It!w2yTj
z0{c%6)Lc8ad7ip^^{fw~y-?d>(d1OA6RGfis_lc2_<2B(RD`5LC6k+rt5<(ik&>r8BpU2z(
zDp=8{H<`x%ystBlK((U!PTi6@1s7Qr<=o&Q51$X@!_Qf_lsj
zm*9;iq{d}m!&cu-iTYMbd2x%-(%(fkfEpTE)aW|{8LnYcpWs+%{3Z<{4lR9Lyx24~
zm4`1L`7kX@j-0HxMbM_Q;XGsVnE(Y%fh5!<$gktW^h~P3ap)MASUJ;nRP}(XDctqv
z5!z!53Ao=P{gTlEP+bUC;5p9>LtIIWelc?c^I52SlcVY%w8trO=1`?>Ew(3)O_mjVB-XSm@Hyxo~I+_e#FbQimD21gYrZFE*};