diff --git a/bliss/cached_dataset.py b/bliss/cached_dataset.py index d9c58f6b9..87863d92f 100644 --- a/bliss/cached_dataset.py +++ b/bliss/cached_dataset.py @@ -69,7 +69,7 @@ def __call__(self, datum_in): class ChunkingSampler(Sampler): def __init__(self, dataset: Dataset) -> None: - super().__init__() + super().__init__(dataset) assert isinstance(dataset, ChunkingDataset), "dataset should be ChunkingDataset" self.dataset = dataset diff --git a/case_studies/galaxy_clustering/README.md b/case_studies/galaxy_clustering/README.md index 5df880db5..e39975475 100644 --- a/case_studies/galaxy_clustering/README.md +++ b/case_studies/galaxy_clustering/README.md @@ -10,28 +10,17 @@ Built on work done in Winter 2024 by [Li Shihang](https://www.linkedin.com/in/sh ## Generation of Data -Data can be generated by running the bash scipt `data-gen.sh`. The bash script has 3 options: -1. `-n`: the number of files to be generated (defaults to 100) -2. `-s`: the size of the image (defaults to 4800) -3. `-t`: the tile size (defaults to 4) - -As an example, the command - -``` -bash data-gen.sh -n 10 -s 2400 -t 8 -``` - -generates 10 images of size 2400 x 2400 tiled with a tile size of 8 x 8. - -The bash command `data-gen.sh` runs three scripts (located under the data_generation directory) for data generation: -1. `catalog_gen.py` which generates catalogs of images and stores them in the data/catalogs subdirectory. Keyword arguments: `image_size` and `nfiles`. -2. `galsim-des.yaml` then reads in these catalogs and uses GalSim to generate corresponding images, which are stored as .fits files (one for each band) in the data/images subdirectory. Keyword arguments: `image_size` and `nfiles`. -3. `file_datum_generation.py` reads in the catalogs and images and saves them as *FileDatum* objects which contain the tile catalog and images in a dictionary. Keyword arguments: `image_size` and `tile_size`. - -Often, after image data has been generated, we would want to retile it with a different tile size. This can be done by just running the file `file_datum_generation.py` with appropriate arguments (you would have to pass in the image size as well since it defaults to 4800). For example, if we have 80 x 80 images that we want to tile with tile size 8, we may run - -``` -python data_generation/file_datum_generation.py image_size=80 tile_size=8 -``` - -Note that you must run the file from the galaxy_clustering directory (since the script takes in the current working directory for the data paths). +The data generation routine proceeds through phases. The entire routine is conveniently wrapped into a single python script `data_gen.py` that draws its parameters from the Hydra configuration, located under `conf/config.yaml` under the `data_gen` key. These phases proceed as follows. + +1. **Catalog Generation.** First, we sample semi-synthetic source catalogs with their relevant properties, which are stored as `.dat` files in the `data_dir/catalogs` subdirectory. +2. **Image Generation.** Then, we take in the aforementioned source catalogs and use GalSim to render them as images, which are stored as `.fits` files (one for each band) in the `data_dir/images` subdirectory. +3. **File Datum Generation.** Finally, we convert the full source catalogs generated in phase 1 into tile catalogs, stack them up with their corresponding images, and store these objects as `.pt` files (which is what the encoder ultimately uses) in the `data_dir/file_data` subdirectory. + +The following parameters can be set within the configuration file `config.yaml`. +1. `data_dir`: the path of the directory where generated data will be stored. +2. `image_size`: size of the image (pixels). +3. `tile_size`: size of tile to be used (pixels). +4. `nfiles`: number of files to be generated. +5. `n_catalogs_per_file`: number of catalogs to be stored in each file datum object. +6. `bands`: survey bands to be used (`["g", "r", "i", "z"]` for DES). +7. `min_flux_for_loss`: minimum flux for filtering. diff --git a/case_studies/galaxy_clustering/config.yaml b/case_studies/galaxy_clustering/conf/config.yaml similarity index 75% rename from case_studies/galaxy_clustering/config.yaml rename to case_studies/galaxy_clustering/conf/config.yaml index 6b1e26260..83c0f0b6c 100644 --- a/case_studies/galaxy_clustering/config.yaml +++ b/case_studies/galaxy_clustering/conf/config.yaml @@ -1,9 +1,17 @@ --- defaults: - - ../../bliss/conf@_here_: base_config + - ../../../bliss/conf@_here_: base_config - _self_ - override hydra/job_logging: stdout +data_gen: + data_dir: /nfs/turbo/lsa-regier/scratch/kapnadak/new_data + image_size: 1280 + tile_size: 128 + nfiles: 5000 + n_catalogs_per_file: 500 + bands: ["g", "r", "i", "z"] + min_flux_for_loss: 0 prior: _target_: case_studies.galaxy_clustering.prior.GalaxyClusterPrior @@ -48,19 +56,15 @@ my_metrics: cluster_membership_acc: _target_: case_studies.galaxy_clustering.encoder.metrics.ClusterMembershipAccuracy +my_image_normalizers: + asinh: + _target_: bliss.encoder.image_normalizer.AsinhQuantileNormalizer + q: [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99, 0.999, 0.9999, 0.99999] + encoder: _target_: case_studies.galaxy_clustering.encoder.encoder.GalaxyClusterEncoder survey_bands: ["g", "r", "i", "z"] - image_normalizer: - _target_: bliss.encoder.image_normalizer.ImageNormalizer - bands: [0, 1, 2, 3] - include_original: true - include_background: false - concat_psf_params: false - num_psf_params: 6 # for SDSS, 4 for DC2 - log_transform_stdevs: null - use_clahe: false - clahe_min_stdev: null + image_normalizers: ${my_image_normalizers} mode_metrics: _target_: torchmetrics.MetricCollection _convert_: "partial" @@ -80,14 +84,21 @@ predict: _target_: case_studies.galaxy_clustering.cached_dataset.CachedDESModule cached_data_path: /nfs/turbo/lsa-regier/scratch/gapatron/desdr-server.ncsa.illinois.edu/despublic/dr2_tiles tiles_per_img: 64 - batch_size: 1 + batch_size: 2 num_workers: 4 trainer: _target_: pytorch_lightning.Trainer accelerator: "gpu" - devices: "6,7" + devices: [6,5] strategy: "ddp" precision: ${train.trainer.precision} + callbacks: + - ${predict.callbacks.writer} + callbacks: + writer: + _target_: case_studies.galaxy_clustering.inference.inference_callbacks.DESPredictionsWriter + output_dir: "/data/scratch/des/dr2_detection_output/run_1" + write_interval: "batch" encoder: ${encoder} weight_save_path: /nfs/turbo/lsa-regier/scratch/gapatron/best_encoder.ckpt device: "cuda:0" diff --git a/case_studies/galaxy_clustering/data_generation/DES_data_extraction.py b/case_studies/galaxy_clustering/data_generation/DES_data_extraction.py deleted file mode 100644 index a70415a0b..000000000 --- a/case_studies/galaxy_clustering/data_generation/DES_data_extraction.py +++ /dev/null @@ -1,48 +0,0 @@ -# flake8: noqa -import os - -import numpy as np -import requests -from astropy import units -from numpy.core.defchararray import startswith -from pyvo.dal import sia - -from case_studies.galaxy_clustering.utils import cluster_utils as utils - -DES_DATAPATH = os.environ["HOME"] + "/bliss/case_studies/galaxy_clustering/data/DES_images" -if not os.path.exists(DES_DATAPATH): - os.makedirs(DES_DATAPATH) - -DATASET = "des_dr1" -DEF_ACCESS_URL = "https://datalab.noirlab.edu/sia/" + DATASET -svc = sia.SIAService(DEF_ACCESS_URL) - - -def compute_fov(m500, z): - m200 = utils.m500_to_m200(m500, z) * units.solMass - r200 = utils.m200_to_r200(m200, z) - da = utils.angular_diameter_distance(z) - fov = (r200 / da) * (360 / (2 * np.pi)) - return fov.value - - -def download_image(m500, z, ra, dec, band, fov_scale=4): - fov = fov_scale * compute_fov(m500, z) - dec_radian = dec * np.pi / 180 - img_table = svc.search((ra, dec), (fov / np.cos(dec_radian), fov), verbosity=2).to_table() - sel = ( - (img_table["proctype"] == "Stack") - & (img_table["prodtype"] == "image") - & (startswith(img_table["obs_bandpass"].astype(str), band)) - ) - row = img_table[sel][0] - url = row["access_url"] # get the download URL - filename = DES_DATAPATH + "/" + str(ra) + "_" + str(dec) + "_" + band + ".fits" - if not url.lower().startswith("http"): - raise ValueError("URL must start with http") - response = requests.get(url, timeout=200) - if response.status_code == 200: - with open(filename, "wb") as file: - file.write(response.content) - else: - raise ValueError("Failed to download file.") diff --git a/case_studies/galaxy_clustering/data_generation/catalog_gen.py b/case_studies/galaxy_clustering/data_generation/catalog_gen.py deleted file mode 100644 index ea5107975..000000000 --- a/case_studies/galaxy_clustering/data_generation/catalog_gen.py +++ /dev/null @@ -1,41 +0,0 @@ -import os -import sys - -import numpy as np -import pandas as pd -from astropy.io import ascii as astro_ascii -from astropy.table import Table - -from case_studies.galaxy_clustering.data_generation.prior import BackgroundPrior, ClusterPrior - -DATA_PATH = "/home/kapnadak/bliss/case_studies/galaxy_clustering/data" -CATALOG_PATH = os.path.join(DATA_PATH, "catalogs") -FILE_PREFIX = "galsim_des" -CLUSTER_PROB = 0.5 - - -def main(**kwargs): - if not os.path.exists(CATALOG_PATH): - os.makedirs(CATALOG_PATH) - - nfiles = int(kwargs.get("nfiles", 100)) - cluster_prior = ClusterPrior(image_size=int(kwargs.get("image_size", 1280))) - background_prior = BackgroundPrior(image_size=int(kwargs.get("image_size", 1280))) - - combined_catalogs = [] - for _ in range(nfiles): - background_catalog = background_prior.sample_background() - if np.random.uniform() < CLUSTER_PROB: - cluster_catalog = cluster_prior.sample_cluster() - combined_catalogs.append(pd.concat([cluster_catalog, background_catalog])) - else: - combined_catalogs.append(background_catalog) - - for i, catalog in enumerate(combined_catalogs): - file_name = f"{CATALOG_PATH}/{FILE_PREFIX}_{i:03}.dat" - catalog_table = Table.from_pandas(catalog) - astro_ascii.write(catalog_table, file_name, format="no_header", overwrite=True) - - -if __name__ == "__main__": - main(**dict(arg.split("=") for arg in sys.argv[1:])) diff --git a/case_studies/galaxy_clustering/data_generation/data-gen.sh b/case_studies/galaxy_clustering/data_generation/data-gen.sh deleted file mode 100644 index e1ca01eba..000000000 --- a/case_studies/galaxy_clustering/data_generation/data-gen.sh +++ /dev/null @@ -1,67 +0,0 @@ -#!/bin/bash - -while getopts ":n:s:t:" opt; do - case $opt in - n) nfiles="$OPTARG" - ;; - s) image_size="$OPTARG" - ;; - t) tile_size="$OPTARG" - ;; - \?) echo "Invalid option -$OPTARG" >&2 - exit 1 - ;; - esac - - case $OPTARG in - -*) echo "Option $opt needs a valid argument" - exit 1 - ;; - esac -done - -echo "Generating Catalogs..." -if [ -z "$nfiles" ]; then - if [ -z "$image_size" ]; then - python3 catalog_gen.py - else - python3 catalog_gen.py image_size="$image_size" - fi -else - if [ -z "$image_size" ]; then - python3 catalog_gen.py nfiles="$nfiles" - else - python3 catalog_gen.py nfiles="$nfiles" image_size="$image_size" - fi -fi -echo "...Done!" -echo "Generating Images..." -if [ -z "$nfiles" ]; then - if [ -z "$image_size" ]; then - galsim galsim-des.yaml - else - galsim galsim-des.yaml variables.image_size="$image_size" - fi -else - if [ -z "$image_size" ]; then - galsim galsim-des.yaml variables.nfiles="$nfiles" - else - galsim galsim-des.yaml variables.nfiles="$nfiles" variables.image_size="$image_size" - fi -fi -echo "...Done!" -echo "Generating File Datums..." -if [ -z "$image_size" ]; then - if [ -z "$tile_size" ]; then - python3 file_datum_generation.py - else - python3 file_datum_generation.py tile_size="$tile_size" - fi -else - if [ -z "$tile_size" ]; then - python3 file_datum_generation.py image_size="$image_size" - else - python3 file_datum_generation.py image_size="$image_size" tile_size="$tile_size" - fi -fi -echo "...Done! Data Generation Complete!" diff --git a/case_studies/galaxy_clustering/data_generation/data_gen.py b/case_studies/galaxy_clustering/data_generation/data_gen.py new file mode 100644 index 000000000..b1c46c26f --- /dev/null +++ b/case_studies/galaxy_clustering/data_generation/data_gen.py @@ -0,0 +1,177 @@ +# flake8: noqa +import os +import subprocess +from pathlib import Path +from typing import Dict, List + +import hydra +import numpy as np +import pandas as pd +import torch +from astropy.io import ascii as astro_ascii +from astropy.io import fits +from astropy.table import Table + +from bliss.catalog import FullCatalog +from case_studies.galaxy_clustering.data_generation.prior import BackgroundPrior, ClusterPrior + +COL_NAMES = ( + "RA", + "DEC", + "X", + "Y", + "MEM", + "FLUX_G", + "FLUX_R", + "FLUX_I", + "FLUX_Z", + "HLR", + "FRACDEV", + "G1", + "G2", + "Z", + "SOURCE_TYPE", +) + +# ============================== Generate Catalogs ============================== + + +def catalog_gen(cfg): + nfiles = int(cfg.nfiles) + image_size = int(cfg.image_size) + data_dir = cfg.data_dir + catalogs_path = f"{data_dir}/catalogs/" + file_prefix = "galsim_des" + if not os.path.exists(catalogs_path): + os.makedirs(catalogs_path) + cluster_prior = ClusterPrior(image_size=image_size) + background_prior = BackgroundPrior(image_size=image_size) + + combined_catalogs = [] + for _ in range(nfiles): + background_catalog = background_prior.sample_background() + if np.random.uniform() < 0.5: + cluster_catalog = cluster_prior.sample_cluster() + combined_catalogs.append(pd.concat([cluster_catalog, background_catalog])) + else: + combined_catalogs.append(background_catalog) + + for i, catalog in enumerate(combined_catalogs): + file_name = f"{catalogs_path}/{file_prefix}_{i:03}.dat" + catalog_table = Table.from_pandas(catalog) + astro_ascii.write(catalog_table, file_name, format="no_header", overwrite=True) + + +# ============================== Generate Images ============================== + + +def image_gen(cfg): + image_size = cfg.image_size + nfiles = cfg.nfiles + data_dir = cfg.data_dir + input_dir = f"{data_dir}/catalogs" + output_dir = f"{data_dir}/images" + args = [] + args.append("galsim") + args.append("galsim-des.yaml") + args.append(f"variables.image_size={image_size}") + args.append(f"variables.nfiles={nfiles}") + args.append(f"variables.input_dir={input_dir}") + args.append(f"variables.output_dir={output_dir}") + subprocess.run(args, shell=False, check=False) + + +# ============================== Generate File Datums ============================== + + +def file_data_gen(cfg): + image_size = int(cfg.image_size) + tile_size = int(cfg.tile_size) + data_dir = cfg.data_dir + n_catalogs_per_file = int(cfg.n_catalogs_per_file) + bands = cfg.bands + min_flux_for_loss = float(cfg.min_flux_for_loss) + catalogs_path = Path(f"{data_dir}/catalogs/") + images_path = f"{data_dir}/images/" + file_path = f"{data_dir}/file_data/" + if not os.path.exists(file_path): + os.makedirs(file_path) + n_tiles = int(image_size / tile_size) + data: List[Dict] = [] + catalog_counter = 0 + file_counter = 0 + + for catalog_path in catalogs_path.glob("*.dat"): + catalog = pd.read_csv(catalog_path, sep=" ", header=None, names=COL_NAMES) + + catalog_dict = {} + catalog_dict["plocs"] = torch.tensor([catalog[["X", "Y"]].to_numpy()]) + catalog_dict["n_sources"] = torch.sum(catalog_dict["plocs"][:, :, 0] != 0, axis=1) + catalog_dict["galaxy_fluxes"] = torch.tensor( + [catalog[["FLUX_G", "FLUX_R", "FLUX_I", "FLUX_Z"]].to_numpy()] + ) + catalog_dict["star_fluxes"] = torch.zeros_like(catalog_dict["galaxy_fluxes"]) + catalog_dict["membership"] = torch.tensor([catalog[["MEM"]].to_numpy()]) + catalog_dict["redshift"] = torch.tensor([catalog[["Z"]].to_numpy()]) + catalog_dict["galaxy_params"] = torch.tensor( + [catalog[["HLR", "G1", "G2", "FRACDEV"]].to_numpy()] + ) + catalog_dict["source_type"] = torch.ones_like(catalog_dict["membership"]) + full_catalog = FullCatalog(height=image_size, width=image_size, d=catalog_dict) + tile_catalog = full_catalog.to_tile_catalog( + tile_slen=tile_size, + max_sources_per_tile=12 * tile_size, + ) + tile_catalog = tile_catalog.filter_by_flux(min_flux=min_flux_for_loss) + tile_catalog = tile_catalog.get_brightest_sources_per_tile(band=2, exclude_num=0) + + membership_array = np.zeros((n_tiles, n_tiles), dtype=bool) + for i, coords in enumerate(full_catalog["plocs"].squeeze()): + if full_catalog["membership"][0, i, 0] > 0: + tile_coord_y, tile_coord_x = ( + torch.div(coords, tile_size, rounding_mode="trunc").to(torch.int).tolist() + ) + membership_array[tile_coord_x, tile_coord_y] = True + + tile_catalog["membership"] = ( + torch.tensor(membership_array).unsqueeze(0).unsqueeze(3).unsqueeze(4) + ) + + tile_catalog_dict = {} + for key, value in tile_catalog.items(): + tile_catalog_dict[key] = torch.squeeze(value, 0) + + image_bands = [] + for band in bands: + fits_filepath = images_path / Path(f"{catalog_path.stem}_{band}.fits") + with fits.open(fits_filepath) as hdul: + image_data = hdul[0].data.astype(np.float32) + image_bands.append(torch.from_numpy(image_data)) + stacked_image = torch.stack(image_bands, dim=0) + + data.append( + { + "tile_catalog": tile_catalog_dict, + "images": stacked_image, + } + ) + catalog_counter += 1 + if catalog_counter == n_catalogs_per_file: + stackname = f"{file_path}/file_data_{file_counter}_size_{n_catalogs_per_file}.pt" + torch.save(data, stackname) + file_counter += 1 + data, catalog_counter = [], 0 + + +# ============================== CLI ============================== + + +@hydra.main(config_path="../conf", config_name="config", version_base=None) +def main(cfg): + catalog_gen(cfg.data_gen) + image_gen(cfg.data_gen) + file_data_gen(cfg.data_gen) + + +if __name__ == "__main__": + main() diff --git a/case_studies/galaxy_clustering/data_generation/dict_creation.py b/case_studies/galaxy_clustering/data_generation/dict_creation.py new file mode 100644 index 000000000..663f333ba --- /dev/null +++ b/case_studies/galaxy_clustering/data_generation/dict_creation.py @@ -0,0 +1,28 @@ +import os +import pickle +from pathlib import Path + +import pandas as pd +from astropy.io import fits + + +def main(): + des_dir = Path( + "/nfs/turbo/lsa-regier/scratch/gapatron/desdr-server.ncsa.illinois.edu/despublic/dr2_tiles/" + ) + des_subdirs = [d for d in os.listdir(des_dir) if d.startswith("DES")] + obj_tile_mapping = {} + output_filename = "/data/scratch/des/obj_to_tile.pickle" + for des_subdir in des_subdirs: + catalog_path = des_dir / Path(des_subdir) / Path(f"{des_subdir}_dr2_main.fits") + catalog_data = fits.getdata(catalog_path) + source_df = pd.DataFrame(catalog_data) + for obj_id in source_df["COADD_OBJECT_ID"]: + obj_tile_mapping[obj_id] = des_subdir + + with open(output_filename, "wb") as handle: + pickle.dump(obj_tile_mapping, handle, protocol=pickle.HIGHEST_PROTOCOL) + + +if __name__ == "__main__": + main() diff --git a/case_studies/galaxy_clustering/data_generation/file_datum_generation.py b/case_studies/galaxy_clustering/data_generation/file_datum_generation.py deleted file mode 100644 index 50e6e1471..000000000 --- a/case_studies/galaxy_clustering/data_generation/file_datum_generation.py +++ /dev/null @@ -1,110 +0,0 @@ -import os -import sys -from pathlib import Path - -import numpy as np -import pandas as pd -import torch -from astropy.io import fits - -from bliss.catalog import FullCatalog - -min_flux_for_loss = 0 -DATA_PATH = "/home/kapnadak/bliss/case_studies/galaxy_clustering/data" -CATALOGS_PATH = DATA_PATH / Path("catalogs") -IMAGES_PATH = DATA_PATH / Path("images") -FILE_DATA_PATH = DATA_PATH / Path("file_data") -if not os.path.exists(FILE_DATA_PATH): - os.makedirs(FILE_DATA_PATH) -COL_NAMES = ( - "RA", - "DEC", - "X", - "Y", - "MEM", - "FLUX_R", - "FLUX_G", - "FLUX_I", - "FLUX_Z", - "HLR", - "FRACDEV", - "G1", - "G2", - "Z", - "SOURCE_TYPE", -) -BANDS = ("g", "r", "i", "z") -N_CATALOGS_PER_FILE = 50 - - -def main(**kwargs): - image_size = int(kwargs.get("image_size", 1280)) - tile_size = int(kwargs.get("tile_size", 128)) - n_tiles = int(image_size / tile_size) - data = [] - - for catalog_path in CATALOGS_PATH.glob("*.dat"): - catalog = pd.read_csv(catalog_path, sep=" ", header=None, names=COL_NAMES) - - catalog_dict = {} - catalog_dict["plocs"] = torch.tensor([catalog[["X", "Y"]].to_numpy()]) - n_sources = torch.sum(catalog_dict["plocs"][:, :, 0] != 0, axis=1) - catalog_dict["n_sources"] = n_sources - catalog_dict["galaxy_fluxes"] = torch.tensor( - [catalog[["FLUX_R", "FLUX_G", "FLUX_I", "FLUX_Z"]].to_numpy()] - ) - catalog_dict["star_fluxes"] = torch.zeros_like(catalog_dict["galaxy_fluxes"]) - catalog_dict["membership"] = torch.tensor([catalog[["MEM"]].to_numpy()]) - catalog_dict["redshift"] = torch.tensor([catalog[["Z"]].to_numpy()]) - catalog_dict["galaxy_params"] = torch.tensor( - [catalog[["HLR", "G1", "G2", "FRACDEV"]].to_numpy()] - ) - catalog_dict["source_type"] = torch.ones_like(catalog_dict["membership"]) - full_catalog = FullCatalog(height=image_size, width=image_size, d=catalog_dict) - tile_catalog = full_catalog.to_tile_catalog( - tile_slen=tile_size, - max_sources_per_tile=12 * tile_size, - ) - tile_catalog = tile_catalog.filter_by_flux(min_flux=min_flux_for_loss) - tile_catalog = tile_catalog.get_brightest_sources_per_tile(band=2, exclude_num=0) - - membership_array = np.zeros((n_tiles, n_tiles), dtype=bool) - for i, coords in enumerate(full_catalog["plocs"].squeeze()): - if full_catalog["membership"][0, i, 0] > 0: - tile_coord_y, tile_coord_x = ( - torch.div(coords, tile_size, rounding_mode="trunc").to(torch.int).tolist() - ) - membership_array[tile_coord_x, tile_coord_y] = True - - tile_catalog["membership"] = ( - torch.tensor(membership_array).unsqueeze(0).unsqueeze(3).unsqueeze(4) - ) - - tile_catalog_dict = {} - for key, value in tile_catalog.items(): - tile_catalog_dict[key] = torch.squeeze(value, 0) - - filename = catalog_path.stem - image_bands = [] - for band in BANDS: - fits_filepath = IMAGES_PATH / Path(f"{filename}_{band}.fits") - # Should the ordering in the bands matter? It does here. - with fits.open(fits_filepath) as hdul: - image_data = hdul[0].data.astype(np.float32) - image_bands.append(torch.from_numpy(image_data)) - stacked_image = torch.stack(image_bands, dim=0) - - file_datum = { - "tile_catalog": tile_catalog_dict, - "images": stacked_image, - "background": stacked_image, - } - data.append(file_datum) - - chunks = [data[i : i + N_CATALOGS_PER_FILE] for i in range(0, len(data), N_CATALOGS_PER_FILE)] - for i, chunk in enumerate(chunks): - torch.save(chunk, f"{FILE_DATA_PATH}/file_data_{i}_size_{N_CATALOGS_PER_FILE}.pt") - - -if __name__ == "__main__": - main(**dict(arg.split("=") for arg in sys.argv[1:])) diff --git a/case_studies/galaxy_clustering/data_generation/galsim-des.yaml b/case_studies/galaxy_clustering/data_generation/galsim-des.yaml index 0d8a1b564..d189a4798 100644 --- a/case_studies/galaxy_clustering/data_generation/galsim-des.yaml +++ b/case_studies/galaxy_clustering/data_generation/galsim-des.yaml @@ -21,9 +21,9 @@ variables : nfiles : 100 - image_size: 4800 - input_dir : /home/kapnadak/bliss/case_studies/galaxy_clustering/data/catalogs - output_dir : /home/kapnadak/bliss/case_studies/galaxy_clustering/data/images + image_size: 1280 + input_dir : /nfs/turbo/lsa-regier/scratch/kapnadak/new_data_2/catalogs + output_dir : /nfs/turbo/lsa-regier/scratch/kapnadak/new_data_2/images psf: type: Convolve diff --git a/case_studies/galaxy_clustering/data_generation/prior.py b/case_studies/galaxy_clustering/data_generation/prior.py index 1f243c52b..6df7ba9bd 100644 --- a/case_studies/galaxy_clustering/data_generation/prior.py +++ b/case_studies/galaxy_clustering/data_generation/prior.py @@ -18,7 +18,7 @@ class Prior: - def __init__(self, image_size=4800): + def __init__(self, image_size=1280): super().__init__() self.width = image_size self.height = image_size @@ -96,15 +96,15 @@ def make_catalog( source_types, membership, ): - """Makes a background catalog from generated samples. + """Makes a catalog from generated samples. Args: flux_samples: flux samples in all bands hlr_samples: samples of HLR g1_size_samples: samples of G1 g2_size_samples: samples of G2 - gal_locs: samples of background locations in galactic coordinates - cartesian_locs: samples of background locations in cartesian coordinates + gal_locs: samples of locations in galactic coordinates + cartesian_locs: samples of locations in cartesian coordinates source_types: source types for each source (0 for star, 1 for galaxy) membership: background (0) or cluster (1) @@ -131,7 +131,7 @@ def make_catalog( class ClusterPrior(Prior): - def __init__(self, image_size=4800): + def __init__(self, image_size=1280): super().__init__(image_size) self.full_cluster_df = Table.read(CLUSTER_CATALOG_PATH).to_pandas() @@ -246,7 +246,7 @@ def sample_cluster(self): class BackgroundPrior(Prior): - def __init__(self, image_size=4800): + def __init__(self, image_size=1280): super().__init__(image_size) self.pixel_scale = 0.263 @@ -266,15 +266,9 @@ def sample_n_sources(self): def sample_des_catalog(self): """Sample a random DES dataframe.""" tile_choice = random.choice(DES_SUBDIRS) - main_path = DES_DIR / Path(tile_choice) / Path(f"{tile_choice}_dr2_main.fits") - flux_path = DES_DIR / Path(tile_choice) / Path(f"{tile_choice}_dr2_flux.fits") - main_data = fits.getdata(main_path) - main_df = pd.DataFrame(main_data) - flux_data = fits.getdata(flux_path) - flux_df = pd.DataFrame(flux_data) - self.source_df = pd.merge( - main_df, flux_df, left_on="COADD_OBJECT_ID", right_on="COADD_OBJECT_ID", how="left" - ) + catalog_path = DES_DIR / Path(tile_choice) / Path(f"{tile_choice}_dr2_main.fits") + catalog_data = fits.getdata(catalog_path) + self.source_df = pd.DataFrame(catalog_data) def sample_sources(self, n_sources): """Samples random sources from the current DES catalog. @@ -326,6 +320,23 @@ def sample_hlr(self, sources): hlr_samples = self.pixel_scale * np.array(sources["FLUX_RADIUS_R"]) return 1e-4 + (hlr_samples * (hlr_samples > 0)) + def sample_shapes(self, sources): + """Samples shapes for each source in the catalog. + Shapes are from DES Table in the form of (a, b) + Converted to (g1, g2) + + Args: + sources: Dataframe of DES sources + + Returns: + samples for g1, g2 for each source + """ + a = np.array(sources["A_IMAGE"]) + b = np.array(sources["B_IMAGE"]) + g = (a - b) / (a + b) + angle = np.arctan(b / a) + return g * np.cos(angle), g * np.sin(angle) + def sample_fluxes(self, sources): """Samples fluxes for all bands for each source. @@ -335,17 +346,18 @@ def sample_fluxes(self, sources): Returns: 5-band array containing fluxes (clamped at 1 from below) """ - fluxes = np.array( + mags = np.array( sources[ [ - "FLUX_AUTO_G_x", - "FLUX_AUTO_R_x", - "FLUX_AUTO_I_x", - "FLUX_AUTO_Z_x", + "MAG_AUTO_G", + "MAG_AUTO_R", + "MAG_AUTO_I", + "MAG_AUTO_Z", ] ] ) + fluxes = 1000 * convert_mag_to_nmgy(mags) return fluxes * (fluxes > 0) def sample_background(self): @@ -364,7 +376,7 @@ def sample_background(self): gal_source_locs = self.cartesian_to_gal(cartesian_source_locs) source_types = self.sample_source_types(des_sources) flux_samples = self.sample_fluxes(des_sources) - g1_size_samples, g2_size_samples = self.sample_shape(n_sources) + g1_size_samples, g2_size_samples = self.sample_shapes(n_sources) hlr_samples = self.sample_hlr(des_sources) return self.make_catalog( flux_samples, diff --git a/case_studies/galaxy_clustering/encoder/encoder.py b/case_studies/galaxy_clustering/encoder/encoder.py index cb755171d..06a39b95d 100644 --- a/case_studies/galaxy_clustering/encoder/encoder.py +++ b/case_studies/galaxy_clustering/encoder/encoder.py @@ -18,10 +18,10 @@ def initialize_networks(self): """ power_of_two = (self.tile_slen != 0) & (self.tile_slen & (self.tile_slen - 1) == 0) assert power_of_two, "tile_slen must be a power of two" - ch_per_band = self.image_normalizer.num_channels_per_band() + ch_per_band = sum(inorm.num_channels_per_band() for inorm in self.image_normalizers) num_features = 256 self.features_net = GalaxyClusterFeaturesNet( - len(self.image_normalizer.bands), + len(self.survey_bands), ch_per_band, num_features, tile_slen=self.tile_slen, @@ -39,7 +39,8 @@ def get_features_and_parameters(self, batch): batch_size, _n_bands, h, w = batch["images"].shape[0:4] ht, wt = h // self.tile_slen, w // self.tile_slen - x = self.image_normalizer.get_input_tensor(batch) + input_lst = [inorm.get_input_tensor(batch) for inorm in self.image_normalizers] + x = torch.cat(input_lst, dim=2) x_features = self.features_net(x) mask = torch.zeros([batch_size, ht, wt]) context = self.make_context(None, mask).to("cuda") @@ -74,8 +75,7 @@ def predict_step(self, batch, batch_idx, dataloader_idx=0): } def _compute_loss(self, batch, logging_name): - batch_size, _n_bands, h, w = batch["images"].shape[0:4] - ht, wt = h // self.tile_slen, w // self.tile_slen + batch_size = batch["images"].shape[0] target_cat = TileCatalog(self.tile_slen, batch["tile_catalog"]) @@ -91,11 +91,8 @@ def _compute_loss(self, batch, logging_name): # make predictions/inferences pred = {} - x = self.image_normalizer.get_input_tensor(batch) - x_features = self.features_net(x) - mask = torch.zeros([batch_size, ht, wt]) - context = self.make_context(None, mask).to("cuda") - pred["x_cat_marginal"] = self.catalog_net(x_features, context) + x_features, x_cat_marginal = self.get_features_and_parameters(batch) + pred["x_cat_marginal"] = x_cat_marginal x_features = x_features.detach() # is this helpful? doing it here to match old code loss = self.var_dist.compute_nll(pred["x_cat_marginal"], target_cat1) diff --git a/case_studies/galaxy_clustering/DES_inference.py b/case_studies/galaxy_clustering/inference/DES_inference.py similarity index 95% rename from case_studies/galaxy_clustering/DES_inference.py rename to case_studies/galaxy_clustering/inference/DES_inference.py index 5d313a167..e68e56fa7 100644 --- a/case_studies/galaxy_clustering/DES_inference.py +++ b/case_studies/galaxy_clustering/inference/DES_inference.py @@ -24,6 +24,7 @@ def inference(predict_cfg): encoder = load_encoder(predict_cfg) trainer = instantiate(predict_cfg.trainer) dataset = instantiate(predict_cfg.cached_dataset) + enc_output = trainer.predict(encoder, datamodule=dataset) gpu_rank = ( distributed.get_rank() if distributed.is_available() and distributed.is_initialized() else 0 @@ -33,7 +34,7 @@ def inference(predict_cfg): def main(): - with initialize(config_path=".", version_base=None): + with initialize(config_path="../conf", version_base=None): cfg = compose("config") predict_cfg = cfg.predict diff --git a/case_studies/galaxy_clustering/cached_dataset.py b/case_studies/galaxy_clustering/inference/cached_dataset.py similarity index 99% rename from case_studies/galaxy_clustering/cached_dataset.py rename to case_studies/galaxy_clustering/inference/cached_dataset.py index 6d37df5bc..cf43758aa 100644 --- a/case_studies/galaxy_clustering/cached_dataset.py +++ b/case_studies/galaxy_clustering/inference/cached_dataset.py @@ -130,7 +130,7 @@ def _get_dataloader(self, dataset): seed=42, ) else: - sampler = DESSampler() + sampler = DESSampler(dataset) return DataLoader( dataset, diff --git a/case_studies/galaxy_clustering/inference/inference_callbacks.py b/case_studies/galaxy_clustering/inference/inference_callbacks.py new file mode 100644 index 000000000..cd885c89e --- /dev/null +++ b/case_studies/galaxy_clustering/inference/inference_callbacks.py @@ -0,0 +1,30 @@ +import os + +import torch +from pytorch_lightning.callbacks import BasePredictionWriter + + +class DESPredictionsWriter(BasePredictionWriter): + def __init__(self, output_dir, write_interval): + super().__init__(write_interval) + self.output_dir = output_dir + + def write_on_batch_end( + self, trainer, pl_module, prediction, batch_indices, batch, batch_idx, dataloader_idx + ): + # this will create N (num processes) files in `output_dir` each containing + # the predictions of it's respective rank + name = f"rank_{trainer.global_rank}_batchIdx_{batch_idx}_dataloaderIdx_{dataloader_idx}.pt" + torch.save( + prediction, + os.path.join( + self.output_dir, + name, + ), + ) + + # optionally, you can also save `batch_indices` to get the information about the data index + # from your prediction data + torch.save( + batch_indices, os.path.join(self.output_dir, f"batch_indices_{trainer.global_rank}.pt") + ) diff --git a/case_studies/galaxy_clustering/notebooks/DESInference.ipynb b/case_studies/galaxy_clustering/notebooks/DESInference.ipynb new file mode 100644 index 000000000..19049bafe --- /dev/null +++ b/case_studies/galaxy_clustering/notebooks/DESInference.ipynb @@ -0,0 +1,255 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "import numpy as np\n", + "import os\n", + "import random" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "DES_DIR = \"/nfs/turbo/lsa-regier/scratch/gapatron/desdr-server.ncsa.illinois.edu/despublic/dr2_tiles\"\n", + "DES_BANDS = (\"g\", \"r\", \"i\", \"z\")\n", + "DES_SUBDIRS = [d for d in os.listdir(DES_DIR) if d.startswith(\"DES\")]\n", + "tiles_per_img = 64" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def convert_to_global_idx(tile_idx, gpu_idx):\n", + " num_gpus = 2\n", + " tiles_per_img = 64\n", + " batch_size = 2\n", + " dir_idx = int(num_gpus * (tile_idx // (tiles_per_img / batch_size)) + gpu_idx)\n", + " subimage_idx = [(batch_size * tile_idx + i) % tiles_per_img for i in range(batch_size)]\n", + " return dir_idx, subimage_idx\n", + "\n", + "def convert_to_tile_idx(dir_idx):\n", + " num_gpus = 2\n", + " tiles_per_img = 64\n", + " batch_size = 2\n", + " gpu_idx = dir_idx % num_gpus\n", + " tile_starting_idx = (tiles_per_img / batch_size) * (dir_idx // num_gpus)\n", + " return int(tile_starting_idx), int(gpu_idx)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Directory: DES0053-2041\n", + "Subimage: [0, 1]\n" + ] + } + ], + "source": [ + "tile_idx = 0\n", + "gpu_idx = 0\n", + "\n", + "dir_idx, subimage_idx = convert_to_global_idx(tile_idx, gpu_idx)\n", + "print(f\"Directory: {DES_SUBDIRS[dir_idx]}\")\n", + "print(f\"Subimage: {subimage_idx}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "torch.Size([64, 1280, 1280])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "memberships = torch.empty((0,10,10))\n", + "output_dir = \"/data/scratch/des/dr2_detection_output/run_0\"\n", + "dir_idx = 10167\n", + "tile_starting_idx, gpu_idx = convert_to_tile_idx(dir_idx)\n", + "for tile in range(tile_starting_idx, tile_starting_idx + 32):\n", + " file = torch.load(f\"{output_dir}/tile_{tile}_gpu_{gpu_idx}.pt\")\n", + " memberships = torch.cat((memberships, file[\"mode_cat\"][\"membership\"].squeeze()), dim=0)\n", + "\n", + "expanded_memberships = torch.repeat_interleave(memberships, repeats=128, dim=1)\n", + "expanded_memberships = torch.repeat_interleave(expanded_memberships, repeats=128, dim=2)\n", + "expanded_memberships.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "500" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# gpu 0 : 0, 2, 4, 6, 8, ...\n", + "# gpu 1 : 1, 3, 5, 7, 9, ...\n", + "\n", + "# tile_0 : 0, 1 -- image 0\n", + "# tile_1 : 2, 3 -- image 0\n", + "# tile_2 : 4, 5 -- image 0\n", + "# ...\n", + "# tile_31 : 62, 63 -- image 0\n", + "# tile_32 : 64, 65 -- image 1 (image 2 overall)\n", + "\n", + "# tile_t_gpu_g --> dir_id = 2 * (t // 32) + g\n", + "# --> sub_id = (2*t % 64), (2*t + 1) % 64" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def count_num_clusters(dir_idx):\n", + " memberships = torch.empty((0,10,10))\n", + " tile_starting_idx, gpu_idx = convert_to_tile_idx(dir_idx)\n", + " for tile in range(tile_starting_idx, tile_starting_idx + 32):\n", + " file = torch.load(f\"{output_dir}/tile_{tile}_gpu_{gpu_idx}.pt\")\n", + " memberships = torch.cat((memberships, file[\"mode_cat\"][\"membership\"].squeeze()), dim=0)\n", + " memberships = torch.repeat_interleave(memberships, repeats=128, dim=1)\n", + " memberships = torch.repeat_interleave(memberships, repeats=128, dim=2)\n", + " return torch.any(memberships.view(memberships.shape[0], -1), dim=1).sum()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "random_sample = random.sample(list(enumerate(DES_SUBDIRS)), 100)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "num_clusters = {}\n", + "for dir_idx , dir in random_sample:\n", + " num_clusters[dir] = count_num_clusters(dir_idx)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "num_clusters_vals = list(num_clusters.values())" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAK9CAYAAADG5r/mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjkElEQVR4nO3deVwV1eP/8fcVEUQENJXFENx30Vxwyz2RzLTMlI8mmll91NRcKiv3StPcKtNWMbNMrdSvW+5aSrlFLqmpoWjuG7jkBuf3Rz/uxysXFETB4fV8POahc2bmzJnD3OHN3HPn2owxRgAAAIBF5crqBgAAAAB3E4EXAAAAlkbgBQAAgKUReAEAAGBpBF4AAABYGoEXAAAAlkbgBQAAgKUReAEAAGBpBF4AAABYGoEXSIfg4GB16dIlq5theWPHjlWJEiXk4uKiqlWrZnVz0rRmzRrZbDatWbMmq5tyX4uKipLNZtOBAwfu+b6PHz+up556Sg888IBsNpsmTpx4z9uQ03Tp0kXBwcEOZTabTcOGDcuS9sD6CLzIsZJ/wW7evNnp8kaNGqlSpUp3vJ/FixdzEU+HZcuW6ZVXXlG9evU0bdo0vfPOO1ndJFjcyy+/rB9//FGDBg3SjBkz1KJFi1TX/fbbb9WpUyeVLl1aNptNjRo1SnXdK1eu6NVXX1VAQIDy5s2r0NBQLV++3Om6GzZsUP369eXh4SE/Pz/17t1bFy5cuGXbDxw4IJvNZp9cXV1VqFAh1a1bV6+//rri4uJSbJP8R1pq06xZs+zrXr16VZMmTVK1atXk5eUlHx8fVaxYUc8//7x2796dZtuOHDmiYcOGKSYm5pbHAdxtubO6AcD9ZM+ePcqVK31/Jy5evFiTJ08m9N6mVatWKVeuXPr888+VJ0+erG4OcoBVq1apdevWGjBgwC3XnTJlirZs2aKaNWvq9OnTaa7bpUsXzZ07V3379lXp0qUVFRWlRx99VKtXr1b9+vXt68XExKhp06YqX768xo8fr8OHD+u9997T3r17tWTJkts6hoiICD366KNKSkrS2bNntWnTJk2cOFGTJk3S559/rg4dOqTYpnfv3qpZs2aK8jp16tj/37ZtWy1ZskQRERHq3r27rl27pt27d2vhwoWqW7euypUrl2qbjhw5ouHDhys4ODjFOzWffvqpkpKSbuvYgMxA4AXSwc3NLaubkG4XL15Uvnz5sroZt+3EiRPKmzdvhsPupUuX5OHhkcmtgpWdOHFCPj4+t7XujBkzVLRoUeXKlSvNd4A2btyoWbNmaezYsfYg3blzZ1WqVEmvvPKKNmzYYF/39ddfV4ECBbRmzRp5eXlJ+nf4VPfu3bVs2TI1b978lu166KGH1KlTJ4eygwcPqnnz5oqMjFT58uUVEhLisPzhhx/WU089lWqdmzZt0sKFC/X222/r9ddfd1j24Ycf6ty5c7dsV2pcXV0zvC2QEQxpANLh5jG8165d0/Dhw1W6dGm5u7vrgQceUP369e1vW3bp0kWTJ0+WJIe3DJNdvHhR/fv3V2BgoNzc3FS2bFm99957MsY47Peff/5R7969VahQIeXPn1+PP/64/v777xRj3oYNGyabzaY//vhD//nPf1SgQAH7naRt27apS5cuKlGihNzd3eXn56dnn302xV2q5Dr+/PNPderUSd7e3ipcuLAGDx4sY4wOHTqk1q1by8vLS35+fho3btxt9d3169c1cuRIlSxZUm5ubgoODtbrr7+uK1eu2Nex2WyaNm2aLl68aO+rqKioVOtMHnayZcsWNWjQQB4eHvZfzPPnz1fLli0VEBAgNzc3lSxZUiNHjlRiYqLTOv744w81btxYHh4eKlq0qMaMGZNif4cPH1abNm2UL18+FSlSRC+//LJD+280Z84cVa9eXXnz5lWhQoXUqVMn/f333w7rdOnSRZ6enoqLi9Njjz0mT09PFS1a1H7ObN++XU2aNFG+fPkUFBSkr7/++rb6etasWapevbry588vLy8vVa5cWZMmTbIvP3PmjAYMGKDKlSvL09NTXl5eCg8P1++//+5QT/Jb37Nnz9bw4cNVtGhR5c+fX0899ZTi4+N15coV9e3bV0WKFJGnp6e6du2aoj9sNpt69eqlmTNnqmzZsnJ3d1f16tW1bt262zqWJUuW6OGHH1a+fPmUP39+tWzZUjt37rytbf/66y+1a9dOBQsWlIeHh2rXrq1FixbZlycPazLGaPLkySlen84EBgbe1rs8c+fOlYuLi55//nl7mbu7u7p166bo6GgdOnRIkpSQkKDly5erU6dO9rAr/RuOPT09NXv27Ns6VmeCgoIUFRWlq1evOj2fb2X//v2SpHr16qVY5uLiogceeCDVbdesWWO/e9y1a9cUr2dnY3id+fvvv/Xss8/K19dXbm5uqlixor744osU633wwQeqWLGiPDw8VKBAAdWoUeO2Xy/IGbjDixwvPj5ep06dSlF+7dq1W247bNgwjRo1Ss8995xq1aqlhIQEbd68WVu3btUjjzyiF154QUeOHNHy5cs1Y8YMh22NMXr88ce1evVqdevWTVWrVtWPP/6ogQMH6u+//9aECRPs63bp0kWzZ8/WM888o9q1a2vt2rVq2bJlqu1q166dSpcurXfeeccenpcvX66//vpLXbt2lZ+fn3bu3KlPPvlEO3fu1C+//JLiF3379u1Vvnx5jR49WosWLdJbb72lggUL6uOPP1aTJk307rvvaubMmRowYIBq1qypBg0apNlXzz33nKZPn66nnnpK/fv316+//qpRo0Zp165d+uGHHyT9e/fsk08+0caNG/XZZ59JkurWrZtmvadPn1Z4eLg6dOigTp06ydfXV9K/YcbT01P9+vWTp6enVq1apSFDhighIUFjx451qOPs2bNq0aKFnnzyST399NOaO3euXn31VVWuXFnh4eGS/v2jo2nTpoqLi1Pv3r0VEBCgGTNmaNWqVSnaFBUVpa5du6pmzZoaNWqUjh8/rkmTJmn9+vX67bffHO4mJiYmKjw8XA0aNNCYMWM0c+ZM9erVS/ny5dMbb7yhjh076sknn9TUqVPVuXNn1alTR8WLF0+1P5YvX66IiAg1bdpU7777riRp165dWr9+vfr06SPp3yA4b948tWvXTsWLF9fx48f18ccfq2HDhvrjjz8UEBDgUOeoUaOUN29evfbaa9q3b58++OADubq6KleuXDp79qyGDRumX375RVFRUSpevLiGDBnisP3atWv17bffqnfv3nJzc9NHH32kFi1aaOPGjWneJZ0xY4YiIyMVFhamd999V5cuXdKUKVNUv359/fbbb2kGpuPHj6tu3bq6dOmSevfurQceeEDTp0/X448/rrlz5+qJJ55QgwYNNGPGDD3zzDN65JFH1Llz51TrS6/ffvtNZcqUcQixklSrVi1J/w5jCAwM1Pbt23X9+nXVqFHDYb08efKoatWq+u233+6oHXXq1FHJkiWdjh0+f/6802tf8of3goKCJEkzZ85UvXr1lDv37UeG8uXLa8SIERoyZIief/55Pfzww5Ju/Xq+0fHjx1W7dm37H02FCxfWkiVL1K1bNyUkJKhv376S/h0e0bt3bz311FPq06ePLl++rG3btunXX3/Vf/7zn9veHyzOADnUtGnTjKQ0p4oVKzpsExQUZCIjI+3zISEhpmXLlmnup2fPnsbZS23evHlGknnrrbccyp966iljs9nMvn37jDHGbNmyxUgyffv2dVivS5cuRpIZOnSovWzo0KFGkomIiEixv0uXLqUo++abb4wks27duhR1PP/88/ay69evmwcffNDYbDYzevRoe/nZs2dN3rx5HfrEmZiYGCPJPPfccw7lAwYMMJLMqlWr7GWRkZEmX758adaXrGHDhkaSmTp1aoplzo73hRdeMB4eHuby5csp6vjyyy/tZVeuXDF+fn6mbdu29rKJEycaSWb27Nn2sosXL5pSpUoZSWb16tXGGGOuXr1qihQpYipVqmT++ecf+7oLFy40ksyQIUMcjlWSeeedd+xlyX1qs9nMrFmz7OW7d+9O8fN2pk+fPsbLy8tcv3491XUuX75sEhMTHcpiY2ONm5ubGTFihL1s9erVRpKpVKmSuXr1qr08IiLC2Gw2Ex4e7lBHnTp1TFBQkENZ8mtp8+bN9rKDBw8ad3d388QTT9jLkl+PsbGxxhhjzp8/b3x8fEz37t0d6jt27Jjx9vZOUX6zvn37Gknmp59+spedP3/eFC9e3AQHBzscvyTTs2fPNOtzpmLFiqZhw4apLmvSpEmK8p07dzqcs3PmzEnxGkzWrl074+fnl2YbYmNjjSQzduzYVNdp3bq1kWTi4+ONMf/7uaY2HT161BhjTFJSkv314evrayIiIszkyZPNwYMH02xTsk2bNhlJZtq0aSmWRUZGOj1Xbjy/u3XrZvz9/c2pU6cc1uvQoYPx9va2v8Zbt26d4loN3IwhDcjxJk+erOXLl6eYqlSpcsttfXx8tHPnTu3duzfd+128eLFcXFzUu3dvh/L+/fvLGGP/sMrSpUslST169HBY76WXXkq17hdffDFFWd68ee3/v3z5sk6dOqXatWtLkrZu3Zpi/eeee87+fxcXF9WoUUPGGHXr1s1e7uPjo7Jly+qvv/5KtS3Sv8cqSf369XMo79+/vyQ5vM2cXm5uburatWuK8huPN/lO1sMPP6xLly6l+HS5p6enw/jHPHnyqFatWg7HtXjxYvn7+zuMefTw8HB4y1qSNm/erBMnTqhHjx5yd3e3l7ds2VLlypVzeqw39nVyn+bLl09PP/20vbxs2bLy8fG5ZV/7+Pjo4sWLqT4NQPq3z5Lflk9MTNTp06fl6empsmXLOj0XOnfu7DDmMjQ0VMYYPfvssw7rhYaG6tChQ7p+/bpDeZ06dVS9enX7fLFixdS6dWv9+OOPKYaYJFu+fLnOnTuniIgInTp1yj65uLgoNDRUq1evTrMfFi9erFq1ajl8OMzT01PPP/+8Dhw4oD/++CPN7e/UP//843TMf/I58c8//zj8m9q6ycvvhKenp6R/Xwc3GjJkiNNrX8GCBSX9Oxzlxx9/1FtvvaUCBQrom2++Uc+ePRUUFKT27dvf0RjeWzHG6LvvvlOrVq1kjHE4B8LCwhQfH28/V318fHT48GFt2rTprrUH9z+GNCDHq1WrVoq3EyWpQIECTt/uu9GIESPUunVrlSlTRpUqVVKLFi30zDPP3FZYPnjwoAICApQ/f36H8vLly9uXJ/+bK1euFG9jlypVKtW6nb3lfebMGQ0fPlyzZs3SiRMnHJbFx8enWL9YsWIO897e3nJ3d1ehQoVSlN/q0+rJx3Bzm/38/OTj42M/1owoWrSo0w+47dy5U2+++aZWrVqlhIQEh2U3H++DDz6YYkhHgQIFtG3bNodjKFWqVIr1ypYt6zCffCw3l0tSuXLl9PPPPzuUubu7q3Dhwg5l3t7eTtvk7e2ts2fPpqj3Rj169NDs2bMVHh6uokWLqnnz5nr66acdHrWVlJSkSZMm6aOPPlJsbKxD6HQ2LtPZuSD9O5715vKkpCTFx8c71FO6dOkUdZYpU0aXLl3SyZMn5efnl2J58h+RTZo0cXqcNw8VuNnBgwcVGhqaovzG11dmPHYwNXnz5nU6vvvy5cv25Tf+m9q6N/7hllHJjze7+VpTuXJlNWvWLM1t3dzc9MYbb+iNN97Q0aNHtXbtWk2aNEmzZ8+Wq6urvvrqqztunzMnT57UuXPn9Mknn+iTTz5xuk7ydezVV1/VihUrVKtWLZUqVUrNmzfXf/7zH6djj5FzEXiBO9CgQQPt379f8+fP17Jly/TZZ59pwoQJmjp1qsNdu3vN2S/Jp59+Whs2bNDAgQNVtWpVeXp6KikpSS1atHD6eCAXF5fbKpOU4kN2qbnVB4Iywtmxnjt3Tg0bNpSXl5dGjBihkiVLyt3dXVu3btWrr76a4njv9LjuRGr7zmibihQpopiYGP34449asmSJlixZomnTpqlz586aPn26JOmdd97R4MGD9eyzz2rkyJEqWLCgcuXKpb59+972uXAnbbwdye2YMWOG00CcnvGkWcHf3z/FhxQl6ejRo5JkHyft7+/vUH7zujePp86IHTt2qEiRIrf8I+FW/P391aFDB7Vt21YVK1bU7NmzFRUVdVd+Fsk//06dOikyMtLpOsk3FsqXL689e/Zo4cKFWrp0qb777jt99NFHGjJkiIYPH57pbcP9KXtfMYD7QMGCBdW1a1d17dpVFy5cUIMGDTRs2DB74E0t5AUFBWnFihU6f/68w52X5Lfbkz8wEhQUpKSkJMXGxjrcKdu3b99tt/Hs2bNauXKlhg8f7vCBoowMxciI5GPYu3ev/Q6b9O+HUs6dO2c/1syyZs0anT59Wt9//73Dh+liY2MzXGdQUJB27NghY4zDz3TPnj0p1ksuv/nu5J49ezL9WJ3JkyePWrVqpVatWikpKUk9evTQxx9/rMGDB6tUqVKaO3euGjdurM8//9xhu3PnzqW4g58ZnJ1nf/75pzw8PFLc3U5WsmRJSf8G+FvdhXQmKCgoxc9GSvn6uluqVq2q1atXKyEhwSFo/vrrr/blklSpUiXlzp1bmzdvdhjCcvXqVcXExDiUZUR0dLT279+f4pFld8LV1VVVqlTR3r17derUKad/kEh39gdu4cKFlT9/fiUmJt7Wzz9fvnxq37692rdvr6tXr+rJJ5/U22+/rUGDBjkMLULOxRhe4A7c/Fa+p6enSpUq5fD2ZPIzcG8e7/boo48qMTFRH374oUP5hAkTZLPZ7E8HCAsLkyR99NFHDut98MEHt93O5DtxN995u1dfofroo4863d/48eMlKc0nTmSEs+O9evVqij5Mj0cffVRHjhzR3Llz7WWXLl1K8XZrjRo1VKRIEU2dOtXhPFiyZIl27dqV6cd6s5vPyVy5ctnvhCW3x8XFJcW5MGfOHKd3JDNDdHS0w9jgQ4cOaf78+WrevHmqd4nDwsLk5eWld955x+kTU06ePJnmPh999FFt3LhR0dHR9rKLFy/qk08+UXBwsCpUqJDBo7k9Tz31lBITEx3OjytXrmjatGkKDQ21Dwfx9vZWs2bN9NVXXzmMsZ0xY4YuXLigdu3aZbgNBw8eVJcuXZQnTx4NHDgw3dvv3bvX6Te1nTt3TtHR0SpQoECqf7BIqV/7boeLi4vatm2r7777Tjt27Eix/Maf/83nfJ48eVShQgUZY27raTvIGbjDC9yBChUqqFGjRqpevboKFiyozZs3a+7cuerVq5d9neQP6/Tu3VthYWFycXFRhw4d1KpVKzVu3FhvvPGGDhw4oJCQEC1btkzz589X37597Xe4qlevrrZt22rixIk6ffq0/bFkf/75p6Tbu4vi5eVlf+zVtWvXVLRoUS1btuyO7nimR0hIiCIjI/XJJ5/Yhxts3LhR06dPV5s2bdS4ceNM3V/dunVVoEABRUZGqnfv3rLZbJoxY8YdvdXevXt3ffjhh+rcubO2bNkif39/zZgxI8WXXLi6uurdd99V165d1bBhQ0VERNgfSxYcHKyXX375Tg8vTc8995zOnDmjJk2a6MEHH9TBgwf1wQcfqGrVqva764899phGjBihrl27qm7dutq+fbtmzpypEiVK3JU2VapUSWFhYQ6PJZOU5tvNXl5emjJlip555hk99NBD6tChgwoXLqy4uDgtWrRI9erVS/HH4o1ee+01ffPNNwoPD1fv3r1VsGBBTZ8+XbGxsfruu+/S/Y2JydatW2d/hvDJkyd18eJFvfXWW5L+HeKU/I5CaGio2rVrp0GDBunEiRMqVaqUpk+frgMHDqS4s/7222+rbt26atiwoZ5//nkdPnxY48aNU/PmzdP8muMbbd26VV999ZWSkpJ07tw5bdq0Sd9995393Hf2uYKffvrJPqb4RlWqVFGVKlX0+++/6z//+Y/Cw8P18MMPq2DBgvr77781ffp0HTlyRBMnTkz1Dxbp37v0Pj4+mjp1qvLnz698+fIpNDQ0zcfq3Wj06NFavXq1QkND1b17d1WoUEFnzpzR1q1btWLFCp05c0aS1Lx5c/n5+alevXry9fXVrl279OGHH6ply5Ypxi0jB8uCJ0MA2ULyY5A2bdrkdHnDhg1v+Viyt956y9SqVcv4+PiYvHnzmnLlypm3337b4RFO169fNy+99JIpXLiwsdlsDo8oO3/+vHn55ZdNQECAcXV1NaVLlzZjx441SUlJDvu9ePGi6dmzpylYsKDx9PQ0bdq0MXv27DGSHB4TlvxIsZMnT6Y4nsOHD5snnnjC+Pj4GG9vb9OuXTtz5MiRVB9tdnMdqT0uzFk/OXPt2jUzfPhwU7x4cePq6moCAwPNoEGDHB4RltZ+nElr3+vXrze1a9c2efPmNQEBAeaVV14xP/74o8MjxNKqw9ljkw4ePGgef/xx4+HhYQoVKmT69Oljli5dmqJOY4z59ttvTbVq1Yybm5spWLCg6dixozl8+PBtHWtqbQoKCrrlY/Dmzp1rmjdvbooUKWLy5MljihUrZl544QX7o6aM+fexZP379zf+/v4mb968pl69eiY6Oto0bNjQ4TFbyY+vmjNnjsM+UnvtODt39P8f+fXVV1+Z0qVLGzc3N1OtWrUU/XXzY8lubENYWJjx9vY27u7upmTJkqZLly4OjzlLzf79+81TTz1lfHx8jLu7u6lVq5ZZuHBhivWS23g7ko/R2XTzI+P++ecfM2DAAOPn52fc3NxMzZo1zdKlS53W+9NPP5m6desad3d3U7hwYdOzZ0+TkJBwy/YkP5YsecqdO7cpWLCgCQ0NNYMGDXL6CLFbPZYs+TiOHz9uRo8ebRo2bGj8/f1N7ty5TYECBUyTJk3M3Llzb6u/5s+fbypUqGBy587t8Iiy23ksWXIbevbsaQIDA42rq6vx8/MzTZs2NZ988ol9nY8//tg0aNDAPPDAA8bNzc2ULFnSDBw40P4YNsAYY2zG3INPZQDIdDExMapWrZq++uordezYMaubAzhls9nUs2fPNO/GAsDdxhhe4D7g7FmcEydOVK5cuW75DWcAAOR0jOEF7gNjxozRli1b1LhxY+XOndv+uKnnn38+xbNQAQCAIwIvcB+oW7euli9frpEjR+rChQsqVqyYhg0bpjfeeCOrmwYAQLbHGF4AAABYGmN4AQAAYGkEXgAAAFgaY3idSEpK0pEjR5Q/f/47+mpEAAAA3B3GGJ0/f14BAQG3/DIZAq8TR44c4ZPvAAAA94FDhw7pwQcfTHMdAq8TyV9FeOjQIXl5eWVxawAAAHCzhIQEBQYG3tZXSBN4nUgexuDl5UXgBQAAyMZuZ/gpH1oDAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGm5s7oBAABkhuDXFmV1E+66A6NbZnUTgPsSd3gBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaVkaeNetW6dWrVopICBANptN8+bNc1hus9mcTmPHjk21zmHDhqVYv1y5cnf5SAAAAJBdZWngvXjxokJCQjR58mSny48ePeowffHFF7LZbGrbtm2a9VasWNFhu59//vluNB8AAAD3gdxZufPw8HCFh4enutzPz89hfv78+WrcuLFKlCiRZr25c+dOsS0AAABypvtmDO/x48e1aNEidevW7Zbr7t27VwEBASpRooQ6duyouLi4NNe/cuWKEhISHCYAAABYw30TeKdPn678+fPrySefTHO90NBQRUVFaenSpZoyZYpiY2P18MMP6/z586luM2rUKHl7e9unwMDAzG4+AAAAssh9E3i/+OILdezYUe7u7mmuFx4ernbt2qlKlSoKCwvT4sWLde7cOc2ePTvVbQYNGqT4+Hj7dOjQocxuPgAAALJIlo7hvV0//fST9uzZo2+//Tbd2/r4+KhMmTLat29fquu4ubnJzc3tTpoIAACAbOq+uMP7+eefq3r16goJCUn3thcuXND+/fvl7+9/F1oGAACA7C5LA++FCxcUExOjmJgYSVJsbKxiYmIcPmSWkJCgOXPm6LnnnnNaR9OmTfXhhx/a5wcMGKC1a9fqwIED2rBhg5544gm5uLgoIiLirh4LAAAAsqcsHdKwefNmNW7c2D7fr18/SVJkZKSioqIkSbNmzZIxJtXAun//fp06dco+f/jwYUVEROj06dMqXLiw6tevr19++UWFCxe+ewcCAACAbMtmjDFZ3YjsJiEhQd7e3oqPj5eXl1dWNwcAcBuCX1uU1U246w6MbpnVTQCyjfTktftiDC8AAACQUQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJaWpYF33bp1atWqlQICAmSz2TRv3jyH5V26dJHNZnOYWrRocct6J0+erODgYLm7uys0NFQbN268S0cAAACA7C5LA+/FixcVEhKiyZMnp7pOixYtdPToUfv0zTffpFnnt99+q379+mno0KHaunWrQkJCFBYWphMnTmR28wEAAHAfyJ2VOw8PD1d4eHia67i5ucnPz++26xw/fry6d++url27SpKmTp2qRYsW6YsvvtBrr712R+0FAADA/Sfbj+Fds2aNihQporJly+q///2vTp8+neq6V69e1ZYtW9SsWTN7Wa5cudSsWTNFR0enut2VK1eUkJDgMAEAAMAasvQO7620aNFCTz75pIoXL679+/fr9ddfV3h4uKKjo+Xi4pJi/VOnTikxMVG+vr4O5b6+vtq9e3eq+xk1apSGDx+e6e0HrCT4tUVZ3YS77sDollndBADAXZCtA2+HDh3s/69cubKqVKmikiVLas2aNWratGmm7WfQoEHq16+ffT4hIUGBgYGZVj8AAACyTrYf0nCjEiVKqFChQtq3b5/T5YUKFZKLi4uOHz/uUH78+PE0xwG7ubnJy8vLYQIAAIA13FeB9/Dhwzp9+rT8/f2dLs+TJ4+qV6+ulStX2suSkpK0cuVK1alT5141EwAAANlIlgbeCxcuKCYmRjExMZKk2NhYxcTEKC4uThcuXNDAgQP1yy+/6MCBA1q5cqVat26tUqVKKSwszF5H06ZN9eGHH9rn+/Xrp08//VTTp0/Xrl279N///lcXL160P7UBAAAAOUuWjuHdvHmzGjdubJ9PHkcbGRmpKVOmaNu2bZo+fbrOnTungIAANW/eXCNHjpSbm5t9m/379+vUqVP2+fbt2+vkyZMaMmSIjh07pqpVq2rp0qUpPsgGAACAnCFLA2+jRo1kjEl1+Y8//njLOg4cOJCirFevXurVq9edNA0AAAAWcV+N4QUAAADSi8ALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAAS8ud1Q0AANwbwa8tyuom4A5Z/Wd4YHTLrG4CLIo7vAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAAS8vSwLtu3Tq1atVKAQEBstlsmjdvnn3ZtWvX9Oqrr6py5crKly+fAgIC1LlzZx05ciTNOocNGyabzeYwlStX7i4fCQAAALKrLA28Fy9eVEhIiCZPnpxi2aVLl7R161YNHjxYW7du1ffff689e/bo8ccfv2W9FStW1NGjR+3Tzz//fDeaDwAAgPtA7qzceXh4uMLDw50u8/b21vLlyx3KPvzwQ9WqVUtxcXEqVqxYqvXmzp1bfn5+mdpWAAAA3J/uqzG88fHxstls8vHxSXO9vXv3KiAgQCVKlFDHjh0VFxeX5vpXrlxRQkKCwwQAAABruG8C7+XLl/Xqq68qIiJCXl5eqa4XGhqqqKgoLV26VFOmTFFsbKwefvhhnT9/PtVtRo0aJW9vb/sUGBh4Nw4BAAAAWeC+CLzXrl3T008/LWOMpkyZkua64eHhateunapUqaKwsDAtXrxY586d0+zZs1PdZtCgQYqPj7dPhw4dyuxDAAAAQBbJ0jG8tyM57B48eFCrVq1K8+6uMz4+PipTpoz27duX6jpubm5yc3O706YCAAAgG8rWd3iTw+7evXu1YsUKPfDAA+mu48KFC9q/f7/8/f3vQgsBAACQ3WVp4L1w4YJiYmIUExMjSYqNjVVMTIzi4uJ07do1PfXUU9q8ebNmzpypxMREHTt2TMeOHdPVq1ftdTRt2lQffvihfX7AgAFau3atDhw4oA0bNuiJJ56Qi4uLIiIi7vXhAQAAIBvI0iENmzdvVuPGje3z/fr1kyRFRkZq2LBhWrBggSSpatWqDtutXr1ajRo1kiTt379fp06dsi87fPiwIiIidPr0aRUuXFj169fXL7/8osKFC9/dgwEAAEC2lKWBt1GjRjLGpLo8rWXJDhw44DA/a9asO20WAAAALCRbj+EFAAAA7hSBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJaWocC7detWbd++3T4/f/58tWnTRq+//rquXr2aaY0DAAAA7lSGAu8LL7ygP//8U5L0119/qUOHDvLw8NCcOXP0yiuvZGoDAQAAgDuRocD7559/qmrVqpKkOXPmqEGDBvr6668VFRWl7777LjPbBwAAANyRDAVeY4ySkpIkSStWrNCjjz4qSQoMDNSpU6cyr3UAAADAHcpQ4K1Ro4beeustzZgxQ2vXrlXLli0lSbGxsfL19c3UBgIAAAB3IkOBd8KECdq6dat69eqlN954Q6VKlZIkzZ07V3Xr1s3UBgIAAAB3IndGNgoJCXF4SkOysWPHKnfuDFUJAAAA3BUZusNbokQJnT59OkX55cuXVaZMmTtuFAAAAJBZMhR4Dxw4oMTExBTlV65c0eHDh++4UQAAAEBmSdf4gwULFtj//+OPP8rb29s+n5iYqJUrV6p48eKZ1zoAAADgDqUr8LZp00aSZLPZFBkZ6bDM1dVVwcHBGjduXKY1DgAAALhT6Qq8yc/eLV68uDZt2qRChQrdlUYBAAAAmSVDj1SIjY3N7HYAAAAAd0WGnyG2cuVKrVy5UidOnLDf+U32xRdf3HHDAAAAgMyQocA7fPhwjRgxQjVq1JC/v79sNltmtwsAAADIFBkKvFOnTlVUVJSeeeaZzG4PAAAAkKky9Bzeq1evZspXCK9bt06tWrVSQECAbDab5s2b57DcGKMhQ4bI399fefPmVbNmzbR3795b1jt58mQFBwfL3d1doaGh2rhx4x23FQAAAPenDAXe5557Tl9//fUd7/zixYsKCQnR5MmTnS4fM2aM3n//fU2dOlW//vqr8uXLp7CwMF2+fDnVOr/99lv169dPQ4cO1datWxUSEqKwsDCdOHHijtsLAACA+0+GhjRcvnxZn3zyiVasWKEqVarI1dXVYfn48eNvq57w8HCFh4c7XWaM0cSJE/Xmm2+qdevWkqQvv/xSvr6+mjdvnjp06OB0u/Hjx6t79+7q2rWrpH+HXyxatEhffPGFXnvttds9RAAAAFhEhgLvtm3bVLVqVUnSjh07HJZl1gfYYmNjdezYMTVr1sxe5u3trdDQUEVHRzsNvFevXtWWLVs0aNAge1muXLnUrFkzRUdHp7qvK1eu6MqVK/b5hISETDkGAAAAZL0MBd7Vq1dndjtSOHbsmCTJ19fXodzX19e+7GanTp1SYmKi0212796d6r5GjRql4cOH32GLAQAAUhf82qKsbsJdd2B0y6xuglMZGsNrNYMGDVJ8fLx9OnToUFY3CQAAAJkkQ3d4GzdunObQhVWrVmW4Qcn8/PwkScePH5e/v7+9/Pjx4/bhFDcrVKiQXFxcdPz4cYfy48eP2+tzxs3NTW5ubnfcZgAAAGQ/GbrDW7VqVYWEhNinChUq6OrVq9q6dasqV66cKQ0rXry4/Pz8tHLlSntZQkKCfv31V9WpU8fpNnny5FH16tUdtklKStLKlStT3QYAAADWlqE7vBMmTHBaPmzYMF24cOG267lw4YL27dtnn4+NjVVMTIwKFiyoYsWKqW/fvnrrrbdUunRpFS9eXIMHD1ZAQIDatGlj36Zp06Z64okn1KtXL0lSv379FBkZqRo1aqhWrVqaOHGiLl68aH9qAwAAAHKWDAXe1HTq1Em1atXSe++9d1vrb968WY0bN7bP9+vXT5IUGRmpqKgovfLKK7p48aKef/55nTt3TvXr19fSpUvl7u5u32b//v06deqUfb59+/Y6efKkhgwZomPHjqlq1apaunRpig+yAQAAIGfI1MAbHR3tEEZvpVGjRjLGpLrcZrNpxIgRGjFiRKrrHDhwIEVZr1697Hd8AQAAkLNlKPA++eSTDvPGGB09elSbN2/W4MGDM6VhAAAAQGbIUOD19vZ2mM+VK5fKli2rESNGqHnz5pnSMAAAACAzZCjwTps2LbPbAQAAANwVdzSGd8uWLdq1a5ckqWLFiqpWrVqmNAoAAADILBkKvCdOnFCHDh20Zs0a+fj4SJLOnTunxo0ba9asWSpcuHBmthEAAADIsAx98cRLL72k8+fPa+fOnTpz5ozOnDmjHTt2KCEhQb17987sNgIAAAAZlqE7vEuXLtWKFStUvnx5e1mFChU0efJkPrQGAACAbCVDd3iTkpLk6uqaotzV1VVJSUl33CgAAAAgs2Qo8DZp0kR9+vTRkSNH7GV///23Xn75ZTVt2jTTGgcAAADcqQwF3g8//FAJCQkKDg5WyZIlVbJkSRUvXlwJCQn64IMPMruNAAAAQIZlaAxvYGCgtm7dqhUrVmj37t2SpPLly6tZs2aZ2jgAAADgTqXrDu+qVatUoUIFJSQkyGaz6ZFHHtFLL72kl156STVr1lTFihX1008/3a22AgAAAOmWrsA7ceJEde/eXV5eXimWeXt764UXXtD48eMzrXEAAADAnUpX4P3999/VokWLVJc3b95cW7ZsueNGAQAAAJklXYH3+PHjTh9Hlix37tw6efLkHTcKAAAAyCzpCrxFixbVjh07Ul2+bds2+fv733GjAAAAgMySrsD76KOPavDgwbp8+XKKZf/884+GDh2qxx57LNMaBwAAANypdD2W7M0339T333+vMmXKqFevXipbtqwkaffu3Zo8ebISExP1xhtv3JWGAgAAABmRrsDr6+urDRs26L///a8GDRokY4wkyWazKSwsTJMnT5avr+9daSgAAACQEen+4omgoCAtXrxYZ8+e1b59+2SMUenSpVWgQIG70T4AAADgjmTom9YkqUCBAqpZs2ZmtgUAAADIdOn60BoAAABwvyHwAgAAwNIIvAAAALC0DI/hBQCrCX5tUVY3AQBwF3CHFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAICl5c7qBgBWEfzaoqxuAgDc17iO4m7hDi8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAAS8v2gTc4OFg2my3F1LNnT6frR0VFpVjX3d39HrcaAAAA2UXurG7ArWzatEmJiYn2+R07duiRRx5Ru3btUt3Gy8tLe/bssc/bbLa72kYAAABkX9k+8BYuXNhhfvTo0SpZsqQaNmyY6jY2m01+fn53u2kAAAC4D2T7IQ03unr1qr766is9++yzad61vXDhgoKCghQYGKjWrVtr586dadZ75coVJSQkOEwAAACwhvsq8M6bN0/nzp1Tly5dUl2nbNmy+uKLLzR//nx99dVXSkpKUt26dXX48OFUtxk1apS8vb3tU2Bg4F1oPQAAALKCzRhjsroRtyssLEx58uTR//3f/932NteuXVP58uUVERGhkSNHOl3nypUrunLlin0+ISFBgYGBio+Pl5eX1x23GzlD8GuLsroJAABkqQOjW96zfSUkJMjb2/u28lq2H8Ob7ODBg1qxYoW+//77dG3n6uqqatWqad++famu4+bmJjc3tzttIgAAALKh+2ZIw7Rp01SkSBG1bJm+vxwSExO1fft2+fv736WWAQAAIDu7LwJvUlKSpk2bpsjISOXO7XhTunPnzho0aJB9fsSIEVq2bJn++usvbd26VZ06ddLBgwf13HPP3etmAwAAIBu4L4Y0rFixQnFxcXr22WdTLIuLi1OuXP/L7WfPnlX37t117NgxFShQQNWrV9eGDRtUoUKFe9lkAAAAZBP31YfW7pX0DIIGkvGhNQBATpddP7R2XwxpAAAAADKKwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLy9aBd9iwYbLZbA5TuXLl0txmzpw5KleunNzd3VW5cmUtXrz4HrUWAAAA2VG2DrySVLFiRR09etQ+/fzzz6muu2HDBkVERKhbt2767bff1KZNG7Vp00Y7duy4hy0GAABAdpLtA2/u3Lnl5+dnnwoVKpTqupMmTVKLFi00cOBAlS9fXiNHjtRDDz2kDz/88B62GAAAANlJtg+8e/fuVUBAgEqUKKGOHTsqLi4u1XWjo6PVrFkzh7KwsDBFR0enuY8rV64oISHBYQIAAIA1ZOvAGxoaqqioKC1dulRTpkxRbGysHn74YZ0/f97p+seOHZOvr69Dma+vr44dO5bmfkaNGiVvb2/7FBgYmGnHAAAAgKyVrQNveHi42rVrpypVqigsLEyLFy/WuXPnNHv27Ezdz6BBgxQfH2+fDh06lKn1AwAAIOvkzuoGpIePj4/KlCmjffv2OV3u5+en48ePO5QdP35cfn5+adbr5uYmNze3TGsnAAAAso9sfYf3ZhcuXND+/fvl7+/vdHmdOnW0cuVKh7Lly5erTp0696J5AAAAyIaydeAdMGCA1q5dqwMHDmjDhg164okn5OLiooiICElS586dNWjQIPv6ffr00dKlSzVu3Djt3r1bw4YN0+bNm9WrV6+sOgQAAABksWw9pOHw4cOKiIjQ6dOnVbhwYdWvX1+//PKLChcuLEmKi4tTrlz/y+x169bV119/rTfffFOvv/66SpcurXnz5qlSpUpZdQgAAADIYjZjjMnqRmQ3CQkJ8vb2Vnx8vLy8vLK6ObhPBL+2KKubAABAljowuuU921d68lq2HtIAAAAA3CkCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLy53VDcC/gl9blNVNuOsOjG6Z1U0AAAA5EHd4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWRuAFAACApRF4AQAAYGkEXgAAAFgagRcAAACWlq0D76hRo1SzZk3lz59fRYoUUZs2bbRnz540t4mKipLNZnOY3N3d71GLAQAAkN1k68C7du1a9ezZU7/88ouWL1+ua9euqXnz5rp48WKa23l5eeno0aP26eDBg/eoxQAAAMhucmd1A9KydOlSh/moqCgVKVJEW7ZsUYMGDVLdzmazyc/P7243DwAAAPeBbH2H92bx8fGSpIIFC6a53oULFxQUFKTAwEC1bt1aO3fuTHP9K1euKCEhwWECAACANdw3gTcpKUl9+/ZVvXr1VKlSpVTXK1u2rL744gvNnz9fX331lZKSklS3bl0dPnw41W1GjRolb29v+xQYGHg3DgEAAABZ4L4JvD179tSOHTs0a9asNNerU6eOOnfurKpVq6phw4b6/vvvVbhwYX388cepbjNo0CDFx8fbp0OHDmV28wEAAJBFsvUY3mS9evXSwoULtW7dOj344IPp2tbV1VXVqlXTvn37Ul3Hzc1Nbm5ud9pMAAAAZEPZ+g6vMUa9evXSDz/8oFWrVql48eLpriMxMVHbt2+Xv7//XWghAAAAsrtsfYe3Z8+e+vrrrzV//nzlz59fx44dkyR5e3srb968kqTOnTuraNGiGjVqlCRpxIgRql27tkqVKqVz585p7NixOnjwoJ577rksOw4AAABknWwdeKdMmSJJatSokUP5tGnT1KVLF0lSXFyccuX6343qs2fPqnv37jp27JgKFCig6tWra8OGDapQocK9ajYAAACykWwdeI0xt1xnzZo1DvMTJkzQhAkT7lKLAAAAcL/J1mN4AQAAgDtF4AUAAIClEXgBAABgaQReAAAAWBqBFwAAAJZG4AUAAIClEXgBAABgaQReAAAAWFq2/uIJWEvwa4uyugkAACAH4g4vAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALI3ACwAAAEsj8AIAAMDSCLwAAACwNAIvAAAALO2+CLyTJ09WcHCw3N3dFRoaqo0bN6a5/pw5c1SuXDm5u7urcuXKWrx48T1qKQAAALKbbB94v/32W/Xr109Dhw7V1q1bFRISorCwMJ04ccLp+hs2bFBERIS6deum3377TW3atFGbNm20Y8eOe9xyAAAAZAc2Y4zJ6kakJTQ0VDVr1tSHH34oSUpKSlJgYKBeeuklvfbaaynWb9++vS5evKiFCxfay2rXrq2qVatq6tSpt7XPhIQEeXt7Kz4+Xl5eXplzILcQ/Nqie7IfAACAu+XA6Jb3bF/pyWu571GbMuTq1avasmWLBg0aZC/LlSuXmjVrpujoaKfbREdHq1+/fg5lYWFhmjdvXqr7uXLliq5cuWKfj4+Pl/RvR94rSVcu3bN9AQAA3A33Mjsl7+t27t1m68B76tQpJSYmytfX16Hc19dXu3fvdrrNsWPHnK5/7NixVPczatQoDR8+PEV5YGBgBloNAACQM3lPvPf7PH/+vLy9vdNcJ1sH3ntl0KBBDneFk5KSdObMGT3wwAOy2Wx3ff8JCQkKDAzUoUOH7tkQivsFfeMc/ZI6+sY5+sU5+iV19I1z9Evq7nXfGGN0/vx5BQQE3HLdbB14CxUqJBcXFx0/ftyh/Pjx4/Lz83O6jZ+fX7rWlyQ3Nze5ubk5lPn4+GSs0XfAy8uLF08q6Bvn6JfU0TfO0S/O0S+po2+co19Sdy/75lZ3dpNl66c05MmTR9WrV9fKlSvtZUlJSVq5cqXq1KnjdJs6deo4rC9Jy5cvT3V9AAAAWFu2vsMrSf369VNkZKRq1KihWrVqaeLEibp48aK6du0qSercubOKFi2qUaNGSZL69Omjhg0baty4cWrZsqVmzZqlzZs365NPPsnKwwAAAEAWyfaBt3379jp58qSGDBmiY8eOqWrVqlq6dKn9g2lxcXHKlet/N6rr1q2rr7/+Wm+++aZef/11lS5dWvPmzVOlSpWy6hBuyc3NTUOHDk0xrAL0TWrol9TRN87RL87RL6mjb5yjX1KXnfsm2z+HFwAAALgT2XoMLwAAAHCnCLwAAACwNAIvAAAALI3ACwAAAEsj8GayKVOmqEqVKvaHLtepU0dLliyxL798+bJ69uypBx54QJ6enmrbtm2KL8q4mTFGQ4YMkb+/v/LmzatmzZpp7969d/tQMlVa/XLmzBm99NJLKlu2rPLmzatixYqpd+/eio+PT7POLl26yGazOUwtWrS4F4eTqW51zjRq1CjFcb744otp1mn1c+bAgQMp+iR5mjNnTqp1WuWcudHo0aNls9nUt29fe1lOvc7c7Oa+yenXmmTOzpmcep252c19k5OvNcOGDUtxDOXKlbMvv++uMwaZasGCBWbRokXmzz//NHv27DGvv/66cXV1NTt27DDGGPPiiy+awMBAs3LlSrN582ZTu3ZtU7du3TTrHD16tPH29jbz5s0zv//+u3n88cdN8eLFzT///HMvDilTpNUv27dvN08++aRZsGCB2bdvn1m5cqUpXbq0adu2bZp1RkZGmhYtWpijR4/apzNnztyjI8o8tzpnGjZsaLp37+5wnPHx8WnWafVz5vr16w79cfToUTN8+HDj6elpzp8/n2qdVjlnkm3cuNEEBwebKlWqmD59+tjLc+p15kbO+ianX2uMSf2cyanXmRs565ucfK0ZOnSoqVixosMxnDx50r78frvOEHjvgQIFCpjPPvvMnDt3zri6upo5c+bYl+3atctIMtHR0U63TUpKMn5+fmbs2LH2snPnzhk3NzfzzTff3PW2303J/eLM7NmzTZ48ecy1a9dS3T4yMtK0bt36LrUua93YNw0bNnT4xXQrOfWcqVq1qnn22WfT3N5K58z58+dN6dKlzfLlyx3OEa4zqfeNMznpWpNWv+T060x6zpmccq0ZOnSoCQkJcbrsfrzOMKThLkpMTNSsWbN08eJF1alTR1u2bNG1a9fUrFkz+zrlypVTsWLFFB0d7bSO2NhYHTt2zGEbb29vhYaGprpNdndzvzgTHx8vLy8v5c6d9nejrFmzRkWKFFHZsmX13//+V6dPn74bTb5nUuubmTNnqlChQqpUqZIGDRqkS5cupVpHTjxntmzZopiYGHXr1u2WdVnlnOnZs6datmzp8HOWxHVGqfeNMznpWnOrfsnJ15nbPWdy2rVm7969CggIUIkSJdSxY0fFxcVJuj+vM9n+m9buR9u3b1edOnV0+fJleXp66ocfflCFChUUExOjPHnyyMfHx2F9X19fHTt2zGldyeXJ3yx3O9tkV6n1y81OnTqlkSNH6vnnn0+zvhYtWujJJ59U8eLFtX//fr3++usKDw9XdHS0XFxc7tZh3BVp9c1//vMfBQUFKSAgQNu2bdOrr76qPXv26Pvvv3daV048Zz7//HOVL19edevWTbM+q5wzs2bN0tatW7Vp06YUy44dO5ajrzNp9c3NctK15lb9kpOvM+k5Z3LStSY0NFRRUVEqW7asjh49quHDh+vhhx/Wjh077svrDIH3LihbtqxiYmIUHx+vuXPnKjIyUmvXrs3qZmW51PrlxgCTkJCgli1bqkKFCho2bFia9XXo0MH+/8qVK6tKlSoqWbKk1qxZo6ZNm96tw7gr0uqbG38ZV65cWf7+/mratKn279+vkiVLZmGr777bOWf++ecfff311xo8ePAt67PCOXPo0CH16dNHy5cvl7u7e1Y3J1tJT9/kpGvN7fRLTr3OpOecyWnXmvDwcPv/q1SpotDQUAUFBWn27NnKmzdvFrYsYxjScBfkyZNHpUqVUvXq1TVq1CiFhIRo0qRJ8vPz09WrV3Xu3DmH9Y8fPy4/Pz+ndSWX3/zJx7S2ya5S65dk58+fV4sWLZQ/f3798MMPcnV1TVf9JUqUUKFChbRv377Mbvpdd6u+uVFoaKgkpXqcOemckaS5c+fq0qVL6ty5c7rrvx/PmS1btujEiRN66KGHlDt3buXOnVtr167V+++/r9y5c8vX1zfHXmdu1TeJiYmSct615nb75UY55TqTnr7Jadeam/n4+KhMmTLat2/ffZlnCLz3QFJSkq5cuaLq1avL1dVVK1eutC/bs2eP4uLiUh3LWrx4cfn5+Tlsk5CQoF9//TXVbe4Xyf0i/XtMzZs3V548ebRgwYIM3bk6fPiwTp8+LX9//8xu6j13Y9/cLCYmRpJSPc6ccs4k+/zzz/X444+rcOHC6a7vfjxnmjZtqu3btysmJsY+1ahRQx07drT/P6deZ27VNy4uLjnyWnM7/XKznHKdSU/f5LRrzc0uXLig/fv3y9/f//7MM3f9Y3E5zGuvvWbWrl1rYmNjzbZt28xrr71mbDabWbZsmTHm38d4FCtWzKxatcps3rzZ1KlTx9SpU8ehjrJly5rvv//ePj969Gjj4+Nj5s+fb7Zt22Zat2593z36Ja1+iY+PN6GhoaZy5cpm3759Do9AuX79ur2OG/vl/PnzZsCAASY6OtrExsaaFStWmIceesiULl3aXL58OasOM0PS6pt9+/aZESNGmM2bN5vY2Fgzf/58U6JECdOgQQOHOnLaOZNs7969xmazmSVLljitw6rnzM1u/lR5Tr3OOHNj3+T0a82NbuyXnHydccbZUxpy4rWmf//+Zs2aNSY2NtasX7/eNGvWzBQqVMicOHHCGHP/XWcIvJns2WefNUFBQSZPnjymcOHCpmnTpg6/oP/55x/To0cPU6BAAePh4WGeeOIJc/ToUYc6JJlp06bZ55OSkszgwYONr6+vcXNzM02bNjV79uy5V4eUKdLql9WrVxtJTqfY2Fh7HTf2y6VLl0zz5s1N4cKFjaurqwkKCjLdu3c3x44dy4KjuzNp9U1cXJxp0KCBKViwoHFzczOlSpUyAwcOTPF8zJx2ziQbNGiQCQwMNImJiU7rsOo5c7Obf0Hn1OuMMzf2TU6/1tzoxn7JydcZZ5wF3px4rWnfvr3x9/c3efLkMUWLFjXt27c3+/btsy+/364ztv/fIAAAAMCSGMMLAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALAAAASyPwAgAAwNIIvAAAALA0Ai8AAAAsjcALOHHgwAHZbDb798lnB7t371bt2rXl7u6uqlWrOl2nUaNG6tu37z1tlxVERUXJx8cnq5uRbXTp0kVt2rRJ1zY2m03z5s27K+25XcOGDUv1tXG/uV+OxUrXnJv7PCOvA2RfBF5kS126dJHNZtPo0aMdyufNmyebzZZFrcpaQ4cOVb58+bRnzx6tXLnS6Trff/+9Ro4ceY9blr3cL0HhXshoGJk0aZKioqLStc3Ro0cVHh6e7n1lZ1kZ5gYMGJDq6/x+N2zYMNlsNtlsNuXOnVuFChVSgwYNNHHiRF25csVh3UaNGtnXvXF68cUX7eusXbtWTZo0UcGCBeXh4aHSpUsrMjJSV69eTbUNzv5As3Kfg8CLbMzd3V3vvvuuzp49m9VNyTRpXYBvZf/+/apfv76CgoL0wAMPOF2nYMGCyp8/f4b3AUiSt7d3uu94+/n5yc3N7e406D6Xkde9p6dnqq9zK6hYsaKOHj2quLg4rV69Wu3atdOoUaNUt25dnT9/3mHd7t276+jRow7TmDFjJEl//PGHWrRooRo1amjdunXavn27PvjgA+XJk0eJiYnpapPV+zynI/Ai22rWrJn8/Pw0atSoVNdxdjdv4sSJCg4Ots8nvy31zjvvyNfXVz4+PhoxYoSuX7+ugQMHqmDBgnrwwQc1bdq0FPXv3r1bdevWlbu7uypVqqS1a9c6LN+xY4fCw8Pl6ekpX19fPfPMMzp16pR9eaNGjdSrVy/17dtXhQoVUlhYmNPjSEpK0ogRI/Tggw/Kzc1NVatW1dKlS+3LbTabtmzZohEjRshms2nYsGFO67n5jlRwcLDeeustde7cWZ6engoKCtKCBQt08uRJtW7dWp6enqpSpYo2b95s3+b06dOKiIhQ0aJF5eHhocqVK+ubb75x2M/58+fVsWNH5cuXT/7+/powYUKKfV+5ckUDBgxQ0aJFlS9fPoWGhmrNmjVO233jcU6ZMkXh4eHKmzevSpQooblz5zqs8+qrr6pMmTLy8PBQiRIlNHjwYF27dk3Sv0MThg8frt9//91+Jyj5TuW5c+f0wgsvyNfX1/7zXLhwoUPdP/74o8qXLy9PT0+1aNFCR48edVj+2WefqXz58nJ3d1e5cuX00Ucf2ZddvXpVvXr1kr+/v9zd3RUUFJTmuStJX3zxhSpWrCg3Nzf5+/urV69e9mVxcXH2n5GXl5eefvppHT9+3L48+dyfMWOGgoOD5e3trQ4dOtjDQpcuXbR27VpNmjTJ3hcHDhxQYmKiunXrpuLFiytv3rwqW7asJk2a5NCum9/KbdSokXr37q1XXnlFBQsWlJ+fX4pz8MY7ZslDgr7//ns1btxYHh4eCgkJUXR0tMM2n376qQIDA+Xh4aEnnnhC48ePv2XQPnz4sCIiIlSwYEHly5dPNWrU0K+//up0XWd3aNu0aaMuXbrY5z/66COVLl1a7u7u8vX11VNPPZVm/0kZe90bYzRs2DAVK1ZMbm5uCggIUO/evVM9ztTeXn/vvffk7++vBx54QD179rSf+6mZP3++HnroIbm7u6tEiRIaPny4rl+/bl8+fvx4Va5cWfny5VNgYKB69OihCxcuONSxfv16NWrUSB4eHipQoIDCwsIcbkQkJSWleW44kzt3bvn5+SkgIECVK1fWSy+9pLVr12rHjh169913Hdb18PCQn5+fw+Tl5SVJWrZsmfz8/DRmzBhVqlRJJUuWVIsWLfTpp58qb968Tved/PvhiSeekM1ms8/f6t2hpKQkjRo1yv7aCQkJSXF9QvZF4EW25eLionfeeUcffPCBDh8+fEd1rVq1SkeOHNG6des0fvx4DR06VI899pgKFCigX3/9VS+++KJeeOGFFPsZOHCg+vfvr99++0116tRRq1atdPr0aUn/BqgmTZqoWrVq2rx5s5YuXarjx4/r6aefdqhj+vTpypMnj9avX6+pU6c6bd+kSZM0btw4vffee9q2bZvCwsL0+OOPa+/evZL+fbu4YsWK6t+/v44ePaoBAwbc9rFPmDBB9erV02+//aaWLVvqmWeeUefOndWpUydt3bpVJUuWVOfOnWWMkSRdvnxZ1atX16JFi7Rjxw49//zzeuaZZ7Rx40Z7nf369dP69eu1YMECLV++XD/99JO2bt3qsN9evXopOjpas2bN0rZt29SuXTu1aNHCfkypGTx4sNq2bavff/9dHTt2VIcOHbRr1y778vz58ysqKkp//PGHJk2apE8//VQTJkyQJLVv3179+/e33z06evSo2rdvr6SkJIWHh2v9+vX66quv9Mcff2j06NFycXGx13vp0iW99957mjFjhtatW6e4uDiHfp45c6aGDBmit99+W7t27dI777yjwYMHa/r06ZKk999/XwsWLNDs2bO1Z88ezZw50+EPr5tNmTJFPXv21PPPP6/t27drwYIFKlWqlKR/f7G2bt1aZ86c0dq1a7V8+XL99ddfat++vUMd+/fv17x587Rw4UItXLhQa9eutQ8DmjRpkurUqeNwdywwMFBJSUl68MEHNWfOHP3xxx8aMmSIXn/9dc2ePTvNn8v06dOVL18+/frrrxozZoxGjBih5cuXp7nNG2+8oQEDBigmJkZlypRRRESEPWytX79eL774ovr06aOYmBg98sgjevvtt9Os78KFC2rYsKH+/vtvLViwQL///rteeeUVJSUlpbldajZv3qzevXtrxIgR2rNnj5YuXaoGDRpISr3/Mvq6/+677zRhwgR9/PHH2rt3r+bNm6fKlSunq72rV6/W/v37tXr1ak2fPl1RUVFpDj356aef1LlzZ/Xp00d//PGHPv74Y0VFRTn0c65cufT+++9r586dmj59ulatWqVXXnnFvjwmJkZNmzZVhQoVFB0drZ9//lmtWrVyuHuakXPDmXLlyik8PFzff//9bW/j5+eno0ePat26dbe9zaZNmyRJ06ZN09GjR+3ztzJq1Ch9+eWXmjp1qnbu3KmXX35ZnTp1SnEjBNmUAbKhyMhI07p1a2OMMbVr1zbPPvusMcaYH374wdx42g4dOtSEhIQ4bDthwgQTFBTkUFdQUJBJTEy0l5UtW9Y8/PDD9vnr16+bfPnymW+++cYYY0xsbKyRZEaPHm1f59q1a+bBBx807777rjHGmJEjR5rmzZs77PvQoUNGktmzZ48xxpiGDRuaatWq3fJ4AwICzNtvv+1QVrNmTdOjRw/7fEhIiBk6dGia9TRs2ND06dPHPh8UFGQ6depknz969KiRZAYPHmwvi46ONpLM0aNHU623ZcuWpn///sYYYxISEoyrq6uZM2eOffm5c+eMh4eHfd8HDx40Li4u5u+//3aop2nTpmbQoEGp7keSefHFFx3KQkNDzX//+99Utxk7dqypXr26fd7ZOfHjjz+aXLly2X8uN5s2bZqRZPbt22cvmzx5svH19bXPlyxZ0nz99dcO240cOdLUqVPHGGPMSy+9ZJo0aWKSkpJSbeuNAgICzBtvvOF02bJly4yLi4uJi4uzl+3cudNIMhs3brQfp4eHh0lISLCvM3DgQBMaGmqfv/l8SE3Pnj1N27Zt7fM3vv6S66lfv77DNjVr1jSvvvqqfV6S+eGHH4wx/3v9fPbZZynav2vXLmOMMe3btzctW7Z0qLNjx47G29s71XZ+/PHHJn/+/Ob06dNOl9/8s3d2/K1btzaRkZHGGGO+++474+Xl5dCHN3K2fUZf9+PGjTNlypQxV69eTfX40jqW5OvY9evX7WXt2rUz7du3T7WOpk2bmnfeecehbMaMGcbf3z/VbebMmWMeeOAB+3xERISpV69equvfzrlxM2ev0WSvvvqqyZs3r0P9rq6uJl++fA7TV199ZYz599rdpUsXI8n4+fmZNm3amA8++MDEx8enun9jHM/X1Np14+vg8uXLxsPDw2zYsMFhm27dupmIiIg094XsgTu8yPbeffddTZ8+3eEuX3pVrFhRuXL973T39fV1uLvi4uKiBx54QCdOnHDYrk6dOvb/586dWzVq1LC34/fff9fq1avl6elpn8qVKyfp3ztvyapXr55m2xISEnTkyBHVq1fPobxevXp3dMzJqlSpYv+/r6+vJDkce3JZ8rEnJiZq5MiRqly5sgoWLChPT0/9+OOPiouLkyT99ddfunbtmmrVqmWvw9vbW2XLlrXPb9++XYmJiSpTpoxD/6xdu9ahb5y5sc+T52/sh2+//Vb16tWTn5+fPD099eabb9rblpqYmBg9+OCDKlOmTKrreHh4qGTJkvZ5f39/e59cvHhR+/fvV7du3RyO56233rIfT5cuXRQTE6OyZcuqd+/eWrZsWar7OnHihI4cOaKmTZs6Xb5r1y4FBgYqMDDQXlahQgX5+Pg49EVwcLDDmO0b25yWyZMnq3r16ipcuLA8PT31ySef3LIPbzyPbndfN27j7+8v6X/n2Z49exzOIUkp5m8WExOjatWqqWDBgmmud7seeeQRBQUFqUSJEnrmmWc0c+ZMXbp0Kc1tMvq6b9eunf755x+VKFFC3bt31w8//OAwtOB2VKxY0eFdiVv9DH7//XeNGDHCoa3Jd6yTj3PFihVq2rSpihYtqvz58+uZZ57R6dOn7cuT7/CmJSPnRmqMMSk+mNyxY0fFxMQ4TI8//rikf6/d06ZN0+HDhzVmzBgVLVpU77zzjv1dnsyyb98+Xbp0SY888ohDf3755Ze3vKYhe8id1Q0AbqVBgwYKCwvToEGDHMbeSf++HWf+/1vxyZyNaXN1dXWYt9lsTsvS89bohQsX1KpVqxTjzaT//XKXpHz58t12nXfDjceZ/IvEWVnysY8dO1aTJk3SxIkT7WP7+vbtm64P3ly4cEEuLi7asmWLwy9o6d8PhmRUdHS0OnbsqOHDhyssLEze3t6aNWuWxo0bl+Z2qY3lu5Gz8yH53Eoe0/jpp58qNDTUYb3k43vooYcUGxurJUuWaMWKFXr66afVrFkzp2P8bqc9tyMj5/CsWbM0YMAAjRs3TnXq1FH+/Pk1duzYVMfB3sm+0jrPMiK9/Xar60P+/Pm1detWrVmzRsuWLdOQIUM0bNgwbdq0KdWxxBl93QcGBmrPnj1asWKFli9frh49emjs2LFau3Ztir5NTXp/BhcuXNDw4cP15JNPpljm7u6uAwcO6LHHHtN///tfvf322ypYsKB+/vlndevWTVevXpWHh0eGXzsZ/Tnv2rVLxYsXdyjz9va2D/dJTdGiRfXMM8/omWee0ciRI1WmTBlNnTpVw4cPz1A7bpZ8DVi0aJGKFi3qsIwPa94fuMOL+8Lo0aP1f//3fyk+9FK4cGEdO3bM4ZdaZj4795dffrH///r169qyZYvKly8v6d+As3PnTgUHB6tUqVIOU3pCrpeXlwICArR+/XqH8vXr16tChQqZcyDpsH79erVu3VqdOnVSSEiISpQooT///NO+vESJEnJ1dXUY9xYfH++wTrVq1ZSYmKgTJ06k6Bs/P780939jnyfPJ/f5hg0bFBQUpDfeeEM1atRQ6dKldfDgQYf1nX06u0qVKjp8+LBDG9PD19dXAQEB+uuvv1Icz42/nL28vNS+fXt9+umn+vbbb/Xdd9/pzJkzKerLnz+/goODU30EUvny5XXo0CEdOnTIXvbHH3/o3Llz6TonnPXF+vXrVbduXfXo0UPVqlVTqVKlsuQOVdmyZVOMnbzVWMoqVaooJibGaZ86U7hwYYe7fImJidqxY4fDOrlz51azZs00ZswYbdu2TQcOHNCqVaskOe+/O3nd582bV61atdL777+vNWvWKDo6Wtu3b7+tY8mIhx56SHv27EnRzlKlSilXrlzasmWLkpKSNG7cONWuXVtlypTRkSNHHOqoUqXKPXtU1+7du7V06VK1bdv2juopUKCA/P39dfHixVTXcXV1TddTHCpUqCA3NzfFxcWl6Msb34lB9sUdXtwXKleurI4dO+r99993KG/UqJFOnjypMWPG6KmnntLSpUu1ZMkS+yd479TkyZNVunRplS9fXhMmTNDZs2f17LPPSpJ69uypTz/9VBEREfZPKO/bt0+zZs3SZ599luLOZloGDhyooUOHqmTJkqpataqmTZummJgYzZw5M1OOIz1Kly6tuXPnasOGDSpQoIDGjx+v48eP24NW/vz5FRkZaX/CRZEiRTR06FDlypXLfhevTJky6tixozp37qxx48apWrVqOnnypFauXKkqVaqoZcuWqe5/zpw5qlGjhurXr6+ZM2dq48aN+vzzz+1ti4uL06xZs1SzZk0tWrRIP/zwg8P2wcHBio2NtQ9jyJ8/vxo2bKgGDRqobdu2Gj9+vEqVKqXdu3fLZrOpRYsWt9Uvw4cPV+/eveXt7a0WLVroypUr2rx5s86ePat+/fpp/Pjx8vf3V7Vq1ZQrVy7NmTNHfn5+qd4pHDZsmF588UUVKVJE4eHhOn/+vNavX6+XXnpJzZo1s5/zEydO1PXr19WjRw81bNhQNWrUuK32JvfFr7/+qgMHDsjT01MFCxZU6dKl9eWXX+rHH39U8eLFNWPGDG3atCnFXbW77aWXXlKDBg00fvx4tWrVSqtWrdKSJUvSfM52RESE3nnnHbVp00ajRo2Sv7+/fvvtNwUEBKQYCiNJTZo0Ub9+/bRo0SKVLFlS48eP17lz5+zLFy5cqL/++ksNGjRQgQIFtHjxYiUlJdmH5zjrv4y+7qOiopSYmKjQ0FB5eHjoq6++Ut68eRUUFHRnHZmGIUOG6LHHHlOxYsX01FNPKVeuXPr999+1Y8cOvfXWWypVqpSuXbumDz74QK1atXL6wdpBgwapcuXK6tGjh1588UXlyZPH/hixQoUKZbht169f17Fjx5SUlKTTp09rzZo1euutt1S1alUNHDjQYd1Lly7p2LFjDmVubm4qUKCAPv74Y8XExOiJJ55QyZIldfnyZX355ZfauXOnPvjgg1T3n/wHZ7169ex1pSV//vwaMGCAXn75ZSUlJal+/fqKj4/X+vXr5eXlpcjIyAz3Be4N7vDivjFixIgUb5OVL19eH330kSZPnqyQkBBt3LgxXU8wuJXRo0dr9OjRCgkJ0c8//6wFCxbYL/LJd2UTExPVvHlzVa5cWX379pWPj4/DeOHb0bt3b/Xr10/9+/dX5cqVtXTpUi1YsEClS5fOtGO5XW+++aYeeughhYWFqVGjRvLz80vxbUPjx49XnTp19Nhjj6lZs2aqV6+e/XFdyaZNm6bOnTurf//+Klu2rNq0aaNNmzapWLFiae5/+PDhmjVrlqpUqaIvv/xS33zzjT1sP/7443r55ZfVq1cvVa1aVRs2bNDgwYMdtm/btq1atGihxo0bq3DhwvZHqn333XeqWbOmIiIiVKFCBb3yyivpusPz3HPP6bPPPtO0adNUuXJlNWzYUFFRUfagmD9/fo0ZM0Y1atRQzZo1deDAAS1evDjVcyEyMlITJ07URx99pIoVK+qxxx6zP8HCZrNp/vz5KlCggBo0aKBmzZqpRIkS+vbbb2+7vdK/D9J3cXFRhQoVVLhwYcXFxemFF17Qk08+qfbt2ys0NFSnT59Wjx490lVvZqhXr56mTp2q8ePHKyQkREuXLtXLL7/scA7dLE+ePFq2bJmKFCmiRx99VJUrV07xtI0bPfvss4qMjFTnzp3VsGFDlShRQo0bN7Yv9/Hx0ffff68mTZqofPnymjp1qr755htVrFhRkvP+y+jr3sfHR59++qnq1aunKlWqaMWKFfq///u/u/rc17CwMC1cuFDLli1TzZo1Vbt2bU2YMMEeskNCQjR+/Hi9++67qlSpkmbOnJniUXplypTRsmXL9Pvvv6tWrVqqU6eO5s+fr9y57+x+2c6dO+Xv769ixYqpUaNGmj17tgYNGqSffvopxbCnTz/9VP7+/g5TRESEpH/HfV+4cEEvvviiKlasqIYNG+qXX37RvHnz1LBhw1T3P27cOC1fvlyBgYGqVq3abbV55MiRGjx4sEaNGqXy5curRYsWWrRo0T3/YxEZYzM3D3ACgHS6ePGiihYtqnHjxqlbt24Zrsdms+mHH37g6zxzqO7du2v37t366aefsropACyGIQ0A0u23337T7t27VatWLcXHx2vEiBGSpNatW2dxy3A/ee+99/TII48oX758WrJkiaZPn+7wZR4AkFkIvAAy5L333tOePXuUJ08eVa9eXT/99NMdjelDzrNx40aNGTNG58+fV4kSJfT+++/rueeey+pmAbAghjQAAADA0vjQGgAAACyNwAsAAABLI/ACAADA0gi8AAAAsDQCLwAAACyNwAsAAABLI/ACAADA0gi8AAAAsLT/B3w6rFZXrTMPAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "plt.figure(figsize=(8,8))\n", + "plt.xlabel(\"Number of image patches containing clusters in each DES tile\")\n", + "plt.ylabel(\"Counts\")\n", + "plt.title(\"Histogram of random sample of 100 DES tiles\")\n", + "plt.hist(num_clusters_vals)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "41.04" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.array(num_clusters_vals).mean()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}