From 25fb244ea5d96631acafe98bf12733b9e47e4fb8 Mon Sep 17 00:00:00 2001 From: Olivier Labayle Date: Fri, 6 Sep 2024 15:07:31 +0100 Subject: [PATCH 1/2] clean repo from sieve and try link docker build to tagbot --- .../{Release.yml => PublishDockerImage.yml} | 6 +- Project.toml | 3 - docs/make.jl | 2 +- docs/src/cli.md | 2 +- docs/src/sieve_variance.md | 13 - docs/src/tmle_estimation.md | 84 +--- experiments/parameters.binary.gwas.yaml | 118 ------ experiments/parameters.continuous.gwas.yaml | 118 ------ experiments/parameters.phewas.yaml | 391 ------------------ experiments/runtime.jl | 89 ---- src/TMLECLI.jl | 8 +- src/cli.jl | 59 +-- src/sieve_variance.jl | 256 ------------ test/data/grm/test.grm.bin | Bin 75660 -> 0 bytes test/data/grm/test.grm.id | 194 --------- test/runner.jl | 9 +- test/runtests.jl | 1 - test/sieve_variance.jl | 368 ----------------- test/summary.jl | 3 +- tmle.jl | 2 +- 20 files changed, 21 insertions(+), 1705 deletions(-) rename .github/workflows/{Release.yml => PublishDockerImage.yml} (95%) delete mode 100644 docs/src/sieve_variance.md delete mode 100644 experiments/parameters.binary.gwas.yaml delete mode 100644 experiments/parameters.continuous.gwas.yaml delete mode 100644 experiments/parameters.phewas.yaml delete mode 100644 experiments/runtime.jl delete mode 100644 src/sieve_variance.jl delete mode 100644 test/data/grm/test.grm.bin delete mode 100644 test/data/grm/test.grm.id delete mode 100644 test/sieve_variance.jl diff --git a/.github/workflows/Release.yml b/.github/workflows/PublishDockerImage.yml similarity index 95% rename from .github/workflows/Release.yml rename to .github/workflows/PublishDockerImage.yml index 56b39ad..1c0429c 100644 --- a/.github/workflows/Release.yml +++ b/.github/workflows/PublishDockerImage.yml @@ -2,8 +2,10 @@ name: Publish Docker Image on: workflow_dispatch: - release: - types: [published] + workflow_run: + workflows: [TagBot] + types: + - completed jobs: push-to-registries: diff --git a/Project.toml b/Project.toml index fcb66b7..7a3b515 100644 --- a/Project.toml +++ b/Project.toml @@ -15,14 +15,12 @@ EvoTrees = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5" GLMNet = "8d5ece8b-de18-5317-b113-243142960cc6" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7" MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d" MLJLinearModels = "6ee0df7b-362f-4a72-a706-9e79364fb692" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7" MLJXGBoostInterface = "54119dfa-1dab-4055-a167-80440f4f7a91" -Mmap = "a63ad114-7e13-5084-954f-fe012c677804" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -42,7 +40,6 @@ EvoTrees = "0.16.5" GLMNet = "0.7" JLD2 = "0.4.22" JSON = "0.21.4" -MKL = "0.6, 0.7" MLJ = "0.20.0" MLJBase = "1.0.1" MLJLinearModels = "0.10.0" diff --git a/docs/make.jl b/docs/make.jl index 3dd0d16..fe5f8e3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,7 +15,7 @@ makedocs( modules = [TMLECLI], pages=[ "Home" => "index.md", - "Command Line Interface" => ["cli.md", "tmle_estimation.md", "sieve_variance.md", "make_summary.md"], + "Command Line Interface" => ["cli.md", "tmle_estimation.md", "make_summary.md"], "MLJ Extensions" => ["models.md", "resampling.md"], ], pagesonly=true, diff --git a/docs/src/cli.md b/docs/src/cli.md index 17dc2d6..e492f05 100644 --- a/docs/src/cli.md +++ b/docs/src/cli.md @@ -25,6 +25,6 @@ Bellow is a description of the functionalities offered by the CLI. ## CLI Description ```@contents -Pages = ["tmle_estimation.md", "sieve_variance.md", "make_summary.md"] +Pages = ["tmle_estimation.md", "make_summary.md"] Depth = 5 ``` diff --git a/docs/src/sieve_variance.md b/docs/src/sieve_variance.md deleted file mode 100644 index 780812e..0000000 --- a/docs/src/sieve_variance.md +++ /dev/null @@ -1,13 +0,0 @@ -# Sieve Variance Plateau Estimation - -If the i.i.d. (independent and identically distributed) hypothesis is not satisfied, most of the traditional statistical inference theory falls apart. This is typically possible in population genetics where a study may contain related individuals. Here we leverage a non-parametric method called [Sieve Variance Plateau](https://biostats.bepress.com/ucbbiostat/paper322/) (SVP) estimation. The hypothesis is that the dependence between individuals is sufficiently small, so that our targeted estimator will still be asymptotically unbiased, but its variance will be under estimated. In brief, the SVP estimator computes a variance estimate for a range of thresholds 𝜏, by considering individuals to be independent if their distance exceeds 𝜏. As the distance threshold 𝜏 increases, fewer individuals are assumed to be independent. The maximum of this curve is the most conservative estimate of the variance of the target parameter estimator and constitutes our SVP corrected variance estimator. - -## [Usage](@id svp_command) - -```bash -tmle sieve-variance-plateau --help -``` - -```@docs -sieve_variance_plateau -``` diff --git a/docs/src/tmle_estimation.md b/docs/src/tmle_estimation.md index 980ce73..19caff9 100644 --- a/docs/src/tmle_estimation.md +++ b/docs/src/tmle_estimation.md @@ -133,86 +133,4 @@ For further information, we recommend you have a look at both: ## Note on Outputs -We can output results in three different formats: HDF5, JSON and JLS. By default no output is written, so you need to specify at least one. An output can be generated by specifying an output filename for it. For instance `--outputs.json.filename=output.json` will output a JSON file. Note that you can generate multiple formats at once, e.g. `--outputs.json.filename=output.json --outputs.hdf5.filename=output.hdf5` will output both JSON and HDF5 result files. Another important output option is the `pval_threshold`. Each estimation result is accompanied by an influence curve vector and by default these vectors are erased before saving the results because they typically take up too much space and are not usually needed. In some occasions you might want to keep them and this can be achieved by specifiying the output's `pval_threhsold`. For instance `--outputs.hdf5.pval_threshold=1.` will keep all such vectors because all p-values lie in between 0 and 1. - -In order to run sieve variance plateau correction after a TMLE run you need to save the results in HDF5 format with influence curve vectors. Furthermore, you will need to save the sample-ids associated with each result. A complete option set for this could be: `--outputs.hdf5.filename=output.hdf5 --outputs.hdf5.pval_threshold=0.05 --sample_ids=true`. In this case, only those results with an individual p-value of less than ``0.05`` will keep track of their influence curves and be considered for sieve variance correction. - -## Runtime - -Targeted Learning can quickly become computationally intensive compared to traditional parametric inference. Here, we illustrate typical runtimes using examples from population genetics. This is because population genetics is currently the main use case for this package, but it shouldn't be understood as the only scope. In fact, the two most prominent study designs in population genetics are perfect illustrations of the computational complexity associated with Targeted Learning. - -### Preliminary - -Remember that for each estimand of interest, Targeted Learning requires 3 main ingredients that drive computational complexity: - -- An estimator for the propensity score: `G(T, W) = P(T|W)`. -- An estimator for the outcome's mean: `Q(T, W) = E[Y|T, W]`. -- A targeting step towards the estimand of interest. - -While the targeting step has a fixed form, both `G` and `Q` require specification of learning algorithms that can range from simple generalized linear models to complex Super Learners. In general, one doesn't know how the data has been generated and the model space should be kept as large as possible in order to provide valid inference. This means we recommend the use Super Learning for both `G` and `Q` as it comes with asymptotic theoretical guarantees. However, Super Learning is an expensive procedure and, depending on the context, might become infeasible. Also, notice that while the targeting step is specific to a given estimand, `G` and `Q` are only specific to the variables occuring in the causal graph. This means that they can potentially be cleverly reused across the estimation of multiple estimands. Note that this clever reuse, is already baked into this package, and nothing needs to be done beside specifying the learning algorithms for `G` and `Q`. The goal of the subsequent sections is to provide some examples, guiding the choice of those learning algorithms. - -In what follows, `Y` is an outcome of interest, `W` a set of confounding variables and `T` a genetic variation. Genetic variations are usually represented as a pair of alleles corresponding to an individual's genotype. We will further restrict the scope to bi-allelic single nucleotide variations. This means that, at a given locus where the two alleles are `A` and `C`, an individual could have any of the following genotype: `AA`, `AC`, `CC`. Those will be our treatment values. - -For all the following experiments: - -- The Julia script can be found at [experiments/runtime.jl](https://github.com/TARGENE/TMLECLI.jl/tree/main/experiments/runtime.jl). -- The various estimators used below are further described in the[estimators-configs](https://github.com/TARGENE/TMLE.jl/tree/main/estimators-configs) folder. - -### Multiple treatment contrasts - -In a classic randomized control trial, the treatment variable can only take one of two levels: `treated` or `not treated`. In out example however, any genetic variation takes its values from three different levels. As such, the `treated` and `not treated` levels need to be defined and any of the following contrasts can be of interest: - -- `AA` -> `AC` -- `AC` -> `CC` -- `AA` -> `CC` - -For a given outcome and genetic variation, for each contrast, both `G` and `Q` are actually invariant. This shows a first level of reduction in computational complexity. **Both `G` and `Q` need to be fitted only once across multiple treatment contrasts and only the targeting step needs to be carried out again.** - -### The PheWAS study design - -In a PheWAS, one is interested in the effect of a genetic variation across many outcomes (typically around 1000). Because the treatment variable is always the same, the propensity score `G` can be reused across all parameters, which drastically reduces computational complexity. - -```@raw html -
-PheWAS -
-``` - -With this setup in mind, the computational complexity is mostly driven by the specification of the learning algorithms for `Q`, which will have to be fitted for each outcome. For 10 outcomes, we estimate the 3 Average Treatment Effects corresponding to the 3 possible treatment contrasts defined in the previous section. There are thus two levels of reuse of `G` and `Q` in this study design. In the table below are presented some runtimes for various specifications of `G` and `Q` using a single cpu. The "Unit runtime" is the average runtime across all estimands and can roughly be extrapolated to bigger studies. - -| Estimator | Unit runtime (s) | Extrapolated runtime to 1000 outcomes | -| --- | :---: | :---: | -| `glm.` | 4.65 | ≈ 1h20 | -| `glmnet` | 7.19 | ≈ 2h | -| `G-superlearning-Q-glmnet` | 50.05| ≈ 13h45 | -| `superlearning` | 168.98 | ≈ 46h | - -Depending on the exact setup, this means one can probably afford to use Super Learning for at least the estimation of `G` (and potentially also for `Q` for a single PheWAS). This turns out to be a great news because TMLE is a double robust estimator. As a reminder, it means that only one of the estimators for `G` or `Q` needs to converge sufficiently fast to the ground truth to guarantee that our estimates will be asymptotically unbiased. - -Finally, note that those runtime estimates should be interpreted as worse cases, this is because: - -- Only 1 cpu is used. -- Most modern high performance computing platform will allow further parallelization. -- In the case where `G` only is a Super Learner, since the number of parameters is still relatively low in this example, it is possible that the time to fit `G` still dominates the runtime. -- Runtimes include precompilation which becomes negligible with the size of the study. - -### The GWAS study design - -In a GWAS, the outcome variable is held fixed and we are interested in the effects of very many genetic variations on this outcome (typically 800 000 for a genotyping array). The propensity score cannot be reused across parameters resulting in a more expensive run. - -```@raw html -
-GWAS -
-``` - -Again, we estimate the 3 Average Treatment Effects corresponding to the 3 possible treatment contrasts. However we now look at 3 different genetic variations and only one outcome. In the table below are presented some runtimes for various specifications of `G` and `Q` using a single cpu. The "Unit runtime" is the average runtime across all estimands and can roughly be extrapolated to bigger studies. - -| Estimator file | Continuous outcome unit runtime (s) | Binary outcome unit runtime (s) | Projected Time on HPC (200 folds //) | -| --- | :---: | :---: | :---: | -| `glm` | 5.64 | 6.14 | ≈ 6h30 | -| `glmnet` | 17.46 | 22.24 | ≈ 22h | -| `G-superlearning-Q-glmnet` | 430.54 | 438.67 | ≈ 20 days | -| `superlearning` | 511.26 | 567.72 | ≈ 24 days | - -We can see that modern high performance computing platforms definitely enable this study design when using GLMs or GLMNets. It is unlikely however, that you will be able to use Super Learning for any of `P(V|W)` or `E[Y|V, W]` if you don't have privileged access to such platform. While the double robustness guarantees will generally not be satisfied, our estimate will still be targeted, which means that its bias will be reduced compared to classic inference using a parametric model. +We can output results in three different formats: HDF5, JSON and JLS. By default no output is written, so you need to specify at least one. An output can be generated by specifying an output filename for it. For instance `--outputs.json.filename=output.json` will output a JSON file. Note that you can generate multiple formats at once, e.g. `--outputs.json.filename=output.json --outputs.hdf5.filename=output.hdf5` will output both JSON and HDF5 result files. Another important output option is the `pval_threshold`. Each estimation result is accompanied by an influence curve vector and by default these vectors are erased before saving the results because they typically take up too much space and are not usually needed. In some occasions you might want to keep them and this can be achieved by specifiying the output's `pval_threhsold`. For instance `--outputs.hdf5.pval_threshold=1.` will keep all such vectors because all p-values lie in between 0 and 1. \ No newline at end of file diff --git a/experiments/parameters.binary.gwas.yaml b/experiments/parameters.binary.gwas.yaml deleted file mode 100644 index 54d93a8..0000000 --- a/experiments/parameters.binary.gwas.yaml +++ /dev/null @@ -1,118 +0,0 @@ -Parameters: - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs3755967 = (case = "CT", control = "CC"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs3755967 = (case = "TT", control = "CT"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs3755967 = (case = "TT", control = "CC"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs1045570 = (case = "TG", control = "GG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs1045570 = (case = "TT", control = "TG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs1045570 = (case = "TT", control = "GG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE diff --git a/experiments/parameters.continuous.gwas.yaml b/experiments/parameters.continuous.gwas.yaml deleted file mode 100644 index 3b88f9b..0000000 --- a/experiments/parameters.continuous.gwas.yaml +++ /dev/null @@ -1,118 +0,0 @@ -Parameters: - - target: Basophill percentage - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs3755967 = (case = "CT", control = "CC"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs3755967 = (case = "TT", control = "CT"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs3755967 = (case = "TT", control = "CC"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs1045570 = (case = "TG", control = "GG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs1045570 = (case = "TT", control = "TG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs1045570 = (case = "TT", control = "GG"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE diff --git a/experiments/parameters.phewas.yaml b/experiments/parameters.phewas.yaml deleted file mode 100644 index c7053c3..0000000 --- a/experiments/parameters.phewas.yaml +++ /dev/null @@ -1,391 +0,0 @@ -Parameters: - - target: ankylosing spondylitis - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: ankylosing spondylitis - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: ankylosing spondylitis - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: J32 Chronic sinusitis - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: J32 Chronic sinusitis - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: J32 Chronic sinusitis - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Basophill percentage - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Non-oily fish intake - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Non-oily fish intake - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: Non-oily fish intake - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: H01 Other inflammation of eyelid - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: H01 Other inflammation of eyelid - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: H01 Other inflammation of eyelid - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: L55-L59 Radiation-related disorders of the skin and subcutaneous tissue - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: L55-L59 Radiation-related disorders of the skin and subcutaneous tissue - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: L55-L59 Radiation-related disorders of the skin and subcutaneous tissue - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E70-E90 Metabolic disorders - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: D00-D09 In situ neoplasms - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: D00-D09 In situ neoplasms - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: D00-D09 In situ neoplasms - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: I95-I99 Other and unspecified disorders of the circulatory system - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: I95-I99 Other and unspecified disorders of the circulatory system - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: I95-I99 Other and unspecified disorders of the circulatory system - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E50-E64 Other nutritional deficiencies - treatment: (rs7971418 = (case = "CA", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E50-E64 Other nutritional deficiencies - treatment: (rs7971418 = (case = "CC", control = "CA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE - - target: E50-E64 Other nutritional deficiencies - treatment: (rs7971418 = (case = "CC", control = "AA"),) - confounders: - - PC1 - - PC2 - - PC3 - - PC4 - - PC5 - - PC6 - covariates: - - Age-Assessment - - Genetic-Sex - type: ATE diff --git a/experiments/runtime.jl b/experiments/runtime.jl deleted file mode 100644 index 953539a..0000000 --- a/experiments/runtime.jl +++ /dev/null @@ -1,89 +0,0 @@ -using ArgParse -using TMLECLI - -const ESTIMATORS = [ - "glm", - "glmnet", - "G-superlearning-Q-glmnet", - "superlearning" -] -const PARAMETERS = [ - "experiments/parameters.phewas.yaml", - "experiments/parameters.continuous.gwas.yaml", - "experiments/parameters.binary.gwas.yaml" -] - -function parse_commandline() - s = ArgParseSettings( - description = "Runs TMLE for 100 SNPs, 1 binary and one continuous trait. Runtime will depend on the platform.", - commands_are_required = false) - - @add_arg_table s begin - "data" - help = string("Path to the dataset, a copy is stored on datastore at: ", - "/exports/igmm/datastore/ponting-lab/olivier/misc_datasets/sample_ukb_data.csv") - required = true - default = "/exports/igmm/datastore/ponting-lab/olivier/misc_datasets/sample_ukb_data.csv" - "--estimator-file" - arg_type = String - help = "Any estimator file from: docs/src/estimators/" - "--param-file" - arg_type = String - help = """ - Any parameter file from: - experiments/parameters.phewas.yaml - experiments/parameters.binary.gwas.yaml - experiments/parameters.continuous.gwas.yaml - """ - "--verbosity" - arg_type = Int - help = "Verbosity level" - default = 0 - required = false - end - - return parse_args(s) -end - - -function main(parsed_args) - param_files = parsed_args["param-file"] isa Nothing ? PARAMETERS : [parsed_args["param-file"]] - estimator_files = parsed_args["estimator-file"] isa Nothing ? ESTIMATORS : [parsed_args["estimator-file"]] - for (paramfile, estimatorfile) in Iterators.product(param_files, estimator_files) - tmle_args = Dict( - "data" => parsed_args["data"], - "param-file" => paramfile, - "estimator-file" => estimatorfile, - "csv-out" => "runtime_output.csv", - "hdf5-out" => nothing, - "pval-threshold" => 0.05, - "chunksize" => 100, - "verbosity" => parsed_args["verbosity"], - ) - nparams = length(TMLECLI.parameters_from_yaml(paramfile)) - - # Time it: this will include precompilation time - t_start = time() - tmle_estimation(tmle_args) - t_end = time() - - totaltime = t_end - t_start - unittime = totaltime / nparams - - summary_info = - """ - TMLE was run with: - - parameter file: $paramfile. - - estimatorfile: $estimatorfile. - - Total time was: $totaltime seconds. - Unit time was: $unittime seconds. - """ - @info summary_info - end -end - -parsed_args = parse_commandline() - -main(parsed_args) - diff --git a/src/TMLECLI.jl b/src/TMLECLI.jl index e2d6a8a..75ff8a0 100644 --- a/src/TMLECLI.jl +++ b/src/TMLECLI.jl @@ -1,9 +1,5 @@ module TMLECLI -if occursin("Intel", Sys.cpu_info()[1].model) - using MKL -end - using ArgParse using DataFrames using MLJBase @@ -18,7 +14,6 @@ using JLD2 using CategoricalArrays using GLMNet using MLJModels -using Mmap using Serialization using Combinatorics using Tables @@ -35,7 +30,6 @@ include("cache_managers.jl") include("outputs.jl") include("runner.jl") include("utils.jl") -include("sieve_variance.jl") include("summary.jl") include("resampling.jl") include(joinpath("models", "glmnet.jl")) @@ -44,7 +38,7 @@ include(joinpath("models", "biallelic_snp_encoder.jl")) include(joinpath("models", "registry.jl")) include("cli.jl") -export Runner, tmle, sieve_variance_plateau, make_summary, main +export Runner, tmle, make_summary export GLMNetRegressor, GLMNetClassifier export RestrictedInteractionTransformer, BiAllelicSNPEncoder export AdaptiveCV, AdaptiveStratifiedCV, JointStratifiedCV diff --git a/src/cli.jl b/src/cli.jl index 78b484b..6e4eda9 100644 --- a/src/cli.jl +++ b/src/cli.jl @@ -5,10 +5,6 @@ function cli_settings() "tmle" action = :command help = "Run TMLE." - - "svp" - action = :command - help = "Run Sieve Variance Plateau." "merge" action = :command @@ -77,42 +73,6 @@ function cli_settings() end - @add_arg_table! s["svp"] begin - "input-prefix" - arg_type = String - help = "Input prefix to HDF5 files generated by the tmle CLI." - - "--out" - arg_type = String - help = "Output filename." - default = "svp.hdf5" - - "--grm-prefix" - arg_type = String - help = "Prefix to the aggregated GRM." - default = "GRM" - - "--verbosity" - arg_type = Int - default = 0 - help = "Verbosity level" - - "--n-estimators" - arg_type = Int - default = 10 - help = "Number of variance estimators to build for each estimate." - - "--max-tau" - arg_type = Float64 - default = 0.8 - help = "Maximum distance between any two individuals." - - "--estimator-key" - arg_type = String - help = "Estimator to use to proceed with sieve variance correction." - default = "1" - end - @add_arg_table! s["merge"] begin "prefix" arg_type = String @@ -134,8 +94,8 @@ function cli_settings() return s end -function main(args=ARGS) - settings = parse_args(args, cli_settings()) +function julia_main()::Cint + settings = parse_args(ARGS, cli_settings()) cmd = settings["%COMMAND%"] cmd_settings = settings[cmd] if cmd ∈ ("tmle", "merge") @@ -163,18 +123,7 @@ function main(args=ARGS) ) end else - sieve_variance_plateau(cmd_settings["input-prefix"]; - out=cmd_settings["out"], - grm_prefix=cmd_settings["grm-prefix"], - verbosity=cmd_settings["verbosity"], - n_estimators=cmd_settings["n-estimators"], - max_tau=cmd_settings["max-tau"], - estimator_key=cmd_settings["estimator-key"] - ) + throw(ArgumentError(string("Unknown command: ", cmd))) end -end - -function julia_main()::Cint - main() return 0 -end \ No newline at end of file +end diff --git a/src/sieve_variance.jl b/src/sieve_variance.jl deleted file mode 100644 index 628da92..0000000 --- a/src/sieve_variance.jl +++ /dev/null @@ -1,256 +0,0 @@ -GRMIDs(file) = CSV.File(file, - header=["FAMILY_ID", "SAMPLE_ID"], - select=["SAMPLE_ID"], - types=String) |> DataFrame - -function readGRM(prefix) - ids = GRMIDs(string(prefix, ".id")) - n = size(ids, 1) - grm_size = n*(n + 1) ÷ 2 - GRM = mmap(string(prefix, ".bin"), Vector{Float32}, grm_size) - - return GRM, ids -end - -function align_ic(ic, sample_ids, grm_ids) - leftjoin!( - grm_ids, - DataFrame(IC=ic, SAMPLE_ID=sample_ids), - on=:SAMPLE_ID - ) - aligned_ic = grm_ids.IC - select!(grm_ids, Not(:IC)) - return coalesce.(aligned_ic, 0) -end - -""" - bit_distances(sample_grm, nτs) - -Returns a matrix of shape (n_samples,nτs) where n_samples is the -size of sample_grm. -The sample_grm comes from the gcta software. -The process is as follows: -- Round between -1 and 1 as some elements may be beyond that limit -- Take 1 - this value to convert this quantity to a distance -- For each τ return if the distance between individuals is less or equal than τ -""" -function bit_distances(sample_grm, τs) - distances = 1 .- max.(min.(sample_grm, 1), -1) - return convert(Matrix{Float32}, permutedims(distances) .<= τs) -end - - -default_τs(nτs;max_τ=2) = Float32[max_τ*(i-1)/(nτs-1) for i in 1:nτs] - -retrieve_sample_ids(sample_ids::AbstractVector, batch_results) = sample_ids - -retrieve_sample_ids(index::Int, batch_results) = batch_results[index].SAMPLE_IDS - -function update_work_lists_with!(result::TMLE.JointEstimate, sample_ids, batch_results, grm_ids, results, influence_curves, n_obs) - for estimate in result.estimates - update_work_lists_with!(estimate, sample_ids, batch_results, grm_ids, results, influence_curves, n_obs) - end -end - -function update_work_lists_with!(result, sample_ids, batch_results, grm_ids, results, influence_curves, n_obs) - if length(result.IC) > 0 - sample_ids = string.(retrieve_sample_ids(sample_ids, batch_results)) - push!(influence_curves, align_ic(result.IC, sample_ids, grm_ids)) - push!(n_obs, size(sample_ids, 1)) - push!(results, result) - end -end - -function build_work_list(prefix, grm_ids; estimator_key=1) - dirname_, prefix_ = splitdir(prefix) - dirname__ = dirname_ == "" ? "." : dirname_ - hdf5files = filter( - x -> startswith(x, prefix_) && endswith(x, ".hdf5"), - readdir(dirname__) - ) - hdf5files = sort([joinpath(dirname_, x) for x in hdf5files]) - - influence_curves = Vector{Float32}[] - n_obs = Int[] - results = [] - for hdf5file in hdf5files - jldopen(hdf5file) do io - for key in keys(io) - batch_results = io[key] - for nt_result in batch_results - result = nt_result[estimator_key] - result isa FailedEstimate && continue - sample_ids = nt_result.SAMPLE_IDS - update_work_lists_with!( - result, - sample_ids, - batch_results, - grm_ids, results, - influence_curves, - n_obs - ) - end - end - end - end - influence_curves = length(influence_curves) > 0 ? reduce(vcat, transpose(influence_curves)) : Matrix{Float32}(undef, 0, 0) - return results, influence_curves, n_obs -end - - -""" - normalize(variances, n_observations) - -Divides the variance estimates by the effective number of observations -used for each phenotype at estimation time. -""" -normalize!(variances, n_observations) = - variances ./= permutedims(n_observations) - - -""" - aggregate_variances(influence_curves, indicator, sample) - -This function computes the sum for a single index i, see also `compute_variances`. -As the GRM is symetric it is performed as : - 2 times off-diagonal elements with j < i + diagonal term -and this for all τs. Some diagonal elements are not equal to 1 while in theory they should be, -this is an artefact of the GRM. We thus always include diagonal elements in the aggregate -whatever the value of the indicator. -""" -function aggregate_variances(influence_curves, indicator, sample) - @views begin - D_off_diag = transpose(influence_curves[:, 1:sample-1]) - D_diag = transpose(influence_curves[:, sample]) - return D_diag .* (2indicator[:, 1:sample-1] * D_off_diag .+ D_diag) - end -end - - -""" - compute_variances(influence_curves, nτs, grm_files) - -An overall variance estimate for a distance function d, a threshold τ -and influence curve D is given by: - σ̂ = 1/n ∑ᵢ∑ⱼ1(dᵢⱼ ≤ τ)DᵢDⱼ - = 1/n 2∑ᵢDᵢ∑ⱼ<ᵢ1(dᵢⱼ ≤ τ)DᵢDⱼ + 1(dᵢᵢ ≤ τ)DᵢDᵢ - -This function computes those variance estimates at each τ for all phenotypes -and queries. - -# Arguments: -- influence_curves: Array of size (n_samples, n_queries, n_phenotypes) -- grm: Vector containing the lower elements of the GRM ie: n(n+1)/2 elements -- τs: 1 row matrix containing the distance thresholds between individuals -- n_obs: Vector of size (n_phenotypes,), containing the number of effective -observations used during estimation - -# Returns: -- variances: An Array of size (nτs, n_curves) where n_curves is the number of influence curves. -""" -function compute_variances(influence_curves, grm, τs, n_obs) - n_curves, n_samples = size(influence_curves) - variances = zeros(Float32, length(τs), n_curves) - start_idx = 1 - for sample in 1:n_samples - # lower diagonal of the GRM are stored in a single vector - # that are accessed one row at a time - end_idx = start_idx + sample - 1 - sample_grm = view(grm, start_idx:end_idx) - indicator = bit_distances(sample_grm, τs) - variances .+= aggregate_variances(influence_curves, indicator, sample) - start_idx = end_idx + 1 - end - normalize!(variances, n_obs) - return variances -end - -function grm_rows_bounds(n_samples) - bounds = Pair{Int, Int}[] - start_idx = 1 - for sample in 1:n_samples - # lower diagonal of the GRM are stored in a single vector - # that are accessed one row at a time - end_idx = start_idx + sample - 1 - push!(bounds, start_idx => end_idx) - start_idx = end_idx + 1 - end - return bounds -end - -function save_results(filename, results, τs, variances) - jldopen(filename, "w") do io - io["taus"] = τs - io["variances"] = variances - io["results"] = results - end -end - -corrected_stderrors(variances) = - sqrt.(view(maximum(variances, dims=1), 1, :)) - -with_updated_std(estimate::T, std) where T <: TMLE.Estimate = T( - estimate.estimand, - estimate.estimate, - convert(Float64, std), - estimate.n, - Float64[] -) - -with_updated_std(results::AbstractVector, stds) = - [NamedTuple{(:SVP,)}([with_updated_std(result, std)]) for (result, std) in zip(results, stds)] - - -""" - sieve_variance_plateau(input_prefix; - out="svp.hdf5", - grm_prefix="GRM", - verbosity=0, - n_estimators=10, - max_tau=0.8, - estimator_key=1 - ) - -Sieve Variance Plateau CLI. - -# Args - -- `input-prefix`: Input prefix to HDF5 files generated by the tmle CLI. - -# Options - -- `-o, --out`: Output filename. -- `-g, --grm-prefix`: Prefix to the aggregated GRM. -- `-v, --verbosity`: Verbosity level. -- `-n, --n_estimators`: Number of variance estimators to build for each estimate. -- `-m, --max_tau`: Maximum distance between any two individuals. -- `-e, --estimator-key`: Estimator to use to proceed with sieve variance correction. -""" -function sieve_variance_plateau(input_prefix::String; - out::String="svp.hdf5", - grm_prefix::String="GRM", - verbosity::Int=0, - n_estimators::Int=10, - max_tau::Float64=0.8, - estimator_key::String="1" - ) - estimator_key_ = tryparse(Int, estimator_key) - estimator_key = estimator_key_ === nothing ? Symbol(estimator_key) : estimator_key_ - τs = default_τs(n_estimators;max_τ=max_tau) - grm, grm_ids = readGRM(grm_prefix) - verbosity > 0 && @info "Preparing work list." - results, influence_curves, n_obs = build_work_list(input_prefix, grm_ids, estimator_key=estimator_key) - - if length(influence_curves) > 0 - verbosity > 0 && @info "Computing variance estimates." - variances = compute_variances(influence_curves, grm, τs, n_obs) - std_errors = corrected_stderrors(variances) - results = with_updated_std(results, std_errors) - else - variances = Float32[] - end - save_results(out, results, τs, variances) - - verbosity > 0 && @info "Done." - return 0 -end diff --git a/test/data/grm/test.grm.bin b/test/data/grm/test.grm.bin deleted file mode 100644 index 1a1ba976923b6a4e135560b8de9197400d591140..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 75660 zcmb4scRi_vd(BV0x=IH_dY;72^Zc$_(` z)L2gAeHtfF+SPko@pOR-{k}ej|H_xk{ocGA#ckT5E7*=632U|&)Y4IoXTC{?;`RD9 z1~MdH$ou!69Is9A7g?C$I^6IB+o@g@H|VJcBxWt-c#Q&1@i5{)7*n6-PtZ}KXHDsz zEz_m&@4Fh)c(;O;)JGbNK3D&Sm6YoB?n=XNB5u~=D{y1}ByybV=9_3ED<37LZIe%m zcAeMceD*wOOjOS1!=(9_DIb65P>Sbe>_y8{oaYV|!!p4nvb{JjVGUVm-+$S;G!MN#J?D8Duk(Wx*Jt|G z3@D>`F7tq3ZUVS^>0&J z=_c0ps>k%9@Y?(_Jgvic3fDK*8{1A6GyW;)qv2qXzbjtRypS zepA}xGirF*M@(u>VMoX15Q*|L_#NC-CXba>=!L_SSN*{z)L+<8iTc?#4s~P*`k=4u zSO&Y!L{OZKrcFhcCN!cl-X!EuzuD(Nn)Y+S8kx^zJLDpAc`l#;gzfw&r`Ae&$2)!o?8*W`J# z9W*UJhmA$G<<~p8Ml{v7121PGZTpXUkUk+FqOGbZEy72wVf2B{qp0kc9lr8Bx~6u6 zGh?^Y@;W0*r8Jo&5~dfM@-jpoAi2<=Z<$bCcHS=Sq%&_vDfI%m{_mP?qHQm??L;m|kT3XC3LUtd zaevfcSG?;H=_R%m0 zHO>QR#5x0sw`&H4EwyV-`F)%j04~u{M7dslj_W<_2(;_Jhq&jkJoO&Zi?^TmgTr`P z^mDvTX?AUS4U>EVX}vJeQI^428MXXW&Ux;4C?+o!IZ>Poowvgrms3P?$d2mpXwncE zz2_3RbPuYNMwZ!E$ZoDM-ejgi+IshWG_9p?6pZLR=0Ejq zUw)cYIm!P&)az68B+6&pOqNHo_XolG0NO%q9FR()4*E*A#XIuyMLMU8>_3dv-OKeT z(F%qn&8Kj+{|2(ZYjaATJ}5xeV`l zJ&3-_WvWXbQ?x1Bf7)i&Espai_#6$p96`H}xe<@|@N@z8Klu^5&-gfg*2?HV>8>s; zrMyGt7(;BK6mA{dO7V_0ZbE(QGnu0M9hAgJ5=v#d%G+j1|7B1)MvLO)ANel&Z8;E* zY&s4u*4D7h@CW7T7b|c7r7E6qH(YLOFU&3^&n(n=y?<_dl|0x-Hw8u9$M9Bb&$d%pD|?EtMFkvi|A&u1S=3?t??1;Vdu+xwD(AA{ zBVs(B*#*}=uA;CV0n8piysTyyDSpkf*CfVt0oB#@);o$f&6L^0h=*>PDeDU+FJ^6U zvHN}y7k4D3{=pQVVJRs+xGNv$dxo-k%GjKB(5dJIm6cl=KyW@F;A0zj{wb1=3*YOV z<@vdE|12pU%;tjvK1q5lG8?^9boagv)pxMpBgz-kJNOpRI`(ZrJ^Eix)n)TAfp$z! z9rEpg2A3_sPe!hzG7s&Y0R6`WgYw4J#N*-yaP2#VsCwR@_>=wSfm+Nb1=@rLJ;;$b z`B;biVv<=sT-{TX>xFQ;;Qogg+khP{cs)hF=z^VH#mVZd&YdmRgQI5~aejNRD%zNgLjJ4z&XD7d zY{BfeyPGpPY@WLKmRN0@_pHw9RdY;M>$!j6e&Idi2d@$ca$l<>;w0?|-0U)YCKrA;QAtmxL(j3O1D4Ij^i~Nb6s@OfsKW`AKP*` z$`68S9yi!rEs87!=sn^VmrZdF_JrLR+w%VtVIECZQykmf^}+YO4{ROL8nO?0fnC^E zVpNd<&qldG*lYRv$;&0ozC&7cO=b2q`likm)c38v_gAKT+_{ATn46<{OqW$$j zFQ4Lvg(|jX%NGqO&3t8MC;GKbrhV4yY!y(^jfWRaU-5Ag<&RZkbLNQwcgZ`gozO*x zwVCOjo|MnCcRfYf+s6_wvyRZTCmYk}b<^YJbLT`y@zfsnJdQj~LNO@V`Tw+aZ#q;AWNVd)*8$F$ zrxUAYXB6@g?PQ`5PmLQ)`EIXRs{lt`*90rh$0CR76F-NSRNi^gkJ3FV$%D<&F@JTF z_Yp!l1Yf>t^R&Tl(U5Sa7s6kDm5;M*e?Y)yKH5S!4qzLw&heILmz@#K78_F;-;-Mt z@ul$!x<&`CqJ2lB>5BgWa?`J1a_k0@?*`+t7AFHGclY&XP2s|*8<4c( zFxLh7pmS{cQwEbtUgYOBP}{?18#F;=zJl?eJhod)u;gkjIJb)W=~Y$t7LeP9RF+yK4DFus)LR! zO)A-L(3IkIE|@_5E?qy6oTkTNg$tV}qKxPvB{c2F?mKY&ZV{ES<5xC?VR}(ZHqLe# z6UJc(M|b+fN4Wh0?Pa;!U=R%+mz80HnGc13{7~DrV|de73F{Vey8qWu_2qZ8-T&Aj zuWyW2UWa2|sq^y1v^Hs76m)#4iEHUfN?J>jp~Kqt@OIr6GBs5(e@RPO{FPd{=!dy8Fb4!5rE`!Z>iFv$rJQ`mu~oDls3n z5gLYP@-hrd*-Y!j$!B|L86x$JUF&(cv+EquJ!27#yF^Pl9P!*H^*Hh=H1nC%xd z>B+{oZ~faqQFRE#KatMrK)HCZXl01}e?cD8-hCqpGfR1%f(Nm+2uu^4txp%oQ<)39 zkF#};ef>wlZ_ zDqGZnobPjy&qs=!r@^rIwp_ngUS>4x2ZoS2qYVv{<3k{G*8wtXQW070*%vDFb3{?E z{V2ale&N)AzQmOCKzYmRb%iHpO%z}qsyFgD?h_*(?hoMFD6R`(T$m5AaX>0rk>M?| ze&+xK0_M7MtW~tCTEooy6(`T5_@DAe^ zqMC&Ia6rQB{^PoA?pAv7tEepe94`Z;PgXuk`97V{ncEI=vrmbxI4DtFhowfze5qvR zFE)lF&CTv2@p-d%ohHUk{R_|bk^c|M~Z_n+ey(vBe0w3!^>RGiS-LH_v`k-s{YIdNz%O_s{AvP z=W(q6YszPM5wo+ASDV)>xPCb=K5*RRsvz1noRj8L*v~V6=Rh4Pi0gw=$?wkk;^XsJ8}2ZDFPGm;_P=+zt_*6^nT;aw+M&VhYoley z$@LjxiXUq}o%%uQ>6||6v?rB+K8F2Y!#}e!P5f{g{OSk7ve7o8>CQed;=(|#%kfJI@Vd_`QcC3YEyS0#J(T(6$u6o(Ng~#J zsl?Ggg>0=E3mv9g^SU+QYAenUX;g<8!zjaSg*x`LjU%kK)+ZJBSY1NA+HU3WO(eFj z;v5arwT2vkBYzrG9dtXhzPUEt_zB-=T6I3FGhyw!5$n<~6jnG-aZD2Ej*tt)<5VEk z`RO&bFD;~fF#Q1M+uR~QHfqy!^LAN0A77_EqI^BPK2!g#(J=DU8`rX=lGQ`m8XEFh z*q5#M6?Ejs2U$3Hl@1apNmY%VE3J(h!i(sUC2vy=13hW@*#|cM9yq)MMm?0T zfgw$Yt~$`~Htv(wwKsw*+@W6fR0x|`$azj&$?T!*)#0F=#cZjwv3J45dK(!q{UI59 zxrey&S{hsoeox%aPXXW3Jus;)>n}quvhuw(jrA$M1Mi8PKPhwGsE^=G!$KH|tUpc=g`)A13`8~j^!i$!Z?n*nF_TY>W87+At zifQZ%>u$5QhjLdv(1IRQ#)DhqCEz~78+u2LCH^tb{{BZ&$?JQppYCh42j;EZ&*dad zXh6gLiqys1wl$;SPa1kq96Xfcw9CliVW$*t%Cqm>m;C=zJAaUzHXTB7`sP$}T%kPv z%fI0FikvoP`?;9jp!hx6KL2(d+W2E;P;!H{{T{3D5Z_HDaLUpXt`@L$(JQuVxE{6T zwApu2fmfAlTNvFEReX z`aicJ_E0%l-i9WgN+5@a20=;x!IW3+{5|}*2zF2G0=Dm&tuW1%(YY;RYYnw=T?Q1B zF=8iB`D^c&?`QTJ;#m64pu7)C27-1)6pZrQrXb6|#adg$YwfaJwC!PsF zExT8P_kwXSC&w1DS`8x|%$m}2(l;1M>2q!|zDV0?w7N)Tsigv6n^WB>zOAuU?O+1jg5zb+Kl-*q1K zUFxxQj@`i~G!EO7dAk#}uB)Hf(b*KQ@mO=J2g(=R+?9&uFTl|+5}zYCA2gw~y*t#W z{=ohI)OUQofcl8H@&0Jax3xjCNdL%LGE(|Y)VfW3QW((^hR$MRfybJqbz~#0hwFOS zctgH6HMoZb<$F$BRoqQ|4}{EM{ouSEZ~v1Xd2Dr#p=ajsybzw1Kp5_kp^r{j7sB}f z+Cm6pZ9pm^OW6DsVdzGLv9&Y|ld3f&Y5ZM^f3BV-3@u}J-O6FrusM7a<-cBbcJpzk z=Hdx+w~6MKwi0*p&a9(nNBzMRRvzKRbwd9Cx<;ClNv~4`>fa56)t}|_RwFxR{~>oeIhWwZ%MA5(Iw;8D(;V(%^_TdseOAGdg?BzKV=q`zdvIq=sjEmmsj1SIHCrO zcKXGuBuAyvHq?!cHG*D856N|Sdh5{Ne96`sI`rohpjekg8 z4Zq*kUN>f(XUWFU-V54^jr>v+%BmD3+>SuZ7u_Q$pM%BnGbhdZ7>nLdEd%@4_sOmc zY@G>lJ%*luh3k#M-)T5Y&<_m`E3!c0sA>@Eu&=T8VdW^ zud9MC$kVi^7EQx2I`GW|FN3!Vu&<+7JGZ*-B%Wt6AHrgo&P#UmB8^P5$)zIU3^wcU zxAs|1u6Er>LiaR-mHp-GQB&X8Q+ufVs1oJ>DR~juCA~|-Uh^B1tL-&N+c~TsD1LE- zT+bgu+kAVwda^QJO7ZaQH0p_G=J31_o|V9TVubzc*nEINu??(h;VJ2t&TRd8=R3oW z05{U6XFRd0I7NBRNIfZ1@l1#3nR1(Y+o4SOGW8m~anpnuy^oM|Q8F*XW5tDb_qVSk zuj8kI>IBw@1o*CjQRibw>7h1MUuE5Gc0aF2(eR?PJEr;DhY(k1KX~^`e!f+sE3>}@Ug-bH zTF-S;-^|t`e*H87qZ8^>j-9BQgor~qPNyx+X}G7?4w`O#Wd}_5?+7FN094-HL{=we zK={dV6lZK_W%wjpgR4o-=kz`SY%E*1+8k0Bej>N0hY^*xwf81uds4S`Kxbz*VhSFTf*4Ap&Sbv?cMt+{vKt&%e=#&$`3;N)G;w&+_XUKKJJ2ZHI05!sMqj+Wx z&kNyM3EU?Z_Q-^_(K(xqscf``kOtd;p<$$CL0k>l*wPhd`sn}7i`?e_xRR~uqzyku zHlJMr3-8PK+k=zOa9xm};F61U$)EddyhhkO=UJen*Pi3HwK5k)5AF*U-N=4fmTSfnUN0X_dOtH6RjcS z*JzPw=|EAzh{ylgPfU1ZEZ%!zA>~7ES@QoCc^$vAhr(??KBu%cBZrXW!*<% zk?yedbaLpb2B|i#QqT$0(p3LYn2N><@VO=e&-Kh!ndzB<EP6ft(zkry0wE$MV;@di|^Eof-}xK(C5+}GAy7SEbQD3ybizOeCr#t|8n{7`ZRBs zY;MTXq>|tawg-44xg9B(txT2_yKox3L&K!7Uv_3rNr|16kf-33jFYS)x`X5Xv)kE{%nrZdpk4VU^LGq%M$; zDUVjPg2^VOFz8VxY3q5LY*vj@_26V(+=(lg=CiGyVK2oI{?))qE|6de17E-Y zSDp}m8NvErA&!1ga6P!wp4k<#Nl-`L)tcvZU^QJ_MFU>S&$BRoQpv?c*2V??c4llW z3$@(0jxtj`c|YWv?Ggpnh<&FKFb?$MA`YIO(8g=~!3VN(T~20dxomPjep>Eq!< zt{dK^#Cv45`7M9W>W#^$o??1dLbgwgdt|sjfNP_**M-m)!f^oU?GH;O;bVolR&9kEFEyOPzeck@?%-qAKOFsan&YeHvU!-Z^f@`Mx`pd|(7iMGUT0@`pwb_{-6+ln>ygfe7Ux1Z{X%swqACojYklec1-PHc?wKKvN+ zKeq*2jg3^l@%3}4zJ7Ks#X}{BcwPlr>wV4;(}m35wz<=U>btiWD_6vsUAc}FkL?D- zn(XK0GHkpt#GNUoZS;vNvuWQ}G8^~*SS z$l?94az2}vc6clAJCXmio8g=u!LT)@MpElJRSKO-iqaI$ZI{?3R21Bi`R&I zmKahO_2cx2n-c5+n>H&``oqtbDe#e08j0<E@Zcr0gJ#P1b!{fxfxEp86!gQplz5L! z;DtV(8$}<_3!#tu#KImKu8krL=fvovErdRnv)<;p+(tskK~MQQpi-v_SQ0mc$dtoZLCC07Jhh zvkdOoAB92Fy2I%QihEZA-gyHy@56Z4Nimd9%!IyFm#+4X(6;v&8ov3sE!=A($4%XN z4<34E)3neFZ@CS0a&U8!bW)!rB-icx_e?D%=RZb}Yvtp~lRqm+^*}9Bsp$YaoaJlC z6Rc7x&&GKwaAJ1@8WwaF`VDzp9%|ZsGSo=fITbmSg)-CGcyBfMW1VvPIlEM3ot&j0 z$M*9%*u6yqYUaz^-S{VLFDU-z3>3c(|rkP;(9DsP>GKRMsM zKMl_~l1b?&%pWiMn%o9PT((irDPLzE=bx4P1lBAcr;v{F4%9d0{TjyAMs%e#x+5EL z+{)^G&^7WCtuH&b?kE-jm&3zBNp&BE551r`=fv~|ZZ$lg2OP`6V|X6u!-~H;D=7{6!rF0*!*Y8B=?ahE z0M)GoPS<28=scozB&9w2%^DgM*ejGR(xKaZ(vji}OL|CYr<$w4q7rsCywjAM(+q;gA?%Bep*(}wy^>(9c(Y&I?-Uv!hA+AHw8RpTPo%+!aIscvGo-p5I&QTF0r zOAW=^V>QKP(G56m$55%l|5*P78^;*e|N50mg!{yJj|}h7{Cj_(HWluX;r@UC!#OdA z1zQN)!0cBGBwAsEIL?{ViuQy+i~iJJ4v@N?Y{DZA=dW`0M|Gj448&ntlioW0)AYX48sw=N6BCDqq=YoeT zsdubB>6yoDu|8K*K{;&-Xln1JJOv#br@MpUOSZOwupal>985^xI*a9BWr7u_d*>Fx z!!N zN!F&+55CCStX5KMc-)GOtACnSb2^V`W^=_KYy^S3*}TH*XFSY{-woA=Cz12}{mA7F zwlI12V_G&hd!3=QAzBSZ-&=1XN882``?78{t(KZJXg;{Lx_3K+7+gq z8`YGQUaDn>{3W>1&67-S@KbcUbrjD(mcQ`65WXuR@J3&_M~3$YgtOBK$MZtMSqa=D z!~Fq-yB#Mvt(x1-tK)Z*c6-P|!RQlV_h&g(8GEa zuS*5#+6s7xJ9p(%(Y$R>D7@LSWsv?niH5sN?7?J|A$&fz9mYS3Af^{r!IKyFd48f^ zjOo4p)qRnf>rgn_WdZoTTPX5rwTkoqvonh0q8$7Wm@n%KqDO$i;ieRKPvmNd>#q!{ znpZ?NDXr{=Ph3l8&#NyUD|_$bc>jImY+p6lpSm7wjNd@xJwth3kXKsc9^hpWMzUu! z8zQh&RqW9Y*A1CFy_@xkKF@@AD3p}^-;E{PLf1jnX*Taf%=#v*{TX?Tr~05KTRZ%#2KAzsv%Wj;r9b3l zvhllU{y~bu<6(lGlfrul$wq6ZmB6XoJV~3$zTSSee*FOvW`AUGY6n zbusANfVkQkvii`8$@tdUojmK5N&=tfLAMvI{%8&S0*1b<98ib<&20_Y4g0S)p!^$6 z*-84g9MAt}P1*gyarR9}u0?NN7HjUavFOBb7mDX}Xc#mcFW;{kv|iDMGjsQX#~KGI zy-j=8Uj*L`2KUVgWJ0@HB6~kWI9kw-rYk>`?{hw%0ifE=najmIxD9Bo5dYuumC$qo z#f?|m3G>=H^E?SO>|F_YABAwdLxcAR@a#068^tqo2*t6zhgh%^1KkZZjop~R21k!`+gjV66&n8ks{)HOL3iS{M9BD_v~@`-i%Nh@6Sdjn^aS?Xk|E*}$irTh~?5v{ZErq;B-+u((hc}_Hz&33t-wEZJVBYvQ z4Zm?%OnGi<%-+EhWC=dLD}i_F@jkHttL?LMqf&`*UI_PzF&*~@gtbwO<9tA{g)sib zZ?R-~GY80Y{7U8fn_f}K8{(jIesoRrX`~4Zd%YGqExS#`I<}O*%Cfr@-`g>l`p2>e z)o5&-Q9vxU68=ydGhV61Nl8wW*}(w?sH6UdnMX zZC9T#D0%%Anr*g#9`>t=eH^pf&2KjZo$qq{eoQYmRwK>7uE8kK^Tx;0dYAMxc1Q((u^XE2y8iyaSxglFw&5o3in2`KgBRqKdUcl)umz z+Yal;G7^(SZide}PpAFN&Y7rnk!banul;Bxv-QxD9u=VZp#{a?zUUR4l%2QE>9ZET zdX0pMq&^Tx;Ji@Q73WoOMmL3)n{}Wm^*!CSN$Y~<5VfO_tnM;Gp-hpc(0?YBr1piJ zPE#qbmKrbW;2@shT1=DQH4$Lh5sjg|=`mi82%oE5OfX}m495LbXNBf zocOyHTz{9QfFry&D!em??+YO=zAJ%u>hV4?-lhEa9vPk+MLcwPR)U58{>s)yg>@kX zxa|KKY#A+S-Nua8-8iqkPD39M2@4?}=uV|Pjy8|FCte**s-bL+=ciHR7=*la^u)h;8oBY?WXg=QOWJ4QkHjkL+*;JgQF$HWr3L(Gs zQ_6GRZdTvC?atdavHSp&zA4t-xlb#J^}5rd(6z((|K(ki2^ZZu^El@5w_Zca&)wfs z0p3Z8)z|pZ%x)eqFAPTgxJ2yDw!?;S)}K8b*j+ri#)Z-&FCBGeS7SJ&DZ2xN;hk$# zI1a|qsU6ZH1&#(H$0?CC?~~J4!2KtGxDFo=-X$M1>`1P@oDZX9zZLGC{^2Q=vEi}h61^yFea!&sZ%eJ>^(){TYp*tKbW7Y~kA~_^t%rsmJ@oLb~AN*=amCif86(pBKVC zGTa{!@M_PAaU8%l@SvYJw~_AdyhuvJ+CY<-FW}U_{eO->a|>8MHu5~a?VxD0p`4Tp zp3us41l7^>z)rY6sR4PQuP54nW&qV=Zs%Co>oktEHEc(5njiAz`B}Rz3i>{0ec(tB zcW9m##?yp63~T-#PHp%|ejN#-{H*tCiANb|@jTvAJwdi?DXv3zChHc*IU-5`ULTJUfl&Mg=_dai19X$Z&r^fdAK=c*4~P zNwkt3Y9*q-Z+^rC-%qNTBsoaize#yM6)3uU&><^Am&$=;GnX+F-q0y+C{*V0WY zdDL?NsYx*-Puj}+1nvCxTpsG-*KP!~>h{=nR`hJjXWWUbT0CUuVfT)KJ{H-4Z?7Xh zI^n%hd}j{d*uwXP@NE^AD*wfMWI_zjP77i5g|ia4M^+n-b7Gtipe=-OYy%tfW=SUO zv=a{=T1b{(K1}t?xqXwS_x^bt%BvoM-ug{k?~Ar4Xxbj-Zk+d&m@jrg%}dBr?Ra?n zmW{p0|HQ4YoNlx1+_LB84wUBSg@+uDcvT0Zh~%ath578<2ECo1*nUvD0tugY!@|bu zATE{|o+}?Lk zke}ri7Qb*Je=Bjt7az@yDvhnHjPuJH;wp+6>dAIHV`lI`qnqt=qg zq;tD1nG@6Waapn$(x>FtO7AHZYGPWlAyBE`MpZ6`& zCoTo{jTi)3{;b~-=-hg<_Hi`eEEtYh1?CRlNn*FtT5L z8h=)$Lj1(VWMry4r;GDy1!3=2(s)V_Ru0HBsPqlN?+@VH>w*mQ@vTLCXAb{6;e8=| zR|3=UJ~7^<#Cv3ThX(Ht;MwWg=SK0Ygs@MHG_@Vh2hbM6aRA#ud{h**k&YT3f*l(3 z$y5Vxs^`46?Ecn~IMxT;pRge$qM zt-N?4w^1hpJI~T-HnTg_n{5Q&yT9PAjRWO{`kf8F0vQ49oq-J>yNH`vui|O7^~dn= zo!iK~0f{6vn5{?uo8Iv^o4ccubgHExr7G>b~2k`Cn+Vb$N zMInyw3kh$l2zTl+4eycFzC$CN8x`o#$2~IiacvZRbT}VCAIAaoU$tn>`ZcZTENxWwj?)9YzaSi!oV(8YR%$S8u1d0*D96t#b( z1)YXOQ=Z2P9ie!Qy13@!PcrN>%kP_?Y@Sdqv86f=b8&`BwdrKG$v{fed};~xcOK~t z<&~A9HFuYB{aP1pfiFK-L&tnJ4-zlb60aybPIdA7!pf(qV-%Ei+DK^{Yeta~L9My| z@@z^aL07+1xc{#hShbd&aeQ8yOny!h!L@tL9?6S)M;09z2rJe-pt!B>vO6ozcdpR5 zk^2n0y)oxRRZ;q|`tSc+>N}I`E0@PYU*^phIaITGgw;y+cIs5)jj&>fJ3V7$BpbWT zvvMf!xv|+I{0cvb@U zi4h0)2hhiLA)FIq7;PaO2QWPD`*c2c#&|bN_Fs;xFSok_^!md)<)fLo)U{S(W3PeI@i*NVdWp~2uCz;5y_)QVhyPk z*Kb705{i5Fju!Rjg^6j}x&DngkD-&J>zq+EAGs2u51B)$>@6dduwqgS$6wS=-e>>X z{)BWepUBfPD`f8ktFZC&v91S%{yi7aM5m1M2`)1Se`R)dB}L~hrHkIwo~Bi6O`$%^ zk&=?9lN)UD*{J~22$}@50_1DPh%4E$73@`W`TiurJ`cCxI$<1LulTMMwjv{q7^J#@ z)+Ah`Vts7Styv^x#YJ9+*>`9d9lpJeA$)Te-y6j+zB4Dhu_e5%f^fV~jK1LT{y-nk zKU^0)FN9|$aGx0W$OJrG7s5F)#?cnSaR9?{TNm)QP~pbrU0b7dZ8)a>_V$W3Ca zDQ|;9+Jh<0$e`b>&LeDTf-N+hSWo=1WQsI$pmd{(5H@H`4$*yy}o_Q~q&zT=@Wn&-XApgDk zjM)dm_sE3r&o^*OmE~n0sOop*7D-iC@9HQY>yqDUV~P;y?DOr4M+jM zi_FGGmfQMf(6lj@>P(Nm9@O7?eki9AL_ePQmj1WDx5)c^Ax-er-$YP4?`fu--eiwF z--gYtDEzeL6Q0Mz8#lnO4{U4}_@I9=YdOh#EC04j#Kc6R>CW12TDc~K6mEwdDr{~w zGeQ3SEYG#94QUid65Cc~RNio_pQM{A>zA<2n%3(G`#o7de&Kgdier>k0J1OTk?|9^ zK%*&{Tz}N1Y9jk~K&4NdU7rDLKOf=f{B3?6RQa;8 zEdA+13UApfndZBlng!4Qt)Txtk0X4CM)>{!$`an(#kUsmow?d?Y+*XysYf3j-Xp_1 zGedUrzq~waS_}CG)=Xt=RH^*x^&_Jxa<`G%s zX$S`wHsWL6mM6nFPHkTJ|1Ns6hn7Q^xb5&-cJ@g~k6L?+Z0)!kLhqg@UuOIvNtW3> z4RJg-E1n+^=mh`4C2u(MJ`GN-Ud+oG(-s@9fJ^CLDKD=`wlhx``hksTWJCNtc5Ats%rXN-V zT@OE4+f#$|Sl>sNrysm)OMb{f}Z!P^RG`8yFlY+!Bi`{q=dHt8jsciZ@9!pMxr zb!3kYI?BMeOc3yR?b`Gx(8rD*t$Li&;s2t~8 zmpE)qzI@H|<`ql&9}C~)!@SjXQc3MMcQGB`*b?yYT?xEXFT~NWeTPPX<9Q)^R^snI zG42oG+NiKDg!2Kx7Q!}=)-ym-SYl84j(aVySIAHML;{un;i&A}m!@;LE+$y)7k z@Y|w@jLF+bqGjihZU)QKk+vfDFvWK)?o0img4s~oUK>uFj)lkNjBoroTX;L<90|{t zg3;hl3VA~wWfzr*?~p1{%B*)h?Z=>9T#sGB;XI7^=$6XvxTC*T%bur=oRSZXU9M6d z3B3kUc>MH5WYyK$*>>tL+@A=x)bCUV^nr}h0nNKbDnJ`r0D4o5?HT95Iiu=9Oi*muze zW&{+LUeTrW(s=&TNpT9lA`>su4`3J_ zo|(h*LM+C9aett8TnOVhfNdb%Hk{8FFm38uHm7xcSIp@QHs#asA*Wr$=6WPeG3>(H~*zp>=TFVH-TS zd=d;-Dy|pJ+Od@XFGw>%_NG?j-K;&RmTHsf<3~VPn-wDCrE=TreE*REEpws$r3x@F zXaF5m+5B-r60@tG4mb&)v;!fszSNdD6EY~@P~7Vzd&8&m`Amp!s77+@$=edjHNBg~ z^&l%g@_aO|+%C#{$;L*+^IW`!r(qo31pUXNvvb)TcxA^S#Q7>~!x@&=qOaYuDE_u< zM@Zs;bWrWL3~DAMf@%0!64a)#xYuPD*b^MR0ygz_vr}5kVN1B^NaZh7uo~;tbg45scD0JAZ zN;W;Rg|W#}j*C3azs^&LZ#1`slat?2T@N|QZT670-6_wUwe0_QpTCUKP5XOy+9Wa> znx6hlZuBt|t5v6So>}upfu>1&@E)B@_LhATeSQ-|(`I#GYbvM%x~tvQ#mi=~zVMY@ z-TtV#(I)bwBU>Z*@%KFH=&88ZC2!-5lk7VKx0AVi0cKRr=JWmL^n+l>ljL0wHL>aq z9g$5utAmJl(Yl8Mo#W0_l4ZJr=DS{i45S#yIUX!m?`n|{@mSK@(v6g&px@nprSUYKv-M4m- z_71F_Kb@ZkPd3Iw=e=ytTHxJhIa?Emob5>EK8ig;^C6^L#VIMo5iWd}62C_#;Gr*k zf1vie)A-gRzB7k!Y}LlaJN0;<7-7OaGQ2;4aXdGQXXb?SLI^`A><{3&5XNynP}>$F zd#CcYV7ft=`_=v5;|i4E#Ps8&+(GU#N7;J!MlwKr8cZ+sB{Rmd{h%=)HN=)tlVC{NN-ATfCd)tO z>~8cfU0brYfe0QPSx#m2NLDV*QHxeZ9WK3*i!cPi=prM`IjravU1 zr{aF_fKf(#E+q12?~iv$Nrs7ze-!Ei=3_|gILf=1cb;f!pJt?d{c|o0)Bnpgd}TuU zr~x~h`f0RmJ%0eR(;v@e<(@F>0$dMZbGF;=Y<~%PHM_AAnzY{s@1t?f$o_lfdW)&d zWh);4r|d=Qzq#Iq32gq|T$6yomjaHTZ?+hE8K(e#r@r<#y79Y|_)Wfl^C5hHKzMgr zNXK{P@Qp2eUr2ac1@F`g_lfZynQ(_jpsDTSKCuwTwb9z^LINE9BV$5nTZp#4OlAe% z1MlUX#N_u}9heUv%C{-l$4Gav!w%ANv9D;XT^rK>G_y;u+@A>+n*+J-Gp_FCyi0UsZ&JkH zA-6|bP?)_x?hDG@%{%cN7olsnu`wINOH!j~nn4rRcOXpeymkYax1c(Svp zsq(t9KTuiRvO%!Odr*JMd;X-~G>>Zto3F08YEO(4)WtbYCl&M}-DG@tW|aa#%1N5{b^6Gz_81YNT%d%^J5g(Tbg8|NXU z?U_^v1FPmii?$lDUTY@j*WmBl1p=N+NfD?uU}f5juy*BrFA2WXXZOwoy#MlhRQ>nf zJmTVaDFvEOFT&v8@6h18)A;5tzBfwWT9myrC%my$`)w7xPmFgdkr&>f!TSSvcAA|V zmDMlVc_BP2f$6wMhWi5;7Mw62U~#D=saqHy7Z9fV)d!kC*MVkcJHfzjChSt8?q%x&g6@`HY#myN*Y?fVz7!4jWwz8A!}( z`vVC3*Wnvm|HAM-G2W#V;CO$)@=QM>oEsI+3*kPout!#VZ4~Ly7Q%4=!^`gn@wRYy z=xWM$&1NB@NVC?7iuzBQE8n>(RQOL8KH={dU zbf3caYSt6CxOdla({Fiz+wIGhr_<>(IoV%F@ks&|6I&-Q0!SDW7_tJ7r5w1zZ& zhg0i1;{)QV^$n(VKYKEp6I`3XWLjO{2X2I|32eDo1Oa*tD4$t=ji_H+&YQ>T;`a;W zYf2724@HYKq&$!1$;?i;*iKFCy5t9?AFuO+e5~-I;o9Hm#_v)J-{iw@F$v$H!S4?U zGV#q_AuPNzhwlpsFuYTb_lfZ?CEgNxm`sloS z1WKH);kiGlr2lwJC~xHtmo35|WEvY28J29HC6K)bm|*aemg~W+FhUGQk?9vFK^6(* z{N_Br1!AWMT<2%6)nJKL7bqPr?@Lz(y(FHajlgkiJM#W~6sTQSBL|WeavnQ$>h|qr zGhdOemn=EotRo-E^)p%UMOnU=kGk#*mO=|Vru(ZdDnvc=EVt`C-pR{a$1#-4YaaFY z+}u3b8?VDmA@(bqcO$>HJ~7a+XCET+*+=o8kO&AFmII$BzX7Mx6I34$i(-=ZM2%e6 zy8&5A)5!i43*gRb`5Mov!JWyod#z=2n|biZJB9L5_iYSq_w*ng&9aEU>@A6yA@V&9 zi#2=5mEv?*{DO_?zfZ{93+jZiM$>$Y}p;8@xQ>}?Kdc@ zw4gL!vejt2SvLC(FH`AofS0CW(7*RYsF2`zB9-KRDmrEPZ^SP6vW=}N zVH{oj0cPV2zO|d`xc|U8GHI$5!tZ6mP`fm;roFs9FbdhX_$~6 zrt1;X3=D3kK{MGm%@J4Svxxe6*|K?}bs35El7h~bOC+SHJ}kN}--o+DKbYbn&4Jys zVfJmj(V}FWuMs^Y+Cnb;_;t%Dbz*%det;vpANG< z!evi8!|4G42eR250Od5_TM66poG4!JD>`DwIVlSH+~s|YB$~#;yqlX~S#A`h4tYyu zC)^JM(*x1aFlhyuw2IkUS4Qn2_N&-97P4$A6j#Q>6>G)uarWtT(ByjpOfS_GJPVee}qVd@`29d8HK0v|e z!Ti6GS{#w3$?b~W?`k02Qw(oHITIE4~%ui;)XR1ri))^4|pa+F_ z95I=vBhH+gVmm{RvkEYT2b4NNjHraSxvjTZz8y8?xD=sxK*!~ ziryZLm7VjEzd;-9%+B}Od(H>PF{yR-jFI+TQU%u~BY^2F@Lap_AOuIVvG8`+`cQsB zO0Jn+Bn4FjQi|BP_Ul)DPA|~a_Lq*WC-zArcAMmymtmMfO5ONkS485g{r2mMuyll0Eyr?@RXIea@KIy_fgr`}oZtcV^Bx zb7o$z`+Ci5=FC~Ll>c3~Ka|UC8ubtEXY2FSDSeEFufNOhIhmZlZk&3=W<^`X3% z#9~4m8o&^L0KV7ZcN#u-0sF6{>sLbJVMBwm)ScGVnZr4;JR<`Qa-@uM zBlv!yu_g|TS*X50C_!Oh5Vd&}TfVHd_;;W<37)kS=}dJMlqPo{~p3G15av@G%p{cJ0~QLRrGTR>jqBP=HQh~|ZB0#dXG0sUqqY~dvVV*(f5&FlXxU>n$*G`vK@lN`eB1cC>eDV2 z|Gr=BB%Ya%k$LP?fH;S45HrbXW;2|L=qFnjjeXBD?OsS1d({1}cpZ6N{)|_1{yJid zQf{_X?5EMzdyIA+RI-}bqwDgNWyGmRjBdmyMr>9IM?5mbVv=>lAAs+5>33Sf;b#tI z>0bzbC16tzdt%s9!j23!G?MBVvi;J^cq`Pea~3z%bsl;MnAjHrDqdz zdEcdKc`4U!=)9N1&eVKDSMCYh0Leq-8XuSm`#HVQ;dQL<%a5#8dJwOBL~Wt!y%zHF zBew_tT=%fKO~^VrhL$}Pd%?O>O8i5EJ_;6zqtcVJXw{{c&O4`w@w{CSLHCv#lH)Rb z?=+q!G&XjtHMdhf0WZ67Wq} zeOsYNT{qeMnosLqt)gn@*3dJ9349)@&^<|yzQp?e>~7?CO1CqonEqvHw6VswPK_^= zQrnKIj*jBAadAtC){bDk{A4Qqt~+7 zODm~;s`w9Zea|w*bBwZXd}74p(~U=l7#x4cAArA6_*&HUGY8))@GF5i1@^?SrGykUm0u)!ut^lD)~Tep58wMm%BZHUlDrbZcI>}7lwdooa~V9j5v!d}VxM#HB>vh)0H8Cd40*W%%5cW#rPwmh`Ox zn|jz2!3MApsCFlzFVU^(P^ArM{!X(VS7Vdf zFJ%Y%iha|B`(YHjG=)a=e8lff4J!xg?_G>sJCvk+jZHJ!=geDws-RyueE>+FAx-5cOV#eM$Op>02|mO3a4Ht}NblU&w? z&SY$*eG?AyKG$pP-AyePt5^#$4_+yq$**gX|6`uFF%aL!1*7ijU9b5-u^vy_(T&5q z7~Q1s$pI8IV=Y_!v#j3#A0GIYGMW z^K=`gT?0`4D45Y*yrm$%lAkG8U6I28V%xuHR|UB1s=x@U0>s z8edQvdt%s9!j4SBpgRq{QCWws1kQ<3$N7MS|DEFP$x+2Hvjw$5Z3f0?qu`a2e9mn~d&P`Iq925J-Yu0<{ zqHTBHPNQ1N@rf52I9>W2#a^rnXKz2GtA_j3pfBY@sQn`$Gi)$zuZ5wlx$#`hzP2_EJQdBaW5u+hv@+A2hU^yK zmG)m6|1*x}1*%a(W+`>0<6g9DOg=MeJb}|J8~BD2#)xlNw+;6>tZd8{C+qRzzND`d zIwIM(Zm>N1!u#7We*ud%c}lfzZKwRY@w`oyz#RT=9JMM-T{ZqoFwR@e=8@CbFcwk5 zg(544&>W3E5#rPXC+8f*r4)7jUm1S@ey8Dc7g*_Q5oP!nf^QX+kzh{@TS{q12Kxc% zPD5`Lc<2j3R|4n6I3oiLNy>$=9}p1n+ahSHi^T)6F1cFWfXcbIqS_JtI8BvdBYFFJ zZQiiV%cCgEXOyaPQZPHUzKwIu+v1#f#k)vqZC!vKZ`AnOs@|Fnth9p7b}Yf$J4ffT ze5dmK+1YC|l=Hf{a--g}q&t>2rs1r-X%{x>%SPV*_EB->9n>j;4V*lMeB&6|wRxe8 zKm7P^&a+H$@$Js-J)1&jw|0(845QgSYAJM=3pmkQO-z~eLLT(yaT5hTETjn){ydJ? z2aaFNd2Mo#|C>g4Qlp3F@*#C(cma*n|+LzuWUD> z8j+=VxmiSW&il;5##AQToLcmDp{8aHX?$&!YG4H+OL{+AtPlNd>48 zTyoaPq_pxn^n35@C)N7{(VV|cvm`p@dx|_&A9x$!&3`0N-}Yi%V%77R0tfE@N)IyH zbJ~tqD^c{D)y&^BoQT{8AKU;K&x=n2p-hWCXtqcadi#GVl;qI4Tl248i zFAp77gXf6Rjrhcf%_`zjYR(}3#3TEAT<|+BVeqvGKXdT01^+_uts?zOz@8X%(vA%F z1JIp@-l%L3T?w2M%eu_vo)~K(dp(qoI8hG?&Ca~wDxlLe_=~`Dka62wWge3o!CHC{gHF{Seuj*GG*}Z|6-pq{?v$bO5f$x*)*tFAWiTRdzO69wH(*$ zgP3nFPB=opySO0Fw@vOxo3~A+b+en$xl;@2OgHguPA-+n z@gEj9CeK$Jcv)Vj3%4oK_|yC)o>d7){_&>1{Ci{nVR~(YGn!JfR0m?m8?p5Tr#UZP zBK*Vtx?`qXANY88bL91JA6M}9PW}#*c|q*CCq)aHuU9=G_W?c9=|C^aoG;E8hCNzA z&xSSEoVV{|?me~sumq*}t!Xaqnd#`_Y+5n#{7byKg)Lum`)wFa^B>RP&n05$*f_&O;>r>)8~yi^|~~&jCf=+787D{ApU^#y$+wd(%-1`Gbh_2hfTe-Cx#sv zY-nIVAa$oXjYeMxx)L}i#u*vT2XuL$jJ1%2?RM>}ybE=6Ea;ozG7Zk(rEqi{{phxb z(;W?*PZc_DSK7$7!|$!6vmT2T{=QT6_p6d_c8vbLk?i zG@*O@RjOMrZFyZa!ipVDPvvxV9$He#tVnt`^CnHb?#c1*J6z{^Y@jvGD6x~8l@hu@ z7!%Tv>^{o41VsGic!<5qy4z-|dRgL~I>Sn-TLd3ujvC)y!)j>PgdIml@IKJ;-sU*= zW2QDH*sNnhUee`Z8Kocl8U9LnNpumVtlR&#>Sg}?NB%B|Q!jCd&5F2`h{-2mx;cEW z!|ya;@HeW9`@0N#V%So`j!fFnKzAB?qtKa?`a(D-#u*vfLpF*Wdt&Sduoglc=~Ko$ zm!-W1EA_n5H+XqMKwlPR7q5DgBi6h-AGTI#&SVs%No`tC-j1nM_RMZ_KkUPtyWCdp zr{JuKWm{k*zoujlq4R!7eTCW1j(2)rx&_S(GG(iZoS`oT#X3i}NtiNI>1RT#ag<#t zm6d+ch4U`jPpo57t7ni!gYycVY+Guq5$~_bidQTrK-ee7#WW@7yD3!f_Eb7Fc^O36-Y&;+r8#J>wZZ*!*ZRJ zlP>lbGZL!OD8Ez6^CjE0vpq}2a+3cUTV~7`4jVOZD$iFODo?%UpJ4CpkEklHyu$m> zc|ZT3c5$ZR99CQ7&#Zow#XorTn-^_Dg*x@7w0*(!;e8mV`MFZe3BYf;zn&pZJz{jD zj@Yb-ONlZPVll}$G^k78>+rb?f1~iV2tRZ1u_gTr!LNj_O+DJdh6eTnC_`@)xzrcJ zIWf-2a6SOpDC9!e6JtLh*Fty)D!$Z@zwn`38}`9nd=D;q?04!L8N&Hqt|rbF96yQm zc154F{2n$rqg~sf|Lwze{qx@*zkGtsjf89zG)Q&@#F|UgG_u`ZJ$XN;d{(MFg0HIX z-K|2YsX-hjX^^+yoy1P24JON6aW)XWKE9o)UJZk(PPZKtF$m`; z{x06@E&19aWx7@u<%DLx4c^T)&rE#VcdCgwtDTSK-=%<)cmC-gbah@{<=5CJ$BQn# z^kpfx#PORSsr zSJBEU+b$Wd@J73h1-C1FE;}w_^}h60>UWlk@5qedV$U7gs4(0AvM8}#=IXWQo^n1; zW9@0Li|8}(oD*4|5poX1sYi@%^znCWR={+V^t}$ByMRkyi@G$rW!O@}j!fFnKzCYK zZxpmRC&n2W%8-phE+qHFSPKD;^tD-Um%+wboBK+4-P!!o+PS5+L85X$^1~id%uA_f zq=`R!`7yJoF`o09Fn1mE&8S56>kVOnp#z+&9=NLXQTar!>fS3c@8tGa!{PT5?bYdf z&QsLbaMjB@LMDXqBzb&QC0!J2bkvc?T9l{d+h?evS6Am_DVLDKjx3q4`Z`XW_dZ^m z$krar=6o+6(C9R~HB#4Yya%&UC%icQx$TVyx zeH$X~!`?weKaMn3FD%xR_g|=Q73$hqtfdWq?_m};yHqx3r*OV6&K+cpmj?23Er%E~ z%}!M@gGIbv$|ass%AZoX+8yTgA~sIcKDjtO=$pmgn>rmsX}n8o8sZT~7rxcuG&??u zy=(T}dflVuS^2&Zk=!Dof@1UE*-lJL%N^rT8X)?;6Agj*q5n zj|w!`us$^m6Yt}S9uqj<)h?Das9qK4icRx*+fzrb(y`b5sqv4VWLtJMhc)Xi&N|We z`JI_G#ov^iPBr1}Q;fxP031@A1>y`GWyGmRE)(KXqKtTCG8U7>!S_1mEa`Js`dWmK zE$LrK(k>{cT|evE&_Hig(jdn;&_NZBa%#9}`HjOJ`z^A6bf2$^c%%nO|UcSxw- zxQpyeA>U=%HRWQi?ZaA zVhz;GPRz;U!UAc{(rD^)BSC2cUVocwf({`-8)^UiCpU6)(f`6z`~UjOuF)4m%}&>L}mlFoPOWIf&z zzH76Q`uY7*%-pg*$2S@iK^Oe~_x>2|Za?2d6^))L>*|+G%(10fi*>hI#VTsY95eC<0krc1WD(p)7SJ za7KnYWTTJ^p^W_ia(oBy4md{yDsm*ib(4Ql!8_t%Kgu~}s9sk&nmxXGpB3%ztc(Nn zKTWiA%IX;H`Ni<5M|fXXzx|}^uA7uL*L;U6W!?V%E&ID1dKyl?!Qwf*Ft-$cZi;LZ z_u=(dA;XsKvdqNVJb!F=MaC#`%i1OBJ#%g5)n*n^4VlyY7JMwvhG^ws_am<<&vomy zUFqfl?bwSBz2|(=?>P0y6LOX0s#;m=+$9NYVEa<4qs4?>--g=yd#UE@UbEI){Me@+ zJNfvsR_~)^Hyhs9oyRAYJ~8IFrAd1CXx{FV z;t!xN_??zMcTtC*Ir!Lue-; z{Qzi1@^`6}x4^}~}IRT@(0Klw?lGv>non?)Y)R5==j}G9aM5LIhcrg zrNmR5*8s=Y?BM9;&j`a9QI|Xl{$LL1YV9z(az=e+kENlhy83-+wJuF${tGZ zo1^lL?VyZXqp7Hr<#oxtYJwr>x97`rI%l<%Jk3AydaAr@JhKcq4x_N!rQ*$Awdpb+V`Ifq`;&gdUGdQo@3c=1EC6k%g z9B~ijSdb6g*ND>|sM_v7$9#UFct6{HG@|(OUhGnaCFiNjLzYF1ZtXe^G5HXWOyXq> z4qe~t5(htXC`w63^ehSVgC2i+K>dKLojR_`-SG*qpk3%5t${i>XV- z359>nJh9&azLVb_4qrDjg%y9>minF$vP0n7n$=L+pf1yB?@rD=?7wi_nfBuQ0^Gw_ z>HMCCjaRE*ZmvdVpJq_ggiRcFJ}!r~HeSz6r!S%^g?jV0zA?Ay@Zp8LobznGPG8z* zPWCy=b~lUC_|1`Z-Q2b33|c(hNnL5qXBNFKjwS|HV{@i2VL_Hp`1fm%pEEf`=ji?4 zXus4~_|UFauO#cUthT~?$-NTHe5)}fUlI3Lw=Y>nOg_XTLo6o15Q76be6PdrG<@#D z-za=7!p|IVGQqEeu1!7c$Y4Vw@z5KUW#~%aoET?hI3Iv)RN`bVzXSd1db^w`RfNBH zBMY3ReUC!ue4&*JuUyL%KE}_g<{Wlxk4ArAvv&iZ2E()HNoO%nX3p3`%?CL(^Q*U4 zwfWLE-o{K5Uk-FRaj}%SKs>*w|D7uQ@x`eNpPkl^t_FYP|3$_YF=y-I2Ob>47MYLb z{K?;ljy(~3PnA~-4l@ot$MZc^*D3$u$VaPaVQDYw*Y}MA)8#44nJ+f7W14ucCC`05{`*W3c;7$YRZ07IUgK*qg*W2V%NX5=&5F2`pqKH;WE>jZ z_ydw3d@aJy9P02d1m7wW2Ycdvu%*=4kx7`;orcaF^o0PI=fv`SK+@{wcn69V?dLLY z@^)6^$t1dxrG4hYH8yqLAq{yS4|Z6pFI%7B-tx>6=?b zoftiZOlEZb!`I(s;a08vewAT<3g7!5s zSPE8`rAt*yDfhJ0NA0>5ymZs1Hk)W&)@Zh1+Xi;qWFhPR@`ZA5dhHbJe~X~Uw7ykg z1>Ssykcs+)bfOEjE$Pf1aSi5+uPA-?4cA{QZu~NkZI-L{Z4pZp#rPMgWn4(AOo`37leQ0ZU``96%xO7d(g4Sp6Ag5!UIm~yRkV!P%+>n}9{XkCBm$K|jzxg;U z4~n@$@|3yJbTLO&RtX+!N|j@mCt53GpFOlAFQ>l@qLgPDtj*O4cd)ez`Z`D&PY-f-`bk0w|~UKqoOD%uNAxITR`CvvgHEV{_^AZ({y$s z^IkKDo?ghaYNp(4$-mKQj3iqA8!!{JI%=#A->(=Qp)yc%gU&}kSBs*Xwdwv^xmaLU49r_a^_jIT4 z{9eTTzRzMbI>yz{ucqafFVf&f_c{Dr{i2k(;3i$j7VCLmpT4Y13GLqCS&i1p*cbjR zuV$;7^15l_D{^l7l#g%ZRq<@CnQzN+MZCqj>c-8XtXA>}THI?khi~iaL)Qjr-&=W| z$d4@)`*B^EEISPr&+w^%73p=`T+VB?#!o_88+X;bLVKxv>j8hBeeoPtzT1p{w-T-o zqk3chU9(+STa>aYp5%N~O&;q#E9kSICO)x@&5D?OfFTwW;?SV{UkSbzrJp(Y7Xl6Z zO2DRG(!-8Smj74I&&V|A15!2$I_wAJS_tpJ^}6HqWUA(kY}l1Cr#bJeL1lSf?^Abj zbXrA+&m|~*uBjS8Q~Tam>T(Qw=MCleuTV*GPrA+)^3YmutrS|Xp-0)f^ixWGlt&9* z4sT|rDz?*vUt{ExJTi_wtm;;4)1UuC^x!nJv~mvs`Mt6RC71fDN3*=4cATKeFSmq! zJ#t5Bx8J` zy}l_Gm_17wKl+gARxNY&Qm;-5F7&2Yn*kqm*;tj-HiP4OZPTu&W!ooj%5mH&djs}i zh1mD6IS#p&{&%@~LwinV*VT^ayJ8RW>kpb6O3mX$*97Y{9<}VDKw}9QMSprGy7wxphMU2jZ1yyhwWOi@Mia^K-X2dPPD z3VR(Q;(Y4zLAhqaziU51Bc8`+L$b-yGmjO|QBnUB4d{$rV|C=3!BoD7my0LLI(U;8y}R z^{^*~EhX&8U_%2ubf=*=in3N;NE35Rb4CW3Zo+;*u7&Up?0%$)?XmqL+nD)*(^hwC z#mA%D2IU6NF0n?fwfnyn2ggw1?Jg|($Oif_w<47-eVE;Dl|Z@gDlvx}QB<;bCpPoK zS-t+5y}9%1ZHGDk!M35apnEy0o>EwkM)rH&FH~*%+f8;jic48h|(oXyz1&snIe*X@= z>yy`wyzPxkPRS;3=<=79?Ak{0tX=LX?2m7J-bfiXZ`nTg&J^S-)`VztrN&$G`zgL* zsIO0RphLBO^YNfQdU`SyI($jBDt_w!?02EBHI04VMQMX}vl5%r?3Ny6*TY;@%T>D{ z%R8}O;p5P62`?KSwP6#tJfixmH}ZRu{_2ZzZ@>%jiIK}BIEzKeeg~ZvJpx9cfPjQ;6L}v-{p+E&;{tQ{0Q~G9Icw!{25F0YpnEBbe510 zNgijjHS5f-{%rP-r&OtrkwPc&l|4n<=JPE%uI0gJ$~_<2LMw0B>?GFf*2Oz;I{PX6 z|KYXo>0OoS++YrW_EW4U^6PC-_VcpOVP0+;-Enai1=z&AbG*-8?4H8&(PS~V0$=S= zb&j7JyjdCB%%kEQ0C=SBb$yh!?vc54tVI=yo%MZR5+pz)8zJw=@Q{}rPfv6v8t1~E9$PS^Lk^fwAPd~6|? zN!rxImQvc0!G1vE<(vs!37iw-j11=kkc~nvggr6#16T{m?|@I9i0?f1%5`>m%4hQE ze1T10Exs|5@0)SU)DstJ?K_qyL@0fx#kN(K=;o-rcMXaaPzU~;sK^FwwzlH5F|$hP zUCYY(N0j>G)UOIZ^!fZoK?R3;r3A6BM7iakwPIa2bKb{sXEAS995s>MH-EwTtol-# z-iH=e-}UjOw~n4f_1`GZ4R|2c=qq?07(9~qHEeYmO86f(E;nAp5s#%l8Is)x{$316vyM%K3;zmcQ<${VzKfogepO{Ve zdB~=z&hAw)Sbfy>TTYRnFw&RCa9% zbxe9mO0h|DE?#0a$1Ss+&+}_7*HFL5;j}hf{LkQF;-!5KQ9pHWg(~%8sxk(^T3W6C z=a~oG=*bUNoX1 zXoE*J-p791(o~_n_BjMi3)lTLc26kpvrw+~TBc9Pyb+5@mJxpdW%!+j&s~&teJ!Gn1iuo}rXIGG(vA%F z1JIoYjjqld&WUkG2Dq+l6#D_Jh438!ZhYWyJvkC!NFyWrs!FuTq(gz}e5@9i#q(O@ zPA;$SaI3?Dj|;!PHMaKF=ua$S?Wg%OgE8f)>7z+xFkINogV)?+HI(aJdf0)=cW=zr zRgNITX3r>j@eBnA+EVt;s^mcNjCZ#!qJH;!8^>1~CFU*AAQ>5obKtUe<*0I&E!lYP zRlY^Q&n{|4W)6jT+j&9dRa=IKv2MQADRuV{)vp?2JO$1yVAcrEe8on|ev{xBUS0g@BXi#3(~H3b_!4x%f@YuS=A8>r?ialN%lShGD>G`c-s4sklq z0uOmUw}SBRHlVDlx~~0ER)LnX4y!yku0x-0JjeCxrU$tdG;n{H3(d7U0J?ZtK9iix z$2fc6233z7I~w|Eu5vA3zZRg3ChlG~w<|k#qzb7=-{Ww% z8)Nj^bliVQul~7^*sp`%v}A8;H0}_E9q?CZsvGIXT+^IagneIUHuj*^V9RG!iP8D zx@j5xRvn1f%JqReqtU65G&Ls(`;P}UBU$aN=e#ego=>UrodDI-FrmwQ8a|myPelqJor&zmAtx7E^^0f9l^}$SEe)rb8 zkTSl~BljtMEPOh7ApJBcaeiyKeF@3x|?<6ltS^v24yTwd0jK9*dq@K`@IgerfiOBdaKahl#mgs#9C zjqF+Xne|8%^SGhL=R3yP=&q2(-kxSm1swxf-CET++&1ch(%1e_#%jx}=2|i7Io2g#+31KsT z&L!)s+Id#CNB-IQK5ygSM8p#+aC9!4v86DF0bjzi9xtERxQcFmj$y-1rZSI>%M?6l zK7K#VY`wH+VY@yFc~8YX)7gu216b|T73hoYVOB!teU0mFq;AqOiSG3@Rj(Y>jh3XG zRsOGsekW+r8RC7MFRYz6&=*ph7Gcb3lm*9kXg7iOb$zJ%@>S3Sf9+uv@>(nQfuD|U zrIlTVDdRoYxhCsgQS<(nX~ZgP{2WefE#|DxCLfjl!2|J$ks~f8V)7v#8S02bgBTo$ zKOlXt!{_eb{zl%r_5%_pbLdLoocQlEGKt4p2;Tv`1DBG;oO#+w z$k##h=6y1m6uijCy<*~QDs{r0&0eehu3K!6IbOnKeqz53T~g)gwSPK(vnGaLA|J!&&)>*wZyP;)_JziD|CW$5NQR@rl+uOiE+j;r+ylV=~J$@P6wd*{48z-LijNfe)oNHk*_l`d)<~n2ZV|2*q4NGp_ zN}-80d##kq`&MNSZk$%e(D%#|_IC7Y>XI$K4G+f$Qr#nCSUq)RPSZYNHAR_yQ7vtK zm%8U_&#Y>WDC?Z=DBeX}P7if)&*$T7*9|4-*82$|w_mH0N_3aUe$&-(nsAN;~C z7&9v>_ECpNmiy0lSe35b6Brf~^1`F$`hLc$yWZC0@4HRJb!g+D|BPO%qf!TqGp2S` zYw6KrEFFz1suQn^IRW*-^^M7WXn(48Zm~k!v+)+C44O-u#hS=z&0SU1#+TWwmHKN> z(8&12C?h7Hj7NqzG=S;GAArwY>2Fl{nbTS4>iCsFJDH^Yz~6eKI471gkc~nvggW*E z$nhP(JFx6Vs3J$Qu6UY-=1-#mZko03@%7|WqX8{RT+Ifah+~yw#aa$+k<6Q{R~5JR zQtpR&su=~jk5u~fi|tH}s#m0=hf-MiHd-0nsm5NkVA)OfZ0uUP`uh@H?6aA=JUzku zd<(E2&zmazHlKdX%dwS&4dIssJ^A?2_p6;16gt$MPm8@N$^~l)o!jZ9W^3Ydi1WAV zR}Qhro!+!@TxGp>`JYFVO}XaOVXt^kc9%D%mNiZIH3SS1{)gi4ucDt%VmYi=Oct&4 zyUEK%ns?#PYnMNBN;vz1DSSwU_0Re|WU6D&|kejQ{PcGppWPh0jxi3z|4%+WqMK+QP1B)lV_+O;2yH zvVJjz9zDCub~$#YvA4AI*>#f>?DqE||BS(`uc7+c$?kM%qKW1|(>}w9PprKfQTjU; z6QJ-rja=8?sPwS~I1+5?r9Cn1$N+==0CL?VbtP~{1~_D+$gwBJegI{eBDO7Ysl8`_ z0$Wf|=mG*?B;z^t(a1niwwgJGjE4%{HPn&1oxIL*9%Z(%gM0FMd1AL0Y=+4oYBgWT z=IX{ea=2)r8B0NHO5C-B-aR#>IaRB3oXr989hEf5S02)SXHY*=;j3yuCmYIl=t>qX zT2tliY2>vnn3CLEQ2xPV%J@MmQ_U!`r$D_%CGB(5>y>uS18!R3TY6&z?C(SlF68;B zI30Q+_PAd)KK^uRMj2bE-X4EXZa4CsiyvT=ThNV+zG?Tk=MN^RjytFEamHTtBEKQ} zYbEgNY`2)d@0;wBIGsZkv3{>>q204}JRsKMHFLI7OQTIH`~I=a#Y)`cny)WX@RONz zx0?1`G<$Vhy7T%aAFJi|F!~&$zYa|>zU@@V;mkkxZ(SnpiS{`T`CU}G-d|Q`6dQYt z>O2?ru<}0mM(ts*S2$97a&?;iXb1b|A=VR!Q;!(ky77r+OgY*pb1825{(3LvK`9XHK3I%ku%y>F$Y9$2)N0(Ng8R0N5(SbXERPp&S3> z)j)oaqUM%i-t#nP2AX7d=PfPwxk_Va8&Jsz{kKi{gqeg|yfk_VL1)ru}ZoT2cp@ZkXY`-(9>KWs`9t=6$R>lbl+ zcBWYGrS31H?s#sI(l2-+-AxyA2b3q(a8Y;~@62N@rbjSq--dkt7_qh*eaw%b&5h5n zLgiL)x~$q_jW9%`2Yr72Th=(WG{tod;C*ycYp=t8XjQVRA!LpBUT9>-npA$~VfOIW zC;qJeC|!U)U#U+Kt!#K-H_HE1E$&l<>Zbj>w>t8yHESMGwkdf0uXG_IkyxG?S(BS3TvQeI=+So8ZSoh|;b*B1?5th$%xdYT%-Qz`Y$m%w zLT3(rT2CCJ+S*5~o&VPu{tNEiOnftBoO;A3{+kZ5m;jS8IR3`N*CPDP!N(T-3;m6U zEv2+01Ma_crz`zAC&u}Jl#OCfjQs%CLXu}48?C)_l}h4UBXKhK8`FnLmGkE zi0$E`@VR+yH!BdTiBFlXm3@HE|4Li;cBbr;sjQWIH-!giAN0A&>%POpwVXQKhrSjU zx=SsN=*JTU{g^dE#yDOzLXVHzdhs5EMy8Hk-!q>_M%3zkJk{A9p>w_dE~mZg(Tb5R zlxG!fJR7LkliA`d74^+IF-jjv=Y;KEfuGZL_<-u)Y*tuLSIgfk%Q3jVNhkv>LrptA?ybv z40(5_V3$H39%NWEm;4u3;keZy!berZk6#tq4^x7u$CLXii;W%Vokaj`&%R4zA5N!l z+n*^k$Bla`Nv%VD{3UMCiXT74KE)CeLRLr#+``JE6 zl>ROSJ|x56Vyy_=w#VX`9AkIWdE;f>J%PW=(ODarX|%Ymnh!EqM@Kd1R@ zm{CXCd3Gv2{PS%X)L{Yzo_(X}Ek?Y&tBR?4;GZ^yLN}-?mS3x~{Z-;0+$&9-+3j-{ za(Ky88hzfv$&_wUM2`-!SrL~KG5HXW46&FHhenJ@`wO4D+M0d|ekEX24|`(RQo@c* zQb2baWnG;))Nw{8&j%nE!k$=nKY(}OdOdxAszwj4kx7pBzW|^AmG1bN^Dz|aWk8-| z)YQei6Z<+tyH-HEMTeg8vE6K*$#ygs|HrwFV&1z~Qp_QpY(yNT$bCZY?zG#B9Z4QV zFALdG(u{+=kLI+H9Xx%4^D68UsEJQhO8MXa9$)We`g8P)vJZH>fNopq?qpvYG|_<8 zTWfvNq5bi6Q{{f5zQ5-NUcdUnL>+ijE5Dh{zv=!iHQk0Viz-z)&D0ON^Eg^-&UVzM z^dw;=y@?N_l$uqjTdvkd8~9r4+WdL0dR@D3@c!739NYTR)*mzJNN8haJRjba=jCVD z&DE({=h)Zn#d-Z{scqCmBbx%hX!G0r8VihB$L=qzpK7^Y&S)#s^9LCUt%M=3_E>v=OwJT#SDV-4e9A5tYu;0T;?yHP zF=DgIwunb2V=)1S_yh304!_f=OMjy%!^f8NFNC^Gu%(n`*bmIxsMQsb`a<%Y80Q0^ zgp--}8L`~HQgJ9)ju??^W1tKolM|3Mv}`nv=56l;(A$39ZZ>_N(v#gVI}x{};4>rwv8- zlKbQs_9MoJri{x}9XS+6zvCl_U2n)1-_eg}CEFgcpF+KwiF28Q&g*zT8S&SdYthGS zo!w>C@GJ}3d)ZvMSLH3W@9fY%U8&&nINs()w+5tYE3O4`>cJoJi4mI>aVcegh{c35 zVsIdb?{!_j)9|$@{mj9?ko2trn|jz2qYdoHU_&GA2c+I8^o4*&!WkKI$VQRNJ+WL1 z;T^DPEBxjS7(9T}fNre48x8y{=DMi{jtW2EvXYy#(1r^sF7p;mn~QX^T*(4 z#3F&O}FgTqjM{Gnx6C- z%VCf68PAiZ)Fb;}LVo!6))_r|$q#wjw^n>?wfAXaeHyvbt^^CJzoj_c?pcIRAG0E# z62k7uudKMP9K%7Je$&rhH2lhN7U(48b&_YTTJM#<(mYEG+L0b?!rQ5O@132&T2435 z;ww3>*hp8V{~LEQc$qjGSWsymZ~LUKF>g~b-h*})HdH%gH)M!YFZ)JpR+JHw5Anzl ziwSXP0F&_t;CEX3+yxwd=HO!sW%yQsUkTXM!=4znl&~WM4)z1cp*IShIh3I*fpcP< zk>Pv*vQfbQEA=o4b=jw}r|7j~A)7vRwyMNF@n6TE=ucH$#XGKX5YO934SGYVt%ZEO z<=I&Z&#TLKiL!lt)oM5Wd0zIpY4)FI{14!s{|)8-02CA{Y++9Q66fEe9$Bl$Tf|YH zgkls?N>@H5p7BKk{(G(mul%$b6#L4LwF?k?-J?(Ai4MG^jYEXYqT%iqe7rV8i}F1F z=X~W_(7w-Dao^|Pe?w+BOLE-Dn3a6IRlb$t^${oX_1b)SaglAyw@{yb6~*D>oVM^h zd)gz`a?%xr)@6Q0UbgiaOM#}P$oNS&R_3+~joR3P8g|H}Vb6{EI{$UCQO=K}W0-fx zZ7TH=u};L;bkh+fX0+zqXwD)&zQ}+ecJx6Rsy0+>Q=;3Ze#~d-z82wU4nDTvUkJWc;8y}T*b~E+ z5_V+Lh6cLRQg0MAGU1F2=L1qUiajy*16T{;JAij!b@NgS z#9XEeLplGMnCDThdic+~qbVOgEo4l9BXzWR!F<+^=KUOX$>!J6qlYiA@7(xa^=0xV z4NNdv3t3uGgo&eo2}*dwZQX71n;w zfkQGbYR6JT3RAf68|TckV(t?d1B2o%n=Ajtuj_58`D%0esfktQ5IT)k{fVjh=!%d7 zeaaC3aj|cw$g=ArHsYXI$6q}%k42R7<6|o@po_xKc$=`(^srbl=tEGB;{SY%>vw-(=Vd`G~tXH1(gRP9d$6?0aO_cG;d$7ER#;1vD zr2@~Mx00UR8%h&~_8{|=IA+*7h5~P9u~~m&c*g$^e|)VIop~RUFX9s;HY?&%A|@Z= zks%h7#3TLye6PdrG<@#D-zadpe&*m?Mf#P1Juz%4L8EI!qpLRxU5UTXi9xTsC&qUG z?|@~^DSGD%p!+-BYF(5fJcWD!xX96d?0eN)w8pnA9WG;|T#Ia%=^=ETW&KwcC421GCLhnUPG;)VG8U=;9X}&21Rh!C1~(n6iyNXVEZE%LsrzOf7b(^>1m z=IYT=U--E0#apU3yB%l83aw&?ha1wnTBB&-HXFKTzn_&ZzLZUQU6=Raw^d(<@JFAQ z9IneB<in*n!spZ0vI?KI4K=pHMO-PXRW};dBbhB2qMakBZJ0}! zGvYbz#(kf7`|7B& z*b@T}J2KeNpbXvV5r6cRQHQPs&WUkGhPq7H6JtLh*Fty)REvh_eHTtnC`6Bzs%dVU zSxVb&R@&=z8dpdYQx|JTf0u-`2MW!qd78KYuKIhcnf-gvQ;hPD`3aRQ~E|W4Y2T`&AmRU@-pox_H7`iO|=rw!qIReR%=l|TGP-<8NcN7I;aNy z7}1haLdAMc!jZS@Cf2eO>k4^MvE#dFPklRP^VpP%?lqv%b(2-GPKK&UeovISqn!Og z-p5nFTP((NDJ6VDzwI2v%M;3nP|_Z; z_U`^In{_-MtK5$Zi-i1O|C@^%zbT?D`=YYRwDbH&Rcn{cs%v@6$*azA-tM!JR+jy+ zsaP8v$nQdkQ!it5g9qYLB1b$j03qG)VJ5zQo@c**Mh$}i8oxO=X>pSFy+s>SRWoY%wg_39Yxn12#l#w< zZG-#lytCLV`+1$A6PoyGul;|h{4_DoEXI`NV+gT%$6hrPdt~tIQ|BoktL0ih$_OgQ z)Pj?}To2Yp02;ot0V)OO{%7oX_N{H+eNK*vh-M91a~^fGyyou?1) z&XQqM2i|V5Q#G~P!i~+g{>J-vUa*OI*UP2I^VjG_;gafzCCkaDU8mnD1 zNS>xU(v&es*x$Kpy{-Qz&6XD0@1RZg75emAtmFUYgL=ZIVf69mMfS+;0p}U*Dc(WE z=tg{E#AXE@V)7v#8DcS^{=X7@E$aH2lXTpsUSm%z>&Rh00NrWmjiL^HA>=qGMlR0> zAQwU%Nv?%ZPHjHf#d63-j{9}B3-fJYsk!g^d)N=B`rjjyadfW7!uk7gb?Ltah2sD-j%LBnfTAU z`#xv4YW8@sW)4oZ{s(_%qK`7R`ZdKIfU*Cb5^HYbbfY?GpTWvkHc(b(uxdg(u~**k z@hWW!c+JPO;)b|);Nz@%&YE@W$SfNiWd-Vq{qc^$XP8Cr5T$R>*{v_9-Z#{QA{PeJ zsw@5quAX}nUjFR1o1CV&DE~#bLA$8_z5=|y-O`EEO#4(0=#iwZ3Y5rnNW1XxXBg^tMF`ZJeLYzJ;G8^8-zl|7FkjETs%yhs_e% z)N`aRrjcU}Hao#_%Q}m9?%sw%${2=Ji&DyHH)i<-W|4H4HZ4BSvfYAIPGMsI(=KN| zr|Viq+b8H{`q;It(~JdI*<+)93cnN4{(57Qw7~%p3JqXKRK?i!Ut%vQX?pDz`qi(4 z#QA`PPtf?wK>11c_SCtNJPZ5AreB;(@na{kV`GGTz$SPZEogg(dNwqt z0c#A^#oOmHTa#UEWVQ2bSpk*8d*cQoB zo_rB|9iQv%^~QT#W3Ml1kVm_%Q`);H)Tb%O>#L?5^`-ZnwR7R{;3W#p=*dT!`MUvR zx2Qekt{$w%Yi-(QwyN?$b~!@-UFx%?8#Qwc=hybRttKA$g~>GUke$-DnnvDs!mJki zaHS8;vzw&AkMGMSulQY*{YJ>?>)sOQ8M@=I;aHxJt-Dze9az~#X{$?fdZqUH0(_Qp zY2|tryo%>-^5@*4YRhbB`7|Nx>~5{q>A2KJ{O4oG=-=y07xedw;C*g$3%0ZTEV}Pj zN1eGUobx`oeYdK)XB zmugya`1dEb$=7{4AHQ9j7wp_AvDcI_x-kaCr9>_h;?U^E;DFz0_}qoRQTSSvbTWru n36x<^j9k}_OzKWUXAU?e{W&9(b?k{z##-p?lN?Ri-{t=S8vFXf diff --git a/test/data/grm/test.grm.id b/test/data/grm/test.grm.id deleted file mode 100644 index 734f11d..0000000 --- a/test/data/grm/test.grm.id +++ /dev/null @@ -1,194 +0,0 @@ -1 11 -2 21 -3 31 -4 41 -5 51 -6 61 -7 71 -8 81 -9 91 -10 101 -11 111 -12 121 -13 131 -14 141 -15 151 -16 161 -17 171 -18 181 -19 191 -20 201 -21 211 -22 221 -23 231 -24 241 -25 251 -26 261 -27 271 -28 281 -29 291 -30 301 -31 311 -32 321 -33 331 -34 341 -35 351 -36 361 -37 371 -38 381 -39 391 -40 401 -41 411 -42 421 -43 431 -44 441 -45 451 -46 461 -47 471 -48 481 -49 491 -50 501 -51 511 -52 521 -53 531 -54 541 -55 551 -56 561 -57 571 -58 581 -59 591 -60 601 -61 611 -62 621 -63 631 -64 641 -65 651 -66 661 -67 671 -68 681 -69 691 -70 701 -71 711 -72 721 -73 731 -74 741 -75 751 -76 761 -77 771 -78 781 -79 791 -80 801 -81 811 -82 821 -83 831 -84 841 -85 851 -86 861 -87 871 -88 881 -89 891 -90 901 -91 911 -92 921 -93 931 -94 941 -95 951 -96 961 -97 971 -98 981 -99 991 -100 1001 -101 1011 -102 1021 -103 1031 -104 1041 -105 1051 -106 1061 -107 1071 -108 1081 -109 1091 -110 1101 -111 1111 -112 1121 -113 1131 -114 1141 -115 1151 -116 1161 -117 1171 -118 1181 -119 1191 -120 1201 -121 1211 -122 1221 -123 1231 -124 1241 -125 1251 -126 1261 -127 1271 -128 1281 -129 1291 -130 1301 -131 1311 -132 1321 -133 1331 -134 1341 -135 1351 -136 1361 -137 1371 -138 1381 -139 1391 -140 1401 -141 1411 -142 1421 -143 1431 -144 1441 -145 1451 -146 1461 -147 1471 -148 1481 -149 1491 -150 1501 -151 1511 -152 1521 -153 1531 -154 1541 -155 1551 -156 1561 -157 1571 -158 1581 -159 1591 -160 1601 -161 1611 -162 1621 -163 1631 -164 1641 -165 1651 -166 1661 -167 1671 -168 1681 -169 1691 -170 1701 -171 1711 -172 1721 -173 1731 -174 1741 -175 1751 -176 1761 -177 1771 -178 1781 -179 1791 -180 1801 -181 1811 -182 1821 -183 1831 -184 1841 -185 1851 -186 1861 -187 1871 -188 1881 -189 1891 -190 1901 -191 1911 -192 1921 -193 1931 -194 1941 diff --git a/test/runner.jl b/test/runner.jl index 7d9b96e..00387ec 100644 --- a/test/runner.jl +++ b/test/runner.jl @@ -155,14 +155,16 @@ end TMLE.write_json(estimandsfile, configuration) output = joinpath(tmpdir, "output.json") # Using the main entry point - main([ + copy!(ARGS, [ "tmle", datafile, "--estimands", estimandsfile, "--estimators=ose--glm", "--pvalue-threshold=1e-15", - "--json-output", output] + "--json-output", output + ] ) + TMLECLI.julia_main() # Essential results results_from_json = TMLE.read_json(output, use_mmap=false) @@ -241,7 +243,7 @@ end hdf5_output = joinpath(tmpdir, "output.hdf5") json_output = joinpath(tmpdir, "output.json") # Using the main entry point - main([ + copy!(ARGS, [ "tmle", datafile, "--estimands", estimandsfile, @@ -251,6 +253,7 @@ end "--hdf5-output", hdf5_output, "--jls-output", jls_output ]) + TMLECLI.julia_main() # JLS Output results = [] diff --git a/test/runtests.jl b/test/runtests.jl index 890f7d0..8f76868 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,6 @@ TESTDIR = joinpath(pkgdir(TMLECLI), "test") @test include(joinpath(TESTDIR, "outputs.jl")) @test include(joinpath(TESTDIR, "cache_managers.jl")) @test include(joinpath(TESTDIR, "utils.jl")) - @test include(joinpath(TESTDIR, "sieve_variance.jl")) @test include(joinpath(TESTDIR, "runner.jl")) @test include(joinpath(TESTDIR, "summary.jl")) @test include(joinpath(TESTDIR, "resampling.jl")) diff --git a/test/sieve_variance.jl b/test/sieve_variance.jl deleted file mode 100644 index 4677c66..0000000 --- a/test/sieve_variance.jl +++ /dev/null @@ -1,368 +0,0 @@ -module TestSievePlateau - -using Test -using DataFrames -using CSV -using JLD2 -using TMLE -using CategoricalArrays -using TMLECLI -using StableRNGs -using Distributions -using LogExpFunctions - -TESTDIR = joinpath(pkgdir(TMLECLI), "test") - -include(joinpath(TESTDIR, "testutils.jl")) - -function write_sieve_dataset(datafile, sample_ids) - rng = StableRNG(123) - n = size(sample_ids, 1) - # Confounders - W₁ = rand(rng, Uniform(), n) - W₂ = rand(rng, Uniform(), n) - # Covariates - C₁ = rand(rng, n) - # Treatment | Confounders - T₁ = rand(rng, Uniform(), n) .< logistic.(0.5sin.(W₁) .- 1.5W₂) - T₂ = rand(rng, Uniform(), n) .< logistic.(-3W₁ - 1.5W₂) - # target | Confounders, Covariates, Treatments - μ = 1 .+ 2W₁ .+ 3W₂ .- 4C₁.*T₁ .+ T₁ + T₂.*W₂.*T₁ - y₁ = μ .+ rand(rng, Normal(0, 0.01), n) - y₂ = rand(rng, Uniform(), n) .< logistic.(μ) - # Add some missingness - y₂ = vcat(missing, y₂[2:end]) - - dataset = DataFrame( - SAMPLE_ID = string.(sample_ids), - T1 = categorical(T₁), - T2 = categorical(T₂), - W1 = W₁, - W2 = W₂, - C1 = C₁, - ) - - dataset[!, "CONTINUOUS, OUTCOME"] = y₁ - dataset[!, "BINARY/OUTCOME"] = categorical(y₂) - dataset[!, "COUNT_OUTCOME"] = rand(rng, [1, 2, 3, 4], n) - - CSV.write(datafile, dataset) -end - -function build_tmle_output_file(dir, sample_ids, estimandfile, outprefix; - pvalue_threshold=1., - estimatorfile=joinpath(TESTDIR, "config", "tmle_ose_config.jl") - ) - datafile = joinpath(dir, "data.csv") - write_sieve_dataset(datafile, sample_ids) - outputs = TMLECLI.Outputs( - hdf5=joinpath(dir, string(outprefix, ".hdf5")), - ) - tmle(datafile; - estimands=estimandfile, - estimators=estimatorfile, - outputs=outputs, - pvalue_threshold=pvalue_threshold, - save_sample_ids=true - ) -end - -function basic_variance_implementation(matrix_distance, influence_curve, n_obs) - variance = 0.f0 - n_samples = size(influence_curve, 1) - for i in 1:n_samples - for j in 1:n_samples - variance += matrix_distance[i, j]*influence_curve[i]* influence_curve[j] - end - end - variance/n_obs -end - -function distance_vector_to_matrix!(matrix_distance, vector_distance, n_samples) - index = 1 - for i in 1:n_samples - for j in 1:i - # enforce indicator = 1 when i =j - if i == j - matrix_distance[i, j] = 1 - else - matrix_distance[i, j] = vector_distance[index] - matrix_distance[j, i] = vector_distance[index] - end - index += 1 - end - end -end - -function test_initial_output(output, expected_output) - # Metadata columns - for col in [:PARAMETER_TYPE, :TREATMENTS, :CASE, :CONTROL, :OUTCOME, :CONFOUNDERS, :COVARIATES] - for index in eachindex(output[!, col]) - if expected_output[index, col] === missing - @test expected_output[index, col] === output[index, col] - else - @test expected_output[index, col] == output[index, col] - end - end - end -end -@testset "Test readGRM" begin - prefix = joinpath(TESTDIR, "data", "grm", "test.grm") - GRM, ids = TMLECLI.readGRM(prefix) - @test eltype(ids.SAMPLE_ID) == String - @test size(GRM, 1) == 18915 - @test size(ids, 1) == 194 -end - -@testset "Test build_work_list" begin - grm_ids = TMLECLI.GRMIDs(joinpath(TESTDIR, "data", "grm", "test.grm.id")) - tmpdir = mktempdir() - configuration = statistical_estimands_only_config() - - # CASE_1: pval = 1. - # Simulate multiple runs that occured - config_1 = TMLE.Configuration(estimands=configuration.estimands[1:3]) - estimandsfile_1 = joinpath(tmpdir, "configuration_1.json") - TMLE.write_json(estimandsfile_1, config_1) - build_tmle_output_file(tmpdir, grm_ids.SAMPLE_ID, estimandsfile_1, "tmle_output_1") - - config_2 = TMLE.Configuration(estimands=configuration.estimands[4:end]) - estimandsfile_2 = joinpath(tmpdir, "configuration_2.json") - TMLE.write_json(estimandsfile_2, config_2) - build_tmle_output_file(tmpdir, grm_ids.SAMPLE_ID, estimandsfile_2, "tmle_output_2") - - results, influence_curves, n_obs = TMLECLI.build_work_list(joinpath(tmpdir, "tmle_output"), grm_ids) - # Check n_obs - @test n_obs == [194, 194, 194, 193, 193, 194] - # Check influence curves - expected_influence_curves = [size(r.IC, 1) == 194 ? r.IC : vcat(0, r.IC) for r in results] - for rowindex in 1:6 - @test convert(Vector{Float32}, expected_influence_curves[rowindex]) == influence_curves[rowindex, :] - end - # Check results - all(x isa TMLE.TMLEstimate for x in results) - all(size(x.IC, 1) > 0 for x in results) - - # CASE_2: pval = 0.1 - tmpdir = mktempdir() - pvalue_threshold = 0.1 - estimandsfile = joinpath(tmpdir, "configuration.json") - TMLE.write_json(estimandsfile, configuration) - build_tmle_output_file(tmpdir, grm_ids.SAMPLE_ID, estimandsfile, "tmle_output"; pvalue_threshold=pvalue_threshold) - results, influence_curves, n_obs = TMLECLI.build_work_list(joinpath(tmpdir, "tmle_output"), grm_ids) - # Check n_obs - @test n_obs == [194, 193, 193, 194] - # Check influence curves - expected_influence_curves = [size(r.IC, 1) == 194 ? r.IC : vcat(0, r.IC) for r in results] - for rowindex in 1:4 - @test convert(Vector{Float32}, expected_influence_curves[rowindex]) == influence_curves[rowindex, :] - end - # Check results - @test all(x isa TMLE.TMLEstimate for x in results) - @test all(size(x.IC, 1) > 0 for x in results) -end - -@testset "Test bit_distance" begin - sample_grm = Float32[-0.6, -0.8, -0.25, -0.3, -0.1, 0.1, 0.7, 0.5, 0.2, 1.] - nτs = 6 - τs = TMLECLI.default_τs(nτs, max_τ=0.75) - @test τs == Float32[0.0, 0.15, 0.3, 0.45, 0.6, 0.75] - τs = TMLECLI.default_τs(nτs) - @test τs == Float32[0., 0.4, 0.8, 1.2, 1.6, 2.0] - d = TMLECLI.bit_distances(sample_grm, τs) - @test d == [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 - 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 - 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 - 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 - 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0] -end - -@testset "Test aggregate_variances" begin - # 2 influence curves containing 5 individuals - influence_curves = [1. 2. 3. 4. 5. - 6. 7. 8. 9. 10.] - # distance indicator with 3 τs and corresponding to row 4 - indicator = [1. 0. 0. 0.2 - 0. 0. 1. 1. - 1. 0. 1. 1.] - sample = 4 - var_ = TMLECLI.aggregate_variances(influence_curves, indicator, sample) - @test var_ == [24.0 189.0 - 40.0 225.0 - 48.0 333.0] -end - -@testset "Test normalize!" begin - # 2 τs and 3 curves - n_obs = [10, 10, 100] - variances = [1. 2. 3. - 4. 5. 6.] - TMLECLI.normalize!(variances, n_obs) - @test variances == [0.1 0.2 0.03 - 0.4 0.5 0.06] -end - -@testset "Test compute_variances" begin - n_curves = 3 - n_samples = 5 - nτs = 5 - n_obs = [3, 4, 4] - τs = TMLECLI.default_τs(nτs) - # The GRM has 15 lower triangular elements - grm = Float32[0.4, 0.1, 0.5, 0.2, -0.2, 0.6, 0.3, -0.6, - 0.4, 0.3, 0.6, 0.3, 0.7, 0.3, 0.1] - influence_curves = Float32[0.1 0. 0.1 0.3 0. - 0.1 0.2 0.1 0.0 0.2 - 0.0 0. 0.1 0.3 0.2] - - - variances = TMLECLI.compute_variances(influence_curves, grm, τs, n_obs) - @test size(variances) == (nτs, n_curves) - - # when τ=2, all elements are used - for curve_id in 1:n_curves - s = sum(influence_curves[curve_id, :]) - var = sum(s*influence_curves[curve_id, i] for i in 1:n_samples)/n_obs[curve_id] - @test variances[end, curve_id] ≈ var - end - - # Decreasing variances with τ as all inf curves are positives - for nτ in 1:nτs-1 - @test all(variances[nτ, :] .<= variances[nτ+1, :]) - end - - # Check against basic_variance_implementation - matrix_distance = zeros(Float32, n_samples, n_samples) - for τ_id in 1:nτs - vector_distance = TMLECLI.bit_distances(grm, [τs[τ_id]]) - distance_vector_to_matrix!(matrix_distance, vector_distance, n_samples) - for curve_id in 1:n_curves - influence_curve = influence_curves[curve_id, :] - var_ = basic_variance_implementation(matrix_distance, influence_curve, n_obs[curve_id]) - @test variances[τ_id, curve_id] ≈ var_ - end - end - - # Check by hand for a single τ=0.5 - @test variances[2, :] ≈ Float32[0.03666667, 0.045, 0.045] -end - -@testset "Test grm_rows_bounds" begin - n_samples = 5 - grm_bounds = TMLECLI.grm_rows_bounds(n_samples) - @test grm_bounds == [1 => 1 - 2 => 3 - 4 => 6 - 7 => 10 - 11 => 15] -end - -@testset "Test corrected_stderrors" begin - variances = [ - 1. 2. 6. - 4. 5. 3. - ] - stderrors = TMLECLI.corrected_stderrors(variances) - # sanity check - @test stderrors == sqrt.([4., 5., 6.]) -end - -@testset "Test SVP" begin - # Generate data - grm_ids = TMLECLI.GRMIDs(joinpath(TESTDIR, "data", "grm", "test.grm.id")) - tmpdir = mktempdir() - configuration = statistical_estimands_only_config() - pvalue_threshold = 0.1 - config_1 = TMLE.Configuration(estimands=configuration.estimands[1:3]) - estimandsfile_1 = joinpath(tmpdir, "configuration_1.json") - TMLE.write_json(estimandsfile_1, config_1) - build_tmle_output_file(tmpdir, grm_ids.SAMPLE_ID, estimandsfile_1, "tmle_output_1"; pvalue_threshold=pvalue_threshold) - - config_2 = TMLE.Configuration(estimands=configuration.estimands[4:end]) - estimandsfile_2 = joinpath(tmpdir, "configuration_2.json") - TMLE.write_json(estimandsfile_2, config_2) - build_tmle_output_file(tmpdir, grm_ids.SAMPLE_ID, estimandsfile_2, "tmle_output_2"; pvalue_threshold=pvalue_threshold) - - # Using the main command - output = joinpath(tmpdir, "svp.hdf5") - main([ - "svp", - joinpath(tmpdir, "tmle_output"), - "--out", output, - "--grm-prefix", joinpath(TESTDIR, "data", "grm", "test.grm"), - "--max-tau", "0.75" - ]) - - io = jldopen(output) - # Check τs - @test io["taus"] == TMLECLI.default_τs(10; max_τ=0.75) - # Check variances - @test size(io["variances"]) == (10, 4) - # Check results - svp_results = io["results"] - - tmleout1 = jldopen(x -> x["Batch_1"], joinpath(tmpdir, "tmle_output_1.hdf5")) - tmleout2 = jldopen(x -> x["Batch_1"], joinpath(tmpdir, "tmle_output_2.hdf5")) - src_results = [tmleout1..., tmleout2...] - - for svp_result in svp_results - src_result_index = findall(x.TMLE.estimand == svp_result.SVP.estimand for x in src_results) - src_result = src_results[only(src_result_index)] - @test src_result.TMLE.std != svp_result.SVP.std - @test src_result.TMLE.estimate == svp_result.SVP.estimate - @test src_result.TMLE.n == svp_result.SVP.n - @test svp_result.SVP.IC == [] - end - - close(io) -end - -@testset "Test SVP: causal and composed estimands" begin - # Generate data - grm_ids = TMLECLI.GRMIDs(joinpath(TESTDIR, "data", "grm", "test.grm.id")) - tmpdir = mktempdir() - configuration = causal_and_joint_estimands_config() - configfile = joinpath(tmpdir, "configuration.json") - TMLE.write_json(configfile, configuration) - build_tmle_output_file( - tmpdir, - grm_ids.SAMPLE_ID, - configfile, - "tmle_output"; - estimatorfile=joinpath(TESTDIR, "config", "ose_config.jl") - ) - - # Using the main command - output = joinpath(tmpdir, "svp.hdf5") - main([ - "svp", - joinpath(tmpdir, "tmle_output"), - "--out", output, - "--grm-prefix", joinpath(TESTDIR, "data", "grm", "test.grm"), - "--max-tau", "0.75", - "--estimator-key", "OSE" - ]) - - # The JointEstimate std is not updated but each component is. - src_results = jldopen(x -> x["Batch_1"], joinpath(tmpdir, "tmle_output.hdf5")) - io = jldopen(output) - svp_results = io["results"] - standalone_estimates = svp_results[1:2] - from_composite = svp_results[3:4] - @test standalone_estimates[1].SVP.estimand == from_composite[1].SVP.estimand - @test standalone_estimates[2].SVP.estimand == from_composite[2].SVP.estimand - - # Check std has been updated - for i in 1:2 - @test standalone_estimates[i].SVP.estimand == src_results[i].OSE.estimand - @test standalone_estimates[i].SVP.estimate == src_results[i].OSE.estimate - @test standalone_estimates[i].SVP.std != src_results[i].OSE.std - end - - close(io) -end - -end - -true diff --git a/test/summary.jl b/test/summary.jl index c47b86e..44bf8b8 100644 --- a/test/summary.jl +++ b/test/summary.jl @@ -46,13 +46,14 @@ include(joinpath(TESTDIR, "testutils.jl")) json_output = joinpath(tmpdir, "summary.json") jls_output = joinpath(tmpdir, "summary.jls") hdf5_output = joinpath(tmpdir, "summary.hdf5") - main([ + copy!(ARGS, [ "merge", joinpath(tmpdir, "tmle_output"), "--json-output", json_output, "--jls-output", jls_output, "--hdf5-output", hdf5_output ]) + TMLECLI.julia_main() # Test correctness inputs = TMLECLI.read_results_from_files([joinpath(tmpdir, "tmle_output_1.hdf5"), joinpath(tmpdir, "tmle_output_2.hdf5")]) diff --git a/tmle.jl b/tmle.jl index 2acd5ca..7cfd1e3 100644 --- a/tmle.jl +++ b/tmle.jl @@ -1 +1 @@ -using TMLECLI; main() \ No newline at end of file +using TMLECLI; julia_main() \ No newline at end of file From 851ea871d447431896089b216ab9f367263a1cd1 Mon Sep 17 00:00:00 2001 From: Olivier Labayle Date: Fri, 6 Sep 2024 15:19:20 +0100 Subject: [PATCH 2/2] up version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7a3b515..333baa2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TMLECLI" uuid = "2573d147-4098-46ba-9db2-8608d210ccac" authors = ["Olivier Labayle"] -version = "0.9.1" +version = "0.10.0" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"