From ef0c87482160a503ec27b7a36bbea2a480576c31 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Sun, 29 Oct 2023 19:14:33 -0400 Subject: [PATCH 1/5] matlab works now --- data/{650^oC.CSV => 650^oC.csv} | 0 data/{700^oC.CSV => 700^oC.csv} | 0 data/{750^oC.CSV => 750^oC.csv} | 0 data/{800^oC.CSV => 800^oC.csv} | 0 pyproject.toml | 7 +- repro/data_io.py | 7 + repro/main.py | 12 +- repro/mixtures.py | 45 ++++++ repro/simulator.py | 150 ++++-------------- ...Regolith_Spectra_Generator_and_Retrieval.m | 32 ++-- 10 files changed, 111 insertions(+), 142 deletions(-) rename data/{650^oC.CSV => 650^oC.csv} (100%) rename data/{700^oC.CSV => 700^oC.csv} (100%) rename data/{750^oC.CSV => 750^oC.csv} (100%) rename data/{800^oC.CSV => 800^oC.csv} (100%) create mode 100644 repro/mixtures.py diff --git a/data/650^oC.CSV b/data/650^oC.csv similarity index 100% rename from data/650^oC.CSV rename to data/650^oC.csv diff --git a/data/700^oC.CSV b/data/700^oC.csv similarity index 100% rename from data/700^oC.CSV rename to data/700^oC.csv diff --git a/data/750^oC.CSV b/data/750^oC.csv similarity index 100% rename from data/750^oC.CSV rename to data/750^oC.csv diff --git a/data/800^oC.CSV b/data/800^oC.csv similarity index 100% rename from data/800^oC.CSV rename to data/800^oC.csv diff --git a/pyproject.toml b/pyproject.toml index 8858ea4..f057ddb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,11 +16,14 @@ authors = [ license = {text = "CC-BY 4.0"} dependencies = [ "matplotlib", - "numpy", - "pint", + "numpy>=1.20.0", "scipy", ] [build-system] requires = ["setuptools"] build-backend = "setuptools.build_meta" + +[tool.ruff.format] +# use single quotes for non-triple quote strings. +quote-style = "single" diff --git a/repro/data_io.py b/repro/data_io.py index 2c7266d..b49fddb 100644 --- a/repro/data_io.py +++ b/repro/data_io.py @@ -2,3 +2,10 @@ Tools for reading data files included in this repository. ''' +import csv +from numpy.typing import ArrayLike +from pathlib import Path + + +def import_reflectance_measurements(band: float, path: Path) -> ArrayLike: + pass \ No newline at end of file diff --git a/repro/main.py b/repro/main.py index 183ab50..5493f53 100644 --- a/repro/main.py +++ b/repro/main.py @@ -3,7 +3,17 @@ This module executes the simulation and runs the plotting functions. ''' +import numpy as np + + +def main(): + # define wavelength range and spacing for simulations + WLS = np.linspace(1, 4, 601) + + # set the total number of simulations to run + num_spectra = 100 + if __name__ == "__main__": - pass + main() \ No newline at end of file diff --git a/repro/mixtures.py b/repro/mixtures.py new file mode 100644 index 0000000..7345839 --- /dev/null +++ b/repro/mixtures.py @@ -0,0 +1,45 @@ +"""mixtures + +Definitions of end member mixtures. +""" +from collections import namedtuple +from dataclasses import dataclass, field + +Range = namedtuple('Range', 'min max') + + +@dataclass +class Endmember: + name: str + density: float + grain_size: float + water_ppm: float + abundance: Range = field(default_factory=Range) + + def __init__(self, min_abundance: float, max_abundance: float): + self.abundance = Range(min_abundance, max_abundance) + + +class Regolith(Endmember): + density = 1.8 + + +class MatureHighlands(Regolith): + grain_size = 32e-6 + + +class Pyroxene(Endmember): + density = 2.8 + grain_size = 69e-6 + + +class HydratedMorbGlass(Endmember): + """Hydrated mid-ocean-ridge basalt (MORB) glass + + Laboratory reflectance data from step-wise heating experiments are used + to simulate varying levels of hydration in lunar regolith. + """ + spectrum_label = 'MORB D38A' + density = 2.8 + grain_size = 69e-6 + diff --git a/repro/simulator.py b/repro/simulator.py index e9bd706..84cfe8b 100644 --- a/repro/simulator.py +++ b/repro/simulator.py @@ -1,139 +1,57 @@ -'''simulator +"""simulator This module emulates the behavior of the following MATLAB functions: src/Hapke_Inverse_Function_Passive.m src/Hapke_Lidar_R_Function.m src/Hapke_Lidar_SSA_Function.m src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m -''' +""" from collections import namedtuple -from dataclasses import dataclass, field from functools import partial from typing import Any, Callable import numpy as np +from numpy.typing import ArrayLike from scipy.optimize import fmin -DEFAULT_WLS = np.linspace(1, 4, 601) +HFunction = namedtuple('HFunction', 'H H0') -Range = namedtuple('Range', 'min max') - - -@dataclass -class Endmember: - name: str - density: float - grain_size: float - water_ppm: float - abundance: Range = field(default_factory=Range) - - -class HydratedMorbGlass(Endmember): - '''Hydrated mid-ocean-ridge basalt (MORB) glass''' - - spectrum_label = "MORB D38A" - density = 2.8 - grain_size = 69E-6 - abundance = Range(0.0, 0.3) - - -class Regolith(Endmember): - density = 1.8 - grain_size = 32e-6 - - -class Spectrum: - - def __init__(self, reflectance_data, species: Endmember, grid=DEFAULT_WLS): - self.reflectance = reflectance_data - self.species = species - self.grid = grid - - def ssa(self, phasing, - emission_angle=0, - incident_angle=30, - phase_angle=30, - filling_factor=0.41): - return reflectance_to_ssa( - self.reflectance, - self.grid, - phasing, - emission_angle, - incident_angle, - phase_angle, - filling_factor, - ) - - spectrum_label = "MORB D38A" - density = 2.8 - grain_size = 69E-6 - abundance = Range(0.0, 0.3) - - -class Regolith(Endmember): - density = 1.8 - grain_size = 32e-6 - - -# SIMULATION TOOLS -DEFAULT_WLS = np.linspace(1, 4, 601) - - -class Spectrum: - - def __init__(self, reflectance_data, species: Endmember, grid=DEFAULT_WLS): - self.reflectance = reflectance_data - self.species = species - self.grid = grid - - def ssa(self, phasing, - emission_angle=0, - incident_angle=30, - phase_angle=30, - filling_factor=0.41): - return reflectance_to_ssa( - self.reflectance, - self.grid, - phasing, - emission_angle, - incident_angle, - phase_angle, - filling_factor, - ) def ordinary_least_squares(x: Any, y: Callable, yx: Any): - '''Ordinary Least Squares Function. + """Ordinary Least Squares Function. y (Callable): The estimator function. x (Any): The argument to y. yx (Any): The observation. Must be the same type as the result of y. - ''' + """ return sum((y(x) - yx) ** 2) def angular_width(filling_factor): - '''Angular width parameter, see Equation 3. + """Angular width parameter, see Equation 3. Args: filling_factor (float, optional): Filling factor. - ''' + """ return (-3 / 8) * np.log(1 - filling_factor) -def backscattering(h, g): - '''Backscattering function, see Equation 2. +def backscattering(h: float, g: float) -> float: + """Backscattering function, see Equation 2. This describes the opposition effect. Args: h (float): angular width parameter. g (float): phase angle in radians. - ''' + """ return 1 / (1 + (1 / h) * np.tan(g / 2)) -def bidirectional_reflectance(SSA, P, mu, mu0, B): - '''Bidirectional reflectance, see Equation 1. +def bidirectional_reflectance( + SSA: float, P: float, mu: float, mu0: float, B: float +) -> float: + """Bidirectional reflectance, see Equation 1. R = (ω/4) (μ₀ / (μ + μ₀)) {(1 + B)P + H(ω)H₀(ω) - 1} @@ -157,11 +75,11 @@ def bidirectional_reflectance(SSA, P, mu, mu0, B): B (float): backscattering function Returns: - _type_: _description_ - ''' + float: _description_ + """ - def h_func(SSA, mu, mu0): - '''Ambartsumian-Chandrasekhar H functions. + def h_func(SSA: float, mu: float, mu0: float) -> HFunction: + """Ambartsumian-Chandrasekhar H functions. Computed using the approximation from equation 8.57 from Hapke (2012). @@ -169,7 +87,7 @@ def h_func(SSA, mu, mu0): SSA (float): single-scattering albedo, aka ω mu (float): cosine of emission angle mu0 (float): coside of incident angle - ''' + """ gamma = np.sqrt(1 - SSA) r0 = (1 - gamma) / (1 + gamma) @@ -183,7 +101,7 @@ def h_func(SSA, mu, mu0): h5 = r0 + h3 * (h4 / 2) H0 = (1 - SSA * mu * h5) ** -1 - return H, H0 + return HFunction(H, H0) H, H0 = h_func(SSA, mu, mu0) @@ -195,30 +113,30 @@ def h_func(SSA, mu, mu0): def reflectance_to_ssa( - Refl, - WLS, - P=0.15, - emission_angle=0, - incident_angle=30, - phase_angle=30, - filling_factor=0.41, + Refl: ArrayLike, + WLS: ArrayLike, + P: float = 0.15, + emission_angle: float = 0, + incident_angle: float = 30, + phase_angle: float = 30, + filling_factor: float = 0.41, ): - '''Convert reflectance spectrum to single-scattering albedo (SSA) + """Convert reflectance spectrum to single-scattering albedo (SSA) Uses scipy.optimize.fmin (equivalent to MATLAB fminsearch) to minimize ordinary least squares distance between SSA obtained from the supplied reflectance, R, and the SSA from the estimated reflectance at each sample point in WLS. - Hapke, B. (2012). Theory of reflectance and emittance spectroscopy (2nd ed.). - Cambridge University Press. + Hapke, B. (2012). Theory of reflectance and emittance spectroscopy + (2nd ed.). Cambridge University Press. Equivalent to src/Hapke_Lidar_SSA_function.m Default parameter values replicate src/Hapke_Inverse_Function_Passive.m Args: - R (array[float]): Bidirectional reflectance, see Equation 1. - WLS: simulation sample grid? + R (ArrayLike): Bidirectional reflectance, see Equation 1. + WLS (ArrayLike): simulation sample grid? p (float, optional): Scattering phase function. Defaults to 0.15 for ansiotropic scattering on the modeled mean particle phase function for lunar soil (Goguen et al., 2010). @@ -229,7 +147,7 @@ def reflectance_to_ssa( phase_angle (float, optional): Phase angle in degrees. Defaults to 30. filling_factor (float, optional: Particle filling factor. Defaults to 0.41. - ''' + """ mu = np.cos(np.deg2rad(emission_angle)) mu0 = np.cos(np.deg2rad(incident_angle)) g = np.deg2rad(phase_angle) diff --git a/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m b/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m index a14a054..aa5b0c4 100644 --- a/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m +++ b/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m @@ -27,18 +27,15 @@ dMORB=69E-6; %Low wavelength portion of MORB spectrum (<1.5 microns) - -Morb_D38A_LowLam=dir('Morb_D38A_Low_wavelength.txt'); -Morb_D38A_LowLam=Morb_D38A_LowLam.name; -Morb_D38A_LowLam=readtable(Morb_D38A_LowLam); +Morb_D38A_LowLam=readtable('../data/Morb_D38A_Low_wavelength.txt'); Morb_D38A_LowLam=table2array(Morb_D38A_LowLam); WaterPPM=[1522,762,176,22]; %Total water measured from step-heating experiments %Load in MORB step-wise heating spectra -MORB_D38A=dir('*.csv'); +MORB_D38A=dir('../data/*.csv'); %Use the 650C spectrum for the lower wavelength portion of the spectrum -MORB_D38A2=MORB_D38A(1).name; +MORB_D38A2=fullfile(MORB_D38A(1).folder, MORB_D38A(1).name); MORB_D38A2=readtable(MORB_D38A2); MORB_D38A2=table2array(MORB_D38A2); MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns @@ -55,7 +52,7 @@ clear MORB_D38A2 %Perform Normalization of MORB spectra at 650, 700, 750, 800 C for i=1:4 -MORB_D38A2=MORB_D38A(i).name; %start from 650 C spectrum +MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum MORB_D38A2=readtable(MORB_D38A2); MORB_D38A2=table2array(MORB_D38A2); MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns @@ -119,9 +116,7 @@ densReg=1.8; dReg=32E-6; %Mare Mature Endmember -MareMature=dir('Mare_70181_Spectra.txt'); -MareMature=MareMature.name; -MareMature=readtable(MareMature); +MareMature=readtable('../data/Mare_70181_Spectra.txt'); MareMature=table2array(MareMature); MareMatureRaw=MareMature; %Remove Water and organics Signature by fitting a linear continuum between @@ -136,9 +131,7 @@ SSAMareMature=Hapke_Inverse_Function_Passive(PassiveMareMature,0.81,WLS); % Mare Immature Endmember -MareImmature=dir('Mare_71061_Spectra.txt'); -MareImmature=MareImmature.name; -MareImmature=readtable(MareImmature); +MareImmature=readtable('../data/Mare_71061_Spectra.txt'); MareImmature=table2array(MareImmature); MareImmatureRaw=MareImmature; %Remove Water and organics Signature @@ -152,9 +145,7 @@ SSAMareImmature=Hapke_Inverse_Function_Passive(PassiveMareImmature,0.81,WLS); % Highlands Mature Endmember -HighlandsMature=dir('Highlands_62231_Spectra.txt'); -HighlandsMature=HighlandsMature.name; -HighlandsMature=readtable(HighlandsMature); +HighlandsMature=readtable('../data/Highlands_62231_Spectra.txt'); HighlandsMature=table2array(HighlandsMature); HighlandsMatureRaw=HighlandsMature; %Remove Water and organics Signature @@ -168,9 +159,7 @@ SSAHighlandsMature=Hapke_Inverse_Function_Passive(PassiveHighlandsMature,0.81,WLS); % Highlands Immature Endmember -HighlandsImmature=dir('Highlands_61221_Spectra.txt'); -HighlandsImmature=HighlandsImmature.name; -HighlandsImmature=readtable(HighlandsImmature); +HighlandsImmature=readtable('../data/Highlands_61221_Spectra.txt'); HighlandsImmature=table2array(HighlandsImmature); HighlandsImmatureRaw=HighlandsImmature; %Remove Water and organics Signature @@ -186,9 +175,7 @@ %Apollo 15 Pyroxene Endmember densPyrox=3.2; dPyrox=70E-6; -Apollo15_Pyroxene=dir('Apollo15Sample15555ReddishBrownPyroxeneB.txt'); -Apollo15_Pyroxene=Apollo15_Pyroxene.name; -Apollo15_Pyroxene=readtable(Apollo15_Pyroxene); +Apollo15_Pyroxene=readtable('../data/Apollo15Sample15555ReddishBrownPyroxeneB.txt'); Apollo15_Pyroxene=table2array(Apollo15_Pyroxene); Apollo15_PyroxeneRaw=Apollo15_Pyroxene; %Remove Water and organics Signature @@ -322,7 +309,6 @@ MeasuredWater(j,:)=MeasuredAbundances(1)/0.8428*WaterPPM(MorbSelection1)+MeasuredAbundances(2)/0.8428*WaterPPM(MorbSelection2); %retrieve two MORBs %Calculate total water error WatError(j,:) = MeasuredWater(j)-ActWatAmt(j); -j end %% Plot the results showing scatter around 1:1 retrieval and histogram of total water error MareError=WatError(1:end/2); From ac1bef062a24448c5ab172ba2726c31193e5e827 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Sun, 29 Oct 2023 20:44:51 -0400 Subject: [PATCH 2/5] more python --- repro/data_io.py | 11 ------ repro/main.py | 24 ++++++++++++ repro/mixtures.py | 45 --------------------- repro/spectra.py | 99 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 123 insertions(+), 56 deletions(-) delete mode 100644 repro/data_io.py delete mode 100644 repro/mixtures.py create mode 100644 repro/spectra.py diff --git a/repro/data_io.py b/repro/data_io.py deleted file mode 100644 index b49fddb..0000000 --- a/repro/data_io.py +++ /dev/null @@ -1,11 +0,0 @@ -'''data_io - -Tools for reading data files included in this repository. -''' -import csv -from numpy.typing import ArrayLike -from pathlib import Path - - -def import_reflectance_measurements(band: float, path: Path) -> ArrayLike: - pass \ No newline at end of file diff --git a/repro/main.py b/repro/main.py index 5493f53..e574fc6 100644 --- a/repro/main.py +++ b/repro/main.py @@ -4,6 +4,9 @@ functions. ''' import numpy as np +from pathlib import Path + +from . import spectra def main(): @@ -13,6 +16,27 @@ def main(): # set the total number of simulations to run num_spectra = 100 + # Total water measured from step-heating experiments + water_ppm=[1522, 762, 176, 22] + + # Low wavelength portion of MORB spectrum (<1.5 microns) + MORB_D38A_LowLam = spectra.reflectance_from_file( + Path('../data/Morb_D38A_Low_wavelength.txt')) + + # Use the 650C spectrum for the lower wavelength portion of the spectrum + MORB_D38A_MidLam = spectra.reflectance_from_file( + Path('../data/650^oC.csv')) + + # MORB step-wise heating spectra + heated_MORB_spectra = [] + for temperature in [650, 700, 750, 800]: + MORB_spectrum = spectra.normalize_to_wavelength( + spectra.reflectance_from_file( + Path(f'../data/{temperature}^oC.csv')), + norm_wavelength=2.6 # microns + ) + full_spectrum = spectra.concatenate([MORB_D38A_LowLam, MORB_spectrum]) + heated_MORB_spectra.append(full_spectrum) if __name__ == "__main__": diff --git a/repro/mixtures.py b/repro/mixtures.py deleted file mode 100644 index 7345839..0000000 --- a/repro/mixtures.py +++ /dev/null @@ -1,45 +0,0 @@ -"""mixtures - -Definitions of end member mixtures. -""" -from collections import namedtuple -from dataclasses import dataclass, field - -Range = namedtuple('Range', 'min max') - - -@dataclass -class Endmember: - name: str - density: float - grain_size: float - water_ppm: float - abundance: Range = field(default_factory=Range) - - def __init__(self, min_abundance: float, max_abundance: float): - self.abundance = Range(min_abundance, max_abundance) - - -class Regolith(Endmember): - density = 1.8 - - -class MatureHighlands(Regolith): - grain_size = 32e-6 - - -class Pyroxene(Endmember): - density = 2.8 - grain_size = 69e-6 - - -class HydratedMorbGlass(Endmember): - """Hydrated mid-ocean-ridge basalt (MORB) glass - - Laboratory reflectance data from step-wise heating experiments are used - to simulate varying levels of hydration in lunar regolith. - """ - spectrum_label = 'MORB D38A' - density = 2.8 - grain_size = 69e-6 - diff --git a/repro/spectra.py b/repro/spectra.py new file mode 100644 index 0000000..8de1e77 --- /dev/null +++ b/repro/spectra.py @@ -0,0 +1,99 @@ +"""spectra + +Definitions of end member mixtures and their spectra. +""" +from collections import namedtuple +from dataclasses import dataclass, field +import numpy as np +from numpy.typing import ArrayLike +from pathlib import Path +from typing import List + + +Range = namedtuple('Range', ['min', 'max']) +Spectrum = namedtuple('Spectrum', ['wavelength', 'value']) + + +@dataclass +class Endmember: + name: str + density: float + grain_size: float + + def __init__(self, reflectance_data: Path): + self.reflectance = reflectance_from_file(reflectance_data) + + +class Regolith(Endmember): + density = 1.8 + + +class MatureHighlands(Regolith): + grain_size = 32e-6 + + +class Pyroxene(Endmember): + density = 2.8 + grain_size = 69e-6 + + +class HydratedMorbGlass(Endmember): + """Hydrated mid-ocean-ridge basalt (MORB) glass + + Laboratory reflectance data from step-wise heating experiments are used + to simulate varying levels of hydration in lunar regolith. + """ + + spectrum_label = 'MORB D38A' + density = 2.8 + grain_size = 69e-6 + + +@dataclass +class Mixture: + water_ppm: float + abundance: Range = field(default_factory=Range) + reflectance: Spectrum = field(default_factory=Spectrum) + end_members: List[Endmember] = field(default_factory=list) + + +def reflectance_from_file(path: Path) -> Spectrum: + with open(path, 'r') as file: + arr = np.loadtxt(file, delimiter=',') + wavelength = 1e4 / arr[:, 0] # in cm^1, convert to microns + value = arr[:, 1] / 100 # convert percent to decimal + + # since we converted frequency to wavelength, need to resort data + sort_order = np.argsort(wavelength) + return Spectrum(wavelength[sort_order], value[sort_order]) + + +def normalize_to_wavelength( + spect: Spectrum, norm_wavelength: float, tolerance: float = 1e-3 + ) -> Spectrum: + wavelength = spect.wavelength + values = spect.value + + norm = np.where(abs(wavelength - norm_wavelength) <= tolerance)[0][0] + new_values = np.asarray( + [v / norm for v in values] + ) + return Spectrum(wavelength, new_values) + + +def concatenate(spectra: List[Spectrum]) -> Spectrum: + all_wavelengths = [s.wavelength for s in spectra] + all_values = [s.value for s in spectra] + + new_wavelengths = np.concatenate(all_wavelengths) + new_values = np.concatenate(all_values) + + # sort by wavelength + sort_index = np.argsort(new_wavelengths) + + return Spectrum(new_wavelengths[sort_index], new_values[sort_index]) + + +def interpolate(src: Spectrum, new_wavelengths: ArrayLike) -> Spectrum: + w, v = np.interp(new_wavelengths, src.wavelength, src.value) + return Spectrum(w, v) \ No newline at end of file From 7a660ab0fe63ead5c37440eccbb859e5621ae752 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Wed, 1 Nov 2023 03:52:13 -0400 Subject: [PATCH 3/5] notebooks! --- .github/workflows/lint.yml | 4 +- pyproject.toml | 15 +- repro/Hapke_Inverse_Function_Passive.m | 16 + repro/Hapke_Lidar_R_Function.m | 18 + repro/Hapke_Lidar_SSA_function.m | 16 + ...Regolith_Spectra_Generator_and_Retrieval.m | 348 ++++++++++++++++++ repro/main.py | 43 --- repro/results.ipynb | 150 ++++++++ repro/simulator.py | 26 +- repro/spectra.py | 63 ++-- ...Regolith_Spectra_Generator_and_Retrieval.m | 32 +- 11 files changed, 618 insertions(+), 113 deletions(-) create mode 100644 repro/Hapke_Inverse_Function_Passive.m create mode 100644 repro/Hapke_Lidar_R_Function.m create mode 100644 repro/Hapke_Lidar_SSA_function.m create mode 100644 repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m delete mode 100644 repro/main.py create mode 100644 repro/results.ipynb diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index e85501e..6bcbf3f 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -2,9 +2,9 @@ name: Lint on: push: - branches: [ "main" ] + branches: [main] pull_request: - branches: [ "main" ] + branches: [main] jobs: ruff: diff --git a/pyproject.toml b/pyproject.toml index f057ddb..cd3b000 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,20 +10,19 @@ description = """\ requires-python = ">=3.8" keywords = ["python", "nasa", "desci", "simulation", "lunar"] authors = [ - {name = "Philip Linden", email = "phil@moondao.com"}, - {name = "Daniel R. Cremons"}, -] -license = {text = "CC-BY 4.0"} -dependencies = [ - "matplotlib", - "numpy>=1.20.0", - "scipy", + { name = "Philip Linden", email = "phil@moondao.com" }, + { name = "Daniel R. Cremons" }, ] +license = { text = "CC-BY 4.0" } +dependencies = ["matplotlib", "notebook", "numpy>=1.20.0", "pandas", "scipy"] [build-system] requires = ["setuptools"] build-backend = "setuptools.build_meta" +[tool.ruff] +line-length = 80 + [tool.ruff.format] # use single quotes for non-triple quote strings. quote-style = "single" diff --git a/repro/Hapke_Inverse_Function_Passive.m b/repro/Hapke_Inverse_Function_Passive.m new file mode 100644 index 0000000..8b1d4d6 --- /dev/null +++ b/repro/Hapke_Inverse_Function_Passive.m @@ -0,0 +1,16 @@ +function SSA = Hapke(R,p,WLS) +mu=cosd(0); +mu0=cosd(30); +g=30; +h=(-3/8)*log(1-0.41); +B=1/(1+(1/h)*tand(g/2)); +for m=1:numel(WLS) +Refl=R(m); +y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1); +x=WLS(m); +yx=Refl; +w0=0.5; %Initial Guesses here +OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function +opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7); +SSA(m,1) = fminsearch(OLS, w0, opts); % Use ‘fminsearch’ to minimise the ‘OLS’ function +end diff --git a/repro/Hapke_Lidar_R_Function.m b/repro/Hapke_Lidar_R_Function.m new file mode 100644 index 0000000..a2444d9 --- /dev/null +++ b/repro/Hapke_Lidar_R_Function.m @@ -0,0 +1,18 @@ +function R = Hapke(w,p,WLS) +mu=cos(0); +mu0=cos(0); +g=0; +SSATot=w; +for m=1:numel(SSATot) + if SSATot(m)>1 + SSATot(m)=1; + else + end +end +h=(-3/8)*log(1-0.41); %for lunar regolith from Li et al 2011 +B=1/(1+(1/h)*tand(g/2)); +gamma=sqrt(1-SSATot); +r0=(1-gamma)./(1+gamma); +H0=(1-SSATot.*mu0.*(r0+((1-2.*r0.*mu0)/2)*log((1+mu0)./mu0))).^-1; +H=(1-SSATot.*mu.*(r0+((1-2.*r0.*mu)/2)*log((1+mu)./mu))).^-1; +R=(SSATot./(4*(mu0+mu))).*(p.*(1+B)+(H0.*H)-1); %Radiance Coefficeint \ No newline at end of file diff --git a/repro/Hapke_Lidar_SSA_function.m b/repro/Hapke_Lidar_SSA_function.m new file mode 100644 index 0000000..6311e03 --- /dev/null +++ b/repro/Hapke_Lidar_SSA_function.m @@ -0,0 +1,16 @@ +function SSA = Hapke(R,p,WLS) +mu=cosd(0); +mu0=cosd(0); +g=0; +h=(-3/8)*log(1-0.4); +B=1/(1+(1/h)*tand(g/2)); +for m=1:numel(WLS) +Refl=R(m); +y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1); +x=WLS(m); +yx=Refl; +w0=0.5; %Initial Guesses here +OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function +opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7); +SSA(m,1) = fminsearch(OLS, w0, opts); % Use ‘fminsearch’ to minimise the ‘OLS’ function +end diff --git a/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m b/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m new file mode 100644 index 0000000..ae2d9c3 --- /dev/null +++ b/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m @@ -0,0 +1,348 @@ +%% +clear all +close all +clc +%% +%Replace with location of spectra and script files +%cd 'C:\Users...' +%% Define wavelength range and spacing for simulations +WLS=linspace(1,4,601).'; +%% Input the total number of simulations to run. first 1/2 will be mare, second 1/2 will be highlands +NumSpectra=100; +%% Set abundance limits for spectral endmembers in mixture +LowAbunReg=0.0; +HighAbunReg=0.25; + +LowAbunPyr=0.0; +HighAbunPyr=0.25; + +LowAbunMORB=0.0; +HighAbunMORB=0.3; + +rng(0,'twister'); +%% Load in optical constants or reflectance data and calculate SSA using Hapke model + +%Density and grain size of MORB glass +densMORB=2.8; +dMORB=69E-6; + +%Low wavelength portion of MORB spectrum (<1.5 microns) +Morb_D38A_LowLam=readtable('../data/Morb_D38A_Low_wavelength.txt'); +Morb_D38A_LowLam=table2array(Morb_D38A_LowLam); + +WaterPPM=[1522,762,176,22]; %Total water measured from step-heating experiments + +%Load in MORB step-wise heating spectra +MORB_D38A=dir('../data/*.csv'); +%Use the 650C spectrum for the lower wavelength portion of the spectrum +MORB_D38A2=fullfile(MORB_D38A(1).folder, MORB_D38A(1).name); +MORB_D38A2=readtable(MORB_D38A2); +MORB_D38A2=table2array(MORB_D38A2); +MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns +MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal +MORB_D38A2=MORB_D38A2(6481:11648,:); %only use the spectral region of interest +MORB_D38A2=flipud(MORB_D38A2); +MORB_D38A_MidLam=MORB_D38A2; + +%Normalize all MORB spectra to the reflectance at 2.6 microns at the +%lowest water amount and highest temperature (22 ppm, 800 C) +NormalizationValue=0.2363; +NormFactor=MORB_D38A_MidLam(end,2)/NormalizationValue; +MORB_D38A_MidLam(:,2)=MORB_D38A_MidLam(:,2)./NormFactor; +clear MORB_D38A2 +%Perform Normalization of MORB spectra at 650, 700, 750, 800 C +for i=1:4 +MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum +MORB_D38A2=readtable(MORB_D38A2); +MORB_D38A2=table2array(MORB_D38A2); +MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns +MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal +MORB_D38A2=flipud(MORB_D38A2); +%Normalize +NormFactor=MORB_D38A2(12752,2)/NormalizationValue; %normalize at 2.6 micron +MORB_D38A2(:,2)=MORB_D38A2(:,2)./NormFactor; +%Stitch with low wavelength spectrum +Morb_D38A_LowLam2=Morb_D38A_LowLam; +Morb_D38A_LowLam2(:,2)=Morb_D38A_LowLam(:,2)*(MORB_D38A2(6901,2)/0.15); +Morb_D38A_LowLam2=Morb_D38A_LowLam2(1:15,:); +MORB_D38A2=cat(1,Morb_D38A_LowLam2,MORB_D38A2(6902:end,:)); +% Interpolate data to 5 nm spacing between 1 and 4 microns +InterpMorb=interp1(MORB_D38A2(:,1),MORB_D38A2(:,2),WLS); +if i>1 % Use 650 C spectrum below 2.6 microns to isolate 3 micron feature changes + InterpMorb(1:321)=RMORB(1:321,1); +else +end +RMORB(:,i)=InterpMorb; +%Use Hapke model to convert from laboratory reflectance to single +%scattering albedo using reflectance and scattering asymmetry +% factor (p) of 0.81 (see manuscript for details on Hapke model) +SSAMORB(:,i)=Hapke_Inverse_Function_Passive(RMORB(:,i),0.81,WLS); +end +%% Interpolation of MORB Spectra +%Set range and resolution for interpolation +%(standard is 1 ppm spacing from 0 to 1666 ppm) With 30% maximum MORB +%abundance this gives a maximum total water abundance of 500 ppm. +WatInterp=linspace(0,1666,1667); + +for m=1:numel(WLS) + %For each wavelength, fit a line to the ESPAT values for each heating step spectrum + ESPAT(m,:)=polyfit(WaterPPM(1:4),(1-SSAMORB(m,1:4))./SSAMORB(m,1:4),1); %linear function + %For each wavelength, create new synthetic spectra from the linear ESPAT relationship + TotalWatInterpESPAT(:,m)=1./(ESPAT(m,1).*WatInterp+ESPAT(m,2)+1); +end + +% Check match of interpolated spectra with real MORB spectra +figure(1); +RealMorbIndices=[23,177,763,1523]; +for k=1:4 + scatter(WLS,SSAMORB(:,k),20,'filled'); + hold on +end +for k=1:4 +plot(WLS,TotalWatInterpESPAT(RealMorbIndices(k),:),'Color','black','linestyle','-','linewidth',2); +end +xlabel('Wavelength (\mum)','FontSize',16,'FontWeight','bold'); +ylabel('SSA','FontSize',16,'FontWeight','bold'); +ax=gca; +ax.XLim=[2.5,3.5]; +ax.LineWidth = 2; +ax.TickDir='in'; +ax.FontName='arial'; +ax.FontSize=14; +ax.FontWeight='bold'; +box on +grid on +%% Load in other spectral endmembers (regolith and pyroxene) +densReg=1.8; +dReg=32E-6; +%Mare Mature Endmember +MareMature=readtable('../data/Mare_70181_Spectra.txt'); +MareMature=table2array(MareMature); +MareMatureRaw=MareMature; +%Remove Water and organics Signature by fitting a linear continuum between +%2.65 and 3.8 microns +[V ind32]=min(abs(MareMature(:,1)-2.65)); +[V ind36]=min(abs(MareMature(:,1)-3.8)); +SandIQuad=polyfit([MareMature(ind32,1),MareMature(ind36,1)],[MareMature(ind32,2),MareMature(ind36,2)],1); +NewRs=MareMature(ind32:ind36,1)*SandIQuad(1)+SandIQuad(2); +MareMature(ind32:ind36,2)=NewRs; + +PassiveMareMature=interp1(MareMature(:,1),MareMature(:,2),WLS); +SSAMareMature=Hapke_Inverse_Function_Passive(PassiveMareMature,0.81,WLS); + +% Mare Immature Endmember +MareImmature=readtable('../data/Mare_71061_Spectra.txt'); +MareImmature=table2array(MareImmature); +MareImmatureRaw=MareImmature; +%Remove Water and organics Signature +[V ind32]=min(abs(MareImmature(:,1)-2.65)); +[V ind36]=min(abs(MareImmature(:,1)-3.8)); +SandIQuad=polyfit([MareImmature(ind32,1),MareImmature(ind36,1)],[MareImmature(ind32,2),MareImmature(ind36,2)],1); +NewRs=MareImmature(ind32:ind36,1)*SandIQuad(1)+SandIQuad(2); +MareImmature(ind32:ind36,2)=NewRs; + +PassiveMareImmature=interp1(MareImmature(:,1),MareImmature(:,2),WLS); +SSAMareImmature=Hapke_Inverse_Function_Passive(PassiveMareImmature,0.81,WLS); + +% Highlands Mature Endmember +HighlandsMature=readtable('../data/Highlands_62231_Spectra.txt'); +HighlandsMature=table2array(HighlandsMature); +HighlandsMatureRaw=HighlandsMature; +%Remove Water and organics Signature +[V ind32]=min(abs(HighlandsMature(:,1)-2.65)); +[V ind36]=min(abs(HighlandsMature(:,1)-3.8)); +SandIQuad=polyfit([HighlandsMature(ind32,1),HighlandsMature(ind36,1)],[HighlandsMature(ind32,2),HighlandsMature(ind36,2)],1); +NewRs=HighlandsMature(ind32:ind36,1)*SandIQuad(1)+SandIQuad(2); +HighlandsMature(ind32:ind36,2)=NewRs; + +PassiveHighlandsMature=interp1(HighlandsMature(:,1),HighlandsMature(:,2),WLS); +SSAHighlandsMature=Hapke_Inverse_Function_Passive(PassiveHighlandsMature,0.81,WLS); + +% Highlands Immature Endmember +HighlandsImmature=readtable('../data/Highlands_61221_Spectra.txt'); +HighlandsImmature=table2array(HighlandsImmature); +HighlandsImmatureRaw=HighlandsImmature; +%Remove Water and organics Signature +[V ind32]=min(abs(HighlandsImmature(:,1)-2.65)); +[V ind36]=min(abs(HighlandsImmature(:,1)-3.8)); +SandIQuad=polyfit([HighlandsImmature(ind32,1),HighlandsImmature(ind36,1)],[HighlandsImmature(ind32,2),HighlandsImmature(ind36,2)],1); +NewRs=HighlandsImmature(ind32:ind36,1)*SandIQuad(1)+SandIQuad(2); +HighlandsImmature(ind32:ind36,2)=NewRs; + +PassiveHighlandsImmature=interp1(HighlandsImmature(:,1),HighlandsImmature(:,2),WLS); +SSAHighlandsImmature=Hapke_Inverse_Function_Passive(PassiveHighlandsImmature,0.81,WLS); + +%Apollo 15 Pyroxene Endmember +densPyrox=3.2; +dPyrox=70E-6; +Apollo15_Pyroxene=readtable('../data/Apollo15Sample15555ReddishBrownPyroxeneB.txt'); +Apollo15_Pyroxene=table2array(Apollo15_Pyroxene); +Apollo15_PyroxeneRaw=Apollo15_Pyroxene; +%Remove Water and organics Signature +[V, ind32]=min(abs(Apollo15_Pyroxene(:,1)-3.2)); %only need to remove organics here +[V ind36]=min(abs(Apollo15_Pyroxene(:,1)-3.6)); +SandIQuad=polyfit([Apollo15_Pyroxene(ind32,1),Apollo15_Pyroxene(ind36,1)],[Apollo15_Pyroxene(ind32,2),Apollo15_Pyroxene(ind36,2)],1); +NewRs=Apollo15_Pyroxene(ind32:ind36,1)*SandIQuad(1)+SandIQuad(2); +Apollo15_Pyroxene(ind32:ind36,2)=NewRs; + +PassivePyrox=interp1(Apollo15_Pyroxene(:,1),Apollo15_Pyroxene(:,2),WLS); +RPyrox=PassivePyrox; +SSAPyrox=Hapke_Inverse_Function_Passive(RPyrox,0.81,WLS); +%% +%Place densities, grain sizes, and SSAs into arrays +rho=1000*[NaN,densMORB,densMORB,densMORB,densMORB,densReg,densReg,densReg,densReg,densPyrox]; +d=[NaN,dMORB,dMORB,dMORB,dMORB,dReg,dReg,dReg,dReg,dPyrox]; +ConstituentSSAs=[SSAMORB,SSAHighlandsMature,SSAHighlandsImmature,SSAMareMature,SSAMareImmature,SSAPyrox]; +%% Create spectral mixtures +for j=1:NumSpectra + %Use Monte Carlo method to determine abundances in each mixture + PCTsMORB(j)=((HighAbunMORB-LowAbunMORB).*rand(1)+LowAbunMORB); + %First half of spectra will be mare + if j2\u001b[0m \u001b[39mfor\u001b[39;00m spectrum \u001b[39min\u001b[39;00m heated_MORB_spectra:\n\u001b[1;32m 3\u001b[0m \u001b[39m# interpolate along our sampling grid of interest\u001b[39;00m\n\u001b[1;32m 4\u001b[0m interpolated_spectrum \u001b[39m=\u001b[39m simulator\u001b[39m.\u001b[39minterpolate_series(spectrum, WLS)\n\u001b[0;32m----> 5\u001b[0m ssa\u001b[39m.\u001b[39mappend(simulator\u001b[39m.\u001b[39;49mreflectance_to_ssa(interpolated_spectrum))\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/repro/simulator.py:174\u001b[0m, in \u001b[0;36mreflectance_to_ssa\u001b[0;34m(Refl, P, emission_angle, incident_angle, phase_angle, filling_factor)\u001b[0m\n\u001b[1;32m 169\u001b[0m y \u001b[39m=\u001b[39m partial(bidirectional_reflectance, P\u001b[39m=\u001b[39mP, mu\u001b[39m=\u001b[39mmu, mu0\u001b[39m=\u001b[39mmu0, B\u001b[39m=\u001b[39mB)\n\u001b[1;32m 170\u001b[0m OLS \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 171\u001b[0m ordinary_least_squares, y\u001b[39m=\u001b[39my, yx\u001b[39m=\u001b[39mvalue\n\u001b[1;32m 172\u001b[0m ) \u001b[39m# turn least squares into the form y=f(x)\u001b[39;00m\n\u001b[1;32m 173\u001b[0m w\u001b[39m.\u001b[39mappend(\n\u001b[0;32m--> 174\u001b[0m fmin(\n\u001b[1;32m 175\u001b[0m OLS,\n\u001b[1;32m 176\u001b[0m w0,\n\u001b[1;32m 177\u001b[0m args\u001b[39m=\u001b[39;49m(index),\n\u001b[1;32m 178\u001b[0m disp\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m,\n\u001b[1;32m 179\u001b[0m maxiter\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 180\u001b[0m maxfun\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 181\u001b[0m ftol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 182\u001b[0m xtol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 183\u001b[0m )\n\u001b[1;32m 184\u001b[0m )\n\u001b[1;32m 185\u001b[0m \u001b[39mreturn\u001b[39;00m w\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:747\u001b[0m, in \u001b[0;36mfmin\u001b[0;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, initial_simplex)\u001b[0m\n\u001b[1;32m 738\u001b[0m opts \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mxatol\u001b[39m\u001b[39m'\u001b[39m: xtol,\n\u001b[1;32m 739\u001b[0m \u001b[39m'\u001b[39m\u001b[39mfatol\u001b[39m\u001b[39m'\u001b[39m: ftol,\n\u001b[1;32m 740\u001b[0m \u001b[39m'\u001b[39m\u001b[39mmaxiter\u001b[39m\u001b[39m'\u001b[39m: maxiter,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[39m'\u001b[39m\u001b[39mreturn_all\u001b[39m\u001b[39m'\u001b[39m: retall,\n\u001b[1;32m 744\u001b[0m \u001b[39m'\u001b[39m\u001b[39minitial_simplex\u001b[39m\u001b[39m'\u001b[39m: initial_simplex}\n\u001b[1;32m 746\u001b[0m callback \u001b[39m=\u001b[39m _wrap_callback(callback)\n\u001b[0;32m--> 747\u001b[0m res \u001b[39m=\u001b[39m _minimize_neldermead(func, x0, args, callback\u001b[39m=\u001b[39;49mcallback, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mopts)\n\u001b[1;32m 748\u001b[0m \u001b[39mif\u001b[39;00m full_output:\n\u001b[1;32m 749\u001b[0m retlist \u001b[39m=\u001b[39m res[\u001b[39m'\u001b[39m\u001b[39mx\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mfun\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnit\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnfev\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mstatus\u001b[39m\u001b[39m'\u001b[39m]\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:899\u001b[0m, in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)\u001b[0m\n\u001b[1;32m 897\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m 898\u001b[0m \u001b[39mfor\u001b[39;00m k \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(N \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m):\n\u001b[0;32m--> 899\u001b[0m fsim[k] \u001b[39m=\u001b[39m func(sim[k])\n\u001b[1;32m 900\u001b[0m \u001b[39mexcept\u001b[39;00m _MaxFuncCallError:\n\u001b[1;32m 901\u001b[0m \u001b[39mpass\u001b[39;00m\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:620\u001b[0m, in \u001b[0;36m_wrap_scalar_function_maxfun_validation..function_wrapper\u001b[0;34m(x, *wrapper_args)\u001b[0m\n\u001b[1;32m 618\u001b[0m ncalls[\u001b[39m0\u001b[39m] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m1\u001b[39m\n\u001b[1;32m 619\u001b[0m \u001b[39m# A copy of x is sent to the user function (gh13740)\u001b[39;00m\n\u001b[0;32m--> 620\u001b[0m fx \u001b[39m=\u001b[39m function(np\u001b[39m.\u001b[39mcopy(x), \u001b[39m*\u001b[39m(wrapper_args \u001b[39m+\u001b[39;49m args))\n\u001b[1;32m 621\u001b[0m \u001b[39m# Ideally, we'd like to a have a true scalar returned from f(x). For\u001b[39;00m\n\u001b[1;32m 622\u001b[0m \u001b[39m# backwards-compatibility, also allow np.array([1.3]),\u001b[39;00m\n\u001b[1;32m 623\u001b[0m \u001b[39m# np.array([[1.3]]) etc.\u001b[39;00m\n\u001b[1;32m 624\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m np\u001b[39m.\u001b[39misscalar(fx):\n", + "\u001b[0;31mTypeError\u001b[0m: can only concatenate tuple (not \"float\") to tuple" + ] + } + ], + "source": [ + "ssa = []\n", + "for spectrum in heated_MORB_spectra:\n", + " # interpolate along our sampling grid of interest\n", + " interpolated_spectrum = simulator.interpolate_series(spectrum, WLS)\n", + " ssa.append(simulator.reflectance_to_ssa(interpolated_spectrum))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/repro/simulator.py b/repro/simulator.py index 84cfe8b..d36e94e 100644 --- a/repro/simulator.py +++ b/repro/simulator.py @@ -9,14 +9,24 @@ from collections import namedtuple from functools import partial from typing import Any, Callable + import numpy as np from numpy.typing import ArrayLike +from pandas import Series from scipy.optimize import fmin - +from scipy.interpolate import interp1d HFunction = namedtuple('HFunction', 'H H0') +def interpolate_series(data: Series, new_index: ArrayLike) -> Series: + x = data.index.values + y = data.values + f = interp1d(x, y, fill_value='extrapolate', bounds_error=False) + new_y = f(new_index) + return Series(data=new_y, index=new_index) + + def ordinary_least_squares(x: Any, y: Callable, yx: Any): """Ordinary Least Squares Function. @@ -113,8 +123,7 @@ def h_func(SSA: float, mu: float, mu0: float) -> HFunction: def reflectance_to_ssa( - Refl: ArrayLike, - WLS: ArrayLike, + Refl: Series, P: float = 0.15, emission_angle: float = 0, incident_angle: float = 30, @@ -135,8 +144,7 @@ def reflectance_to_ssa( Default parameter values replicate src/Hapke_Inverse_Function_Passive.m Args: - R (ArrayLike): Bidirectional reflectance, see Equation 1. - WLS (ArrayLike): simulation sample grid? + R (Series): Bidirectional reflectance, see Equation 1. p (float, optional): Scattering phase function. Defaults to 0.15 for ansiotropic scattering on the modeled mean particle phase function for lunar soil (Goguen et al., 2010). @@ -155,18 +163,18 @@ def reflectance_to_ssa( B = backscattering(h, g) w = [] - for m, x in zip(Refl, WLS): + for index, value in Refl.items(): w0 = 0.5 # initial guess, ω₀ # turn bidrectional_reflectance() into the form y=f(x) y = partial(bidirectional_reflectance, P=P, mu=mu, mu0=mu0, B=B) OLS = partial( - ordinary_least_squares, y=y, yx=m + ordinary_least_squares, y=y, yx=value ) # turn least squares into the form y=f(x) w.append( fmin( OLS, w0, - args=x, + args=index, disp=False, maxiter=10_000, maxfun=10_000, @@ -174,4 +182,4 @@ def reflectance_to_ssa( xtol=1e-7, ) ) - return np.concatenate(w) + return w diff --git a/repro/spectra.py b/repro/spectra.py index 8de1e77..a13c961 100644 --- a/repro/spectra.py +++ b/repro/spectra.py @@ -4,14 +4,14 @@ """ from collections import namedtuple from dataclasses import dataclass, field -import numpy as np -from numpy.typing import ArrayLike from pathlib import Path from typing import List +import numpy as np +import pandas as pd +from pandas import Series Range = namedtuple('Range', ['min', 'max']) -Spectrum = namedtuple('Spectrum', ['wavelength', 'value']) @dataclass @@ -21,7 +21,7 @@ class Endmember: grain_size: float def __init__(self, reflectance_data: Path): - self.reflectance = reflectance_from_file(reflectance_data) + self.reflectance = spectrum_from_file(reflectance_data) class Regolith(Endmember): @@ -53,47 +53,26 @@ class HydratedMorbGlass(Endmember): class Mixture: water_ppm: float abundance: Range = field(default_factory=Range) - reflectance: Spectrum = field(default_factory=Spectrum) + reflectance: Series = field(default_factory=Series) end_members: List[Endmember] = field(default_factory=list) -def reflectance_from_file(path: Path) -> Spectrum: - with open(path, 'r') as file: - arr = np.loadtxt(file, delimiter=',') - wavelength = 1e4 / arr[:, 0] # in cm^1, convert to microns - value = arr[:, 1] / 100 # convert percent to decimal - - # since we converted frequency to wavelength, need to resort data - sort_order = np.argsort(wavelength) - return Spectrum(wavelength[sort_order], value[sort_order]) - - -def normalize_to_wavelength( - spect: Spectrum, norm_wavelength: float, tolerance: float = 1e-3 - ) -> Spectrum: - wavelength = spect.wavelength - values = spect.value - - norm = np.where(abs(wavelength - norm_wavelength) <= tolerance)[0][0] - new_values = np.asarray( - [v / norm for v in values] +def spectrum_from_file(path: Path, name='response') -> Series: + spectrum = pd.read_csv( + path, + index_col=0, + delimiter=',', + header=None, + names=['wavelength', name], ) - return Spectrum(wavelength, new_values) - - -def concatenate(spectra: List[Spectrum]) -> Spectrum: - all_wavelengths = [s.wavelength for s in spectra] - all_values = [s.value for s in spectra] - - new_wavelengths = np.concatenate(all_wavelengths) - new_values = np.concatenate(all_values) - - # sort by wavelength - sort_index = np.argsort(new_wavelengths) - - return Spectrum(new_wavelengths[sort_index], new_values[sort_index]) + # since we converted frequency to wavelength, need to resort data + spectrum.sort_index() + return spectrum.squeeze() -def interpolate(src: Spectrum, new_wavelengths: ArrayLike) -> Spectrum: - w, v = np.interp(new_wavelengths, src.wavelength, src.value) - return Spectrum(w, v) \ No newline at end of file +def get_normalization_factor(data: Series, normalization_index: float) -> float: + """get value at the index closest to the desired normalization index + """ + wavelength = data.index.values + norm_index = np.argmin(abs(wavelength - normalization_index)) + return data.iloc[norm_index] diff --git a/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m b/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m index aa5b0c4..a14a054 100644 --- a/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m +++ b/src/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m @@ -27,15 +27,18 @@ dMORB=69E-6; %Low wavelength portion of MORB spectrum (<1.5 microns) -Morb_D38A_LowLam=readtable('../data/Morb_D38A_Low_wavelength.txt'); + +Morb_D38A_LowLam=dir('Morb_D38A_Low_wavelength.txt'); +Morb_D38A_LowLam=Morb_D38A_LowLam.name; +Morb_D38A_LowLam=readtable(Morb_D38A_LowLam); Morb_D38A_LowLam=table2array(Morb_D38A_LowLam); WaterPPM=[1522,762,176,22]; %Total water measured from step-heating experiments %Load in MORB step-wise heating spectra -MORB_D38A=dir('../data/*.csv'); +MORB_D38A=dir('*.csv'); %Use the 650C spectrum for the lower wavelength portion of the spectrum -MORB_D38A2=fullfile(MORB_D38A(1).folder, MORB_D38A(1).name); +MORB_D38A2=MORB_D38A(1).name; MORB_D38A2=readtable(MORB_D38A2); MORB_D38A2=table2array(MORB_D38A2); MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns @@ -52,7 +55,7 @@ clear MORB_D38A2 %Perform Normalization of MORB spectra at 650, 700, 750, 800 C for i=1:4 -MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum +MORB_D38A2=MORB_D38A(i).name; %start from 650 C spectrum MORB_D38A2=readtable(MORB_D38A2); MORB_D38A2=table2array(MORB_D38A2); MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns @@ -116,7 +119,9 @@ densReg=1.8; dReg=32E-6; %Mare Mature Endmember -MareMature=readtable('../data/Mare_70181_Spectra.txt'); +MareMature=dir('Mare_70181_Spectra.txt'); +MareMature=MareMature.name; +MareMature=readtable(MareMature); MareMature=table2array(MareMature); MareMatureRaw=MareMature; %Remove Water and organics Signature by fitting a linear continuum between @@ -131,7 +136,9 @@ SSAMareMature=Hapke_Inverse_Function_Passive(PassiveMareMature,0.81,WLS); % Mare Immature Endmember -MareImmature=readtable('../data/Mare_71061_Spectra.txt'); +MareImmature=dir('Mare_71061_Spectra.txt'); +MareImmature=MareImmature.name; +MareImmature=readtable(MareImmature); MareImmature=table2array(MareImmature); MareImmatureRaw=MareImmature; %Remove Water and organics Signature @@ -145,7 +152,9 @@ SSAMareImmature=Hapke_Inverse_Function_Passive(PassiveMareImmature,0.81,WLS); % Highlands Mature Endmember -HighlandsMature=readtable('../data/Highlands_62231_Spectra.txt'); +HighlandsMature=dir('Highlands_62231_Spectra.txt'); +HighlandsMature=HighlandsMature.name; +HighlandsMature=readtable(HighlandsMature); HighlandsMature=table2array(HighlandsMature); HighlandsMatureRaw=HighlandsMature; %Remove Water and organics Signature @@ -159,7 +168,9 @@ SSAHighlandsMature=Hapke_Inverse_Function_Passive(PassiveHighlandsMature,0.81,WLS); % Highlands Immature Endmember -HighlandsImmature=readtable('../data/Highlands_61221_Spectra.txt'); +HighlandsImmature=dir('Highlands_61221_Spectra.txt'); +HighlandsImmature=HighlandsImmature.name; +HighlandsImmature=readtable(HighlandsImmature); HighlandsImmature=table2array(HighlandsImmature); HighlandsImmatureRaw=HighlandsImmature; %Remove Water and organics Signature @@ -175,7 +186,9 @@ %Apollo 15 Pyroxene Endmember densPyrox=3.2; dPyrox=70E-6; -Apollo15_Pyroxene=readtable('../data/Apollo15Sample15555ReddishBrownPyroxeneB.txt'); +Apollo15_Pyroxene=dir('Apollo15Sample15555ReddishBrownPyroxeneB.txt'); +Apollo15_Pyroxene=Apollo15_Pyroxene.name; +Apollo15_Pyroxene=readtable(Apollo15_Pyroxene); Apollo15_Pyroxene=table2array(Apollo15_Pyroxene); Apollo15_PyroxeneRaw=Apollo15_Pyroxene; %Remove Water and organics Signature @@ -309,6 +322,7 @@ MeasuredWater(j,:)=MeasuredAbundances(1)/0.8428*WaterPPM(MorbSelection1)+MeasuredAbundances(2)/0.8428*WaterPPM(MorbSelection2); %retrieve two MORBs %Calculate total water error WatError(j,:) = MeasuredWater(j)-ActWatAmt(j); +j end %% Plot the results showing scatter around 1:1 retrieval and histogram of total water error MareError=WatError(1:end/2); From aac2c60bad70a47e9e7c6b5bcff49ebc16df2da1 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Thu, 2 Nov 2023 23:48:58 -0400 Subject: [PATCH 4/5] some cleanup --- repro/Hapke_Inverse_Function_Passive.m | 16 +- ...Regolith_Spectra_Generator_and_Retrieval.m | 170 +++++++++--------- repro/results.ipynb | 44 ++--- repro/simulator.py | 37 ++-- 4 files changed, 130 insertions(+), 137 deletions(-) diff --git a/repro/Hapke_Inverse_Function_Passive.m b/repro/Hapke_Inverse_Function_Passive.m index 8b1d4d6..543fcc9 100644 --- a/repro/Hapke_Inverse_Function_Passive.m +++ b/repro/Hapke_Inverse_Function_Passive.m @@ -5,12 +5,12 @@ h=(-3/8)*log(1-0.41); B=1/(1+(1/h)*tand(g/2)); for m=1:numel(WLS) -Refl=R(m); -y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1); -x=WLS(m); -yx=Refl; -w0=0.5; %Initial Guesses here -OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function -opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7); -SSA(m,1) = fminsearch(OLS, w0, opts); % Use ‘fminsearch’ to minimise the ‘OLS’ function + Refl=R(m); + y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1); + x=WLS(m); + yx=Refl; + w0=0.5; %Initial Guesses here + OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function + opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7); + SSA(m,1) = fminsearch(OLS, w0, opts); % Use ‘fminsearch’ to minimise the ‘OLS’ function end diff --git a/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m b/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m index ae2d9c3..d5cb330 100644 --- a/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m +++ b/repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m @@ -52,43 +52,43 @@ clear MORB_D38A2 %Perform Normalization of MORB spectra at 650, 700, 750, 800 C for i=1:4 -MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum -MORB_D38A2=readtable(MORB_D38A2); -MORB_D38A2=table2array(MORB_D38A2); -MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns -MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal -MORB_D38A2=flipud(MORB_D38A2); -%Normalize -NormFactor=MORB_D38A2(12752,2)/NormalizationValue; %normalize at 2.6 micron -MORB_D38A2(:,2)=MORB_D38A2(:,2)./NormFactor; -%Stitch with low wavelength spectrum -Morb_D38A_LowLam2=Morb_D38A_LowLam; -Morb_D38A_LowLam2(:,2)=Morb_D38A_LowLam(:,2)*(MORB_D38A2(6901,2)/0.15); -Morb_D38A_LowLam2=Morb_D38A_LowLam2(1:15,:); -MORB_D38A2=cat(1,Morb_D38A_LowLam2,MORB_D38A2(6902:end,:)); -% Interpolate data to 5 nm spacing between 1 and 4 microns -InterpMorb=interp1(MORB_D38A2(:,1),MORB_D38A2(:,2),WLS); -if i>1 % Use 650 C spectrum below 2.6 microns to isolate 3 micron feature changes - InterpMorb(1:321)=RMORB(1:321,1); -else -end -RMORB(:,i)=InterpMorb; -%Use Hapke model to convert from laboratory reflectance to single -%scattering albedo using reflectance and scattering asymmetry -% factor (p) of 0.81 (see manuscript for details on Hapke model) -SSAMORB(:,i)=Hapke_Inverse_Function_Passive(RMORB(:,i),0.81,WLS); + MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum + MORB_D38A2=readtable(MORB_D38A2); + MORB_D38A2=table2array(MORB_D38A2); + MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns + MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal + MORB_D38A2=flipud(MORB_D38A2); + %Normalize + NormFactor=MORB_D38A2(12752,2)/NormalizationValue; %normalize at 2.6 micron + MORB_D38A2(:,2)=MORB_D38A2(:,2)./NormFactor; + %Stitch with low wavelength spectrum + Morb_D38A_LowLam2=Morb_D38A_LowLam; + Morb_D38A_LowLam2(:,2)=Morb_D38A_LowLam(:,2)*(MORB_D38A2(6901,2)/0.15); + Morb_D38A_LowLam2=Morb_D38A_LowLam2(1:15,:); + MORB_D38A2=cat(1,Morb_D38A_LowLam2,MORB_D38A2(6902:end,:)); + % Interpolate data to 5 nm spacing between 1 and 4 microns + InterpMorb=interp1(MORB_D38A2(:,1),MORB_D38A2(:,2),WLS); + if i>1 % Use 650 C spectrum below 2.6 microns to isolate 3 micron feature changes + InterpMorb(1:321)=RMORB(1:321,1); + else + end + RMORB(:,i)=InterpMorb; + %Use Hapke model to convert from laboratory reflectance to single + %scattering albedo using reflectance and scattering asymmetry + % factor (p) of 0.81 (see manuscript for details on Hapke model) + SSAMORB(:,i)=Hapke_Inverse_Function_Passive(RMORB(:,i),0.81,WLS); end %% Interpolation of MORB Spectra -%Set range and resolution for interpolation +%Set range and resolution for interpolation %(standard is 1 ppm spacing from 0 to 1666 ppm) With 30% maximum MORB %abundance this gives a maximum total water abundance of 500 ppm. -WatInterp=linspace(0,1666,1667); +WatInterp=linspace(0,1666,1667); -for m=1:numel(WLS) +for m=1:numel(WLS) %For each wavelength, fit a line to the ESPAT values for each heating step spectrum ESPAT(m,:)=polyfit(WaterPPM(1:4),(1-SSAMORB(m,1:4))./SSAMORB(m,1:4),1); %linear function %For each wavelength, create new synthetic spectra from the linear ESPAT relationship - TotalWatInterpESPAT(:,m)=1./(ESPAT(m,1).*WatInterp+ESPAT(m,2)+1); + TotalWatInterpESPAT(:,m)=1./(ESPAT(m,1).*WatInterp+ESPAT(m,2)+1); end % Check match of interpolated spectra with real MORB spectra @@ -99,7 +99,7 @@ hold on end for k=1:4 -plot(WLS,TotalWatInterpESPAT(RealMorbIndices(k),:),'Color','black','linestyle','-','linewidth',2); + plot(WLS,TotalWatInterpESPAT(RealMorbIndices(k),:),'Color','black','linestyle','-','linewidth',2); end xlabel('Wavelength (\mum)','FontSize',16,'FontWeight','bold'); ylabel('SSA','FontSize',16,'FontWeight','bold'); @@ -203,7 +203,7 @@ PCTsMareImmature(j)=((HighAbunReg-LowAbunReg).*rand(1)+LowAbunReg); PCTsHighlandsImmature(j)=0; PCTsHighlandsMature(j)=0; - %Second half of spectra will be highlands spectra + %Second half of spectra will be highlands spectra else PCTsMareMature(j)=0; PCTsMareImmature(j)=0; @@ -211,35 +211,35 @@ PCTsHighlandsMature(j)=((HighAbunReg-LowAbunReg).*rand(1)+LowAbunReg); end PCTsPyroxene(j)=((HighAbunPyr-LowAbunPyr).*rand(1)+LowAbunPyr); -%Select a random interpolated MORB spectrum to simulate hydration -RandWatAmt(j)=round(rand(1)*max(WatInterp)); -%Round to the nearest 1 ppm -[M I]=min(abs(RandWatAmt(j)-WatInterp)); -%Actual water amount in mixture is morb abundance*interpolated MORB -%spectrum used -ActWatAmt(j)=WatInterp(I)*PCTsMORB(j); -WtPercent(j,:)=[PCTsMORB(j),PCTsHighlandsMature(j),PCTsHighlandsImmature(j),PCTsMareMature(j),PCTsMareImmature(j),PCTsPyroxene(j)]; %Wt% -if j2\u001b[0m \u001b[39mfor\u001b[39;00m spectrum \u001b[39min\u001b[39;00m heated_MORB_spectra:\n\u001b[1;32m 3\u001b[0m \u001b[39m# interpolate along our sampling grid of interest\u001b[39;00m\n\u001b[1;32m 4\u001b[0m interpolated_spectrum \u001b[39m=\u001b[39m simulator\u001b[39m.\u001b[39minterpolate_series(spectrum, WLS)\n\u001b[0;32m----> 5\u001b[0m ssa\u001b[39m.\u001b[39mappend(simulator\u001b[39m.\u001b[39;49mreflectance_to_ssa(interpolated_spectrum))\n", - "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/repro/simulator.py:174\u001b[0m, in \u001b[0;36mreflectance_to_ssa\u001b[0;34m(Refl, P, emission_angle, incident_angle, phase_angle, filling_factor)\u001b[0m\n\u001b[1;32m 169\u001b[0m y \u001b[39m=\u001b[39m partial(bidirectional_reflectance, P\u001b[39m=\u001b[39mP, mu\u001b[39m=\u001b[39mmu, mu0\u001b[39m=\u001b[39mmu0, B\u001b[39m=\u001b[39mB)\n\u001b[1;32m 170\u001b[0m OLS \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 171\u001b[0m ordinary_least_squares, y\u001b[39m=\u001b[39my, yx\u001b[39m=\u001b[39mvalue\n\u001b[1;32m 172\u001b[0m ) \u001b[39m# turn least squares into the form y=f(x)\u001b[39;00m\n\u001b[1;32m 173\u001b[0m w\u001b[39m.\u001b[39mappend(\n\u001b[0;32m--> 174\u001b[0m fmin(\n\u001b[1;32m 175\u001b[0m OLS,\n\u001b[1;32m 176\u001b[0m w0,\n\u001b[1;32m 177\u001b[0m args\u001b[39m=\u001b[39;49m(index),\n\u001b[1;32m 178\u001b[0m disp\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m,\n\u001b[1;32m 179\u001b[0m maxiter\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 180\u001b[0m maxfun\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 181\u001b[0m ftol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 182\u001b[0m xtol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 183\u001b[0m )\n\u001b[1;32m 184\u001b[0m )\n\u001b[1;32m 185\u001b[0m \u001b[39mreturn\u001b[39;00m w\n", + "\u001b[1;32m/home/phil/repos/philiplinden/desci/cremons-etal-2022/repro/results.ipynb Cell 6\u001b[0m line \u001b[0;36m1\n\u001b[1;32m 4\u001b[0m interpolated_spectrum \u001b[39m=\u001b[39m simulator\u001b[39m.\u001b[39minterpolate_series(spectrum, WLS)\n\u001b[1;32m 6\u001b[0m \u001b[39m# Use Hapke model to convert from laboratory reflectance to single\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[39m# scattering albedo using reflectance and scattering asymmetry\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[39m# factor (p) of 0.81 (see manuscript for details on Hapke model)\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m ssa\u001b[39m.\u001b[39mappend(simulator\u001b[39m.\u001b[39;49mreflectance_to_ssa(\n\u001b[1;32m 11\u001b[0m reflectance\u001b[39m=\u001b[39;49minterpolated_spectrum, asymmetry_factor\u001b[39m=\u001b[39;49m\u001b[39m0.81\u001b[39;49m))\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/repro/simulator.py:174\u001b[0m, in \u001b[0;36mreflectance_to_ssa\u001b[0;34m(reflectance, asymmetry_factor, emission_angle, incident_angle, phase_angle, filling_factor)\u001b[0m\n\u001b[1;32m 167\u001b[0m y \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 168\u001b[0m bidirectional_reflectance, P\u001b[39m=\u001b[39masymmetry_factor, mu\u001b[39m=\u001b[39mmu, mu0\u001b[39m=\u001b[39mmu0, B\u001b[39m=\u001b[39mB\n\u001b[1;32m 169\u001b[0m )\n\u001b[1;32m 170\u001b[0m OLS \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 171\u001b[0m ordinary_least_squares, y\u001b[39m=\u001b[39my, yx\u001b[39m=\u001b[39mvalue\n\u001b[1;32m 172\u001b[0m ) \u001b[39m# turn least squares into the form y=f(x)\u001b[39;00m\n\u001b[1;32m 173\u001b[0m w\u001b[39m.\u001b[39mappend(\n\u001b[0;32m--> 174\u001b[0m fmin(\n\u001b[1;32m 175\u001b[0m OLS,\n\u001b[1;32m 176\u001b[0m w0,\n\u001b[1;32m 177\u001b[0m args\u001b[39m=\u001b[39;49mindex,\n\u001b[1;32m 178\u001b[0m disp\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m,\n\u001b[1;32m 179\u001b[0m maxiter\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 180\u001b[0m maxfun\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 181\u001b[0m ftol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 182\u001b[0m xtol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 183\u001b[0m )\n\u001b[1;32m 184\u001b[0m )\n\u001b[1;32m 185\u001b[0m \u001b[39mreturn\u001b[39;00m w\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:747\u001b[0m, in \u001b[0;36mfmin\u001b[0;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, initial_simplex)\u001b[0m\n\u001b[1;32m 738\u001b[0m opts \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mxatol\u001b[39m\u001b[39m'\u001b[39m: xtol,\n\u001b[1;32m 739\u001b[0m \u001b[39m'\u001b[39m\u001b[39mfatol\u001b[39m\u001b[39m'\u001b[39m: ftol,\n\u001b[1;32m 740\u001b[0m \u001b[39m'\u001b[39m\u001b[39mmaxiter\u001b[39m\u001b[39m'\u001b[39m: maxiter,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[39m'\u001b[39m\u001b[39mreturn_all\u001b[39m\u001b[39m'\u001b[39m: retall,\n\u001b[1;32m 744\u001b[0m \u001b[39m'\u001b[39m\u001b[39minitial_simplex\u001b[39m\u001b[39m'\u001b[39m: initial_simplex}\n\u001b[1;32m 746\u001b[0m callback \u001b[39m=\u001b[39m _wrap_callback(callback)\n\u001b[0;32m--> 747\u001b[0m res \u001b[39m=\u001b[39m _minimize_neldermead(func, x0, args, callback\u001b[39m=\u001b[39;49mcallback, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mopts)\n\u001b[1;32m 748\u001b[0m \u001b[39mif\u001b[39;00m full_output:\n\u001b[1;32m 749\u001b[0m retlist \u001b[39m=\u001b[39m res[\u001b[39m'\u001b[39m\u001b[39mx\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mfun\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnit\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnfev\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mstatus\u001b[39m\u001b[39m'\u001b[39m]\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:899\u001b[0m, in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)\u001b[0m\n\u001b[1;32m 897\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m 898\u001b[0m \u001b[39mfor\u001b[39;00m k \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(N \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m):\n\u001b[0;32m--> 899\u001b[0m fsim[k] \u001b[39m=\u001b[39m func(sim[k])\n\u001b[1;32m 900\u001b[0m \u001b[39mexcept\u001b[39;00m _MaxFuncCallError:\n\u001b[1;32m 901\u001b[0m \u001b[39mpass\u001b[39;00m\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:620\u001b[0m, in \u001b[0;36m_wrap_scalar_function_maxfun_validation..function_wrapper\u001b[0;34m(x, *wrapper_args)\u001b[0m\n\u001b[1;32m 618\u001b[0m ncalls[\u001b[39m0\u001b[39m] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m1\u001b[39m\n\u001b[1;32m 619\u001b[0m \u001b[39m# A copy of x is sent to the user function (gh13740)\u001b[39;00m\n\u001b[0;32m--> 620\u001b[0m fx \u001b[39m=\u001b[39m function(np\u001b[39m.\u001b[39mcopy(x), \u001b[39m*\u001b[39m(wrapper_args \u001b[39m+\u001b[39;49m args))\n\u001b[1;32m 621\u001b[0m \u001b[39m# Ideally, we'd like to a have a true scalar returned from f(x). For\u001b[39;00m\n\u001b[1;32m 622\u001b[0m \u001b[39m# backwards-compatibility, also allow np.array([1.3]),\u001b[39;00m\n\u001b[1;32m 623\u001b[0m \u001b[39m# np.array([[1.3]]) etc.\u001b[39;00m\n\u001b[1;32m 624\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m np\u001b[39m.\u001b[39misscalar(fx):\n", @@ -122,7 +105,12 @@ "for spectrum in heated_MORB_spectra:\n", " # interpolate along our sampling grid of interest\n", " interpolated_spectrum = simulator.interpolate_series(spectrum, WLS)\n", - " ssa.append(simulator.reflectance_to_ssa(interpolated_spectrum))" + "\n", + " # Use Hapke model to convert from laboratory reflectance to single\n", + " # scattering albedo using reflectance and scattering asymmetry\n", + " # factor (p) of 0.81 (see manuscript for details on Hapke model)\n", + " ssa.append(simulator.reflectance_to_ssa(\n", + " reflectance=interpolated_spectrum, asymmetry_factor=0.81))" ] } ], diff --git a/repro/simulator.py b/repro/simulator.py index d36e94e..aabb709 100644 --- a/repro/simulator.py +++ b/repro/simulator.py @@ -8,7 +8,7 @@ """ from collections import namedtuple from functools import partial -from typing import Any, Callable +from typing import Callable import numpy as np from numpy.typing import ArrayLike @@ -27,12 +27,12 @@ def interpolate_series(data: Series, new_index: ArrayLike) -> Series: return Series(data=new_y, index=new_index) -def ordinary_least_squares(x: Any, y: Callable, yx: Any): +def ordinary_least_squares(x: float, y: Callable, yx: float) -> float: """Ordinary Least Squares Function. y (Callable): The estimator function. - x (Any): The argument to y. - yx (Any): The observation. Must be the same type as the result of y. + x (float): The argument to y. + yx (float): The observation. Must be the same type as the result of y. """ return sum((y(x) - yx) ** 2) @@ -61,7 +61,9 @@ def backscattering(h: float, g: float) -> float: def bidirectional_reflectance( SSA: float, P: float, mu: float, mu0: float, B: float ) -> float: - """Bidirectional reflectance, see Equation 1. + """Estimate bidirectional reflectance from single-scattering albedo. + + see Equation 1 in the manuscript: R = (ω/4) (μ₀ / (μ + μ₀)) {(1 + B)P + H(ω)H₀(ω) - 1} @@ -85,7 +87,7 @@ def bidirectional_reflectance( B (float): backscattering function Returns: - float: _description_ + float: bidrectional reflectance """ def h_func(SSA: float, mu: float, mu0: float) -> HFunction: @@ -123,14 +125,14 @@ def h_func(SSA: float, mu: float, mu0: float) -> HFunction: def reflectance_to_ssa( - Refl: Series, - P: float = 0.15, + reflectance: Series, + asymmetry_factor: float = 0.81, emission_angle: float = 0, incident_angle: float = 30, phase_angle: float = 30, filling_factor: float = 0.41, ): - """Convert reflectance spectrum to single-scattering albedo (SSA) + """Estimate single-scattering albedo (SSA) from reflectance Uses scipy.optimize.fmin (equivalent to MATLAB fminsearch) to minimize ordinary least squares distance between SSA obtained from the supplied @@ -144,10 +146,8 @@ def reflectance_to_ssa( Default parameter values replicate src/Hapke_Inverse_Function_Passive.m Args: - R (Series): Bidirectional reflectance, see Equation 1. - p (float, optional): Scattering phase function. Defaults to 0.15 for - ansiotropic scattering on the modeled mean particle phase function - for lunar soil (Goguen et al., 2010). + reflectance (Series): Bidirectional reflectance, R. see Equation 1. + asymmetry_factor (float, optional): Scattering asymmetry factor, p. emission_angle (float, optional): Emission angle in degrees. Defaults to 0. incident_angle (float, optional): Incident angle in degrees. @@ -163,13 +163,18 @@ def reflectance_to_ssa( B = backscattering(h, g) w = [] - for index, value in Refl.items(): + for index, value in reflectance.items(): w0 = 0.5 # initial guess, ω₀ # turn bidrectional_reflectance() into the form y=f(x) - y = partial(bidirectional_reflectance, P=P, mu=mu, mu0=mu0, B=B) + y = partial( + bidirectional_reflectance, P=asymmetry_factor, mu=mu, mu0=mu0, B=B + ) + + # formulate ordinary least squares between estimated reflectance, y + # and observed reflectance, yx, and arrange into the form y=f(x) OLS = partial( ordinary_least_squares, y=y, yx=value - ) # turn least squares into the form y=f(x) + ) w.append( fmin( OLS, From 6473afbb02a2cec8c6ed70a973904948526efce8 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Sat, 4 Nov 2023 11:48:44 -0400 Subject: [PATCH 5/5] some cleanup --- repro/results.ipynb | 47 +++++++++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/repro/results.ipynb b/repro/results.ipynb index b23baa4..93e2897 100644 --- a/repro/results.ipynb +++ b/repro/results.ipynb @@ -7,7 +7,7 @@ "outputs": [], "source": [ "%load_ext autoreload\n", - "%autoreload 2" + "%autoreload 2\n" ] }, { @@ -22,7 +22,7 @@ "\n", "# local modules\n", "import simulator\n", - "import spectra" + "import spectra\n" ] }, { @@ -39,7 +39,7 @@ "num_spectra = 100\n", "\n", "# Total water measured from step-heating experiments\n", - "water_ppm = [1522, 762, 176, 22]" + "water_ppm = [1522, 762, 176, 22]\n" ] }, { @@ -51,14 +51,25 @@ "# Low wavelength portion of MORB spectrum (<1.5 microns)\n", "MORB_D38A_LowLam = spectra.spectrum_from_file(\n", " Path(\"../data/Morb_D38A_Low_wavelength.txt\")\n", - ")" + ")\n" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAACEr0lEQVR4nOzdd3hUVf7H8fdMep30Bmm0kJDQizRpkVBUEKUp3bK66k+XVVdc6+ourrq76OqKHSxUFVBApFfpPSEEAuk9Icmkl5n7++PChEgogdTh+3qePDp3zr1zDsMwn5x7ikZRFAUhhBBCCDOlbe4KCCGEEEI0Jgk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZs2yuSvQEIxGI+np6Tg5OaHRaJq7OkIIIYS4AYqiUFRUhJ+fH1pt4/W/mEXYSU9Px9/fv7mrIYQQQoibkJKSQtu2bRvt+mYRdpycnAD1D8vZ2bmZayOEEEKIG6HX6/H39zd9jzcWswg7l25dOTs7S9gRQgghWpnGHoIiA5SFEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpg1CTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIs1avsDN//nz69OmDk5MTXl5ejB8/nri4uOuet3LlSjp37oytrS0RERGsX7++1vOKovDqq6/i6+uLnZ0dkZGRnD17tn4tEUIIIYSoQ73Czo4dO3jyySfZt28fmzZtoqqqipEjR1JSUnLVc3777TemTp3Kww8/zNGjRxk/fjzjx48nOjraVOadd97hgw8+YOHChezfvx8HBweioqIoLy+/+ZYJIYQQQgAaRVGUmz05JycHLy8vduzYwZ133llnmcmTJ1NSUsLatWtNx+644w66d+/OwoULURQFPz8//vznP/Pcc88BUFhYiLe3N4sWLWLKlCnXrYder0en01FYWIizs/PNNkcIIYQQTaipvr9vacxOYWEhAG5ublcts3fvXiIjI2sdi4qKYu/evQAkJCSQmZlZq4xOp6Nfv36mMr9XUVGBXq+v9SOEEEIIUZebDjtGo5Fnn32WgQMHEh4eftVymZmZeHt71zrm7e1NZmam6flLx65W5vfmz5+PTqcz/fj7+99sM4QQQghh5m467Dz55JNER0ezbNmyhqzPDZk3bx6FhYWmn5SUlCavgxBCCCFaB8ubOempp55i7dq17Ny5k7Zt216zrI+PD1lZWbWOZWVl4ePjY3r+0jFfX99aZbp3717nNW1sbLCxsbmZqgshhBDiNlOvnh1FUXjqqadYtWoVW7duJTg4+Lrn9O/fny1bttQ6tmnTJvr37w9AcHAwPj4+tcro9Xr2799vKiOEEEIIcbPq1bPz5JNPsmTJEtasWYOTk5NpTI1Op8POzg6AGTNm0KZNG+bPnw/AM888w5AhQ/jXv/7F2LFjWbZsGYcOHeLTTz8FQKPR8Oyzz/LWW2/RsWNHgoODeeWVV/Dz82P8+PEN2FQhhBBC3I7qFXY+/vhjAIYOHVrr+FdffcWsWbMASE5ORqut6TAaMGAAS5Ys4eWXX+all16iY8eOrF69utag5hdeeIGSkhIee+wxCgoKGDRoEBs2bMDW1vYmmyWEEEIIobqldXZaCllnRwghhGh9WsU6O0IIIYQQLZ2EHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpg1CTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCNAtFUZrkdSyb5FWEEEIIIS4qrSrll7NrObh0YZO8noQdIYQQQjSJhMIEVh3+hqplqxh0pJz7Cw38swleV8KOEEIIIRpNtbGaHSk7WLtvMb7rDjHimIJ9pfpcvs6+SeogYUcIIYQQDa64sphV8avYuHMRfbdnMOeEgqVRfa66XVv8n/oTfn37gqdno9dFwo4QQgghGkx6cTrfxX7HiS0ruGtPCS+drRmEbNGzK36PP4nDoEFotFr0en2T1EnCjhBCCCFu2fGc43xz6huS9mzkgZ3V3J1YE3Lshw3B8+FHsO/du1nqJmFHCCGEEDel2ljN1uStfB2zGOXAMcbvVXg4WQ05iqUFLvfdh/vs2di0a9es9ZSwI4QQQoh6Ka4s5sezP7Lk9BIcTqcyaZeRrpd6ciws0N03Hs8nnsCqTZvmrehFEnaEEEIIcUMyijP4NvZbfjjzPcHxxTy8RyH8Yk8OVla4PTgVt9mzsfLxuf7FClPh0MrGrfBFEnaEEEIIcU2xebEsilnEr4m/EppQxYs7jYSkXXzS0hLd+HF4PP441m3bXvtCRZlwciVE/wjpR6BCVlAWQgghRDNRFIXf0n9jUcwiDqbupW+cwuuHakKOxsYGl4kTcX94Dla+vle/UGUJnF4Hx5fB+W2gXJx/jgba9gG2NHZTJOwIIYQQokaVoYpfEn9hUcwizuWdYegJhQ93G3EvUp/XWFnhMnkyHn94DMurrZFjNELiLjXgxP4ElcU1z/n3g66TIPReMNrCU7pGb5OEHSGEEEJQWFHIyjMrWRq7lJILWQw7ofB/RxW889VbTRaeHrhOnITLlMlYeXnVfZHs03BiGZxYAfq0muOuQdB1ihpy3NvXHG+p6+zs3LmTd999l8OHD5ORkcGqVasYP378VcvPmjWLxYsXX3E8LCyMmJgYAF5//XXeeOONWs+HhIRw+vTp+lZPCCGEEPWQUpTCt6e+ZVX8KoylpYw8onD/XgX78oshx9UVj8f/gMuUKWhtbK68QHEORH8Px5dCxvGa47Y66DIBuk1Re3M0miZq0ZXqHXZKSkro1q0bc+bMYcKECdct//777/P222+bHldXV9OtWzcmTpxYq1yXLl3YvHlzTcUspdNJCCGEaCzHso/x9amv2ZK8BYtKA3cfVBh3UIN9qTqmxqZjB1ynT0c3dixaB4faJ1eVQdx6OL4c4jeDYlCPay2hYxR0m6z+18q2iVtVt3onitGjRzN69OgbLq/T6dDpau7HrV69mvz8fGbPnl27IpaW+NzIVDWgoqKCiooK0+OmWm5aCCGEaM0MRgNbU7ayOGYxx3OOY1WtMPiUwvTfrHDOV79Xrfz98Xjyj+juuQeNhUXtC2TFwKGv4MRyqLjsu7dNb7UHp8sEcHBvwhbdmCbvPvniiy+IjIwkMDCw1vGzZ8/i5+eHra0t/fv3Z/78+QQEBNR5jfnz519x20sIIYQQdSutKmV1/Gq+OfUNqcWpOJcoTD6qYewxS2yLKoAKLH198Zr7J5zHjKkdcqrK4dQaOPQlpOyrOa4LUHtwuk4Gj45N3qb60CiKctOT3DUazXXH7FwuPT2dgIAAlixZwqRJk0zHf/nlF4qLiwkJCSEjI4M33niDtLQ0oqOjcXJyuuI6dfXs+Pv7U1hYiLOz8802RwghhDArmSWZLDu9jJVnVqKv1OOmV3jgkBXDjlRhUaXeerL09cXtoQdxnTYNre1lt53yzsHhr+Dod1B2QT2mtYTOY6HXbAgeAlrtLdVPr9ej0+ka/fu7SXt2Fi9ejIuLyxXh6PLbYl27dqVfv34EBgayYsUKHn744SuuY2Njg01dg6SEEEIIwYmcE3x76ls2Jm3EoBjwzVN47LAdvY+Xoq0uB8A2IgL3ObNxuusuNJfGyRqq1LE4h76E89trLujcFnrNgp7TwenGhpy0JE0WdhRF4csvv2T69OlYW1tfs6yLiwudOnUiPj6+iWonhBBCtG5Vxio2J23m29hvOZFzAhSF8CSFKTEudDyZj0ZR17qx79sXj8f/gH3//mguzZAqTIXDi+DI11CcdfGKGug4EnrPgY53gdaiztdtDZos7OzYsYP4+Pg6e2p+r7i4mHPnzjF9+vQmqJkQQgjReuWX5/P9me9ZFreM7NJsLAwKQ+K0TD3igFtKIaDegnIcPhyPxx7Frnt39URFgfM74OBncHp9zYwqBy+1B6fnTHANrPM1W5t6h53i4uJaPS4JCQkcO3YMNzc3AgICmDdvHmlpaXz99de1zvviiy/o168f4eHhV1zzueee45577iEwMJD09HRee+01LCwsmDp16k00SQghhDB/Z/LP8F3sd6w7v44KQwUOZQoPnbBl1FGwyS8BCtHY2+MyfhyuU6di0/HiIOJyvTqb6sBnkBtXc8HAQdDnYeh8N1he+w5Ma1PvsHPo0CGGDRtmejx37lwAZs6cyaJFi8jIyCA5ObnWOYWFhfzwww+8//77dV4zNTWVqVOnkpeXh6enJ4MGDWLfvn14Xm0ZaiGEEOI2ZDAa2JG6g+9iv+NA5gEAXIsUnjjpyh0H9GjLSgCw9PTEZeoUXKdOxdLVVT05+7Tai3N8Wc32DVYO6pTxPo+Ad1hzNKlJ3NJsrJaiqUZzCyGEEM2hqLKIVWdXsfT0UlKLUwFoe0HDH6K96bQ/A021egvKJiQE90cexjkqCo21NRiqIW6d2ouTuKvmgh6doM+jatCxbb7vTbOcjSWEEEKIGxefH8/yuOX8dO4nSqtLAYjItefRY+74HEoERQ0+9r174/7YozgMHqwOOi7Ohr2L1VlVRenqxTRaCBkDfR9Vp4034/YNTU3CjhBCCNGCVBmq2Jy8meVxyzmcdVg9qCiMzPFlyiFrHI+fB9TVix2HD8f9kUew79lDHXCcsl/txTm1BoxV6rn2HtBrpro2jot/8zSqmUnYEUIIIVqA9OJ0Vp5ZyY9nf+RCuTqDylLRMisvlOE7CrA8k6QWtLBAd/fduD/ysDrouLIUDi9Wx+Nknqy5YNu+ai9O2DiwvL3XppOwI4QQQjQTg9HAnvQ9rIhbwc7UnSiow2iDje48ktGJkB2JGJPUncQ1tra4TJyI+6yZWLVpAwUpsOlVdX2c8kL1gpa2EPGAOh7Hr3vzNKoFkrAjhBBCNLG8sjxWxa/i+zPfk1acZjo+xLknM4444fzzHpTyXRgBrU6H20MP4TrtISzd3CDlIKz8K5z6qWZtHNcg6P0w9JgG9m7N0qaWTMKOEEII0QQUReFo9lGWxy1nU9Imqi6OqXG2cmK24Q7uPFyBcdselIoKFMAmLBTXiRNxvudeLOys1XE4P3wMaYdqLhp8J9zxR+gYdcv7VJkzCTtCCCFEIyqpKmHtubUsP7Ocs/lnTcd7OIYyJzmItr+eoCrxFy720WAbFobH/z2N45AhaMoL4PCn6qBj/cUeIAtriJgEdzwOPhFN3p7WSMKOEEII0QjiLsSxIm4Fa8+vNU0bt7WwZYLbMMYds8JizRaMhSepArT29jiPHYvLxAewjYhAk3sG1s1VFwCsUs/FwVNd/K/3HHD0ar6GtUISdoQQQogGUmmoZGPSRlbEreBo9lHT8SDnIGbaDqP39gzK129AqarCCFgFBuA2cyYu48ahtbODc1vhuwcgfnPNRb3D1VtV4feDlW3TN8oMSNgRQgghblFKUQorz6xk9dnV5FfkA2CpsWS4/zAeLInAc81vlOz6jLKL5e169sR9zmwchw1DYyiH40th/yeQe+ZiCY26AOAdT0DQoNtqAcDGIGFHCCGEuAkGo4FdabtYHrecPWl7TNPGve29mdz2Hkadd6bqvz9TceoXSgC0Wpzuugv32bPUnccLU2HrG7Wnjls7qTuO930M3IKbqWXmR8KOEEIIUQ+5ZbmsOruKlWdWklGSYTo+wG8AD9neSYdt8ejnf01xqTrWRmNnh8uECbjNnIF1QACkHobv50DM6sumjgdDv8eh+4PNuleVuZKwI4QQQlyHoigcyDzAirgVbE3eSrVSDYDORseEoHsZn+qDxeJNlB56i4t9NFi3a4fLhPvQ3X8/ljpniFsPXzwKKftqLhw0WB2P0ykKtBZN37DbhIQdIYQQ4ioKygtYc24N35/5nkR9oul4V8+uPOQ6ku57sin+7xoq8vLUJ7RaHIcPw+3BB7Hv3x9NVSkcWwJ7P4L8hItlrCBiojoex7dr0zfqNiRhRwghhLiMoigcyznGirgVbEzcSKWxEgAHKwfuDhrDAxfa4/DTbop3/JNCRR2nY+nlhcvEibhMfAArHx8oyoStb8LBL6C8QL2wrQv0eVgdj+Pk0zyNu01J2BFCCCGAosoi1p5fy4q4FcQXxJuOh7qFMtVnLHccLqHk1VVUpS6l+OJzDgP64zJlCk7DhqGxsoLcs/DT39X1cQxqSMI1GPo/qY7HsXZo+oYJCTtCCCFubzF5MayMW8n6hPWUVauTw20tbBkdNIpJld1w/WU/RRv+RUGVur2D1tkZl/vuw2XKZGyCL86YSjkIexbA6XVwcVYW/v1gwNPqFHIZj9OsJOwIIYS47ZRWlbIhcQMr4lYQkxdjOt5e154pAeMYckpD2burqYj7Hv3F52wjInCdMgXnMaPVBQAVBc5sVENO0p6ai4eMgYHPQMAdTdomcXUSdoQQQtw2YvJi+PHMj6xPWE9xlXozykprxV2BdzHZZiC+G46hf/MjCkpKANDY2uI8dgyuU6ZiFxGuXsRQBceXw573IftiUNJaQdfJMPD/wDOkOZomrkHCjhBCCLNWVFnE+vPr+eHsD8ReiDUd93fyZ2L7+xmV50fVd6sp2fEXCi4+Zx0UhOvUKejGj8dCp1MPVpbAkW9g74dQmHKxoCP0nq1OH3f2a9J2iRsnYUcIIYTZURSFE7knTDOqyg3lgNqLE+k/ggcs+xJ4IAX9f7+jIOPiwoAaDU6RI3B98EHs77gDzaUtGkovqLuO718IZRfUYw6e6tTx3g+DnUvTN1DUi4QdIYQQZqO0qtQ0oyouP850vJ2uHQ8EjWd4nBUV/1pOxdm1XIwtaJ2d0d17L27Tp2EdGFhzMX2G2otzeBFUXpx/5RqsDjru/iBY2TVZu8StkbAjhBCi1UsoTGB53HLWxK8xjcWxsbAhKiiKibph+P5yhML5n6AvVNc31lhb43DnYHRjxuA4fDha28t2E887p47HOb60Zvq4dwQM/hOEjZeZVa2QhB0hhBCtUpWxiq3JW1kZt5L9mftNxwOcApjcYSIj092pXPwzJbv/jwuXLf7nNnMGLpMmYeHkVPuCWTGw698Q8yMoxosXGwCD50KHSNl5vBWTsCOEEKJVSSxMZHX8atacW0NuWS4AWo2WO9vcyVSPUXTYlUDhB4vJz8w0nWN/xx24zZyB4513orH4Xc9M2hHY+a66d9UlHUfC4D/L9HEzIWFHCCFEi6YoClmlWexO281P537iaPZR03Medh5MaH8f4wraoV2ziaItfyHPoO4kbuHqisv9E3CZNEndbfz30o/C9rfhzIaLBzQQNk4NObJnlVmRsCOEEKLFSS1K5fOTn3Ms+xhpxWmm2VSg9uIMajOIB6zvoEtcOUUv/UBJUpLpebtevXCdMgWnqJFora2vvHjGcTXkXOrJ0WghYpIacjw7NXbTRDOQsCOEEKLFyC/P57OTn7Hs9DKqjFWm4xo0hLmHEeU+iBGnraj+6FcqTv2Di3uNo3VwQDfuXlwmT8E25CqBJTMats+H02svXlQL4Q/AkBfAo2PjNkw0Kwk7Qgghml1ZdRnfxX7HFye/MM2m6u/bn4dCHyIYD5zj0ilZv4GizZ9RUqnOkNJYWWHXuxfOo0aju3ssWoerbLKZFaP25MT+dPGABsLvhyF/kZ6c24SEHSGEEM2msKKQpaeXsiR2CfkV+QB0duvM3PaP0nFnAgX/+gclScmUXHaOTUgILvffj/M9d2Pp6nr1i2fHqiHn1OqLBzTQ5T415Hh1bqwmiWtQFIWMwnJOZ+o5l11CbHLm9U9qAPUOOzt37uTdd9/l8OHDZGRksGrVKsaPH3/V8tu3b2fYsGFXHM/IyMDHx8f0+KOPPuLdd98lMzOTbt268d///pe+ffvWt3pCCCFagbTiNJbELmHlmZWmncbbOrblWa9JdN2aTOErfyGnvGacjlVAAI5DhqAbNw7bLmE1qxvXJSdODTkxqzDtQB42Hoa+CF6hjdcoUUtecQWxGUUk5JWQml/KqXQ9Mel6LpRUmsoYK0qbpC71DjslJSV069aNOXPmMGHChBs+Ly4uDmdnZ9NjLy8v0/8vX76cuXPnsnDhQvr168eCBQuIiooiLi6uVjkhhBCtl6IoHMw8yJLTS9iWsg3jxbVsQlw68UeG0mHDKUq2v0vBxTVxbEJDcZ81E8ehQ2v2p7qWnDOw458Q/QOmkBN6rxpyvLs0UqtEtcFIQm4JpzL0xGYUEZuhJzZDT3ZRRZ3lLbQaOng60sHLET97hZcXNH4dNYpy8W/VzZys0dxwz05+fj4uLi51lunXrx99+vThww8/BMBoNOLv78/TTz/Niy++eN166PV6dDodhYWFtQKVEEKI5ldWXcb68+v57vR3nM0/azo+0Ksfc7I74/nTPspjYkzHHYcOxW3WLOz79b12D84lufEXQ873NYsBdr5bDTk+EQ3dnNteaWU1h5PyOZBwgf3nL3A8tYCKauMV5TQaCHJ3oL2nI746W0J9neni50yIjxO2VupaR031/d1kY3a6d+9ORUUF4eHhvP766wwcOBCAyspKDh8+zLx580xltVotkZGR7N27t85rVVRUUFFRkxj1en3jVl4IIUS9nSs4x/dnvufn8z9TWKFu02Bnacc0mzuJirVB+/kuqnP2UA5obGzQ3TcetxkzsWkXfGMvkHcOdrwDJ1fUhJyQsWrIkXVyGoy+vIpDiRfYfzHcRKcVUm2s3U9ib21BZx8nwvycCfVVf0K8nXCwaRlDgxu9Fr6+vixcuJDevXtTUVHB559/ztChQ9m/fz89e/YkNzcXg8GAt7d3rfO8vb05ffp0ndecP38+b7zxRmNXXQghRD2lFaexK3UX6xPW11r8L9DGjyfTOtNxVxJVp9dhBIyAhYcHbg89iMuUKdcebHy53HjY9S84sRwUdQFBOo1WQ45f94Zu0m2nstrIgYQLbIvLZn9CHqfS9fwu29DGxY5+wW70a+dG7yA3gt0d0Gpb7nYajR52QkJCCAkJMT0eMGAA586d4z//+Q/ffPPNTV1z3rx5zJ071/RYr9fj7+9/y3UVQghRf0bFyM/nfmbp6aXE5NXcjrJUtEwp68pdCQ447j6J4cJGqgCsrHAaOgTd+PE4Dh6Mpq6F/+qSHatu6xCzqqYnp2OUGnLa9Gzwdt0uyqsMHE7K52hyPkeSCziYcIGiiupaZYI9HOgbpIabvsFutHW1b6ba3pxm6V/q27cvu3fvBsDDwwMLCwuysrJqlcnKyqo1W+tyNjY22NjYNHo9hRBCXFvchTje2vcWx3KOAWClaBmrb8ewBHvaHE7BmH0IAANg6eeL+6zZ158yfrmiTHURwFM/QcKOmuOdRsOdz0PbXg3boNtEekEZ+xPy2HfuAuujMygqrx1uPBxtGNHZi0EdPegb7Ia3s+1VrtQ6NEvYOXbsGL6+vgBYW1vTq1cvtmzZYhrobDQa2bJlC0899VRzVE8IIcR1lFSV8NGxj1gSuwSDsZpeqTbMSgvG52gqSoE6BMEIaHU6nCJH4BwVhcOAAWgsb/BrR58Bv86DmNWYZlaBOrvqzudlTE49GYwKBxIusDk2i62ns0nILan1vLezDX2C3OgZ4ErvIFfC/XQt+rZUfdU77BQXFxMfH296nJCQwLFjx3BzcyMgIIB58+aRlpbG119/DcCCBQsIDg6mS5culJeX8/nnn7N161Y2btxousbcuXOZOXMmvXv3pm/fvixYsICSkhJmz57dAE0UQgjRkHam7uRve/9GVkkmEYkKjx5wxud8ARCDAli4uOA4fDhOkZE4DBpY9/5UV2M0wuEvYfMbUHFx8knbPhB6jxp03G5w8LKgstrIb+dy2RCdyaZTWeRdtr6NhVZDeBsdfQJdGdbZi/7t3M0q3PxevcPOoUOHai0SeGnszMyZM1m0aBEZGRkkJyebnq+srOTPf/4zaWlp2Nvb07VrVzZv3lzrGpMnTyYnJ4dXX32VzMxMunfvzoYNG64YtCyEEKL5ZJZk8t+j/+Wncz8RlmTk2d+saJ9YARSYZlM5jxqNfe9eN96Dc7mM47DuOUg9oD726wn3vC+9OPVQXmVg55kcNkRnsjk2C/1lt6d0dlZEhnoTGarennKytWrGmjatW1pnp6WQdXaEEKJxKIpCXH4cS08v5adzP9EhqYpJu4yEJ6lfHRpra1wmT8b90UewutlFYEsvwNY34dBXgALWTjDiVejzMGgtGq4xZqqkopqtp7PZEJPJttPZlFYaTM95ONoQ1cWb0eG+9GvnhpWFthlreiWzW2dHCCFE61FpqGTt+bV8F/sdSVlx9IpXmHdMIeJSyLGywmXiRNz/8BhWN9sLb6iGw1/B1regvEA9Fv4AjHwTnP0apiFmqrCsis2nsvglOpOdZ3OovGxRPz+dLVHhPowO96VXoCsWZnx76kZJ2BFCCFHLkawjvLHtJTr8lsL0GCPtM8Dy0neplRUuEybg8YfHsPK7hUCSuBt+eRGyTqqPvcNh9DsQNPCW62+u8ksq2XQqi/XRGeyJz6XKUHNjJsjdnlHhvowO96FrW92NrTx9G5GwI4QQAoBqYzX/2/cf8r9axEuHjDiX1TxnFRiA85gxuNz/ANZt29z8i+Scgc2vQdx69bGtCwx/GXrNBgv5Svo9RVE4klzAV3sS2BCdWWvl4k7ejqaA09nHSQLONcjfLCGEEGSWZPK/z/7A8GVn8C5Qj1m08cNj5iwchw/Dqk2bW/syLc6B7fPh8CJ11WONBfSaCcNeBgf3hmiC2TAaFY4k5/NrTCYbYjJJuVCTOsN8nRkT4cOocF86eDk2Yy1bFwk7Qghxm9sVvY4zb/6VqcfVPQerPVwIePFlnEdF3dysqstVlsK+j2D3+1BZpB4LGQORb4Bnp1usuXlJLyhj2cEUvj+UQnphuem4rZWWu7v6MWdgMGF+MgnnZkjYEUKI21RJYS4b5v+R4PUnGVAJRg3YTBxHpxdewcLR4dYubjTA8aWw9e9QlK4e8+sBI9+CoEG3XnkzUW0wsj0uh2UHk9l6Otu0B5WTjSUjQr2I6uLDkBBP7K3l6/pWyJ+eEELcZgwlJRz+dD58u4qwEnXkcUGQO13+sQDnnr1v7eLVleou5LsXQN5Z9ZguACJfgy4TQNuypj43l/M5xaw8nMqPR1LJ0leYjt/Rzo2pfQOI6uKDrZVMu28oEnaEEOI2UZWVRcKX/6Nk5Y84laqLzWW7W2D/5KPcMfX/bm1MTnUlHPoS9rxf05Nj5wqDn4M+j4BV695bqSGUVFSz/mQGKw6lcDAx33Tc1d6KB3q1ZXKfABmH00gk7AghhBlTjEZK9+0je+l3lG7ZhtaoYAtkumrIeWAwY554B0d73a29SMJO+OlpyE9UHzv5Qv8nodcssHG6xRa0boqicDSlgOUHUlh7Ip2Siwv+aTUwNMSLib3aMjzUCxtL6cVpTBJ2hBDCDCnV1ejXrSNr4UIMCYkAaIFT/pA0uiv3z36bYa63uM+UoRp2vA073wMUcPSGoS9C94fA0uZWm9CqlVcZ+Ol4Ol/vTSQ6TW86HuRuz8Te/jzQq22r30m8NZGwI4QQZqQyKQn9+vVc+P57DGnq7aQSG9gTpuHM0HY8dPdL3N9mwK2/kD4DfngYkvaoj3vOgKj5YHN734ZJKyjj231JLDuQTH5pFQDWllru7urL5N7+9A12k/VwmoGEHSGEaOUURaF0/wEufPsNxZu3mI4X2sPavlpSIrsws+/jPO0/FK2mAQYIx2+BHx+D0lx1H6t7FkDEA7d+3VZKURT2nb/A4t8S2Xgq0zSjqo2LHdPuCGRyH3/cHOqx87tocBJ2hBCilarOz0e/fj2FP/1E+fETgDp9/ESQhj1hGtL7BvHcoL/S369/w/QmGKrVhQF3/QtQwCcCJi4G9/a3fu1WqLSymtVH1VtVpzOLTMf7t3Nn5oAgIkO9sGxhG2/eriTsCCFEK6IoCqUHD1KwbBlFmzajVKm3SiqtNOzoAuv7aKkK8OaxiMeY0HECVhZWDfPCOXHqIOSU/erjXrNh1Nu35SyrmPRCvtiVwIaYTNMO43ZWFtzXsw0z+wcR4nN7D8puiSTsCCFEK1CVnk7h2nUU/vADlUlJpuNpvtZsDa1md5gGCy9PHol4hAc6PYCNRQMNEC5MUwchH/1O3ebB2hHuef+2u21lMCrsOJPNot+S2Hkmx3Q8yN2eaXcEMrGXPzr7BgqWosFJ2BFCiBasLDqG/O++o3DNGjCqCwAa7Ww4GGHLD12KSfQx4mbrwcPhc5gUMgk7S7uGeeHqCtj+Nuz9CAwXF70LGQuj/wku/g3zGq1AYVkV3+1P4rt9yaQVqHtUaTVwd1c/Zg4IomeAiww4bgUk7AghRAujKAolu3aR++mnlB06bDpe0a0j6zqVsCogiwrrEpytdTwTPpsHOz+IvZV9w1WgOBuWT6u5ZRU4EEa8BgH9Gu41WriUC6Us+i2R5QdTKK5QF2DU2amL/83oH0ig+y1upyGalIQdIYRoIZTKSoo2byb3s8+piI1VD1pZUto/gm+6FrDFKQEARysnHg6bwbSwaThZN/D4kPSjsOwh0KeBjQ7GfwSd74bbpPciOU8NOd/uT6KyWu1J6+TtyCOD23FvNz/ZwqGVkrAjhBDNrCorm4Lly8hfshRDQQEAGnt7qu4ZxgedEthXfRIAe0t7Hgx9kFldZqGzucVVj+ty8ntY8yRUl4N7R5i6DDw6NPzrtEAx6YUs2HyWzbFZKBenjt/Rzo0/DGnPkI6eaLW3R9gzVxJ2hBCiGShGIyV79lCwYgVFW7eBQZ3VY+HhgdX9d/NlSAarc36FanCwcmBa6DSmh01vnJBjNMLWN2H3v9XHHUfC/Z+DbSO8VgsTn13MpzvPsfJwqinkDA3xZNaAIIZ08pTxOGZCwo4QQjSh6vx8CpavoGDlSqrS0kzH7Xr2xHHaFL73TuaL2EWU55SjQcOEjhN4qsdTeNh5NE6FygvVBQLPbFAfD3wWRrwKWvO9XaMoCrvO5vLlngS2x9XMrLq7qy9/uqsT7T1v71WgzZGEHSGEaGTGigqKt26lePsO9L/+ilJeDoDW2RnduHE4PDCOn4zH+OLkf8jJUb98e3r15C99/0KYe1jjVSzzJKyYARfOg6Ut3Ptf6Dqp8V6vBTiRWsBLq06a9qvSaCAy1JvHh7SnV6BrM9dONBYJO0II0UjK4+Io+P4HCn/6CWNhoem4bZcuuE6fht3IEfyc+isLjz9DVmkWAH4Ofvyp95+ICoxq3FsoR7+FdX9Wx+fo/GHS19CmZ+O9XjMyGhV2nM3hqz2JpjVyHKwtmNjbn1kDggjykJlV5k7CjhBCNCBDcTH6tesoWLmS8pgY03FLHx+cx4zBadhQbHr1ZGPSRj76dSpJenWBQG97bx7r+hj3dbiv4VY9rkvpBVj/PER/rz7ucBdM+BTs3RrvNZtJRbWBFQdT+HJPIgm5JYDak3NvNz9eHhuGp9PtvTP77UTCjhBC3CLFaKT00CEKV61Gv2EDSpm6+BxWVjgNG4bLxAdwGDAAtFp2pu7kg7WTOJN/BgBXG1ce7fook0ImNdyqx1cTtwHW/gmK0kFjAcPmwaA/g9a89m8qKK3km71JLD2QTHqhesvQydaSyb39mdE/iAD3BlyTSLQKEnaEEOImVaakULh6DYVr1lCVmmo6bt2uHS4TJ6Ibdy+Wbm4YFSPbUrbz6YlPiclTe3scrRyZ1WUW08Km4WDVyLdRSvLgl+ch+gf1sXsHuO9TaNurcV+3iZVXGVj0WyL/2xaPvlxdCNDb2YY/Du3AA73a4mAjX3m3K3nnhRCiHoylpeg3/ErhqlWUHjxoOq51cMB5zGh0992HXY8eaDQaDEYD68+v57OTnxFfEA+ArYUtD4Y+yJzwOY0zjfz3Tv0E6+ZCSQ5otND/KRg6D6zNp3fDYFRYdTSNf2+MM/XkdPZx4rE72zEmwlcWAhQSdoQQ4kZUnD9P/tJlFK5ejbGoSD2o0eDQvz+6+8bjFBmJ1k7dl6rKUMXa82v5IvoL05gcBysHpnaeyrTQabjbuTd+hUvy4JcXasbmeIbC+P+Z1SBko1FhfXQGH26N53Sm+p746WyZOzKE+3q0wUIWAhQXSdgRQoirMBSXUPjD9+jX/0LZ8eOm41b+/rjcfz+6cfdi5etrOl5hqGDV2VV8Gf0lGSUZAOhsdEwLncbUzlObpicHIGYVrHsOSnPV3pyBz8LQF8HSfAbkHky8wFvrYjmeUgCoY3KeHNaBWQOCpCdHXEHCjhBCXEYxGinatJnirVso2rYdo15djwWtFsehQ3GdOhWHgQPQXDaot7SqlBVxK1h8ajG5ZbkAeNh5MKvLLCZ2mtiwm3ReS3G2Op089if1sWeourdVG/MZm3M0OZ//bD5bawr5I4PbMXtgEC721s1cO9FSSdgRQtz2FEWh/ORJijZtpmjjRiqTkkzPWQcH4/rQQzjdFYmVt3et8/SVepbELuHb2G8prFDX0fF18GVO+Bzu63hf48+uqmmAuq/VLy9A2QV1ptXguXDn82bRm6MoCgcSLvC/7efYcTHkWGg1TOrtz5/u6oiXk20z11C0dPUOOzt37uTdd9/l8OHDZGRksGrVKsaPH3/V8j/++CMff/wxx44do6Kigi5duvD6668TFRVlKvP666/zxhtv1DovJCSE06dP17d6Qghxw6rS0tBv2EDh6tVUnI03Hdc6OeEyaSKOg+/Evm+fWr04ABfKL/DNqW9YdnoZxVXFAAQ4BfBIxCPc3e7uxl0n5/eKMtXp5HHr1cfeEWpvjm+3pqtDIzqcdIF//hLHgcQLgBpyJvRow9PDO8oUcnHD6h12SkpK6NatG3PmzGHChAnXLb9z507uuusu/vGPf+Di4sJXX33FPffcw/79++nRo4epXJcuXdi8eXNNxSyl00kI0fCMJSU1s6kOHTId19ja4jR8GI7DR+A4dAgWjlfuj5RenM7imMX8ePZHyg3qrJ8OLh14rOtjjAwciUVT7ielKHDsO/j1JXV/K60VDHkBBv0JmjJsNQKDUWHdyQwW7UngSHIBANaWWib2astjd7Yj0F1WPBb1U+9EMXr0aEaPHn3D5RcsWFDr8T/+8Q/WrFnDzz//XCvsWFpa4uPjc0PXrKiooKKiwvRYf+meuhBC1EFRFEoPHqTwx1XoN25EKS1Vn9BosO/TB+cxY3AeOwYLJ6c6zz+bf5Yvo7/kl4RfMCjq7uRd3LvwWNfHGOo/FK2miRflu3Aefn4WEnaoj327qzOtvLs0bT0amKIobD2dzfxfThOfrfaYaTUwsZc/f7qrEz46uV0lbk6Td58YjUaKiopwc6u9NPnZs2fx8/PD1taW/v37M3/+fAICAuq8xvz586+47SWEEL9nKC6h6NcN5H32OZWJiabj1oGB6CZMUGdTXeOXrKPZR/ni5BfsSN1hOnaH7x08HPEw/Xz6Ne7eVXVRFDiyGDa8BFUl6uadQ+epa+dYtO7e8FPpev6+/hR74vMA0NlZMWdgMFP7+uPlLCFH3Jom/3S89957FBcXM2lSzc66/fr1Y9GiRYSEhJCRkcEbb7zB4MGDiY6OxqmO37TmzZvH3LlzTY/1ej3+/v5NUn8hRMumVFdTtHkL+l83ULxlK0plJQBae3ucx45Bd98E7Hp0v2pQMRgNbEvZxtenvuZo9lEANGi4K/Au5kTMoYt7M/WeFGXBz/8HZzaojwMHwbj/glu75qlPA8ktruBfG+NYdjAFRQFrCy2zBwXxx6Ed0Nm17ttxouVo0rCzZMkS3njjDdasWYOXl5fp+OW3xbp27Uq/fv0IDAxkxYoVPPzww1dcx8bGBhub1j/DQAjRcKrS08lfuZLC73+gOifHdNw6KAiXB+7HdepUtA5XH+tRWlXKqvhVfHPqG9KK0wCw0loxrsM4ZnWZRaBzYKO3oU6msTl/hfICsLCGEa/BHX9s1XtaVVYbWfRbAv/dEk9Rhbq1w9iuvrw4qjP+bjLwWDSsJgs7y5Yt45FHHmHlypVERkZes6yLiwudOnUiPj7+muWEELc3pbqa4p07yV++nJKdu9RgAFi4ueEy4T6cokZhG97lmrebskqyWB63nBVnVpimj7vYuDCx00Smdp6Kp71nk7SlTsU5am/OpZlWPl3hvk/AO6z56nSLFEVh06ks/rE+lsQ8dexURBsdr94TRp8g89t5XbQMTRJ2li5dypw5c1i2bBljx469bvni4mLOnTvH9OnTm6B2QojWpiozk4KV31Pwww9UZ2aajtv364fr5Ek4Rkaitb72AnMnc07yTew3bErcRLWi9iz4O/kzM2wm93a4FztLu0Ztw3WdXg8/Pa2ugmxhDcP+2urH5pzO1PPm2ppxOZ5ONjwfFcIDPduila0dRCOq96emuLi4Vo9LQkICx44dw83NjYCAAObNm0daWhpff/01oN66mjlzJu+//z79+vUj8+I/THZ2duh06tLpzz33HPfccw+BgYGkp6fz2muvYWFhwdSpUxuijUIIM6AoCqX7D3Dh668p3r4djEYALFxc0N13Hy6TJmITHHzNa1Qbq9mavJVvTn3DsZxjpuO9vHsxLXQaw/yHNe308bqU62HjX+GI+m8oXl1gwqfgE9689boFecUV/GfzGZbsT8aoqNPIHxkUzB+HdcBRdiIXTaDef8sOHTrEsGHDTI8vDRSeOXMmixYtIiMjg+TkZNPzn376KdXV1Tz55JM8+eSTpuOXygOkpqYydepU8vLy8PT0ZNCgQezbtw9Pz2bsPhZCtAgV585RtGkTRVu2Un7ypOm4fZ8+uEyahNPIu9BeZwxfcWUxP579kSWnl5jG41hqLRkTPIZpodMIdQ9t1DbcsDMb1QUC9amABgY8DcNfbrWrIJdUVPPtviQ+3BZPUbnaezY63IeXxoTKuBzRpDSKcvEmdyum1+vR6XQUFhbi7Ozc3NURQtwiRVEoP3GC/CVLKPzpZ9NYHI2VFS4TH8D1oYewad/+utdJK05jSewSfjj7AyVVJYA6HmdSyCSmhExp3vE4lyvJg1/nwYnl6mPXIBj3EQQNatZq3ayC0ko+2Xme7/Ylob8YckJ9nXn17jD6t2+CHd9Fq9FU39/SfyiEaDEMRUUUrvmJgu+/p+Ky7WIchtyJ0/AROA0fhuUN9Pgeyz7GN6e+YXPyZoyKersrWBfM9LDp3NPuHmwtW8i6LYoCMT/C+hdqdii/44/q+Bzr1tnzsf98Hk8vPUp2kbrwa7CHA38c2p4JPdtiIeNyRDORsCOEaHZVWdnkL1lC/nffYSxWV87VWFvjPHo0rg9Oxa7b9fd5qjZWsyV5C1+f+poTOSdMx/v79md62HQGthnY9CsdX4s+Xd2h/NJMK68wuPdDaNs6dyg/nJTPx9vPseV0FooC7Twc+MvoztwV6i2Dj0Wzk7AjhGgWxrIyCn78kYLlK6g4c8Z03Lp9e1ynTkV3z91YXJzEcC1FlUXqeJzYJaSXpAPq+jhj241leth0Orl2arQ23BSjUV0FedOrUKFX97S68zkYNBcsrz2DrCXKK67gzbWnWH0s3XRsXHc/3p7QFTvrZh7sLcRFEnaEEE2qKj2dgtWryV+yFENurum4Xc+euM+ZjePw4VfsMl6X1KJUvov9jlXxq0zjcVxtXJnceTKTQybjYefRaG24aXnn4OdnIHGX+rhNbxj3IXi1kAHS9VBYWsVXvyXw1Z5ECsuq0GjggZ5t+cOQdnTwqnuPMSGai4QdIUSjU6eN76dg5fcUbdyIUlUFgFWbNrjNno3u7rFYuLjc0HWO5xzn61NfsyV5i2k8Tntde6aHTWdsu7EtZzzO5QzVsP9j2Pp3qC4DK3sY/gr0+wM091T3eqo2GFl6IJl/bTpDQan6Pnb2ceKf93elm79L81ZOiKuQsCOEaDTGsjIKf/qZC998TWX8OdNxu969cHngAXRjxqC5zuJ/oO5XtTVlK4uiF3Eit2Y8zgC/AcwIm8EAvwFNvynnjcqKgTVPQfoR9XHwELjnfXC79ppALVF6QRlPLTnCkeQCADp5O/J/IzoyOtxXBh+LFk3CjhCiwVVlZZP/3XcULF+OoVDdgkFrb4/zuHtxmTAB2/DwGwon5dXlrIlfw9envia5SF2/y1przd3t72Za6DQ6unZs1HbckuoK2Pke7P43GKvBRgdRf4ce06ClBrOrqKg2sOxACv/ZrPbmONla8nxUCA/2DcDSogUN+hbiKiTsCCEahFJVRcnevRT+9DP6X3+FS7eq2rbFbfo0dBMmYOF0Y2M59JV6lp9ezrex33Kh/AIAztbOTOk8hamdp7bM8TiXS96njs3JuTh9vvPdMOY9cPZt3nrVk6IorD+ZyfxfYknNLwPUfaw+erAnAe6tc2q8uD1J2BFC3JLq/HxyP/wI/dq1pl4cUG9Vuc2YgdOIEWgsbmxcSk5pDt/EfsOKuBWmQcd+Dn7M6DKD+zrch71VC/+CLc6GTa/B8SXqYwdPNeSEjWtVvTmKonA4KZ93fo3jQIIaNr2dbXh6eEcm9fbH2lJ6c0TrImFHCHFTqvPyyF+2jPxvvsVQUACAhbs7zlFR6CZMwC68yw1fK0WfwlcxX7Emfg2VxkoAOrh0YE74HEYFj8JKa9UYTWg4hmo49IU6APnizun0mA53/Q3sW9dO3rEZep7//jjRaXoAbK20PD6kPX+4s71MJRetloQdIcQNUxSFsiNHyF++nKJfNphmVdl06oTX88/h0L8/Gssb/2fl9IXTfHnyS35N+tU0s6qbZzceiXiEO9ve2bIWAbya5H2w7jnIurhvl293GPsvaNu7Wat1M348kspLq05SXmXE2lLLhB5t+L8RHfFzaeYd4IW4RRJ2hBDXZSwtRf/LBi4sXlxrAUDbbl1xmz4D56iRaKxuvPclJjeGj49/zI7UHaZjg9oM4uHwh+nl3avlzqy63O9vWdm6wIhXodesVjedvKLawN9+PsV3+9VB4Hd28uTfk7rh4dg6NyAV4vck7Aghrqo6P58LX35F/tKlNds42NnhPGqUuo1DRMQNX0tRFA5kHmBRzCJ2p+0GQKvREhUYxcMRDxPiFtIobWhwhmo49CVsfavmllXPGTDidXBofZtcxqQX8uIPJzmZVohGA8+M6MjTwzvKVHJhViTsCCGuUBYdQ8Hy5RSuW4dSWgqos6pcp0zGZeLEG9rG4RJFUdiaspVPjn9C7IVYQA05Y4PH8ljXxwjSBTVGExqHGd2yKq2s5p+/nOabfUkYFXCxt+L9KT0Y0qmF7AQvRAOSsCOEMCnZu5fcjxdSeuCA6ZhNWCieTz6J47BhN7SNwyWXenI+OPKBaSFAWwtbxncYz/Sw6QQ4BzR4/RtNcQ5sfg2Ofac+bsW3rABOphbyzLKjnM9VZ7zd3dWXl8eG4aNrgatPC9EAJOwIcZurzs2lcPVqinfspPTgQfWghYV6q2rKZOx69673GJrDWYf58OiHHMo6BICdpR3TQqcxI2wGLrYuDdyCRmRmt6wulFTy0bZ4Fv+WSLVRwdvZhvcmdmNwR+nNEeZNwo4Qt6mqzEzyvvySgpXfo5SpC8ah0eA6dSrujzyMlZ9fva95POc4Hx39iL0ZewF19/GJnSbyaNdHW/5CgL93bitseAly1Ftv+HaDMf8C/z7NW6+boCgKKw+n8vd1sRSWqTPoxkT48I/7InCxb307rQtRXxJ2hLiNKEYjpQcPkb90KUWbNoHBAIBteDi68eNxHHIn1v7+9b5uTF4MHx39iF1p6m7elhpL7ut4H491fQwfB58GbUOjyz0LG1+GMxvUx3auMPxl6DW7Vd6ySswt4aVVJ/ntXB6gbto5b0wod3b0aB2z3oRoABJ2hDBzlUlJlMfGUnrgAEXbtlOdkWF6zr5vX9z/8BgOA25uI82j2Uf59MSnptlVFhoL7ml/D3/o+gfaOrVtsDY0ibJ82PEOHPhU3ctKawl9H4MhL6iBp5WpMhj5dOd5PthylopqI7ZWWv4U2YmHBwXLflbitiNhRwgzZSwrI3/5crLfeReMRtNxraMjzmPH4vrgVGxD6j/dW1EU9mXs47OTn3EwUx3jo9VoGR08mie6PUGgc2CDtaFJGKrUcTnb56uBB6BjlLppp0cL3mj0Go6lFPDiDyc4nVkEwKAOHvzjvgjZz0rctiTsCGFGSo8eJffjj6nOyqYqLc20No6Vnx8OgwbhOGwoDv37o7Wt/6wbo2JkR8oOPj/5uWl2laXWknHtxzEnfE7rml0FoCgQ9wtsegXy4tVjnp3VkNMhsnnrdpOKK6r518Y4Fv2WiKKAq70Vr9wdxn092sgtK3Fbk7AjRAtiKCqiKj0Dm3bBphWJy06coDwmBpf77wcLCzLfeguNRotT5AgASo8dQyktpTzuDCW7dtW6nlXbtrg+9BBuM2fUa9r45aqN1WxO3swnxz8hvkANBTYWNtzf8X5mh89ufWNyAFIPq1PJEy/+edl7wLCXoOdMsGh9/yxWGYx8vTeJ9zefQV9eDcB9Pdrw8thQ3GUVZCEk7AjR3BRFoWjzZqy8vcn+938o3bcPAP/PPqPsxHFy//shAKWHj2DVpg0FS5cBkL9kSZ3Xs+veHaeRI7ENC8W+b9+bDjkGo4H1Cev55MQnJOmTAHC0cmRSyCSmh01vfbOrALJi1GnkcevVxxY20P+PMGgu2Do3b91u0v7zeby8Opqz2WovXjtPB167p4ssDijEZTSKoijNXYlbpdfr0el0FBYW4uzcOv/BEreX6vx8LJyc0Fhaol+/nrS5f77pa9lGRKB1dMDS3QOPJx7Hpn37W6qbwWhgQ+IGvoz+kjP56j5YOhsdD3Z+kGlh03C2boWfsbxz6pick98DCmi00G0qDH0RXFrZ7beLLpRUMn99LCsPpwLg5mDN81EhTOrtL1s9iFajqb6/pWdHiEakVFaCRlNrk8yCHdvJ+ONTYDBQ4WyLjb78mtewcHPDqk0byk+fxtrfH8ehQ7EN7Yyluzt2PXve1PibuhgVI78m/srHxz8moTABACcrJ+ZEzGFq56k4WDk0yOs0qfwk2PkuHFsCijrNnrDxMOyv4NmpWat2sy6tmTN/fSz5peqaOVP7BvCXUSGyZo4QVyFhR4hGYEiJIW32A5Sov3Rj6WRJhZs9FzyccT2aiuXFyVGXgk68L2S4a9G2D+Kh15ZiWVxO6f792Hbpgk27doD6JddYg0y3p2zn/SPvm8bkOFs7MyNsBpNDJreuFY8vKc6Bne/Aoa/AqAYCOo5U18vx7da8dbsFZ7OK+OvqaA4kXADUNXP+fl84vQLdmrlmQrRsEnaEaGiVpeT/bY4p6ABUF1VjUaTHM0kPQJob/BaqodxGwwVH+C1Mg6LRAMkY45cR5h7G2eBc7vF15tLw0oYOOkbFyL70fXwb+61pMUAnaydmhM1gWug0HK0dG/T1mkRlCez7H+x+HyrVade0G6r25Pj3bdaq3YqC0koW7jjPF7vPU2VQsLOy4NnIjswZFIyVrJkjxHXJmB0hGkrqYQq3vsW2fdGEbL3y94gt3TTkO8KpIIjx1+JoVNh0oYL3w4ezNGtvnZf0tPPkh3t/wNW24Ra1UxSFX5N+5b9H/ktyUTKgLgY4o8sMHol4pHWOyTFUq5t0bvsHFGeqx3y7w11/g3ZDmrVqt0JRFL7bn8y7v8aZtnkY0dmLN8Z1oa2rrJkjWr+m+v6WsCPEzagoAo0FWFhTemIFG1cuwHioiNDkmiJFtvDKFBdGHy9kZ7iWs2017J/yG9XZMXx5diUjjv9ERGE2p6ytebCNL/ZWDhRVFdf5cuvuW4eNhQ3eDt43XWV9pZ4VcStYn7Ces/lnAXVMzr0d7mVyyGSCdcE3fe1moyhw5ld1GnnOafWYSwCMeA26TICbnInWEmTry3nxx5NsPZ0NqLesnhsZwohQL1kzR5gNCTv1IGFHNJnKEgxHfiTp2TeoqNayuY8lfvEKYSk1RSosYW1fDb+FaUnxvOxLKWcyb0XOZmyEr7pcf2UpbPwrHPqSREtLnF2DmeLtRkaJup3D6ODR/JLwS62XH9RmEG8MeAMve6/rVlVRFDJKMtBqtCw4soBNiZuoNFYC6mKAj0U8xswuM7G3aoU9BIoCibth+9uQpG5VgZ0r3Pk89HkELFvv2jKFpVV8uO0s3+xLorzKiLWllhdHdWbmgCCZZSXMTosNOzt37uTdd9/l8OHDZGRksGrVKsaPH3/Nc7Zv387cuXOJiYnB39+fl19+mVmzZtUq89FHH/Huu++SmZlJt27d+O9//0vfvjd2j13CjmhI1bm5GPR608BggKpzJ0h+aCqVegMY6/7CSXfVsqODD8sCpuIQsJVq25Om50rOP4uxQl18z8PRhmcjO/Jg3wC0GmD/J7DhLwCcv38hW5Riunp2pa9PXx7f/Di/pf9W63V8HXwZETCCcI9wxrYbq9bPUIWl1hIFhdKqUr4+9TVrz68lpSil1rkdXDrwUOhDDGozqHUuBgiQcRw2vgIJO9THFjZwxxMw6E9g59KsVbtVx1IKeHrpEVIuqLvQd/d34Z/3dyXEx6mZayZE42ixYeeXX35hz5499OrViwkTJlw37CQkJBAeHs7jjz/OI488wpYtW3j22WdZt24dUVFRACxfvpwZM2awcOFC+vXrx4IFC1i5ciVxcXF4eV3/N1gJO6KhKEYj5+8aRmV6NoFffUbWmy9i5WZHRWYGlSmGK8qfCLIg3bYb+YPH0WlgX3oEuNLe05FNyet5afdLpnIHph7lw23xLDuQQl6J2rsys38gb4wLVwv8/Cwc/gq63AcTF5nOK6osYsDSAQCMChrF4azD5JTlmJ53tnamu1d3DmQcwMrCChQoqiq6op52lnb8b8T/6OXdq/XeAilIURcEPLEcUMDCGnpMV0OOS/13am9JsvTlvPtrHD8cSUVRwN/Njr+NC2doJ8/W+34JcQNabNipdbJGc92w85e//IV169YRHR1tOjZlyhQKCgrYsGEDAP369aNPnz58+KG6UqzRaMTf35+nn36aF1988br1kLAjGkr56VgSxk+4ZplkT1gUqSU6SMufwv/JtG6jsLasPTbEYDSwPWU7a86tITIwknvb3wuoy/ov/i2Rt9bFotHAT08OIqKtDtKPwqfDAAWCBsPUpWCj/jb/yfFPOJ5znLfvfJtvT33Lx8c/vm47tBotf+33V4YHDGftubX08elDF48uN/Vn0uzKCmD3v2HfQjBUqMfCH4ARr4BrUHPW7JYpisKi3xJ5Z0McZVVqmJ7Qow2v3hMma+aI24LZLCq4d+9eIiNrb6oXFRXFs88+C0BlZSWHDx9m3rx5pue1Wi2RkZHs3Vv3DJWKigoqKipMj/V6fcNXXNxWSvYfIPPNv1GRlEhdv0cXOMCSoVq2d60JNY5WjkzrHol1HVN/LbQWjAgcwYjAEbWOW1loeWRwO6LTCll9LJ3XformhycGoPHroa7mu32+ul/Toa9g4P8B8IdufzCdP7jNYFPY6ebZjWR9MgbFwDD/YUwPm05KUQq/JPzCUP+h3NP+HgBmhc+6xT+dZlJdCQc/V9fLubQbeeAgGPkmtOnZvHVrACkXSnl1TTTb4tSeul6Brrw8NpQeAQ03804IoWr0sJOZmYm3d+0ZJN7e3uj1esrKysjPz8dgMNRZ5vTp03Vec/78+bzxxhuNVmdxe1CMRvK/W0Le2nVUnjqJtspQK+hUWML/xmrZG1Y7zLxz5zs4WTvhauuKtcXN/fb94uhQNp7K4khyAauPpXFfj7Yw5C9QWQy//RcOfQH9/nDFQNsIzwjeHPgmLjYuDPUfesV1Q9xCiAxsnTt2mygKxKyCLW9AfqJ6zCNEnUbeKQpa+W2dskoDn+06z/+2x6sDkC20vDRGHYAst6yEaBytclHBefPmMXfuXNNjvV6Pv3/rvmcvmoZSWUn2v/5FeXYOFw4cxDYvF4BLceaMH2S5avh2mJZSG6iwVr98dkzewcmckwQ6BxKkC7rlevjobHlyWAfe/TWOt385zcgwHxxsLNXAc/Rb9Us++kfoPvWKc8d3GH/Lr99iJe6BTa9A2mH1saO3uht592mtcjfyyxmNCj8eTePdX0+TpVd7pvsFu/H3+8Lp4CUDkIVoTI3+r4ePjw9ZWVm1jmVlZeHs7IydnR0WFhZYWFjUWcbHp+7ZIjY2NtjYtN6ppaLpVSYlUZGYyOkvPsPxgPpFevmOUl+M1LKlm4Zqy5rfrPv49CE2L5aJIRNxs3VjiH/DLk738KBglh9MIflCKR9ti+eFUZ3VcTp9/wA73oZTq+sMO2Yp+7Tak3NpN3IrB/U2Xv+nwKYVruT8O4eT8vnbzzEcTy0EoI2LHX8Z3Zl7uvpKb44QTaDRw07//v1Zv359rWObNm2if//+AFhbW9OrVy+2bNliGuhsNBrZsmULTz31VGNXT5gxY3k5F44eJ+mb77DfugmAS1+bFxzhrJ+GJUO1ZLjXfNmEuYdxKu8UAJEBkXwZ9WWj1c/WyoKXx4by2DeH+XxXApP7+BPo7gBh49Swc26rOji3lU+nvqbCVHWc0rEloBjVhRp7zYQhL4LTzS+g2FLEpBfywZaz/Bqj/jLnYG3BU8M7MntgELZWFs1cOyFuH/UOO8XFxcTHx5seJyQkcOzYMdzc3AgICGDevHmkpaXx9ddfA/D444/z4Ycf8sILLzBnzhy2bt3KihUrWLdunekac+fOZebMmfTu3Zu+ffuyYMECSkpKmD17dgM0UdxOFIOBkj17SFu2EsO2LWgUhUtL5mW6QJ4zbOilZX/nmnE4fxvwN8a0G0NsXiw6Gx33rlZnTvk7Nf6t0bvCvBnc0YNdZ3P5x/pYPpneG7xCwSsMsk+pa/AM/Uuj16PJlV6A3f9R23dphlXnu2HEq+AZ0rx1awDFFdV8uvM8/9sWT7VRQaOBib3a8lxUCF5ODbNLvRDixtU77Bw6dIhhw4aZHl8aOzNz5kwWLVpERkYGyck1a+YHBwezbt06/vSnP/H+++/Ttm1bPv/8c9MaOwCTJ08mJyeHV199lczMTLp3786GDRuuGLQsxO8pikLl+fMU79rFhR07KT98CMtKdQ+hS/01u8M0nOqkkNbOyCgrLwqdHIh0D+Xvg/5ea/Xg7l7dqTJU4WLjQkFFAaHuoY1ef41Gwyt3hzHyPzvZdCqLwtIqdPZWcOdz8P0c2Psh9HtMXR3YHFSVwf6FatApV2/pEDgQIl9v1Rt1XqIoChuiM3ntpxiyi9QQNzLMm+eiQujkLeNyhGgusl2EaHUMxcWUHT+Ofu06ig4fwphce5VgvR3sCdOwO0JDe/tyBhsVRnd6AKs+j4Bnp+teP1mfTLWxmnYu7a5btqHc9e8dnM0u5oOpPbi3mx8YjbBwEGTHwB1/hFHzm6wujeLSRp3b34aidPWYVxc15HS8q9XPsAI4mHiB+etjOZJcAECQuz3PR3VmTISPjMsR4irMZp0dIW6VoaiIwlWrKT1yhPITJ6hKT6/9vAbi2sKRDlrSA4yE25bSv6Kcx106Yt/3MQi/H6wdbvj1ApwDGroJ1zU81Iuz2cVsjMlUw45WC1FvwTf3qbd6+jwC7u2bvF63TFHg9FrY8jfIPaMe0wXA8L9CxETQtv5xK/HZxfxzw2k2nVLH5dhZWfDo4GD+OKyDjMsRooWQsCNanOr8fCrOnqV03z5K9h+g7MgR9UvzMno72B+i4WSQhuy2RgYbS5laXEqXag2aoPuhz8PQpler6TEY1cWHT3acZ9vpbMqrDOqXZPvh0OEuiN8ES6fCo1tb18ykpN9g06uQelB9bOd2caPOh1v1Rp2XFJZV8Z9NZ/hmXxIGo4JWA5P7BPCnyI54Ocu4HCFaEgk7otkZS0spPXKUos2bKDt8mIqz8VeUSXXXsLczJHtpOBGswdpSQ1RJEc8UldAltxKNW3sY8mfoMQ3s3ZqhFbemW1sX/HS2pBeWsyE6k/E92qhPjJoPnx+E3Dh1KnqPac1azxuScwY2v3bZNHJ76P8kDHgabHXNW7cGYDQq/HQ8nbfWxZJbrI7LiQz15sXRIbJejhAtlIQd0eQUo5HSffso3rmLsuPHKT91CuWy7T8ALjhZc6qNkdMBRo500JDrDK7V1owt1fNUXiFdKyrRai3VGTy950Dwna2mF6cuWq2GKX0D+PfFngJT2PHoCAOeUjfAjP6hZYedwlTY+R4c+RoUgzqNvOcMGDrPLKaRg7rFwzPLjprG5bTzdOBv94YzqKNH81ZMCHFNEnZEozKWlFB2/DgV5xMojz1FRexpqtLSMBQW1ipX7ubM4SBbdgcVcK6NgQJHIwB21bbcW1rBhPRMQiur1BlWLgEwaJa6qq6ZfIkCTOnrz/tbznI4KZ/zOcW087x4y6rLBDXsnN8B+gxw9m3eiv5ecTbsfBcOLwKDuqM7IWMg8o0bGhDeGmTry/lwWzwrD6VSVmXA0caSx4e045HB7WRcjhCtgIQd0aCq0tMp3rOHyvh4ys+cofTAQTAYriincXIkOSKYnd6lHHFNIcWrBDSlANhW2XN3mYb7CpPoVVaGBYClHXS9H3o8pG4Gqb1y883WzsvJloEdPNh5JoeNp7J4fMjFsOPeHvzvgJR9sOOfcM+CZq2nSXkh/PYh7P0IqkrUY4GD1O0dggY2b90aSHmVgW/3JbFg81mKK6oB6Bvkxr8mdcPfzf46ZwshWgoJO+KmVefkUJmairGklPLok5QdP0Hx7t1QVVWrnKWnJzYhIVQEePGbcxnbK+M47pZMlVWsqYx1pQvDqqyZpo+nW3lyzYacAQOg+4PqqsK25r+swF2hXuw8k8PW2GweH3LZ7KvI1+Cr0eotoj6PgE9481WyoggOfAZ73ofyAvWYX0+1ju2GNl+9GlCVwcjKQ6l8sOUsmfpyALr5u/D8yBAGdnCXqeRCtDISdsQNqcrKwlhURMX585QdP07xtu1Unj9fZ1nbiAjse/fG0s+X+AAnlpbHsj97J8XGfbXKWVZ4c6fiyMzic/QoOlkTcHQB6p5Q3aaAW9OtddMSDOvsBWtiOJR0gQsllbg5XNxVPXAAhN4LsT+piw3OXg8OTTxOpLwQDnyq9uSU5avHPDvDsL9C6D2teszUJYqisPV0Nv9YH8u5HLW3yk9nyzORHZnYyx+ttvW3UYjbkYQdcYXqnBwqzidQEReHUl1Nyf59lOzYWWdZC3d3LHQ6bENDsQ0PxyoshMPeVXwb/QvH89ZSkZhvKqsoWpTSIIZauvNoVQrdcg6i4eKUcit7tfem+4Nme5vqRrR1tSfU15nYDD3b47KZ0LNtzZNj3oPEXerMrIWD4elD9Vo/6KaV5atr/ez7X82qx+4d4M4XIOIBs1grByA6rZB/rI/lt3N5ALg5WPP08A482C8AG0vzaKMQtysJO7cxRVGoSkmh4swZqnNyMBQWUnbiJMVbt9Z9gpUVVr6+2IaE4Dx6FPZ9+mDp6UlxZTG703azIWEzu898TMXp0prXMFhjXRXGYFt/ntBkElK4EU1FUc01AwfW3KaykWm7AEM6eRKboWd3fG7tsOPkDVOWwFdj1FWIj34L/f7QeBUpyoJ9H8HBL6CyWD3mEQJDXoAu95lNyInPLuY/m86w7mQGANaWWuYMDOaPw9rjbGvVzLUTQjQECTu3EWNJCSX7D1Cyby+le/dRnZeH4cKFOstaBQZg3dYfrb091sHB6MaPx6ZdsOn5rJIsfknZypajWzmQeQCDUl3zOtWOWJSF09M+lAcNaURWbEOb8VPNxXUBasDpNgXcghG1DezgzsId5/gtPg9FUWqPDwkcAGP/BevmqruFBw5s+PE7Bcmw5wM4+g1Uq+NV8Oqi7tcVNt5set0yC8tZsPkMKw6lYFTUu3D3dPXj+agQGXwshJmRsGOmSg8douDHVWA0olRXUxEfT+W5cyi/GzyssbbGOjgYKx8fLFx0WAcH4zhsOLYhtacMGxUjMbkx7EzbyfbkHZy6EFPreUOFJ9VFYbTVduVJt2LuddyBVeJKqHWbavzF21QDzeYLszH0CXLD2lJLpr6c87kltPf83arJ3R+E3QugMBkWjYFZ68An4tZfOPu0Ouj45AowXgyvbfvA4OegU5RZjMkBSC8o47Nd51l6IJnyKnWJg8hQb56L6kRnH/MfBC/E7UjCTitmKC6hMuE8xrIySg8epCJOvR1VmZyMIS+vznOs2rTBvndvtE5O2PfqieOQIWjt6/4ttqSqhL3pe9mWso2dqbspqKjpBVIUDcYyf6qLuqBTujHLV2Gq025cE99Ek6CvuUjgoIu3qe6V21Q3yNbKgl4Bruw9n8dv8blXhh0rO3hkM6yYDin74evxMOU7CLij/i9mqIIzG+DwYnVbikvaDYXBf4agwWYTcvTlVSzcfo7PdydQWa2GnN6Brswb05lega1v1W0hxI2TsNPCGYqLMRYXU7xrF8biEspPnqAqMwtjWRkV8fFXTPOuxcIC+969se/XF+u2bbGNiMA6KOiq02YVReFM/hl2p+1mT/oejmQdrXV7SjFYU13SEUNJCMaSMP4Q7s4jwQdxPfMvNOfP1VzIJQC6yW2qWzGwgzt7z+exJz6P6f2Drizg5A0ProCv74WM47BoLETNh76PXj+cKApkRcOpNXDkGyjOvPiEBkLvhoF/gra9GrpJzSa7qJzPdp7nu/3JlFaqaz71C3bjyWEdGNzRQ6aRC3EbkLDTAhj0eqoyMik/dYrq7Gwq4uLUlYejo6/aQ3M5jZ0dduHh2Pfti037dlj5B2DVtg2Wrq7XPTejOIO9GXvZl7GPg5kHyS3LrfW8sdKd6qLOVBeH4mnVmf7+Tvwx9DTt05eijd1OzW0qB+gyHrpNldtUDWBABw/YeIbfzuVSbTBiaVHHn6edC8xaDz89BTGr4Jfn4dCX4N9Xva3VKQosrCH9KFSWgGKEhB0QvwWKMmquY++hLtbYc2br3Fn9KhJzS/h6bxLf7U+i4mJPTkcvR14Y1ZnIUC8JOULcRiTsNDKlqoqqjAwMBQWgKFRlZ1O6dx/lsbHqOjXWVhhycq97HQDr4GAc+t+BXbduaJ2csQ4MUHtqLG58VkxxZTEHMw/yW/pv7MvYR6I+sdbzGsWaqpJ2VBeHUF3cES+7NnRto2N8SApjDWshehXEXzabKmiwepsq9N7WtSN3C9etrQtuDtZcKKnkQMIFNfzUxcYRHvhKHVuz5U3IiVV/ANY/d/UXsLKHoEFqOO18N1haN3wjmsn5nGLe2RDHhphM07EeAS7834iODO3kKSFHiNuQhJ0GUHH2LFXZ2WA0UrxjJ4b8fNBoKI+Opioj44pNLq/GQqfDYdAgrAMD0VhbYR0QgH3fvmgdHNDY2t7UP9KFFYUcyTrC0ZyjHM06ysnckxiUmu0bNGjRVARQpm+HobQ9hrJAUCwJdLfniSG2PGC5C8sTyyDh8ttUgTWzqVyD6l0ncX0WWg13hXqz/FAKG2Iyrx52QL1t1f9JNbic26reokrYBWmH1Od1/ur7ZKhUQ1GHEWrvm6VNk7SlqZxK1/PJznP8fDzdNLvqzo6ezBkUzJ1yu0qI25qEnWuoSkvDWFaGIT+fisREjCUlVJ47R3V+PkppKdUX8jEWF1OVmnr9i2m1YDRi07EjVoEBaG3tsPT2wmHAACwcHLAND0djeetvR2lVKcdyjnE0+yh70vYQnRuNculW00XOlr6UFbZDnx9MdUl7MNoB6vouD3TzYJByBNe4xbBnk3rrAy7eprpPDTkB/eU2VRMYFe7D8kMp/BqTyev3dLn+6r32buoifxEPqI+NRvUb38y/5I8m5/PBlrNsi8sxHRvR2Yu/jO5MJ28ZFC+EuE3DjqIoVKWmUnrwEOUxMWjt7dA6OlGZkkxVYhLVBflUpWeglJZe/2IXaezt0Wi1OPS/A5uQzhhLS7Hy9sL+jjuwDg5Wv3iMxqvOfLpZFYYKjmcf50DmAQ5kHuBk7kmqjdW1yvg5BOBnG0ZcojsZWW0oqqqZeTK4owf3RPhwj8t57GK/ho0/QcVlO5IHDoIe09TtAOQ2VZMa0MEdRxtLsvQVHE8toEfA9cdg1WLmgTSvuIL/bD7Dt/uSAbU3bFS4D08MaU94G10z104I0ZKYVdgpO3ECioqoSEiA6moUgxFjSQnVWZlUJiZSnZOLoaAAjZXVFevNXItVmzZYenujtbXFOigQS19fLFxc0Ggt0NjaYBsaik27ptnDqaSqhOM5xzmafZQjWUc4ln2MSmNlrTI+Dj509+yJTVVHMjKD2HKogrjLnvd3s+P+Hm2Y1DYf36Q1aHatUlfkvcS5DURMhB7TwaNDk7RLXMnG0oJhnb34+Xg6G2Iy6x92zJCiKOw7f4Ev9ySwOTYL5WKn5YSebXh6eEeCPZpg+wwhRKtjVmEnefYcHG9gsK5SVQVaLRY6HVonJ7T29li6u2MbEY6lpydaewcsPT2xC++ChYtL41f8GjKKMziSfYSj2Uc5nnOcM/lnMF66tXSRp50nfXz60NW9Nzo68/PhClb+knXxWXW8kFYDD/Rqy6PhWtpnb0R7ch7sOV1zEVudumVDxCSZTdWCjOriw8/H0/k1OpMXR3W+rcedHEy8wLu/xnEgoWa9p1BfZ14eG8rAa41pEkLc9swq7FjonLFyd8eoL8K+Tx+0zk5Y6HRYeXlh1bYtGktLFEXBpkNHLD3c0draNneVaymrLuNEzgkOZR3iZM5J4gviySrNuqJcG8c2dPfqTk+vnvTx6YOj1pfPdyfw6oYEDMZkUzkPRxvu7urLAyFWhBdsg5P/gWUHay5kYQMho9VenI53md2AVXMwNMQTa0stiXmlxGUV3ZYr/CbklvDGzzFsvzgmx9pCy8TebZk9MJgOXnJrVQhxfWYVdjps3Yqzc+v4MiiqLOJM/hni8+OJyYvhcNZh0orTas2UArDQWNDZrTM9vHrQ3as7Pbx64GXvRVJeCb/GZPL81gwOJ58ydecDjO3qy8O93elWvBuL6P/Csp01A43RQPBgCH9A7cmxc2myNov6c7Cx5M6OnmyOzWJDdOZtFXZOphbyzw2n2Xc+j2qjgqVWw8Te/jw9vAN+LnbNXT0hRCtiVmGnpSqrLiPuQhzRudFE50UTkxtzxfo2l3jbe9PLuxc9vXrSwbUDoW6h2Fupg5r3n8/jvfWpxGac4WRaYa3zurbV8X9DgxhhfQrNiY9gxTqoLqsp0LYPdJkA4RPAyaexmioawahwH1PYeTay0/VPaOUKS6t4f8tZFu9NxGBUU/yQTp68fm8XGZMjhLgpEnYaWJWxivj8eFOoic6NJr4g/ooeG1AHErd3aU+oWyg9vHoQ4hqCl33Nyq5Go8L+hAvsPJvM/vN5HEkuMJ1rqdVwRzt3RoZ5MkaXhEfCj7B+NZRetuKye0foOlmdiizbNrRakaFeWGg1nM4sIimvhEB38/zCT8orYcn+ZJYeSEZfrs4oHBvhy59HdqLd7/cHE0KIepCwcwuKK4uJL4jnTP4ZU6iJuxB3xewoAA87D8LdwwnzCCPcPZwuHl1ws71y88H0gjL2J+RxKDGfPfG5JObVnv6us7PixVEhRLln4Xb+J9i3CvSXrfNj7wHh96sL/vn1MPs1Vm4HLvbW3NHOjT3xefwak8ljd5rPlg4ApzP1fLLjPD8dTzf15HT0cuTlu8MY0smzmWsnhDAHEnZuQJWxihR9CvEF8cQXxHMy9ySxebHklde9b5WTtRNd3LsQ7hFuCjbe9t51zqTJL6nkhyOpHE7KJ62gjBOptW9PWVtoGdnFm25tXehincEdpTvQHpgHefE1hWyc1XVwwu+H4CFgIW+ruRnVxYc98XlsiDafsFNZbeTDbfF8tC3eFHLu7OTJtH4BjAj1xuJ6iygKIcQNkm/Fy1Qbq0nWJ3Ou8BzxBfGcLzhPfEE8ifrEKxbqu8TL3otg52C6e3Wng0sHurh3oa1T22vuLH44KZ8fj6YRn13MseQCKg21p5J39nGiras9I0M9Ga5LxSN1HZxcCzmXTRW3tIVOo9RbVB3uAquWNbNMNKyRXXx49acYjiQXcC6nmPat+LZOcUU1b/8Sy4qDqaa/+yPDvHlqeAe6tnVp3soJIczSbRl2yqvLSdInkVCYwPnC88TkxZBYmEh6SfpVQ42dpR3tde3VMTbuoUR4RBCkC8LZ+tqzY4xGhfO5Jfx4JJXodD2nM/RkF9XeK8vdwRqdnRXjurchqpMzncuPQdyPsPMXKK7ZzBCtFbQfpvbgdB4LNrIU/u3C29mW4SFebDmdzTd7k3j93i7NXaV6K62sZsn+ZD7ded70GfBwtOa1e7pwTze/Zq6dEMKcmXXYuVB+gZjcGJKLkkkvTidRn8j5gvOkFaddsV/UJXaWdnRw6UA7XTvau7Q3/fg6+KLVXH+hPUVRSC8sZ+3xdE5nFnEspYCE3JJaZeytLRgd7ku/YDd6BrrS3voCmvjNEPcR7NsJ1eU1ha0d1Y0bO9+jroUjU8VvWzMHBLHldDbfH07luagQHG1ax8c3MbeEr/Yk8OPRNIouDjz2d7PjzXHhDO7oKberhBCNrnX8a3mDPjj8ARnGDFKKUsguzaakquSqZZ2tnWmna0c7l3Z0dOlIiFsIfo5+Nxxqfi85r5TYTD1f7ErgQOKFK553tbdiYm9/BnbwoE8bW+zT9kLcV7BnKxQk/a5ybaHTSAgZq66JI4v9CWBQBw/aeTpwPqeEBZvO8PLdYc1dpWtKyC3hv1vPsuZYzcDjIHd7Hh/Snvt6tsHG8vqrnQshREPQKIpSdxfHNXz00Ue8++67ZGZm0q1bN/773//St2/fOssOHTqUHTt2XHF8zJgxrFu3DoBZs2axePHiWs9HRUWxYcOGG6qPXq9Hp9MR+nEoFna1/wEN1gXTXtceX0dfAp0CaefSjmBdMO627je99H61wUhhWRXR6Xp2ncnhl+hM0grKapXxcLRm9sBg/FxsifBzwr/yPDaJ2+DcVkg5AIbLbmVpLKBNLwgZpY7D8QqTWVSiThtjMnnsm8NYW2o5/upI7KxbXmBIyivh/c1nWX0sjYsZh6EhnsweGMzgDh7X371dCHHbuPT9XVhY2KiLAte7Z2f58uXMnTuXhQsX0q9fPxYsWEBUVBRxcXF4eXldUf7HH3+ksrJmKnZeXh7dunVj4sSJtcqNGjWKr776yvTYxqb+vRmTOk0i1C8Ufyd/3O3caePYBkfrhhnIWVxRTWyGnnUnMlhzLI380tobiVpoNQS629O1jY7p/YPo6qnFKnkPxP4EmzdDSU7tCzq3gU5RargJHCDjb8QNuSvMGz+dLemF5eyJzyUyzLu5q2SSVlDG57vO8+2+JKoMasoZ3tmL/xvRke7+Ls1bOSHEba3eYeff//43jz76KLNnzwZg4cKFrFu3ji+//JIXX3zxivJubrXXklm2bBn29vZXhB0bGxt8fG5tZd+5vec2aDJceyKdT3acB9S1QC79A36JtYWW3kGuTOrtT2RnDxxzjsLZn2HTTkg/CpcPdrZygOA7of1waDcUPDpK742oN41GQ1S4D1/tSeSbfUmMCPVq1s1BDUaFHWey+W5fMtvisk09OYM7evB8VIjMrhJCtAj1CjuVlZUcPnyYefPmmY5ptVoiIyPZu3fvDV3jiy++YMqUKTg41F4Fdvv27Xh5eeHq6srw4cN56623cHd3r/MaFRUVVFTU3AbS6/X1aUadyioNHEq6wLIDKZRXGUi+UMrZ7OJaZXycbQnzc+auMG+6+7vQ1rYcp4x9cOYb2PRL7dWLAVyDoGOUugaOfz+wtL7legrxYN8AvtmbxI4zOXyzL4kZ/YOavA6KonAwMZ+31p2qtTZU/3buPDG0PXfKYoBCiBakXmEnNzcXg8GAt3ftrnNvb29Onz59lbNqHDhwgOjoaL744otax0eNGsWECRMIDg7m3LlzvPTSS4wePZq9e/diYXHlmIT58+fzxhtv1KfqtSiKgqLAxlNZ7InPJaeogp1ncyitvHJLh7audjzYL4AB7T3o7m6EpD2Q+AMc2Q2Z0XD5rC5bHXSIVHtvggaDa+BN11GIq+no7cSLozvz1rpY5q8/TWSod5NujJmaX8rzK0+w97wa7p1sLZnc25+p/QJa9fo/Qgjz1aSzsb744gsiIiKuGMw8ZcoU0/9HRETQtWtX2rdvz/bt2xkxYsQV15k3bx5z5841Pdbr9fj7+1/ztQtLq5j/Syy5xZUcSMijvNpIZXXtxfycbS2pNio8PCiYiDY6erhV4pl3CJLWw897IPvUlRf2CFFvS4XeDQH9wcLqBv4khLg1Dw8K5teYTA4m5vO3n0+xcHqvRn/NM1lFfL030bQYoLWllvt7tuHp4R1lF3IhRItWr7Dj4eGBhYUFWVlZtY5nZWVdd7xNSUkJy5Yt429/+9t1X6ddu3Z4eHgQHx9fZ9ixsbG55gDm05l6otP0HEnOZ0dcDvqyKoorq6lr3pmnkw3392zL6HAfujqXoEn6DZL+B1v3QN7ZOk7oDIED1UHFQYNkB3HRLDQaDW+OD2fsB7vZEJPJqqOp3NejbYO/zqXNaBf9lsCvMTWf+zvaufHP+7ua7aakQgjzUq+wY21tTa9evdiyZQvjx48HwGg0smXLFp566qlrnrty5UoqKiqYNm3adV8nNTWVvLw8fH1961M9Xlh5nOjcKlLzy+p83tXeihAfJyJDvekb7EZn2wKsU3+DxFXw427IT/zdGRrw7qKGm6CB6n8dPOpVJyEaS2cfZ54Y0p4Pt8Uz78eTRLTR0cGrYWb1XRp4/J9NZzmZpo7JsdBqGBbixSODg7mjXd3j6YQQoiWq922suXPnMnPmTHr37k3fvn1ZsGABJSUlptlZM2bMoE2bNsyfP7/WeV988QXjx4+/YtBxcXExb7zxBvfffz8+Pj6cO3eOF154gQ4dOhAVFVWvuq2PzkRrY49Gg/oPv6cjY7v64u1si7ONJf5koEnaA0mL4NAeKEypfQGNFny6qj02gQMh4A6wv3JnciFairl3deJYSgG743OZuHAvH0ztweCONz84uLC0iuWHkvl6b5LplwZHG0vGRPgwZ1AwnX0abx0MIYRoLPUOO5MnTyYnJ4dXX32VzMxMunfvzoYNG0yDlpOTk9Fqa69AHBcXx+7du9m4ceMV17OwsODEiRMsXryYgoIC/Pz8GDlyJG+++Wa919rp4e/C3b3bM6FnG9wdrCH3DCSuhZg9kLin9j5TAFpL8Otx8bbUQAjopw4yFqKV0Go1/HtSN2Z8eYDTmUXM/PIADw8K5k93dcLe+sY+3oVlVeyJz2VzbBa/nMykrEodqO9ib8WEHm3547D2eDjKKt5CiNbrplZQbmlMKzAe/Qnn4nhIPQTJ+6A0t3ZBC2t1peJLt6Xa9gUbmT0iWr+KagOvrI5mxaFUADwcbXh4UDCDOnhgb2OBh6MNOjsrFEUhraCME6mFHE8t4HBiPkdTCkzbOQCE+joza0Ag47q3wdaq5a3QLIQwH021grJ5hZ0XnXC2uWyBNUtbaNun5rZU295gJbNGhPnaejqL136KIeVC7XFrGg1YWWhBgUqD8Yrz2ns6MKSTF6MjfOgd6NqsCxUKIW4fLXa7iBbNNQiCuqu9N237QpuesommuK0M7+zNoA6e/HgklfXRmRxNyqekshqjgmmpBUuthhAfJ7r7u9C1rY4B7T3wd7Nv5poLIUTjMa+enUZOhkK0NkajggLkFFVQUW3AQqvBw9FGbk8JIVoE6dkRQtyySzuM++hsm7kmQgjRfLTXLyKEEEII0XpJ2BFCCCGEWZOwI4QQQgizJmFHCCGEEGZNwo4QQgghzJqEHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpg1CTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCmDUJO0IIIYQwaxJ2hBBCCGHWJOwIIYQQwqxJ2BFCCCGEWZOwI4QQQgizdlNh56OPPiIoKAhbW1v69evHgQMHrlp20aJFaDSaWj+2tra1yiiKwquvvoqvry92dnZERkZy9uzZm6maEEIIIUQt9Q47y5cvZ+7cubz22mscOXKEbt26ERUVRXZ29lXPcXZ2JiMjw/STlJRU6/l33nmHDz74gIULF7J//34cHByIioqivLy8/i0SQgghhLhMvcPOv//9bx599FFmz55NWFgYCxcuxN7eni+//PKq52g0Gnx8fEw/3t7epucURWHBggW8/PLLjBs3jq5du/L111+Tnp7O6tWr67xeRUUFer2+1o8QQgghRF3qFXYqKys5fPgwkZGRNRfQaomMjGTv3r1XPa+4uJjAwED8/f0ZN24cMTExpucSEhLIzMysdU2dTke/fv2ues358+ej0+lMP/7+/vVphhBCCCFuI/UKO7m5uRgMhlo9MwDe3t5kZmbWeU5ISAhffvkla9as4dtvv8VoNDJgwABSU1MBTOfV55rz5s2jsLDQ9JOSklKfZgghhBDiNmLZ2C/Qv39/+vfvb3o8YMAAQkND+eSTT3jzzTdv6po2NjbY2Ng0VBWFEEIIYcbq1bPj4eGBhYUFWVlZtY5nZWXh4+NzQ9ewsrKiR48exMfHA5jOu5VrCiGEEEJcTb3CjrW1Nb169WLLli2mY0ajkS1bttTqvbkWg8HAyZMn8fX1BSA4OBgfH59a19Tr9ezfv/+GrymEEEIIcTX1vo01d+5cZs6cSe/evenbty8LFiygpKSE2bNnAzBjxgzatGnD/PnzAfjb3/7GHXfcQYcOHSgoKODdd98lKSmJRx55BFBnaj377LO89dZbdOzYkeDgYF555RX8/PwYP358w7VUCCGEELeleoedyZMnk5OTw6uvvkpmZibdu3dnw4YNpgHGycnJaLU1HUb5+fk8+uijZGZm4urqSq9evfjtt98ICwszlXnhhRcoKSnhscceo6CggEGDBrFhw4YrFh8UQgghhKgvjaIoSnNX4lbp9Xp0Oh2FhYU4Ozs3d3WEEEIIcQOa6vtb9sYSQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCmDUJO0IIIYQwaxJ2hBBCCGHWJOwIIYQQwqxJ2BFCCCGEWZOwI4QQQgizJmFHCCGEEGZNwo4QQgghzJqEHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpg1CTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCmDUJO0IIIYQwazcVdj766COCgoKwtbWlX79+HDhw4KplP/vsMwYPHoyrqyuurq5ERkZeUX7WrFloNJpaP6NGjbqZqgkhhBBC1FLvsLN8+XLmzp3La6+9xpEjR+jWrRtRUVFkZ2fXWX779u1MnTqVbdu2sXfvXvz9/Rk5ciRpaWm1yo0aNYqMjAzTz9KlS2+uRUIIIYQQl9EoiqLU54R+/frRp08fPvzwQwCMRiP+/v48/fTTvPjii9c932Aw4OrqyocffsiMGTMAtWenoKCA1atX178FgF6vR6fTUVhYiLOz801dQwghhBBNq6m+v+vVs1NZWcnhw4eJjIysuYBWS2RkJHv37r2ha5SWllJVVYWbm1ut49u3b8fLy4uQkBCeeOIJ8vLyrnqNiooK9Hp9rR8hhBBCiLrUK+zk5uZiMBjw9vauddzb25vMzMwbusZf/vIX/Pz8agWmUaNG8fXXX7Nlyxb++c9/smPHDkaPHo3BYKjzGvPnz0en05l+/P3969MMIYQQQtxGLJvyxd5++22WLVvG9u3bsbW1NR2fMmWK6f8jIiLo2rUr7du3Z/v27YwYMeKK68ybN4+5c+eaHuv1egk8QgghhKhTvXp2PDw8sLCwICsrq9bxrKwsfHx8rnnue++9x9tvv83GjRvp2rXrNcu2a9cODw8P4uPj63zexsYGZ2fnWj9CCCGEEHWpV9ixtramV69ebNmyxXTMaDSyZcsW+vfvf9Xz3nnnHd588002bNhA7969r/s6qamp5OXl4evrW5/qCSGEEEJcod5Tz+fOnctnn33G4sWLiY2N5YknnqCkpITZs2cDMGPGDObNm2cq/89//pNXXnmFL7/8kqCgIDIzM8nMzKS4uBiA4uJinn/+efbt20diYiJbtmxh3LhxdOjQgaioqAZqphBCCCFuV/UeszN58mRycnJ49dVXyczMpHv37mzYsME0aDk5ORmttiZDffzxx1RWVvLAAw/Uus5rr73G66+/joWFBSdOnGDx4sUUFBTg5+fHyJEjefPNN7GxsbnF5gkhhBDidlfvdXZaIllnRwghhGh9WuQ6O0IIIYQQrY2EHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpg1CTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCmDUJO0IIIYQwaxJ2hBBCCGHWJOwIIYQQwqxJ2BFCCCGEWZOwI4QQQgizJmFHCCGEEGZNwo4QQgghzJqEHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLN2U2Hno48+IigoCFtbW/r168eBAweuWX7lypV07twZW1tbIiIiWL9+fa3nFUXh1VdfxdfXFzs7OyIjIzl79uzNVE0IIYQQopZ6h53ly5czd+5cXnvtNY4cOUK3bt2IiooiOzu7zvK//fYbU6dO5eGHH+bo0aOMHz+e8ePHEx0dbSrzzjvv8MEHH7Bw4UL279+Pg4MDUVFRlJeX33zLhBBCCCEAjaIoSn1O6NevH3369OHDDz8EwGg04u/vz9NPP82LL754RfnJkydTUlLC2rVrTcfuuOMOunfvzsKFC1EUBT8/P/785z/z3HPPAVBYWIi3tzeLFi1iypQp162TXq9Hp9NRWFiIs7NzfZojhBBCiGbSVN/f9erZqays5PDhw0RGRtZcQKslMjKSvXv31nnO3r17a5UHiIqKMpVPSEggMzOzVhmdTke/fv2ues2Kigr0en2tHyGEEEKIutQr7OTm5mIwGPD29q513Nvbm8zMzDrPyczMvGb5S/+tzzXnz5+PTqcz/fj7+9enGUIIIYS4jbTK2Vjz5s2jsLDQ9JOSktLcVRJCCCFEC1WvsOPh4YGFhQVZWVm1jmdlZeHj41PnOT4+Ptcsf+m/9bmmjY0Nzs7OtX6EEEIIIepSr7BjbW1Nr1692LJli+mY0Whky5Yt9O/fv85z+vfvX6s8wKZNm0zlg4OD8fHxqVVGr9ezf//+q15TCCGEEOJGWdb3hLlz5zJz5kx69+5N3759WbBgASUlJcyePRuAGTNm0KZNG+bPnw/AM888w5AhQ/jXv/7F2LFjWbZsGYcOHeLTTz8FQKPR8Oyzz/LWW2/RsWNHgoODeeWVV/Dz82P8+PEN11IhhBBC3JbqHXYmT55MTk4Or776KpmZmXTv3p0NGzaYBhgnJyej1dZ0GA0YMIAlS5bw8ssv89JLL9GxY0dWr15NeHi4qcwLL7xASUkJjz32GAUFBQwaNIgNGzZga2vbAE0UQgghxO2s3uvstESyzo4QQgjR+rTIdXaEEEIIIVobCTtCCCGEMGsSdoQQQghh1iTsCCGEEMKsSdgRQgghhFmTsCOEEEIIsyZhRwghhBBmTcKOEEIIIcyahB0hhBBCmDUJO0IIIYQwa/XeG6slurTjhV6vb+aaCCGEEOJGXfrebuydq8wi7OTl5QHg7+/fzDURQgghRH3l5eWh0+ka7fpmEXbc3NwAdcf1xvzDamn0ej3+/v6kpKTcVhugSrul3bcDabe0+3ZQWFhIQECA6Xu8sZhF2NFq1aFHOp3utvpLcomzs7O0+zYi7b69SLtvL7druy99jzfa9Rv16kIIIYQQzUzCjhBCCCHMmlmEHRsbG1577TVsbGyauypNStot7b4dSLul3bcDaXfjtlujNPZ8LyGEEEKIZmQWPTtCCCGEEFcjYUcIIYQQZk3CjhBCCCHMmoQdIYQQQpi1Fhd25s+fT58+fXBycsLLy4vx48cTFxd33fNWrlxJ586dsbW1JSIigvXr19d6XlEUXn31VXx9fbGzsyMyMpKzZ882VjPq7Wba/dlnnzF48GBcXV1xdXUlMjKSAwcO1Coza9YsNBpNrZ9Ro0Y1ZlPq5WbavWjRoivaZGtrW6uMOb7fQ4cOvaLdGo2GsWPHmsq09Pf7448/pmvXrqaF0/r3788vv/xyzXNa+2cb6t9uc/hsQ/3bbQ6fbah/u83hs12Xt99+G41Gw7PPPnvNck3yGVdamKioKOWrr75SoqOjlWPHjiljxoxRAgIClOLi4ques2fPHsXCwkJ55513lFOnTikvv/yyYmVlpZw8edJU5u2331Z0Op2yevVq5fjx48q9996rBAcHK2VlZU3RrOu6mXY/+OCDykcffaQcPXpUiY2NVWbNmqXodDolNTXVVGbmzJnKqFGjlIyMDNPPhQsXmqJJN+Rm2v3VV18pzs7OtdqUmZlZq4w5vt95eXm12hwdHa1YWFgoX331lalMS3+/f/rpJ2XdunXKmTNnlLi4OOWll15SrKyslOjo6DrLm8NnW1Hq325z+GwrSv3bbQ6fbUWpf7vN4bP9ewcOHFCCgoKUrl27Ks8888xVyzXVZ7zFhZ3fy87OVgBlx44dVy0zadIkZezYsbWO9evXT/nDH/6gKIqiGI1GxcfHR3n33XdNzxcUFCg2NjbK0qVLG6fit+hG2v171dXVipOTk7J48WLTsZkzZyrjxo1rhBo2jhtp91dffaXodLqrPn+7vN//+c9/FCcnp1oBqbW934qiKK6ursrnn39e53Pm+Nm+5Frt/j1z+Gxfcq12m+Nn+5L6vN+t/bNdVFSkdOzYUdm0aZMyZMiQa4adpvqMt7jbWL9XWFgIcM1Nwvbu3UtkZGStY1FRUezduxeAhIQEMjMza5XR6XT069fPVKaluZF2/15paSlVVVVXnLN9+3a8vLwICQnhiSeeMO0S3xLdaLuLi4sJDAzE39+fcePGERMTY3rudnm/v/jiC6ZMmYKDg0Ot463l/TYYDCxbtoySkhL69+9fZxlz/GzfSLt/zxw+2zfabnP7bN/M+93aP9tPPvkkY8eOveKzW5em+oy36I1AjUYjzz77LAMHDiQ8PPyq5TIzM/H29q51zNvbm8zMTNPzl45drUxLcqPt/r2//OUv+Pn51fpLMWrUKCZMmEBwcDDnzp3jpZdeYvTo0ezduxcLC4vGqP5Nu9F2h4SE8OWXX9K1a1cKCwt57733GDBgADExMbRt2/a2eL8PHDhAdHQ0X3zxRa3jreH9PnnyJP3796e8vBxHR0dWrVpFWFhYnWXN6bNdn3b/Xmv+bNen3eb02b7Z97s1f7YBli1bxpEjRzh48OANlW+qz3iLDjtPPvkk0dHR7N69u7mr0qRupt1vv/02y5YtY/v27bUG9E2ZMsX0/xEREXTt2pX27duzfft2RowY0aD1vlU32u7+/fvX+g1pwIABhIaG8sknn/Dmm282djUb3M2831988QURERH07du31vHW8H6HhIRw7NgxCgsL+f7775k5cyY7duy44S/+1upm293aP9v1abc5fbZv9v1uzZ/tlJQUnnnmGTZt2nTFwPLm1mJvYz311FOsXbuWbdu20bZt22uW9fHxISsrq9axrKwsfHx8TM9fOna1Mi1Ffdp9yXvvvcfbb7/Nxo0b6dq16zXLtmvXDg8PD+Lj4xuiug3mZtp9iZWVFT169DC1ydzf75KSEpYtW8bDDz983bIt8f22tramQ4cO9OrVi/nz59OtWzfef//9Osua02e7Pu2+xBw+2zfT7kta82f7Ztrd2j/bhw8fJjs7m549e2JpaYmlpSU7duzggw8+wNLSEoPBcMU5TfUZb3FhR1EUnnrqKVatWsXWrVsJDg6+7jn9+/dny5YttY5t2rTJ9BtCcHAwPj4+tcro9Xr2799/w/dQG9vNtBvgnXfe4c0332TDhg307t37uuVTU1PJy8vD19f3VqvcIG623ZczGAycPHnS1CZzfr9BnaZZUVHBtGnTrlu2pb3fdTEajVRUVNT5nDl8tq/mWu2G1v/ZvprrtftyrfGzfTU30u7W/tkeMWIEJ0+e5NixY6af3r1789BDD3Hs2LE6b7c12Wf8xsdXN40nnnhC0el0yvbt22tNsSstLTWVmT59uvLiiy+aHu/Zs0extLRU3nvvPSU2NlZ57bXX6py65uLioqxZs0Y5ceKEMm7cuBY1XfFm2v32228r1tbWyvfff1/rnKKiIkVR1BHxzz33nLJ3714lISFB2bx5s9KzZ0+lY8eOSnl5eZO3sS430+433nhD+fXXX5Vz584phw8fVqZMmaLY2toqMTExpjLm+H5fMmjQIGXy5MlXHG8N7/eLL76o7NixQ0lISFBOnDihvPjii4pGo1E2btyoKIp5frYVpf7tNofPtqLUv93m8NlWlPq3+5LW/Nm+mt/Pxmquz3iLCztAnT+XrzcwZMgQZebMmbXOW7FihdKpUyfF2tpa6dKli7Ju3bpazxuNRuWVV15RvL29FRsbG2XEiBFKXFxcE7ToxtxMuwMDA+s857XXXlMURVFKS0uVkSNHKp6enoqVlZUSGBioPProo1esW9Gcbqbdzz77rBIQEKBYW1sr3t7eypgxY5QjR47Uuq45vt+KoiinT59WANM/mpdrDe/3nDlzlMDAQMXa2lrx9PRURowYUast5vjZVpT6t9scPtuKUv92m8NnW1Fu7u95a/9sX83vw05zfcY1iqIoN94PJIQQQgjRurS4MTtCCCGEEA1Jwo4QQgghzJqEHSGEEEKYNQk7QgghhDBrEnaEEEIIYdYk7AghhBDCrEnYEUIIIYRZk7AjhBBCCLMmYUcI0WJpNBpWr17d3NUA4PXXX6d79+7NXQ0hxE2QsCOEEL/TkkKWEOLWSdgRQgghhFmTsCOEAGDt2rW4uLhgMBgAOHbsGBqNhhdffNFU5pFHHmHatGnk5eUxdepU2rRpg729PRERESxdutRU7tNPP8XPzw+j0VjrNcaNG8ecOXNMj9esWUPPnj2xtbWlXbt2vPHGG1RXV1+1jikpKUyaNAkXFxfc3NwYN24ciYmJpudnzZrF+PHjee+99/D19cXd3Z0nn3ySqqoqU5mMjAzGjh2LnZ0dwcHBLFmyhKCgIBYsWABAUFAQAPfddx8ajcb0+JJvvvmGoKAgdDodU6ZMoaio6Ib+fIUQzUfCjhACgMGDB1NUVMTRo0cB2LFjBx4eHmzfvt1UZseOHQwdOpTy8nJ69erFunXriI6O5rHHHmP69OkcOHAAgIkTJ5KXl8e2bdtM5164cIENGzbw0EMPAbBr1y5mzJjBM888w6lTp/jkk09YtGgRf//73+usX1VVFVFRUTg5ObFr1y727NmDo6Mjo0aNorKy0lRu27ZtnDt3jm3btrF48WIWLVrEokWLTM/PmDGD9PR0tm/fzg8//MCnn35Kdna26fmDBw8C8NVXX5GRkWF6DHDu3DlWr17N2rVrWbt2LTt27ODtt9++yT9xIUSTublN24UQ5qhnz57Ku+++qyiKoowfP175+9//rlhbWytFRUVKamqqAihnzpyp89yxY8cqf/7zn02Px40bp8yZM8f0+JNPPlH8/PwUg8GgKIqijBgxQvnHP/5R6xrffPON4uvra3oMKKtWrTI9FxISohiNRtPzFRUVip2dnfLrr78qiqIoM2fOVAIDA5Xq6mpTmYkTJyqTJ09WFEVRYmNjFUA5ePCg6fmzZ88qgPKf//ynzte95LXXXlPs7e0VvV5vOvb8888r/fr1q/PPQwjRckjPjhDCZMiQIWzfvh1FUdi1axcTJkwgNDSU3f/frv27pPbGcQB/e0uw8gQe+mGB1ZAohkGKEA16qMApLKfItrC2CIKGiMBaggjrH3AQDtTQlkNZUXlraEgcQhQriKKaXCSKMr/Dl3tAbOheiLre9wsE9Xn48DkO8uZ5Pj9/4vDwEM3NzTAajcjn81hcXITVaoUoitBqtdje3sb19bVSy+fzYXNzE8/PzwAAWZYxPDyMHz/+/9tJJBJYWFiAVqtVXn6/H3d3d3h8fCzpLZFIIJPJQBAEZb8oinh6esLFxYWyr6OjAxUVFcrnpqYm5eQmlUqhsrISNptNWW9vb4dOp/vQ79PW1gZBEN6tTUTfV+VXN0BE34ckSQiFQkgkElCr1TCbzZAkCQcHB8hms3C5XACA5eVlrK2tYXV1FVarFTU1NZiamiq6ThoYGEChUEAkEoHD4UAsFkMwGFTWc7kcAoEAvF5vSR8ajabku1wuB7vdDlmWS9bq6+uV92q1umhNpVKVzA79qc+sTUSfh2GHiBS/5naCwaASbCRJwtLSErLZLKanpwEAx8fH8Hg8GB0dBQC8vb0hnU7DYrEotTQaDbxeL2RZRiaTgclkKjpRsdlsSKVSaG9v/1BvNpsNGxsbaGhoQG1t7R89n8lkwuvrK+LxOOx2OwAgk8kgm80W7VOr1cqgNhH9/XiNRUQKnU6Hzs5OyLIMSZIAAE6nE2dnZ0in00oAMhqNiEajODk5QTKZxMTEBB4eHkrq+Xw+RCIRhEIhZTD5l/n5eYTDYQQCAZyfnyOZTGJ9fR1zc3Pv9ubz+VBXVwePx4NYLIarqyscHBxgcnISNzc3H3o+s9mM/v5+jI+P4/T0FPF4HOPj46iqqoJKpVL2tbW1YW9vD/f39yVBiIj+Pgw7RFTE5XIhn88rYUcURVgsFuj1ephMJgDA3NwcbDYb3G43JEmCXq/H4OBgSa3e3l6IoohUKoWRkZGiNbfbja2tLezs7MDhcKC7uxvBYBCtra3v9lVdXY2joyO0tLQos0RjY2N4enr6rZOecDiMxsZGOJ1ODA0Nwe/3QxCEoquzlZUVRKNRGAwGdHV1fbg2EX1PqkKhUPjqJoiIvsrNzQ0MBgN2d3fR19f31e0Q0Sdg2CGif8r+/j5yuRysVivu7u4wMzOD29tbpNPpkgFkIioPHFAmon/Ky8sLZmdncXl5CUEQ0NPTA1mWGXSIyhhPdoiIiKiscUCZiIiIyhrDDhEREZU1hh0iIiIqaww7REREVNYYdoiIiKisMewQERFRWWPYISIiorLGsENERERl7T+n7fVSaRelxQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "# MORB step-wise heating spectra\n", "heated_MORB_spectra = []\n", @@ -70,18 +81,20 @@ " MORB_spectrum.index = 1e4 / MORB_spectrum.index # 1/cm to microns\n", " MORB_spectrum = MORB_spectrum / 100 # percent to decimal\n", " MORB_spectrum = MORB_spectrum.sort_index()\n", - " if temperature == 650:\n", - " # Use the 650C spectrum for the portion of the spectrum below 2.6\n", - " # microns to isolate 3 micron features\n", - " norm_factor = spectra.get_normalization_factor(MORB_spectrum, 2.6)\n", - " MORB_spectrum = MORB_spectrum / norm_factor\n", - " MORB_spectrum.squeeze()\n", - " heated_MORB_spectra.append(MORB_spectrum)" + " \n", + " # Normalize all MORB spectra to the reflectance at 2.6 microns at the\n", + " # lowest water amount and highest temperature (22 ppm, 800 C)\n", + " norm_factor = spectra.get_normalization_factor(MORB_spectrum, 2.6)\n", + " normalized_spectrum = MORB_spectrum / norm_factor\n", + " normalized_spectrum.squeeze()\n", + " \n", + " normalized_spectrum.plot(xlim=(2,4))\n", + " heated_MORB_spectra.append(normalized_spectrum)\n" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -91,8 +104,8 @@ "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m/home/phil/repos/philiplinden/desci/cremons-etal-2022/repro/results.ipynb Cell 6\u001b[0m line \u001b[0;36m1\n\u001b[1;32m 4\u001b[0m interpolated_spectrum \u001b[39m=\u001b[39m simulator\u001b[39m.\u001b[39minterpolate_series(spectrum, WLS)\n\u001b[1;32m 6\u001b[0m \u001b[39m# Use Hapke model to convert from laboratory reflectance to single\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[39m# scattering albedo using reflectance and scattering asymmetry\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[39m# factor (p) of 0.81 (see manuscript for details on Hapke model)\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m ssa\u001b[39m.\u001b[39mappend(simulator\u001b[39m.\u001b[39;49mreflectance_to_ssa(\n\u001b[1;32m 11\u001b[0m reflectance\u001b[39m=\u001b[39;49minterpolated_spectrum, asymmetry_factor\u001b[39m=\u001b[39;49m\u001b[39m0.81\u001b[39;49m))\n", - "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/repro/simulator.py:174\u001b[0m, in \u001b[0;36mreflectance_to_ssa\u001b[0;34m(reflectance, asymmetry_factor, emission_angle, incident_angle, phase_angle, filling_factor)\u001b[0m\n\u001b[1;32m 167\u001b[0m y \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 168\u001b[0m bidirectional_reflectance, P\u001b[39m=\u001b[39masymmetry_factor, mu\u001b[39m=\u001b[39mmu, mu0\u001b[39m=\u001b[39mmu0, B\u001b[39m=\u001b[39mB\n\u001b[1;32m 169\u001b[0m )\n\u001b[1;32m 170\u001b[0m OLS \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 171\u001b[0m ordinary_least_squares, y\u001b[39m=\u001b[39my, yx\u001b[39m=\u001b[39mvalue\n\u001b[1;32m 172\u001b[0m ) \u001b[39m# turn least squares into the form y=f(x)\u001b[39;00m\n\u001b[1;32m 173\u001b[0m w\u001b[39m.\u001b[39mappend(\n\u001b[0;32m--> 174\u001b[0m fmin(\n\u001b[1;32m 175\u001b[0m OLS,\n\u001b[1;32m 176\u001b[0m w0,\n\u001b[1;32m 177\u001b[0m args\u001b[39m=\u001b[39;49mindex,\n\u001b[1;32m 178\u001b[0m disp\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m,\n\u001b[1;32m 179\u001b[0m maxiter\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 180\u001b[0m maxfun\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 181\u001b[0m ftol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 182\u001b[0m xtol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 183\u001b[0m )\n\u001b[1;32m 184\u001b[0m )\n\u001b[1;32m 185\u001b[0m \u001b[39mreturn\u001b[39;00m w\n", + "\u001b[1;32m/home/phil/repos/philiplinden/desci/cremons-etal-2022/repro/results.ipynb Cell 6\u001b[0m line \u001b[0;36m9\n\u001b[1;32m 4\u001b[0m interpolated_spectrum \u001b[39m=\u001b[39m simulator\u001b[39m.\u001b[39minterpolate_series(spectrum, WLS)\n\u001b[1;32m 6\u001b[0m \u001b[39m# Use Hapke model to convert from laboratory reflectance to single\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[39m# scattering albedo using reflectance and scattering asymmetry\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[39m# factor (p) of 0.81 (see manuscript for details on Hapke model)\u001b[39;00m\n\u001b[0;32m----> 9\u001b[0m ssa\u001b[39m.\u001b[39mappend(simulator\u001b[39m.\u001b[39;49mreflectance_to_ssa(\n\u001b[1;32m 10\u001b[0m reflectance\u001b[39m=\u001b[39;49minterpolated_spectrum, asymmetry_factor\u001b[39m=\u001b[39;49m\u001b[39m0.81\u001b[39;49m))\n", + "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/repro/simulator.py:179\u001b[0m, in \u001b[0;36mreflectance_to_ssa\u001b[0;34m(reflectance, asymmetry_factor, emission_angle, incident_angle, phase_angle, filling_factor)\u001b[0m\n\u001b[1;32m 173\u001b[0m \u001b[39m# formulate ordinary least squares between estimated reflectance, y\u001b[39;00m\n\u001b[1;32m 174\u001b[0m \u001b[39m# and observed reflectance, yx, and arrange into the form y=f(x)\u001b[39;00m\n\u001b[1;32m 175\u001b[0m OLS \u001b[39m=\u001b[39m partial(\n\u001b[1;32m 176\u001b[0m ordinary_least_squares, y\u001b[39m=\u001b[39my, yx\u001b[39m=\u001b[39mvalue\n\u001b[1;32m 177\u001b[0m )\n\u001b[1;32m 178\u001b[0m w\u001b[39m.\u001b[39mappend(\n\u001b[0;32m--> 179\u001b[0m fmin(\n\u001b[1;32m 180\u001b[0m OLS,\n\u001b[1;32m 181\u001b[0m w0,\n\u001b[1;32m 182\u001b[0m args\u001b[39m=\u001b[39;49mindex,\n\u001b[1;32m 183\u001b[0m disp\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m,\n\u001b[1;32m 184\u001b[0m maxiter\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 185\u001b[0m maxfun\u001b[39m=\u001b[39;49m\u001b[39m10_000\u001b[39;49m,\n\u001b[1;32m 186\u001b[0m ftol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 187\u001b[0m xtol\u001b[39m=\u001b[39;49m\u001b[39m1e-7\u001b[39;49m,\n\u001b[1;32m 188\u001b[0m )\n\u001b[1;32m 189\u001b[0m )\n\u001b[1;32m 190\u001b[0m \u001b[39mreturn\u001b[39;00m w\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:747\u001b[0m, in \u001b[0;36mfmin\u001b[0;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, initial_simplex)\u001b[0m\n\u001b[1;32m 738\u001b[0m opts \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mxatol\u001b[39m\u001b[39m'\u001b[39m: xtol,\n\u001b[1;32m 739\u001b[0m \u001b[39m'\u001b[39m\u001b[39mfatol\u001b[39m\u001b[39m'\u001b[39m: ftol,\n\u001b[1;32m 740\u001b[0m \u001b[39m'\u001b[39m\u001b[39mmaxiter\u001b[39m\u001b[39m'\u001b[39m: maxiter,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[39m'\u001b[39m\u001b[39mreturn_all\u001b[39m\u001b[39m'\u001b[39m: retall,\n\u001b[1;32m 744\u001b[0m \u001b[39m'\u001b[39m\u001b[39minitial_simplex\u001b[39m\u001b[39m'\u001b[39m: initial_simplex}\n\u001b[1;32m 746\u001b[0m callback \u001b[39m=\u001b[39m _wrap_callback(callback)\n\u001b[0;32m--> 747\u001b[0m res \u001b[39m=\u001b[39m _minimize_neldermead(func, x0, args, callback\u001b[39m=\u001b[39;49mcallback, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mopts)\n\u001b[1;32m 748\u001b[0m \u001b[39mif\u001b[39;00m full_output:\n\u001b[1;32m 749\u001b[0m retlist \u001b[39m=\u001b[39m res[\u001b[39m'\u001b[39m\u001b[39mx\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mfun\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnit\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mnfev\u001b[39m\u001b[39m'\u001b[39m], res[\u001b[39m'\u001b[39m\u001b[39mstatus\u001b[39m\u001b[39m'\u001b[39m]\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:899\u001b[0m, in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)\u001b[0m\n\u001b[1;32m 897\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m 898\u001b[0m \u001b[39mfor\u001b[39;00m k \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(N \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m):\n\u001b[0;32m--> 899\u001b[0m fsim[k] \u001b[39m=\u001b[39m func(sim[k])\n\u001b[1;32m 900\u001b[0m \u001b[39mexcept\u001b[39;00m _MaxFuncCallError:\n\u001b[1;32m 901\u001b[0m \u001b[39mpass\u001b[39;00m\n", "File \u001b[0;32m~/repos/philiplinden/desci/cremons-etal-2022/.venv/lib/python3.11/site-packages/scipy/optimize/_optimize.py:620\u001b[0m, in \u001b[0;36m_wrap_scalar_function_maxfun_validation..function_wrapper\u001b[0;34m(x, *wrapper_args)\u001b[0m\n\u001b[1;32m 618\u001b[0m ncalls[\u001b[39m0\u001b[39m] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m1\u001b[39m\n\u001b[1;32m 619\u001b[0m \u001b[39m# A copy of x is sent to the user function (gh13740)\u001b[39;00m\n\u001b[0;32m--> 620\u001b[0m fx \u001b[39m=\u001b[39m function(np\u001b[39m.\u001b[39mcopy(x), \u001b[39m*\u001b[39m(wrapper_args \u001b[39m+\u001b[39;49m args))\n\u001b[1;32m 621\u001b[0m \u001b[39m# Ideally, we'd like to a have a true scalar returned from f(x). For\u001b[39;00m\n\u001b[1;32m 622\u001b[0m \u001b[39m# backwards-compatibility, also allow np.array([1.3]),\u001b[39;00m\n\u001b[1;32m 623\u001b[0m \u001b[39m# np.array([[1.3]]) etc.\u001b[39;00m\n\u001b[1;32m 624\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m np\u001b[39m.\u001b[39misscalar(fx):\n", @@ -110,7 +123,7 @@ " # scattering albedo using reflectance and scattering asymmetry\n", " # factor (p) of 0.81 (see manuscript for details on Hapke model)\n", " ssa.append(simulator.reflectance_to_ssa(\n", - " reflectance=interpolated_spectrum, asymmetry_factor=0.81))" + " reflectance=interpolated_spectrum, asymmetry_factor=0.81))\n" ] } ],