diff --git a/benchmarks/basic_tests.py b/benchmarks/basic_tests.py index 953c1326..9f5453e0 100644 --- a/benchmarks/basic_tests.py +++ b/benchmarks/basic_tests.py @@ -1351,7 +1351,7 @@ def build_and_draw_prior(t_d="ends",num_reals=500): #glm_long_name_test() #sen_plusplus_test() #parchglim_test() - #unc_file_test() + unc_file_test() #cmdline_test() #secondary_marker_test() #basic_test("ies_10par_xsec") @@ -1382,7 +1382,7 @@ def build_and_draw_prior(t_d="ends",num_reals=500): #mf6_v5_sen_test() #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-opt.exe"),os.path.join("..","bin","win","pestpp-opt.exe")) - mf6_v5_opt_stack_test() + #mf6_v5_opt_stack_test() #mf6_v5_glm_test() #mf6_v5_ies_test() #cmdline_test() diff --git a/documentation/pestpp_users_guide_v5.2.3.docx b/documentation/pestpp_users_guide_v5.2.4.docx similarity index 71% rename from documentation/pestpp_users_guide_v5.2.3.docx rename to documentation/pestpp_users_guide_v5.2.4.docx index 30fa3918..6c94028f 100644 Binary files a/documentation/pestpp_users_guide_v5.2.3.docx and b/documentation/pestpp_users_guide_v5.2.4.docx differ diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md index 6aec840e..49dff663 100644 --- a/documentation/pestpp_users_manual.md +++ b/documentation/pestpp_users_manual.md @@ -1,13 +1,13 @@ A close up of a purple sign Description automatically generated -# Version 5.2.3 +# Version 5.2.4 PEST++ Development Team -April 2023 +May 2023 # Acknowledgements @@ -70,7 +70,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI # Table of Contents -- [Version 5.2.3](#s1) +- [Version 5.2.4](#s1) - [Acknowledgements](#s2) - [Preface](#s3) - [License](#s4) @@ -2693,6 +2693,8 @@ If an observation belongs to an observation group whose name begins with the str Similarly, if an observation belongs to an observation group whose name begins with the string “l\_” or “less\_”, then this observation is a “less than” observation. No objective function penalty is incurred if the modelled value of the pertinent quantity is less than the measured value listed in the “observation data” section of the PEST control file. However, if the model-calculated value is greater than the measured value, the objective function penalty is calculated in the usual manner, that is as the squared residual times the squared weight. +Users that wish to (very) strongly enforce inequality conditions are advised to avoid using extremely high weights for inequalities, and instead are encouraged to use the “drop_violations” functionality described later. + ### 9.1.11 Localization Calculating an empirical cross-covariance between large numbers of parameters and observations from a limited number of realizations is likely to result in spurious cross-correlations. Because of this, some parameters will be adjusted when they should not be adjusted. Furthermore, when large numbers of independent observations comprise a calibration dataset, a small ensemble size will almost certainly not provide enough degrees of freedom to reproduce these data. To combat these problems, users can employ localization. The term “localization” comes from ensemble Kalman filter parlance. It refers to a strategy whereby only “local” covariances inform unmeasured states in a spatially distributed filtering problem. @@ -2705,7 +2707,7 @@ Figure 9.1. An example localization matrix. When applying localization in a history-matching problem involving large numbers of parameters and observations, a user may wish to define a “local” neighbourhood around each observation location wherein parameters are expected to influence the simulated counterparts to observations. This, in effect, creates a series of “local” history-matching problems using subsets of adjustable parameters and observations. The number of degrees of freedom featured in each local problem can be relatively high, this allowing a small ensemble size to better reproduce large numbers of independent observations. Localization also provides protection against “spurious” (non-plausible) correlations between parameters and observations arising from the limited size of parameter ensembles. For example, standard methods of covariance calculation may suggest a correlation between a pumping rate parameter and a head that precedes it in time. Spurious correlations of this type can lead to parameter compensation and predictive bias. See Chen and Oliver (2016) for a good description of the theory and practice of localization. -Through localization, a large-scale parameter estimation problem can be turned into a series of independent parameter estimation problems. Suppose, for example, that localization is employed in the most granular manner, so that the localization matrix contains one column for each adjustable parameter and that each column contains a single non-zero value, this pertaining to a single observation which that parameter is presumed to influence. If large numbers of parameters are being adjusted, the parameter upgrade calculation process for a given lambda will require as many truncated SVD solves as there are adjustable parameters. This can require considerable numerical effort. To overcome this problem, the localized upgrade solution process in PESTPP-IES has been multithreaded; this is possible in circumstances such as these where each local solve is independent of every other local solve. The use of multiple threads is invoked through the *ies_num_threads()* control variable. It should be noted that the optimal number of threads to use is problem-specific. Furthermore, it should not exceed the number of physical cores of the host machine on which the PESTPP-IES master instance is running. Note that pestpp-ies does not require that all adjustable parameters and all non-zero weighted observations be included in the localization matrix.. In this case, pestpp-ies implicitly treats any adjustable parameters not in the localizer as “fixed” and any non-zero weighted observations not in the localizer as zero-weighted. However, users are strongly encouraged to explicitly include all adjustable parameters and all non-zero weighted observations in the localization matrix: explicit is better than implicit… +Through localization, a large-scale parameter estimation problem can be turned into a series of independent parameter estimation problems. Suppose, for example, that localization is employed in the most granular manner, so that the localization matrix contains one column for each adjustable parameter and that each column contains a single non-zero value, this pertaining to a single observation which that parameter is presumed to influence. If large numbers of parameters are being adjusted, the parameter upgrade calculation process for a given lambda will require as many truncated SVD solves as there are adjustable parameters. This can require considerable numerical effort. To overcome this problem, the localized upgrade solution process in PESTPP-IES has been multithreaded; this is possible in circumstances such as these where each local solve is independent of every other local solve. The use of multiple threads is invoked through the *ies_num_threads()* control variable. It should be noted that the optimal number of threads to use is problem-specific. Furthermore, it should not exceed the number of physical cores of the host machine on which the PESTPP-IES master instance is running. Note that pestpp-ies does not require that all adjustable parameters and all non-zero weighted observations be included in the localization matrix. In this case, pestpp-ies implicitly treats any adjustable parameters not in the localizer as “fixed” and any non-zero weighted observations not in the localizer as zero-weighted. However, users are strongly encouraged to explicitly include all adjustable parameters and all non-zero weighted observations in the localization matrix: explicit is better than implicit… PESTPP-IES also supports correlation-based, automatic adaptive localization. Its algorithm is based on that of Luo et al (2018). This functionality is activated by setting the *ies_autoadaloc()* control variable to *true*. If this localization scheme is employed, a user does not need to provide a localization matrix to define which parameters can be informed by which observations. Instead, PESTPP-IES uses the parameter ensemble and resulting model-output observation ensemble to calculate a Pearson correlation coefficient between each adjustable parameter and each non-zero-weighted observation. Because the normal usage context of PESTPP-IES is that of numerically efficient parameter estimation and uncertainty quantification, it is expected that the number of realizations that constitute an ensemble will normally be considerably smaller than the number of adjustable parameters. The estimated correlation coefficient between any given parameter and any observation is therefore likely to be somewhat in error. This error decreases (eventually to zero) as the size of an ensemble rises. @@ -2739,6 +2741,8 @@ While detecting prior-data conflict is relatively simple (and PESTPP-IES will do With these options active, PESTPP-IES will remove observations in a prior-data conflict state from the parameter adjustment process, that is, these observations will not feature in the residual matrix used for upgrade calculations. +Note that observations tagged with the “drop_violations” flag will not be “dropped” for prior-data conflict reasons. + ### 9.1.14 Multi-modal solution process The theory that underpins PEST, PESTPP-GLM and PESTPP-IES is designed to seek a single minimum of the objective function (e.g., a single peak of the likelihood function). For many data assimilation problems however, the objective function surface is not convex with a single minimum located neatly at the bottom. Instead, it maybe pitted with local minimum and/or contain a curving high-dimensional trough of nearly equal objective function minimum. In the context of deterministic parameter estimation, the goal is to seek a unique minimum of the objective function, that is seeking a minimum error variance solution. However, in the context of uncertainty quantification, and especially in the context of evaluating the likelihood of an unwanted outcome, exploring local minima and/or this high-dimensional trough is paramount. @@ -2889,6 +2893,8 @@ Where model runs are based on random parameter realizations, the risk of occasio PESTPP-IES provides a mechanism for detection of model run failure that extends those provided by its run manager. If the objective function associated with a particular model run is calculated to be greater than a certain threshold, PESTPP-IES deems the model run to have failed. This threshold is supplied as the value of *ies_bad_phi()* control variable. In addition the absolute filtering/rejecting of realizations provided by *ies_bad_phi*, users can also activate a relative filter with *ies_bad_phi_sigma*. This option accepts a floating-poin number that represents the number of standard deviations above the mean that a realizations phi value must be before it is considered “bad”. This fitler is adaptive in that as the realizations move toward a minimum of the objective function, the measured phi value of each realization changes. A setting of 2.0 for *ies_bad_phi_sigma* results in realizations with a measured phi exceeding the mean phi plus 2 standard deviations being removed. This obviously assumes a normal distribution of measured phi values, which is often violated. In the case where users are concerned about this normality assumption, the value of *ies_bad_phi_sigma* can be supplied as a negative value, in which case the absolute value of *ies_bad_phi_sigma* is treated at the empirical phi quantile value above which realizations are treated as “bad”. For example, an *ies_bad_phi_sigma* value of -95 indicates that realizations above the upper 95th quantile of phi should removed. +PESTPP-IES also supports a more granular approach to identifying and removing “bad” realizations through the use of an additional column in the observation data section (this requires the version 2 control file). This column must be named “drop_violations” and it should contain entries that are either “true” or “false”. If true, then any realization that violates the observation condition will be dropped. This functionality is designed to work with inequality-type observations so that if a realization violates the inequality, it is removed from the ensemble, rather than penalizing the objective function. Note “drop_violations” also works with standard equality type observation, in which case, any realization that does not numerically equal the *obsval* quantity of the observation (to within 1.0e-7) will also be dropped. The “drop_violations” functionality respects the weight values for observations, so users can activate and de-activate the violation enforcement by simply changing the weights from non-zero to zero, respectively. + To forestall excessive PESTPP-IES run times incurred by occasional model failure, it is a good idea to set the *max_run_fail()* model run control variable to 1 (the default value for PESTPP‑IES), and to choose values for the *overdue_giveup_fac()* and/or *overdue_giveup_minutes()* control variables judiciously; see section 5.3. Note also that model run failure does not hurt PESTPP-IES as much as it hurts PESTPP-GLM or PEST. This is because the value of any model run undertaken by PESTPP-IES is lower than that undertaken by PEST or PESTPP-GLM. For the latter programs a failed model run during finite-difference derivatives calculation may lead to an empty column of the Jacobian matrix. In contrast, because PESTPP-IES uses an entire ensemble to fill a Jacobian matrix, a single failed model run does not result in an empty Jacobian matrix column. The outcome of model run failure is that the number of model runs employed in the averaging process through which this column is calculated is reduced by one. ### 9.2.8 Reporting @@ -2983,7 +2989,7 @@ Table 9.4 lists PESTPP-IES control variables. All of these are optional. If a va Note also that the number of control variables may change with time. Refer to the PEST++ web site for variables used by the latest version of PESTPP-IES. -
VariableTypeRole
ies_num_reals(50)integerThe number of realizations to draw in order to form parameter and observation ensembles.
parcov()textThe name of a file containing the prior parameter covariance matrix. This can be a parameter uncertainty file (extension .unc), a covariance matrix file (extension .cov) or a binary JCO or JCB file (extension .jco or .jcb).
par_sigma_range(4.0)realThe difference between a parameter’s upper and lower bounds expressed as standard deviations.
ies_parameter_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied parameter realizations comprising the initial (prior) parameter ensemble. If this keyword is omitted, PESTPP-IES generates the initial parameter ensemble itself.
ies_observation_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied observation plus noise realizations comprising the observation plus noise ensemble. If this keyword is omitted, PESTPP-IES generates the observation plus noise ensemble itself.
ies_add_base(true)BooleanIf set to true, instructs PESTPP-IES to include a “realization” in the initial parameter ensemble comprised of parameter values read from the “parameter data” section of the PEST control file. The corresponding observation ensemble is comprised of measurements read from the “observation data” section of the PEST control file.
ies_restart_observation_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing model outputs calculated using a parameter ensemble. If it reads this file, PESTPP-IES does not calculate these itself, proceeding to upgrade calculations instead.
ies_restart_parameter_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing a parameter ensemble that corresponds to the ies_restart_observation_ensemble(). This option requires that the ies_restart_observation_ensemble() control variable also be supplied. This ensemble is only used in the calculation of the regularization component of the objective function for a restarted PESTPP-IES analysis.
ies_enforce_bounds(true)BooleanIf set to true PESTPP-IES will not transgress bounds supplied in the PEST control file when generating or accepting parameter realizations, and when adjusting these realizations.
ies_initial_lambda()realThe initial Marquardt lambda. The default value is \(10^{\text{floor}\left( \log_{10}\frac{\mu_{Փ}}{2n} \right)}\text{.\ \ }\)If supplied as a negative value, then the abs(ies_initial_lambda) is used as multiplier of the default initial-phi-based value.
ies_lambda_mults(0.1,1.0,10.0)comma-separated realsFactors by which to multiply the best lambda from the previous iteration to yield values for testing parameter upgrades during the current iteration.
lambda_scale_fac(0.75,1.0,1.1)comma-separated realsLine search factors along parameter upgrade directions computed using different Marquardt lambdas.
ies_subset_size(4)integerNumber of realizations used in testing and evaluation of different Marquardt lambdas. If supplied as a negative value, then abs(ies_subset_size) is treated as a percentage of the current ensemble size – this allows the subset size to fluctuate with the size of the ensemble
ies_use_approx(true)BooleanUse complex or simple formula provided by Chen and Oliver (2013) for calculation of parameter upgrades. The more complex formula includes a function which constrains parameter realizations to respect prior means and probabilities.
ies_reg_factor(0.0)realRegularization objective function as a fraction of measurement objective function when constraining parameter realizations to respect initial values.
ies_bad_phi(1.0E300)realIf the objective function calculated as an outcome of a model run is greater than this value, the model run is deemed to have failed.
ies_bad_phi_sigma(1.0E300)realIf the objective function calculated for a given realization is greater than the current mean objective function of the ensemble plus the objective function standard deviation of the ensemble times ies_bad_phi_sigma(), that realization is treated as failed. If negative, its absolute value is treated as the upper quantile for identifying “bad” realizations.
ies_use_prior_scaling(false)BooleanUse a scaling factor based on the prior parameter distribution when evaluating parameter-to-model-output covariance used in calculation of the randomized Jacobian matrix.
ies_use_empirical_prior(false)BooleanUse an empirical, diagonal parameter covariance matrix for certain calculations. This matrix is contained in a file whose name is provided with the ies_parameter_ensemble() keyword.
Ies_save_lambda_ensembles(false)BooleanSave a set of CSV or JCB files that record parameter realizations used when testing different Marquardt lambdas.
ies_verbose_level(1)0, 1 or 2The level of diagnostic output provided by PESTPP-IES. If set to 2, all intermediate matrices are saved to ASCII files. This can require a considerable amount of storage.
ies_accept_phi_fac(1.05)real > 1.0The factor applied to the previous best mean objective function to determine if the current mean objective function is acceptable.
ies_lambda_dec_fac(0.75)real < 1.0The factor by which to decrease the value of the Marquardt lambda during the next IES iteration if the current iteration of the ensemble smoother process was successful in lowering the mean objective function.
ies_lambda_inc_fac(10.0)real > 1.0The factor by which to increase the current value of the Marquardt lambda for further lambda testing if the current lambda testing cycle was unsuccessful.
ies_subset_how(random)“first”,”last”,
”random”,
”phi_based
How to select the subset of realizations for objective function evaluation during upgrade testing. Default is “random”.
ies_num_threads(-1)integer > 1The number of threads to use during the localized upgrade solution process, the automatic adaptive localization process. If the localizer contains many (>10K) rows, then multithreading can substantially speed up the upgrade calculation process. ies_num_threads() should not be greater than the number of physical cores on the host machine.
ies_localizer()textThe name of a matrix to use for localization. The extension of the file is used to determine the type: .mat is an ASCII matrix file, .jcb/.jco signifies use of (enhanced) Jacobian matrix format (a binary format), while .csv signifies a comma-delimited file. Note that adjustable parameters not listed in localization matrix columns are implicitly treated as “fixed” while non-zero weighted observations not listed in rows of this matrix are implicitly treated as zero-weighted.
ies_group_draws(true)BooleanA flag to draw from the (multivariate) Gaussian prior by parameter/observation groups. This is usually a good idea since groups of parameters/observations are likely to have prior correlation.
ies_save_binary(false)BooleanA flag to save parameter and observation ensembles in binary (i.e., JCB) format instead of CSV format.
ies_csv_by_reals(true)BooleanA flag to save parameter and observation ensemble CSV files by realization instead of by variable name. If true, each row of the CSV file is a realization. If false, each column of the CSV file is a realization.
ies_autoadaloc(false)BooleanFlag to activate automatic adaptive localization.
ies_autoadaloc_sigma_dist(1.0)RealReal number representing the factor by which a correlation coefficient must exceed the standard deviation of background correlation coefficients to be considered significant. Default is 1.0
tie_by_group(false)BooleanFlag to tie all adjustable parameters together within each parameter group. Initial parameter ratios are maintained as parameters are adjusted. Parameters that are designated as already tied, or that have parameters tied to them, are not affected.
ies_enforce_chglim(false)BooleanFlag to enforce parameter change limits (via FACPARMAX and RELPARMAX) in a way similar to PEST and PESTPP-GLM (by scaling the entire realization). Default is false.
ies_center_on()StringA realization name that should be used for the ensemble center in calculating the approximate Jacobian matrix. The realization name must be in both the parameter and observation ensembles. If not passed, the mean vector is used as the center. The value “_MEDIAN_” can also be used, which instructs PESTPP-IES to use the median vector for calculating anomalies.
enforce_tied_bounds(false)BooleanFlag to enforce parameter bounds on any tied parameters. Depending on the ration between the tied and free parameters, this option can greatly limit parameter changes.
ies_no_noise(false)BooleanFlag to not generate and use realizations of measurement noise. Default is False (that is, to use measurement noise).
ies_drop_conflicts(false)BooleanFlag to remove non-zero weighted observations that are in a prior-data conflict state from the upgrade calculations. Default is False.
ies_pdc_sigma_distance()Real > 0.0The number of standard deviations from the mean used in checking for prior-data conflict.
ies_save_rescov(False)BooleanFlag to save the iteration-level residual covariance matrix. If ies_save_binary is True, then a binary format file is written, otherwise an ASCII format (.cov) file is written. The file name is case.N.res.cov/.jcb. Note that this functionality does not scale beyond about 20,000 non-zero-weighted observations
obscov()textThe name of a file containing the observation noise covariance matrix. This can be a parameter uncertainty file (extension .unc), a covariance matrix file (extension .cov) or a binary JCO or JCB file (extension .jco or .jcb). Please see the section on this matrix above to understand the implications of using this matrix
rand_seed(358183147)unsigned integerSeed for the random number generator.
Ies_use_mda(false)BooleanFlag to use the (optionally iterative) Kalman update equation – the number of data assimilation iterations is controlled by NOPTMAX; NOPTMAX = 1 and ies_use_mda(true) results in the standard ensemble smoother Kalman update. If False, the GLM iterative ensemble smoother equation is used. Default is False
Ies_mda_init_fac(10.0)doubleThe initial MDA covariance inflation factor. Only used if ies_use_mda is true. Default is 10.0
Ies_mda_decl_fac(0.5)doubleThe final MDA covariance inflation factor. Only used in ies_use_mda is true. Default is 0.5
Ies_localization_type(local)textCan be either “local” for local analysis or “covariance” for covariance-only localization. Default is “local”
Ies_upgrades_in_memory(true)BooleanFlag to hold parameter upgrade ensembles in memory during testing. If False, parameter ensembles are saved to disk during testing and the best-phi ensemble is loaded from disk after testing – this can reduce memory pressure for very high dimensional problems. Default is True but is only activated if number of parameters > 100K.
Ies_ordered_binary(true)BooleanFlag to write control-file-ordered binary ensemble files. Only used if save_binary is true. If false, hash-ordered binary files are written – for very high dimensional problems, writing unordered binary can save lots of time. If not passed and number of parameters > 100K, then ies_ordered_binary is set to false.
ensemble_output_precision(6)intNumber of significant digits to use in ASCII format ensemble files. Default is 6
ies_multimodal_alpha(1.0)doubleThe fraction of the total ensemble size to use as the local neighborhood realizations in the multimodal solution process. Must be greater than zero and less than 1. Values of 0.1 to 0.25 seem to work well. Default is 1.0 (disable multi-modal solution process)
ies_weight_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied weight vectors for each realization. If this keyword is omitted, PESTPP-IES uses the weight vector in the control file for all realizations. Only used with ies_multimodal_alpha
Ies_phi_factor_filetextA two-column ASCII file that contains observation group “tags” and phi factors. Used to internally adjust weights to implement a balanced objective function using the mean residuals from the initial ensemble.
Ies_phi_factors_by_realBoolA flag to use internal weight balancing for each realization. This option should be used in conjunction with the multi-modal solution process.
+
VariableTypeRole
ies_num_reals(50)integerThe number of realizations to draw in order to form parameter and observation ensembles.
parcov()textThe name of a file containing the prior parameter covariance matrix. This can be a parameter uncertainty file (extension .unc), a covariance matrix file (extension .cov) or a binary JCO or JCB file (extension .jco or .jcb).
par_sigma_range(4.0)realThe difference between a parameter’s upper and lower bounds expressed as standard deviations.
ies_parameter_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied parameter realizations comprising the initial (prior) parameter ensemble. If this keyword is omitted, PESTPP-IES generates the initial parameter ensemble itself.
ies_observation_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied observation plus noise realizations comprising the observation plus noise ensemble. If this keyword is omitted, PESTPP-IES generates the observation plus noise ensemble itself.
ies_add_base(true)BooleanIf set to true, instructs PESTPP-IES to include a “realization” in the initial parameter ensemble comprised of parameter values read from the “parameter data” section of the PEST control file. The corresponding observation ensemble is comprised of measurements read from the “observation data” section of the PEST control file.
ies_restart_observation_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing model outputs calculated using a parameter ensemble. If it reads this file, PESTPP-IES does not calculate these itself, proceeding to upgrade calculations instead.
ies_restart_parameter_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing a parameter ensemble that corresponds to the ies_restart_observation_ensemble(). This option requires that the ies_restart_observation_ensemble() control variable also be supplied. This ensemble is only used in the calculation of the regularization component of the objective function for a restarted PESTPP-IES analysis.
ies_enforce_bounds(true)BooleanIf set to true PESTPP-IES will not transgress bounds supplied in the PEST control file when generating or accepting parameter realizations, and when adjusting these realizations.
ies_initial_lambda()realThe initial Marquardt lambda. The default value is \(10^{\text{floor}\left( \log_{10}\frac{\mu_{Փ}}{2n} \right)}\text{.\ \ }\)If supplied as a negative value, then the abs(ies_initial_lambda) is used as multiplier of the default initial-phi-based value.
ies_lambda_mults(0.1,1.0,10.0)comma-separated realsFactors by which to multiply the best lambda from the previous iteration to yield values for testing parameter upgrades during the current iteration.
lambda_scale_fac(0.75,1.0,1.1)comma-separated realsLine search factors along parameter upgrade directions computed using different Marquardt lambdas.
ies_subset_size(4)integerNumber of realizations used in testing and evaluation of different Marquardt lambdas. If supplied as a negative value, then abs(ies_subset_size) is treated as a percentage of the current ensemble size – this allows the subset size to fluctuate with the size of the ensemble
ies_use_approx(true)BooleanUse complex or simple formula provided by Chen and Oliver (2013) for calculation of parameter upgrades. The more complex formula includes a function which constrains parameter realizations to respect prior means and probabilities.
ies_reg_factor(0.0)realRegularization objective function as a fraction of measurement objective function when constraining parameter realizations to respect initial values.
ies_bad_phi(1.0E300)realIf the objective function calculated as an outcome of a model run is greater than this value, the model run is deemed to have failed.
ies_bad_phi_sigma(1.0E300)realIf the objective function calculated for a given realization is greater than the current mean objective function of the ensemble plus the objective function standard deviation of the ensemble times ies_bad_phi_sigma(), that realization is treated as failed. If negative, its absolute value is treated as the upper quantile for identifying “bad” realizations.
ies_use_prior_scaling(false)BooleanUse a scaling factor based on the prior parameter distribution when evaluating parameter-to-model-output covariance used in calculation of the randomized Jacobian matrix.
ies_use_empirical_prior(false)BooleanUse an empirical, diagonal parameter covariance matrix for certain calculations. This matrix is contained in a file whose name is provided with the ies_parameter_ensemble() keyword.
Ies_save_lambda_ensembles(false)BooleanSave a set of CSV or JCB files that record parameter realizations used when testing different Marquardt lambdas.
ies_verbose_level(1)0, 1 or 2The level of diagnostic output provided by PESTPP-IES. If set to 2, all intermediate matrices are saved to ASCII files. This can require a considerable amount of storage.
ies_accept_phi_fac(1.05)real > 1.0The factor applied to the previous best mean objective function to determine if the current mean objective function is acceptable.
ies_lambda_dec_fac(0.75)real < 1.0The factor by which to decrease the value of the Marquardt lambda during the next IES iteration if the current iteration of the ensemble smoother process was successful in lowering the mean objective function.
ies_lambda_inc_fac(10.0)real > 1.0The factor by which to increase the current value of the Marquardt lambda for further lambda testing if the current lambda testing cycle was unsuccessful.
ies_subset_how(random)“first”,”last”,
”random”,
”phi_based
How to select the subset of realizations for objective function evaluation during upgrade testing. Default is “random”.
ies_num_threads(-1)integer > 1The number of threads to use during the localized upgrade solution process, the automatic adaptive localization process. If the localizer contains many (>10K) rows, then multithreading can substantially speed up the upgrade calculation process. ies_num_threads() should not be greater than the number of physical cores on the host machine.
ies_localizer()textThe name of a matrix to use for localization. The extension of the file is used to determine the type: .mat is an ASCII matrix file, .jcb/.jco signifies use of (enhanced) Jacobian matrix format (a binary format), while .csv signifies a comma-delimited file. Note that adjustable parameters not listed in localization matrix columns are implicitly treated as “fixed” while non-zero weighted observations not listed in rows of this matrix are implicitly treated as zero-weighted.
ies_group_draws(true)BooleanA flag to draw from the (multivariate) Gaussian prior by parameter/observation groups. This is usually a good idea since groups of parameters/observations are likely to have prior correlation.
ies_save_binary(false)BooleanA flag to save parameter and observation ensembles in binary (i.e., JCB) format instead of CSV format.
ies_csv_by_reals(true)BooleanA flag to save parameter and observation ensemble CSV files by realization instead of by variable name. If true, each row of the CSV file is a realization. If false, each column of the CSV file is a realization.
ies_autoadaloc(false)BooleanFlag to activate automatic adaptive localization.
ies_autoadaloc_sigma_dist(1.0)RealReal number representing the factor by which a correlation coefficient must exceed the standard deviation of background correlation coefficients to be considered significant. Default is 1.0
tie_by_group(false)BooleanFlag to tie all adjustable parameters together within each parameter group. Initial parameter ratios are maintained as parameters are adjusted. Parameters that are designated as already tied, or that have parameters tied to them, are not affected.
ies_enforce_chglim(false)BooleanFlag to enforce parameter change limits (via FACPARMAX and RELPARMAX) in a way similar to PEST and PESTPP-GLM (by scaling the entire realization). Default is false.
ies_center_on()StringA realization name that should be used for the ensemble center in calculating the approximate Jacobian matrix. The realization name must be in both the parameter and observation ensembles. If not passed, the mean vector is used as the center. The value “_MEDIAN_” can also be used, which instructs PESTPP-IES to use the median vector for calculating anomalies.
enforce_tied_bounds(false)BooleanFlag to enforce parameter bounds on any tied parameters. Depending on the ration between the tied and free parameters, this option can greatly limit parameter changes.
ies_no_noise(false)BooleanFlag to not generate and use realizations of measurement noise. Default is False (that is, to use measurement noise).
ies_drop_conflicts(false)BooleanFlag to remove non-zero weighted observations that are in a prior-data conflict state from the upgrade calculations. Default is False.
ies_pdc_sigma_distance()Real > 0.0The number of standard deviations from the mean used in checking for prior-data conflict.
ies_save_rescov(False)BooleanFlag to save the iteration-level residual covariance matrix. If ies_save_binary is True, then a binary format file is written, otherwise an ASCII format (.cov) file is written. The file name is case.N.res.cov/.jcb. Note that this functionality does not scale beyond about 20,000 non-zero-weighted observations
obscov()textThe name of a file containing the observation noise covariance matrix. This can be a parameter uncertainty file (extension .unc), a covariance matrix file (extension .cov) or a binary JCO or JCB file (extension .jco or .jcb). Please see the section on this matrix above to understand the implications of using this matrix
rand_seed(358183147)unsigned integerSeed for the random number generator.
Ies_use_mda(false)BooleanFlag to use the (optionally iterative) Kalman update equation – the number of data assimilation iterations is controlled by NOPTMAX; NOPTMAX = 1 and ies_use_mda(true) results in the standard ensemble smoother Kalman update. If False, the GLM iterative ensemble smoother equation is used. Default is False
Ies_mda_init_fac(10.0)doubleThe initial MDA covariance inflation factor. Only used if ies_use_mda is true. Default is 10.0
Ies_mda_decl_fac(0.5)doubleThe final MDA covariance inflation factor. Only used in ies_use_mda is true. Default is 0.5
Ies_upgrades_in_memory(true)BooleanFlag to hold parameter upgrade ensembles in memory during testing. If False, parameter ensembles are saved to disk during testing and the best-phi ensemble is loaded from disk after testing – this can reduce memory pressure for very high dimensional problems. Default is True but is only activated if number of parameters > 100K.
Ies_ordered_binary(true)BooleanFlag to write control-file-ordered binary ensemble files. Only used if save_binary is true. If false, hash-ordered binary files are written – for very high dimensional problems, writing unordered binary can save lots of time. If not passed and number of parameters > 100K, then ies_ordered_binary is set to false.
ensemble_output_precision(6)intNumber of significant digits to use in ASCII format ensemble files. Default is 6
ies_multimodal_alpha(1.0)doubleThe fraction of the total ensemble size to use as the local neighborhood realizations in the multimodal solution process. Must be greater than zero and less than 1. Values of 0.1 to 0.25 seem to work well. Default is 1.0 (disable multi-modal solution process)
ies_weight_ensemble()textThe name of a CSV or JCO/JCB file (recognized by its extension) containing user-supplied weight vectors for each realization. If this keyword is omitted, PESTPP-IES uses the weight vector in the control file for all realizations. Only used with ies_multimodal_alpha
Ies_phi_factor_filetextA two-column ASCII file that contains observation group “tags” and phi factors. Used to internally adjust weights to implement a balanced objective function using the mean residuals from the initial ensemble.
Ies_phi_factors_by_realBoolA flag to use internal weight balancing for each realization. This option should be used in conjunction with the multi-modal solution process.
Table 9.4 PESTPP-IES control variables with default values. Parallel run management variables can be supplied in addition to these. See section 5.3.6. diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index 73b3a603..5db265bc 100644 --- a/src/libs/common/config_os.h +++ b/src/libs/common/config_os.h @@ -2,7 +2,7 @@ #define CONFIG_OS_H_ -#define PESTPP_VERSION "5.2.3"; +#define PESTPP_VERSION "5.2.4"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index 0d3422d4..6421af0a 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -1017,8 +1017,8 @@ void EnsembleSolver::nonlocalized_solve(double cur_lam,bool use_glm_form, Parame obs_resid.transposeInPlace(); par_resid.transposeInPlace(); obs_err.transposeInPlace(); - ensemble_solution(iter,verbose_level,maxsing,0,0,use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid, - par_diff,Am,obs_resid,obs_diff,upgrade_1,obs_err,local_weights,parcov_inv); + UpgradeThread::ensemble_solution(iter,verbose_level,maxsing,0,0,use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid, + par_diff,Am,obs_resid,obs_diff,upgrade_1,obs_err,local_weights,parcov_inv, act_obs_names,act_par_names); pe_upgrade.add_2_cols_ip(act_par_names, upgrade_1); @@ -1043,17 +1043,31 @@ void EnsembleSolver::solve(int num_threads, double cur_lam, bool use_glm_form, P Localizer::How _how = localizer.get_how(); Localizer::LocTyp loctyp = localizer.get_loctyp(); bool use_cov_loc = true; - if (loctyp == Localizer::LocTyp::LOCALANALYSIS) - use_cov_loc = false; + if (loctyp == Localizer::LocTyp::COVARIANCE) { + throw runtime_error("EnsembleSolver::solve(): 'covariacne' localization is deprecated"); + } //LocalAnalysisUpgradeThread worker(performance_log, par_resid_map, par_diff_map, obs_resid_map, obs_diff_map,obs_err_map, // localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, _how); UpgradeThread* ut_ptr; - if (!use_cov_loc) - ut_ptr = new LocalAnalysisUpgradeThread(performance_log, par_resid_map, par_diff_map, obs_resid_map, obs_diff_map, obs_err_map, - localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, _how); - else - ut_ptr = new CovLocalizationUpgradeThread(performance_log, par_resid_map, par_diff_map, obs_resid_map, obs_diff_map, obs_err_map, - localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, _how); + ut_ptr = new LocalAnalysisUpgradeThread(performance_log, par_resid_map, par_diff_map, obs_resid_map, + obs_diff_map, obs_err_map, + localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, + _how); + performance_log->log_event("using local analysis upgrade thread"); +// if (!use_cov_loc) { +// ut_ptr = new LocalAnalysisUpgradeThread(performance_log, par_resid_map, par_diff_map, obs_resid_map, +// obs_diff_map, obs_err_map, +// localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, +// _how); +// performance_log->log_event("using local analysis upgrade thread"); +// } +// else { +// ut_ptr = new CovLocalizationUpgradeThread(performance_log, par_resid_map, par_diff_map, obs_resid_map, +// obs_diff_map, obs_err_map, +// localizer, parcov_inv_map, weight_map, pe_upgrade, loc_map, Am_map, +// _how); +// performance_log->log_event("using covariance localization upgrade thread"); +// } if ((num_threads < 1) || (loc_map.size() == 1)) //if (num_threads < 1) { @@ -1186,6 +1200,261 @@ UpgradeThread::UpgradeThread(PerformanceLog* _performance_log, unordered_map& weights, + const Eigen::DiagonalMatrix& parcov_inv, + const vector& act_obs_names,const vector& act_par_names) +{ + class local_utils + { + public: + + static void save_names(int verbose_level, int tid, int iter, int t_count, string prefix, const vector& names) + { + if (verbose_level < 2) + return; + + if (verbose_level < 3) + return; + stringstream ss; + ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; + string fname = ss.str(); + ofstream f(fname); + if (!f.good()) + cout << "error getting ofstream " << fname << endl; + else + { + for (const auto& name : names) + { + f << name << endl; + } + + } + f.close(); + return; + } + + static void save_mat(int verbose_level, int tid, int iter, int t_count, string prefix, const Eigen::MatrixXd& mat) + { + if (verbose_level < 2) + return; + + if (verbose_level < 3) + return; + //cout << "thread: " << tid << ", " << t_count << ", " << prefix << " rows:cols" << mat.rows() << ":" << mat.cols() << endl; + stringstream ss; + + ss << "thread_" << tid << ".count_" << t_count << ".iter_" << iter << "." << prefix << ".dat"; + string fname = ss.str(); + ofstream f(fname); + if (!f.good()) + cout << "error getting ofstream " << fname << endl; + else + { + try + { + f << mat << endl; + f.close(); + } + catch (...) + { + cout << "error saving matrix " << fname << endl; + } + } + } + }; + + local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_obs_names", act_obs_names); + local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_par_names", act_par_names); + + stringstream ss; + if (!use_glm) + { + if (true) + { + // Low rank Cee. Section 14.3.2 Evenson Book + obs_err = obs_err.colwise() - obs_err.rowwise().mean(); + obs_err = sqrt(cur_lam) * obs_err; + Eigen::MatrixXd s0, V0, U0, s0_i; + SVD_REDSVD rsvd; + rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); + s0_i = s0.asDiagonal().inverse(); + Eigen::MatrixXd X0 = U0.transpose() * obs_err; + X0 = s0_i * X0; + Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; + rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); + + s1_2 = s1.cwiseProduct(s1); + s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); + Eigen::MatrixXd X1 = s0_i * U1; + X1 = U0 * X1; + + Eigen::MatrixXd X4 = s1_2i * X1.transpose(); + Eigen::MatrixXd X2 = X4 * obs_resid; + Eigen::MatrixXd X3 = X1 * X2; + + X3 = obs_diff.transpose() * X3; + upgrade_1 = -1 * par_diff * X3; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); + upgrade_1.transposeInPlace(); + + } + else + { + // Use when ensemble size is larger than number of observation. + // This is a good option when number of observation is small + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); + obs_err = obs_err.colwise() - obs_err.rowwise().mean(); + + Eigen::MatrixXd s2_, s, V, U, cum_sum; + SVD_REDSVD rsvd; + Eigen::MatrixXd C; + C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); + Eigen::VectorXd s2; + + + rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); + s2 = s.cwiseProduct(s); + s2_ = s2.asDiagonal().inverse(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); + + + Eigen::MatrixXd X1 = s2_ * U.transpose(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); + + X1 = X1 * obs_resid; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); + + X1 = U * X1; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); + + X1 = obs_diff.transpose() * X1; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); + + upgrade_1 = -1 * par_diff * X1; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); + upgrade_1.transposeInPlace(); + } + + } + else + { + + Eigen::MatrixXd ivec, s, s2, V, Ut, d_dash; + string key; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "weights", weights.toDenseMatrix()); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); + obs_resid = weights * obs_resid; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_resid", obs_resid); + + int num_reals = par_resid.cols(); + double scale = (1.0 / (sqrt(double(num_reals - 1)))); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); + obs_diff = scale * (weights * obs_diff); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_diff", obs_diff); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); + if (verbose_level > 1) { + if (parcov_inv.size() < 10000) { + + Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); + } + if (act_obs_names.size() > 0) { //this works bc the mm solve doesnt pass these names... + ss.str(""); + ss << "solution scaling factor: " << scale << endl; + ss << "eigthresh: " << eigthresh << endl; + ss << "maxsing: " << maxsing << endl; + cout << ss.str() << endl; + } + } + if (use_prior_scaling) + + par_diff = scale * parcov_inv * par_diff; + else + par_diff = scale * par_diff; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); + SVD_REDSVD rsvd; + + rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); + + Ut.transposeInPlace(); + //obs_diff.resize(0, 0); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); + + s2 = s.cwiseProduct(s); + + ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); + + Eigen::MatrixXd X1 = Ut * obs_resid; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); + + Eigen::MatrixXd X2 = ivec * X1; + //X1.resize(0, 0); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); + + Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; + //X2.resize(0, 0); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); + upgrade_1 = -1.0 * par_diff * X3; + + if (use_prior_scaling) { + //upgrade_1 = parcov_inv * upgrade_1; + } + + upgrade_1.transposeInPlace(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); + //X3.resize(0, 0); + + Eigen::MatrixXd upgrade_2; + if ((!use_approx) && (iter > 1)) { + if (use_prior_scaling) { + par_resid = parcov_inv * par_resid; + } + + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); + Eigen::MatrixXd x4 = Am.transpose() * par_resid; + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); + + par_resid.resize(0, 0); + + Eigen::MatrixXd x5 = Am * x4; + //x4.resize(0, 0); + //Am.resize(0, 0); + + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); + Eigen::MatrixXd x6 = par_diff.transpose() * x5; + //x5.resize(0, 0); + + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); + Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; + //x6.resize(0, 0); + + if (use_prior_scaling) { + upgrade_2 = -1.0 * parcov_inv * par_diff * x7; + } else { + upgrade_2 = -1.0 * (par_diff * x7); + } + //x7.resize(0, 0); + + upgrade_1 = upgrade_1 + upgrade_2.transpose(); + local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); + //upgrade_2.resize(0, 0); + + } + } + +} + void MmUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, Eigen::VectorXd parcov_inv_vec, Eigen::MatrixXd Am) { @@ -1363,67 +1632,6 @@ void MmUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm num_reals = pe_real_names.size(); - //now loop until this thread gets access to all the containers it needs to solve with -// while (true) -// { -// //if all the solution pieces are filled, break out and solve! -// if (((use_approx) || (par_resid.rows() > 0)) && -// (weights.size() > 0) && -// (par_diff.rows() > 0) && -// (obs_resid.rows() > 0) && -// (obs_err.rows() > 0) && -// (obs_diff.rows() > 0) && -// ((use_approx) || (Am.rows() > 0))) -// break; -// -// //get access to the obs_diff container -// if ((obs_diff.rows() == 0) && (obs_diff_guard.try_lock())) -// { -// //piggy back here for thread safety -// //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_svd_pack() == PestppOptions::SVD_PACK::PROPACK) -// // use_propack = true; -// obs_diff = local_utils::get_matrix_from_map(oe_real_names, obs_diff_map); -// obs_diff_guard.unlock(); -// } -// -// //get access to the residual container -// if ((obs_resid.rows() == 0) && (obs_resid_guard.try_lock())) -// { -// obs_resid = local_utils::get_matrix_from_map(oe_real_names, obs_resid_map); -// obs_resid_guard.unlock(); -// } -// -// //get access to the obs noise container -// if ((obs_err.rows() == 0) && (obs_err_guard.try_lock())) -// { -// obs_err = local_utils::get_matrix_from_map(oe_real_names, obs_err_map); -// obs_err_guard.unlock(); -// } -// -// //get access to the par diff container -// if ((par_diff.rows() == 0) && (par_diff_guard.try_lock())) -// { -// par_diff = local_utils::get_matrix_from_map(pe_real_names, par_diff_map); -// par_diff_guard.unlock(); -// } -// -// //get access to the par residual container -// if ((par_resid.rows() == 0) && (par_resid_guard.try_lock())) -// { -// par_resid = local_utils::get_matrix_from_map(pe_real_names, par_resid_map); -// par_resid_guard.unlock(); -// } -// -// //get access to the obs weights container -// if ((weights.rows() == 0) && (weight_guard.try_lock())) -// { -// weights = weight_map.at(key).asDiagonal(); -// weight_guard.unlock(); -// } -// -// -// } - obs_diff = local_utils::get_matrix_from_map(oe_real_names, obs_diff_map); obs_resid = local_utils::get_matrix_from_map(oe_real_names, obs_resid_map); obs_err = local_utils::get_matrix_from_map(oe_real_names, obs_err_map); @@ -1439,7 +1647,6 @@ void MmUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm for (int i=0;i empty_obs_names,empty_par_names; + UpgradeThread::ensemble_solution(iter,verbose_level,maxsing,thread_id,t_count, use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid,par_diff,Am,obs_resid, + obs_diff,upgrade_1,obs_err,weights,parcov_inv,empty_obs_names,empty_par_names); - - Eigen::MatrixXd X1 = s2_ * U.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - X1 = X1 * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); - - X1 = U * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); - - X1 = obs_diff.transpose() * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); - - upgrade_1 = -1 * par_diff * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - upgrade_1.transposeInPlace(); - } - - - } - - -// ---------------------------------- -// glm solution -// ---------------------------------- - else - { - - obs_resid = weights * obs_resid; - obs_diff = scale * (weights * obs_diff); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); - if (use_prior_scaling) - par_diff = scale * parcov_inv * par_diff; - else - par_diff = scale * par_diff; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); - - Ut.transposeInPlace(); - obs_diff.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - - Eigen::MatrixXd s2 = s.cwiseProduct(s); - - ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); - - Eigen::MatrixXd X1 = Ut * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - Eigen::MatrixXd X2 = ivec * X1; - X1.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); - - Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; - X2.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); - upgrade_1 = -1.0 * par_diff * X3; - - if (use_prior_scaling) - { - //upgrade_1 = parcov_inv * upgrade_1; - } - - upgrade_1.transposeInPlace(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - X3.resize(0, 0); - - Eigen::MatrixXd upgrade_2; - if ((!use_approx) && (iter > 1)) - { - if (use_prior_scaling) - { - par_resid = parcov_inv * par_resid; - } - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); - Eigen::MatrixXd x4 = Am.transpose() * par_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); - - par_resid.resize(0, 0); - - Eigen::MatrixXd x5 = Am * x4; - x4.resize(0, 0); - //Am.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); - Eigen::MatrixXd x6 = par_diff.transpose() * x5; - x5.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); - Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; - x6.resize(0, 0); - - if (use_prior_scaling) - { - upgrade_2 = -1.0 * parcov_inv * par_diff * x7; - } - else - { - upgrade_2 = -1.0 * (par_diff * x7); - } - x7.resize(0, 0); - - upgrade_1 = upgrade_1 + upgrade_2.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); - upgrade_2.resize(0, 0); - - } - } +// Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, d_dash; +// +// //---------------------------------- +// //es-mda solution +// //---------------------------------- +// +// if (!use_glm_form) +// { +// if (true) +// { +// // Low rank Cee. Section 14.3.2 Evenson Book +// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); +// obs_err = sqrt(cur_lam) * obs_err; +// Eigen::MatrixXd s0, V0, U0, s0_i; +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); +// s0_i = s0.asDiagonal().inverse(); +// Eigen::MatrixXd X0 = U0.transpose() * obs_err; +// X0 = s0_i * X0; +// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; +// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); +// +// s1_2 = s1.cwiseProduct(s1); +// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); +// Eigen::MatrixXd X1 = s0_i * U1; +// X1 = U0 * X1; +// +// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); +// Eigen::MatrixXd X2 = X4 * obs_resid; +// Eigen::MatrixXd X3 = X1 * X2; +// +// X3 = obs_diff.transpose() * X3; +// upgrade_1 = -1 * par_diff * X3; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// upgrade_1.transposeInPlace(); +// +// } +// else +// { +// // Use when ensemble size is larger than number of observation. +// // This is a good option when number of observation is small +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); +// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); +// +// Eigen::MatrixXd s2_, s, V, U, cum_sum; +// SVD_REDSVD rsvd; +// Eigen::MatrixXd C; +// C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); +// Eigen::VectorXd s2; +// +// +// rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); +// s2 = s.cwiseProduct(s); +// s2_ = s2.asDiagonal().inverse(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); +// +// +// Eigen::MatrixXd X1 = s2_ * U.transpose(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); +// +// X1 = X1 * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); +// +// X1 = U * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); +// +// X1 = obs_diff.transpose() * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); +// +// upgrade_1 = -1 * par_diff * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// upgrade_1.transposeInPlace(); +// } +// +// +// } +// +// +//// ---------------------------------- +//// glm solution +//// ---------------------------------- +// else +// { +// +// obs_resid = weights * obs_resid; +// obs_diff = scale * (weights * obs_diff); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); +// if (use_prior_scaling) +// par_diff = scale * parcov_inv * par_diff; +// else +// par_diff = scale * par_diff; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); +// +// Ut.transposeInPlace(); +// obs_diff.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); +// +// Eigen::MatrixXd s2 = s.cwiseProduct(s); +// +// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); +// +// Eigen::MatrixXd X1 = Ut * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); +// +// Eigen::MatrixXd X2 = ivec * X1; +// X1.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); +// +// Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; +// X2.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); +// upgrade_1 = -1.0 * par_diff * X3; +// +// if (use_prior_scaling) +// { +// //upgrade_1 = parcov_inv * upgrade_1; +// } +// +// upgrade_1.transposeInPlace(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// X3.resize(0, 0); +// +// Eigen::MatrixXd upgrade_2; +// if ((!use_approx) && (iter > 1)) +// { +// if (use_prior_scaling) +// { +// par_resid = parcov_inv * par_resid; +// } +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); +// Eigen::MatrixXd x4 = Am.transpose() * par_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); +// +// par_resid.resize(0, 0); +// +// Eigen::MatrixXd x5 = Am * x4; +// x4.resize(0, 0); +// //Am.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); +// Eigen::MatrixXd x6 = par_diff.transpose() * x5; +// x5.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); +// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; +// x6.resize(0, 0); +// +// if (use_prior_scaling) +// { +// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; +// } +// else +// { +// upgrade_2 = -1.0 * (par_diff * x7); +// } +// x7.resize(0, 0); +// +// upgrade_1 = upgrade_1 + upgrade_2.transpose(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); +// upgrade_2.resize(0, 0); +// +// } +// } //assuming that the fist row is the realization we are after... row_vec = upgrade_1.row(0); @@ -1683,567 +1892,510 @@ MmUpgradeThread::MmUpgradeThread(PerformanceLog* _performance_log, unordered_map } -void CovLocalizationUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, - vector obs_names) -{ - - //declare these helpers in here so they are thread safe... - class local_utils - { - public: - static Eigen::DiagonalMatrix get_matrix_from_map(vector& names, unordered_map& dmap) - { - Eigen::VectorXd vec(names.size()); - int i = 0; - for (auto name : names) - { - vec[i] = dmap.at(name); - ++i; - } - Eigen::DiagonalMatrix m = vec.asDiagonal(); - return m; - } - static Eigen::MatrixXd get_matrix_from_map(int num_reals, vector& names, unordered_map& emap) - { - Eigen::MatrixXd mat(num_reals, names.size()); - mat.setZero(); - - for (int j = 0; j < names.size(); j++) - { - mat.col(j) = emap[names[j]]; - } - - return mat; - } - static void save_mat(int verbose_level, int tid, int iter, int t_count, string prefix, Eigen::MatrixXd& mat) - { - if (verbose_level < 2) - return; - - if (verbose_level < 3) - return; - //cout << "thread: " << tid << ", " << t_count << ", " << prefix << " rows:cols" << mat.rows() << ":" << mat.cols() << endl; - stringstream ss; - - ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; - string fname = ss.str(); - ofstream f(fname); - if (!f.good()) - cout << "error getting ofstream " << fname << endl; - else - { - - try - { - f << mat << endl; - f.close(); - } - catch (...) - { - cout << "error saving matrix " << fname << endl; - } - } - } - }; - - stringstream ss; - - unique_lock ctrl_guard(ctrl_lock, defer_lock); - int maxsing, num_reals, verbose_level, pcount = 0, t_count=0; - double eigthresh; - bool use_approx; - bool use_prior_scaling; - bool use_localizer = false; - bool loc_by_obs = true; - - while (true) - { - if (ctrl_guard.try_lock()) - { - - maxsing = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().maxsing; - eigthresh = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().eigthresh; - use_approx = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_approx(); - use_prior_scaling = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_prior_scaling(); - num_reals = pe_upgrade.shape().first; - verbose_level = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_verbose_level(); - ctrl_guard.unlock(); - //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_localize_how()[0] == 'P') - if (how == Localizer::How::PARAMETERS) - loc_by_obs = false; - else - { - throw runtime_error("Covariance localization only supporte for localization by parameters..."); - } - break; - } - } - ofstream f_thread; - if (verbose_level > 2) - { - ss.str(""); - ss << "thread_" << thread_id << "part_map.csv"; - f_thread.open(ss.str()); - ss.str(""); - } - Eigen::MatrixXd par_resid, par_diff, Am; - Eigen::MatrixXd obs_resid, obs_diff, obs_err, loc; - Eigen::DiagonalMatrix weights, parcov_inv; - vector case_par_names, case_obs_names; - string key; - - - //these locks are used to control (thread-safe) access to the fast look up containers - unique_lock next_guard(next_lock, defer_lock); - unique_lock obs_diff_guard(obs_diff_lock, defer_lock); - unique_lock obs_resid_guard(obs_resid_lock, defer_lock); - unique_lock obs_err_guard(obs_err_lock, defer_lock); - unique_lock par_diff_guard(par_diff_lock, defer_lock); - unique_lock par_resid_guard(par_resid_lock, defer_lock); - unique_lock loc_guard(loc_lock, defer_lock); - unique_lock weight_guard(weight_lock, defer_lock); - unique_lock parcov_guard(parcov_lock, defer_lock); - unique_lock am_guard(am_lock, defer_lock); - - //reset all the solution parts - par_resid.resize(0, 0); - par_diff.resize(0, 0); - obs_resid.resize(0, 0); - obs_diff.resize(0, 0); - obs_err.resize(0, 0); - loc.resize(0, 0); - Am.resize(0, 0); - weights.resize(0); - parcov_inv.resize(0); - Am.resize(0, 0); - - - //loop until this thread gets access to all the containers it needs to solve with - //which in this solution, is all active pars and obs... - -// while (true) +//void CovLocalizationUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, +// vector obs_names) +//{ +// +// //declare these helpers in here so they are thread safe... +// class local_utils // { -// //if all the solution pieces are filled, break out and solve! -// if (((use_approx) || (par_resid.rows() > 0)) && -// (weights.size() > 0) && -// (parcov_inv.size() > 0) && -// (par_diff.rows() > 0) && -// (obs_resid.rows() > 0) && -// (obs_err.rows() > 0) && -// (obs_diff.rows() > 0) && -// ((!use_localizer) || (loc.rows() > 0)) && -// ((use_approx) || (Am.rows() > 0))) -// break; +// public: // +// static void save_names(int verbose_level, int tid, int iter, int t_count, string prefix, vector& names) +// { +// if (verbose_level < 2) +// return; // -// //get access to the obs_diff container -// if ((obs_diff.rows() == 0) && (obs_diff_guard.try_lock())) -// { -// //piggy back here for thread safety -// //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_svd_pack() == PestppOptions::SVD_PACK::PROPACK) -// // use_propack = true; -// obs_diff = local_utils::get_matrix_from_map(num_reals, obs_names, obs_diff_map); -// obs_diff_guard.unlock(); -// } +// if (verbose_level < 3) +// return; +// stringstream ss; // -// //get access to the residual container -// if ((obs_resid.rows() == 0) && (obs_resid_guard.try_lock())) -// { -// obs_resid = local_utils::get_matrix_from_map(num_reals, obs_names, obs_resid_map); -// obs_resid_guard.unlock(); -// } +// ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; +// string fname = ss.str(); +// ofstream f(fname); +// if (!f.good()) +// cout << "error getting ofstream " << fname << endl; +// else +// { +// for (const auto& name : names) +// { +// f << name << endl; +// } // -// //get access to the obs noise container -// if ((obs_err.rows() == 0) && (obs_err_guard.try_lock())) +// } +// f.close(); +// return; +// } +// static Eigen::DiagonalMatrix get_matrix_from_map(vector& names, unordered_map& dmap) // { -// obs_err = local_utils::get_matrix_from_map(num_reals, obs_names, obs_err_map); -// obs_err_guard.unlock(); +// Eigen::VectorXd vec(names.size()); +// int i = 0; +// for (auto name : names) +// { +// vec[i] = dmap.at(name); +// ++i; +// } +// Eigen::DiagonalMatrix m = vec.asDiagonal(); +// return m; // } -// -// //get access to the par diff container -// if ((par_diff.rows() == 0) && (par_diff_guard.try_lock())) +// static Eigen::MatrixXd get_matrix_from_map(int num_reals, vector& names, unordered_map& emap) // { -// par_diff = local_utils::get_matrix_from_map(num_reals, par_names, par_diff_map); -// par_diff_guard.unlock(); -// } +// Eigen::MatrixXd mat(num_reals, names.size()); +// mat.setZero(); // -// //get access to the par residual container -// if ((par_resid.rows() == 0) && (par_resid_guard.try_lock())) -// { -// par_resid = local_utils::get_matrix_from_map(num_reals, par_names, par_resid_map); -// par_resid_guard.unlock(); -// } +// for (int j = 0; j < names.size(); j++) +// { +// mat.col(j) = emap[names[j]]; +// } // -// //get access to the obs weights container -// if ((weights.rows() == 0) && (weight_guard.try_lock())) -// { -// weights = local_utils::get_matrix_from_map(obs_names, weight_map); -// weight_guard.unlock(); +// return mat; // } -// -// //get access to the inverse prior parcov -// if ((parcov_inv.rows() == 0) && (parcov_guard.try_lock())) +// static void save_mat(int verbose_level, int tid, int iter, int t_count, string prefix, Eigen::MatrixXd& mat) // { -// parcov_inv = local_utils::get_matrix_from_map(par_names, parcov_inv_map); -// parcov_guard.unlock(); +// if (verbose_level < 2) +// return; +// +// if (verbose_level < 3) +// return; +// //cout << "thread: " << tid << ", " << t_count << ", " << prefix << " rows:cols" << mat.rows() << ":" << mat.cols() << endl; +// stringstream ss; +// +// ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; +// string fname = ss.str(); +// ofstream f(fname); +// if (!f.good()) +// cout << "error getting ofstream " << fname << endl; +// else +// { +// +// try +// { +// f << mat << endl; +// f.close(); +// } +// catch (...) +// { +// cout << "error saving matrix " << fname << endl; +// } +// } // } +// }; +// +// stringstream ss; // -// //if needed, get access to the Am container - needed in the full glm solution -// if ((!use_approx) && (Am.rows() == 0) && (am_guard.try_lock())) +// unique_lock ctrl_guard(ctrl_lock, defer_lock); +// int maxsing, num_reals, verbose_level, pcount = 0, t_count=0; +// double eigthresh; +// bool use_approx; +// bool use_prior_scaling; +// bool use_localizer = false; +// bool loc_by_obs = true; +// +// while (true) +// { +// if (ctrl_guard.try_lock()) // { -// //Am = local_utils::get_matrix_from_map(num_reals, par_names, Am_map).transpose(); -// int am_cols = Am_map[par_names[0]].size(); -// Am.resize(par_names.size(), am_cols); -// Am.setZero(); // -// for (int j = 0; j < par_names.size(); j++) +// maxsing = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().maxsing; +// eigthresh = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().eigthresh; +// use_approx = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_approx(); +// use_prior_scaling = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_prior_scaling(); +// num_reals = pe_upgrade.shape().first; +// verbose_level = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_verbose_level(); +// ctrl_guard.unlock(); +// //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_localize_how()[0] == 'P') +// if (how == Localizer::How::PARAMETERS) +// loc_by_obs = false; +// else // { -// Am.row(j) = Am_map[par_names[j]]; +// throw runtime_error("Covariance localization only supporte for localization by parameters..."); // } -// am_guard.unlock(); +// break; // } // } +// ofstream f_thread; +// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_obs_names",obs_names); +// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_par_names",par_names); +// if (verbose_level > 2) +// { +// ss.str(""); +// ss << "thread_" << thread_id << ".part_map.csv"; +// f_thread.open(ss.str()); +// ss.str(""); +// } +// Eigen::MatrixXd par_resid, par_diff, Am; +// Eigen::MatrixXd obs_resid, obs_diff, obs_err, loc; +// Eigen::DiagonalMatrix weights, parcov_inv; +// vector case_par_names, case_obs_names; +// string key; +// +// +// //these locks are used to control (thread-safe) access to the fast look up containers +// unique_lock next_guard(next_lock, defer_lock); +// unique_lock obs_diff_guard(obs_diff_lock, defer_lock); +// unique_lock obs_resid_guard(obs_resid_lock, defer_lock); +// unique_lock obs_err_guard(obs_err_lock, defer_lock); +// unique_lock par_diff_guard(par_diff_lock, defer_lock); +// unique_lock par_resid_guard(par_resid_lock, defer_lock); +// unique_lock loc_guard(loc_lock, defer_lock); +// unique_lock weight_guard(weight_lock, defer_lock); +// unique_lock parcov_guard(parcov_lock, defer_lock); +// unique_lock am_guard(am_lock, defer_lock); +// +// //reset all the solution parts +// par_resid.resize(0, 0); +// par_diff.resize(0, 0); +// obs_resid.resize(0, 0); +// obs_diff.resize(0, 0); +// obs_err.resize(0, 0); +// loc.resize(0, 0); +// Am.resize(0, 0); +// weights.resize(0); +// parcov_inv.resize(0); +// Am.resize(0, 0); +// +// +// +// obs_diff = local_utils::get_matrix_from_map(num_reals, obs_names, obs_diff_map); +// obs_resid = local_utils::get_matrix_from_map(num_reals, obs_names, obs_resid_map); +// obs_err = local_utils::get_matrix_from_map(num_reals, obs_names, obs_err_map); +// par_diff = local_utils::get_matrix_from_map(num_reals, par_names, par_diff_map); +// par_resid = local_utils::get_matrix_from_map(num_reals, par_names, par_resid_map); +// weights = local_utils::get_matrix_from_map(obs_names, weight_map); +// parcov_inv = local_utils::get_matrix_from_map(par_names, parcov_inv_map); +// +// if ((!use_approx) && (Am.rows() == 0)) +// { +// //Am = local_utils::get_matrix_from_map(num_reals, par_names, Am_map).transpose(); +// int am_cols = Am_map[par_names[0]].size(); +// Am.resize(par_names.size(), am_cols); +// Am.setZero(); +// +// for (int j = 0; j < par_names.size(); j++) +// { +// Am.row(j) = Am_map[par_names[j]]; +// } +// } +// +// par_diff.transposeInPlace(); +// obs_diff.transposeInPlace(); +// obs_resid.transposeInPlace(); +// par_resid.transposeInPlace(); +// obs_err.transposeInPlace(); +// +// +// //container to quickly look indices +// map par2col_map; +// for (int i = 0; i < par_names.size(); i++) +// par2col_map[par_names[i]] = i; +// +// map obs2row_map; +// for (int i = 0; i < obs_names.size(); i++) +// obs2row_map[obs_names[i]] = i; +// +// +// //form the scaled obs resid matrix +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); +// //Eigen::MatrixXd scaled_residual = weights * obs_resid; +// +// //form the (optionally) scaled par resid matrix +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); +// +// ss.str(""); +// double scale = (1.0 / (sqrt(double(num_reals - 1)))); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); +// +// //calculate some full solution components first... +// +// Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, t; +// Eigen::MatrixXd upgrade_2; +// Eigen::VectorXd loc_vec; +// +// upgrade_2.resize(num_reals, pe_upgrade.shape().second); +// upgrade_2.setZero(); +// +// if (!use_glm_form) +// { +// if (true) +// { +// // Low rank Cee. Section 14.3.2 Evenson Book +// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); +// obs_err = sqrt(cur_lam) * obs_err; +// Eigen::MatrixXd s0, V0, U0, s0_i; +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); +// s0_i = s0.asDiagonal().inverse(); +// Eigen::MatrixXd X0 = U0.transpose() * obs_err; +// X0 = s0_i * X0; +// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; +// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); +// +// s1_2 = s1.cwiseProduct(s1); +// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); +// Eigen::MatrixXd X1 = s0_i * U1; +// X1 = U0 * X1; +// +// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); +// Eigen::MatrixXd X2 = X4 * obs_resid; +// Eigen::MatrixXd X3 = X1 * X2; +// +// X3 = obs_diff.transpose() * X3; +// //upgrade_1 = -1 * par_diff * X3; +// //local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// //upgrade_1.transposeInPlace(); +// t = obs_diff * X3; +// Ut.resize(0, 0); +// obs_diff.resize(0, 0); +// X3.resize(0,0); +// X2.resize(0,0); +// X4.resize(0,0); +// X1.resize(0,0); +// +// } +// else { +// obs_diff = scale * obs_diff;// (H-Hm)/sqrt(N-1) +// par_diff = scale * par_diff;// (K-Km)/sqrt(N-1) +// obs_err = scale * obs_err; // (E-Em)/sqrt(N-1) +// +// SVD_REDSVD rsvd; +// Eigen::MatrixXd C = obs_diff + (cur_lam * obs_err); // curr_lam is the inflation factor +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); +// rsvd.solve_ip(C, s, Ut, V, eigthresh, maxsing); +// Ut.transposeInPlace(); +// V.resize(0, 0); +// C.resize(0, 0); +// obs_err.resize(0, 0); +// +// // s2 = s.asDiagonal().inverse(); +// s2 = s.cwiseProduct(s).asDiagonal().inverse(); +// for (int i = 0; i < s.size(); i++) { +// if (s(i) < 1e-50) { +// s2(i, i) = 0; +// } +// } +// +// t = obs_diff.transpose() * Ut.transpose() * s2 * Ut; +// Ut.resize(0, 0); +// obs_diff.resize(0, 0); +// } +// } +// +// +// //---------------------------------- +// //glm solution +// //---------------------------------- +// else +// { +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); +// obs_diff = scale * (weights * obs_diff); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_diff", obs_diff); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); +// if (verbose_level > 1) { +// if (parcov_inv.size() < 10000) { +// Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); +// } +// ss.str(""); +// ss << "solution scaling factor: " << scale; +// cout << ss.str() << endl; +// } +// if (use_prior_scaling) +// par_diff = scale * parcov_inv * par_diff; +// else +// par_diff = scale * par_diff; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); +// +// Ut.transposeInPlace(); +// obs_diff.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); +// +// Eigen::MatrixXd s2 = s.cwiseProduct(s); +// +// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); +// +// obs_resid = weights * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_resid", obs_resid); +// t = V * s.asDiagonal() * ivec * Ut; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "t", t); +// +// if ((!use_approx) && (iter > 1)) +// { +// if (use_prior_scaling) +// { +// par_resid = parcov_inv * par_resid; +// } +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); +// Eigen::MatrixXd x4 = Am.transpose() * par_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); +// +// par_resid.resize(0, 0); +// +// Eigen::MatrixXd x5 = Am * x4; +// x4.resize(0, 0); +// Am.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); +// Eigen::MatrixXd x6 = par_diff.transpose() * x5; +// x5.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); +// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; +// x6.resize(0, 0); +// +// if (use_prior_scaling) +// { +// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; +// } +// else +// { +// upgrade_2 = -1.0 * (par_diff * x7); +// } +// x7.resize(0, 0); +// +// //add upgrade_2 piece +// //upgrade_1 = upgrade_1 + upgrade_2.transpose(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); +// //upgrade_2.resize(0, 0); +// } +// } +// +// //This is the main thread loop - it continues until all upgrade pieces have been completed +// map loc_map; +// //solve for each case par_name +// Eigen::VectorXd pt = parcov_inv.diagonal(); +// string name; +// while (true) +// { +// //clear the pieces used last time... +// key = ""; +// loc_map.clear(); +// +// //get new pieces... +// while (true) +// { +// if ((key != "") && (loc_map.size() > 0)) +// break; +// if (next_guard.try_lock()) +// { +// //if all the pieces have been completed, return +// if (count == keys.size()) +// { +// if (verbose_level > 3) +// { +// cout << "upgrade thread: " << thread_id << " processed " << pcount << " upgrade parts" << endl; +// } +// if (f_thread.good()) +// f_thread.close(); +// return; +// } +// key = keys[count]; +// pair, vector> p = cases.at(key); +// //we can potentially optimize the speed of this loc type by changing how many pars are +// //passed in each "case" so herein, we support a generic number of pars per case... +// case_par_names = p.second; +// //In this solution, we ignore case obs names since we are using the full set of obs for the solution... +// case_obs_names = p.first; +// +// if (count % 1000 == 0) +// { +// ss.str(""); +// ss << "upgrade thread progress: " << count << " of " << total << " parts done"; +// if (verbose_level > 3) +// cout << ss.str() << endl; +// performance_log->log_event(ss.str()); +// } +// count++; +// t_count = count; +// pcount++; +// next_guard.unlock(); +// +// } +// //get access to the localizer +// if ((key != "") && (loc_map.size() == 0) && (loc_guard.try_lock())) +// { +// //get the nobs-length localizing vector for each case par name +// for (auto& par_name : case_par_names) +// { +// loc_map[par_name] = localizer.get_obs_hadamard_vector(par_name, obs_names); +// } +// loc_guard.unlock(); +// } +// } +// +// if (verbose_level > 2) +// { +// f_thread << t_count << "," << iter; +// for (auto name : case_par_names) +// f_thread << "," << name; +// for (auto name : case_obs_names) +// f_thread << "," << name; +// f_thread << endl; +// } +// upgrade_1.resize(num_reals,case_par_names.size()); +// upgrade_1.setZero(); +// +// for (int i = 0; i < case_par_names.size(); i++) +// { +// name = case_par_names[i]; +// loc_vec = loc_map[name]; +// Eigen::VectorXd par_vec = par_diff.row(par2col_map[name]) * t; +// //apply the localizer +// par_vec = par_vec.cwiseProduct(loc_vec); +// if (!use_glm_form) +// par_vec = -1.0 * par_vec.transpose() * obs_resid; +// else +// { +// par_vec = -1.0 * par_vec.transpose() * obs_resid; +// //if (use_prior_scaling) +// // par_vec *= pt[i]; +// } +// upgrade_1.col(i) += par_vec.transpose(); +// //add the par change part for the full glm solution +// if ((use_glm_form) && (!use_approx) && (iter > 1)) +// { +// //update_2 is transposed relative to upgrade_1 +// upgrade_1.col(i) += upgrade_2.row(par2col_map[name]); +// } +// +// } +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// +// //put this piece of the upgrade vector in +// unique_lock put_guard(put_lock, defer_lock); +// while (true) +// { +// if (put_guard.try_lock()) +// { +// pe_upgrade.add_2_cols_ip(case_par_names, upgrade_1); +// put_guard.unlock(); +// break; +// } +// } +// } +//} - obs_diff = local_utils::get_matrix_from_map(num_reals, obs_names, obs_diff_map); - obs_resid = local_utils::get_matrix_from_map(num_reals, obs_names, obs_resid_map); - obs_err = local_utils::get_matrix_from_map(num_reals, obs_names, obs_err_map); - par_diff = local_utils::get_matrix_from_map(num_reals, par_names, par_diff_map); - par_resid = local_utils::get_matrix_from_map(num_reals, par_names, par_resid_map); - weights = local_utils::get_matrix_from_map(obs_names, weight_map); - parcov_inv = local_utils::get_matrix_from_map(par_names, parcov_inv_map); - - if ((!use_approx) && (Am.rows() == 0)) - { - //Am = local_utils::get_matrix_from_map(num_reals, par_names, Am_map).transpose(); - int am_cols = Am_map[par_names[0]].size(); - Am.resize(par_names.size(), am_cols); - Am.setZero(); - - for (int j = 0; j < par_names.size(); j++) - { - Am.row(j) = Am_map[par_names[j]]; - } - } - - par_diff.transposeInPlace(); - obs_diff.transposeInPlace(); - obs_resid.transposeInPlace(); - par_resid.transposeInPlace(); - obs_err.transposeInPlace(); - - - //container to quickly look indices - map par2col_map; - for (int i = 0; i < par_names.size(); i++) - par2col_map[par_names[i]] = i; - - map obs2row_map; - for (int i = 0; i < obs_names.size(); i++) - obs2row_map[obs_names[i]] = i; - - - //form the scaled obs resid matrix - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); - //Eigen::MatrixXd scaled_residual = weights * obs_resid; - - //form the (optionally) scaled par resid matrix - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); - - ss.str(""); - double scale = (1.0 / (sqrt(double(num_reals - 1)))); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); - - //calculate some full solution components first... - - Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, t; - Eigen::MatrixXd upgrade_2; - Eigen::VectorXd loc_vec; - - upgrade_2.resize(num_reals, pe_upgrade.shape().second); - upgrade_2.setZero(); - - if (!use_glm_form) - { - if (true) - { - // Low rank Cee. Section 14.3.2 Evenson Book - obs_err = obs_err.colwise() - obs_err.rowwise().mean(); - obs_err = sqrt(cur_lam) * obs_err; - Eigen::MatrixXd s0, V0, U0, s0_i; - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); - s0_i = s0.asDiagonal().inverse(); - Eigen::MatrixXd X0 = U0.transpose() * obs_err; - X0 = s0_i * X0; - Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; - rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); - - s1_2 = s1.cwiseProduct(s1); - s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); - Eigen::MatrixXd X1 = s0_i * U1; - X1 = U0 * X1; - - Eigen::MatrixXd X4 = s1_2i * X1.transpose(); - Eigen::MatrixXd X2 = X4 * obs_resid; - Eigen::MatrixXd X3 = X1 * X2; - - X3 = obs_diff.transpose() * X3; - //upgrade_1 = -1 * par_diff * X3; - //local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - //upgrade_1.transposeInPlace(); - t = obs_diff * X3; - Ut.resize(0, 0); - obs_diff.resize(0, 0); - X3.resize(0,0); - X2.resize(0,0); - X4.resize(0,0); - X1.resize(0,0); - - } - else { - obs_diff = scale * obs_diff;// (H-Hm)/sqrt(N-1) - par_diff = scale * par_diff;// (K-Km)/sqrt(N-1) - obs_err = scale * obs_err; // (E-Em)/sqrt(N-1) - - SVD_REDSVD rsvd; - Eigen::MatrixXd C = obs_diff + (cur_lam * obs_err); // curr_lam is the inflation factor - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); - rsvd.solve_ip(C, s, Ut, V, eigthresh, maxsing); - Ut.transposeInPlace(); - V.resize(0, 0); - C.resize(0, 0); - obs_err.resize(0, 0); - - // s2 = s.asDiagonal().inverse(); - s2 = s.cwiseProduct(s).asDiagonal().inverse(); - for (int i = 0; i < s.size(); i++) { - if (s(i) < 1e-50) { - s2(i, i) = 0; - } - } - - t = obs_diff.transpose() * Ut.transpose() * s2 * Ut; - Ut.resize(0, 0); - obs_diff.resize(0, 0); - } - } - - - //---------------------------------- - //glm solution - //---------------------------------- - else - { - - obs_diff = scale * (weights * obs_diff); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); - if (verbose_level > 1) { - Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp ); - ss.str(""); - ss << "solution scaling factor: " << scale; - performance_log->log_event(ss.str()); - } - if (use_prior_scaling) - par_diff = scale * parcov_inv * par_diff; - else - par_diff = scale * par_diff; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); - - Ut.transposeInPlace(); - obs_diff.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - - Eigen::MatrixXd s2 = s.cwiseProduct(s); - - ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); - - obs_resid = weights * obs_resid; - t = V * s.asDiagonal() * ivec * Ut; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "t", t); - - if ((!use_approx) && (iter > 1)) - { - if (use_prior_scaling) - { - par_resid = parcov_inv * par_resid; - } - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); - Eigen::MatrixXd x4 = Am.transpose() * par_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); - - par_resid.resize(0, 0); - - Eigen::MatrixXd x5 = Am * x4; - x4.resize(0, 0); - Am.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); - Eigen::MatrixXd x6 = par_diff.transpose() * x5; - x5.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); - Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; - x6.resize(0, 0); - - if (use_prior_scaling) - { - upgrade_2 = -1.0 * parcov_inv * par_diff * x7; - } - else - { - upgrade_2 = -1.0 * (par_diff * x7); - } - x7.resize(0, 0); - - //add upgrade_2 piece - //upgrade_1 = upgrade_1 + upgrade_2.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); - //upgrade_2.resize(0, 0); - } - } - - //This is the main thread loop - it continues until all upgrade pieces have been completed - map loc_map; - //solve for each case par_name - Eigen::VectorXd pt = parcov_inv.diagonal(); - string name; - while (true) - { - //clear the pieces used last time... - key = ""; - loc_map.clear(); - - //get new pieces... - while (true) - { - if ((key != "") && (loc_map.size() > 0)) - break; - if (next_guard.try_lock()) - { - //if all the pieces have been completed, return - if (count == keys.size()) - { - if (verbose_level > 3) - { - cout << "upgrade thread: " << thread_id << " processed " << pcount << " upgrade parts" << endl; - } - if (f_thread.good()) - f_thread.close(); - return; - } - key = keys[count]; - pair, vector> p = cases.at(key); - //we can potentially optimize the speed of this loc type by changing how many pars are - //passed in each "case" so herein, we support a generic number of pars per case... - case_par_names = p.second; - //In this solution, we ignore case obs names since we are using the full set of obs for the solution... - case_obs_names = p.first; - - if (count % 1000 == 0) - { - ss.str(""); - ss << "upgrade thread progress: " << count << " of " << total << " parts done"; - if (verbose_level > 3) - cout << ss.str() << endl; - performance_log->log_event(ss.str()); - } - count++; - t_count = count; - pcount++; - next_guard.unlock(); - - } - //get access to the localizer - if ((key != "") && (loc_map.size() == 0) && (loc_guard.try_lock())) - { - //get the nobs-length localizing vector for each case par name - for (auto& par_name : case_par_names) - { - loc_map[par_name] = localizer.get_obs_hadamard_vector(par_name, obs_names); - } - loc_guard.unlock(); - } - } - - if (verbose_level > 2) - { - f_thread << t_count << "," << iter; - for (auto name : case_par_names) - f_thread << "," << name; - for (auto name : case_obs_names) - f_thread << "," << name; - f_thread << endl; - } - upgrade_1.resize(num_reals,case_par_names.size()); - upgrade_1.setZero(); - - for (int i = 0; i < case_par_names.size(); i++) - { - name = case_par_names[i]; - loc_vec = loc_map[name]; - Eigen::VectorXd par_vec = par_diff.row(par2col_map[name]) * t; - //apply the localizer - par_vec = par_vec.cwiseProduct(loc_vec); - if (!use_glm_form) - par_vec = -1.0 * par_vec.transpose() * obs_resid; - else - { - par_vec = -1.0 * par_vec.transpose() * obs_resid; - //if (use_prior_scaling) - // par_vec *= pt[i]; - } - upgrade_1.col(i) += par_vec.transpose(); - //add the par change part for the full glm solution - if ((use_glm_form) && (!use_approx) && (iter > 1)) - { - //update_2 is transposed relative to upgrade_1 - upgrade_1.col(i) += upgrade_2.row(par2col_map[name]); - } - - } - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - - //put this piece of the upgrade vector in - unique_lock put_guard(put_lock, defer_lock); - while (true) - { - if (put_guard.try_lock()) - { - pe_upgrade.add_2_cols_ip(case_par_names, upgrade_1); - put_guard.unlock(); - break; - } - } - } -} - - -void ensemble_solution(const int iter, const int verbose_level,const int maxsing, const int thread_id, const int t_count, const bool - use_prior_scaling,const bool use_approx, const bool use_glm, const double cur_lam, - const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff, - const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1, - Eigen::MatrixXd& obs_err, const Eigen::DiagonalMatrix& weights, - const Eigen::DiagonalMatrix& parcov_inv) +void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, + vector par_names, vector obs_names) { - class local_utils - { - public: - static void save_mat(int verbose_level, int tid, int iter, int t_count, string prefix, const Eigen::MatrixXd& mat) + + //declare these helpers in here so they are thread safe... + class local_utils + { + public: + static void save_names(int verbose_level, int tid, int iter, int t_count, string prefix, vector& names) { if (verbose_level < 2) return; if (verbose_level < 3) return; - //cout << "thread: " << tid << ", " << t_count << ", " << prefix << " rows:cols" << mat.rows() << ":" << mat.cols() << endl; stringstream ss; ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; @@ -2253,198 +2405,16 @@ void ensemble_solution(const int iter, const int verbose_level,const int maxsing cout << "error getting ofstream " << fname << endl; else { - - try - { - f << mat << endl; - f.close(); - } - catch (...) + for (const auto& name : names) { - cout << "error saving matrix " << fname << endl; + f << name << endl; } - } - } - }; - - stringstream ss; - if (!use_glm) - { - if (true) - { - // Low rank Cee. Section 14.3.2 Evenson Book - obs_err = obs_err.colwise() - obs_err.rowwise().mean(); - obs_err = sqrt(cur_lam) * obs_err; - Eigen::MatrixXd s0, V0, U0, s0_i; - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); - s0_i = s0.asDiagonal().inverse(); - Eigen::MatrixXd X0 = U0.transpose() * obs_err; - X0 = s0_i * X0; - Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; - rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); - - s1_2 = s1.cwiseProduct(s1); - s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); - Eigen::MatrixXd X1 = s0_i * U1; - X1 = U0 * X1; - - Eigen::MatrixXd X4 = s1_2i * X1.transpose(); - Eigen::MatrixXd X2 = X4 * obs_resid; - Eigen::MatrixXd X3 = X1 * X2; - - X3 = obs_diff.transpose() * X3; - upgrade_1 = -1 * par_diff * X3; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - upgrade_1.transposeInPlace(); - - } - else - { - // Use when ensemble size is larger than number of observation. - // This is a good option when number of observation is small - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); - obs_err = obs_err.colwise() - obs_err.rowwise().mean(); - - Eigen::MatrixXd s2_, s, V, U, cum_sum; - SVD_REDSVD rsvd; - Eigen::MatrixXd C; - C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); - Eigen::VectorXd s2; - - - rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - s2 = s.cwiseProduct(s); - s2_ = s2.asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); - - - Eigen::MatrixXd X1 = s2_ * U.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - X1 = X1 * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); - - X1 = U * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); - - X1 = obs_diff.transpose() * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); - - upgrade_1 = -1 * par_diff * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - upgrade_1.transposeInPlace(); - } - - } - else - { - - Eigen::MatrixXd ivec, s, s2, V, Ut, d_dash; - string key; - obs_resid = weights * obs_resid; - int num_reals = par_resid.cols(); - double scale = (1.0 / (sqrt(double(num_reals - 1)))); - obs_diff = scale * (weights * obs_diff); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); - if (verbose_level > 1) { - Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp ); -// ss.str(""); -// ss << "solution scaling factor: " << scale; -// cout << ss.str() << endl; - } - if (use_prior_scaling) - - par_diff = scale * parcov_inv * par_diff; - else - par_diff = scale * par_diff; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); - - Ut.transposeInPlace(); - //obs_diff.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - - s2 = s.cwiseProduct(s); - - ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); - - Eigen::MatrixXd X1 = Ut * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - Eigen::MatrixXd X2 = ivec * X1; - //X1.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); - - Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; - //X2.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); - upgrade_1 = -1.0 * par_diff * X3; - - if (use_prior_scaling) { - //upgrade_1 = parcov_inv * upgrade_1; - } - - upgrade_1.transposeInPlace(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - //X3.resize(0, 0); - - Eigen::MatrixXd upgrade_2; - if ((!use_approx) && (iter > 1)) { - if (use_prior_scaling) { - par_resid = parcov_inv * par_resid; - } - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); - Eigen::MatrixXd x4 = Am.transpose() * par_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); - - par_resid.resize(0, 0); - - Eigen::MatrixXd x5 = Am * x4; - //x4.resize(0, 0); - //Am.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); - Eigen::MatrixXd x6 = par_diff.transpose() * x5; - //x5.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); - Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; - //x6.resize(0, 0); - - if (use_prior_scaling) { - upgrade_2 = -1.0 * parcov_inv * par_diff * x7; - } else { - upgrade_2 = -1.0 * (par_diff * x7); } - //x7.resize(0, 0); - - upgrade_1 = upgrade_1 + upgrade_2.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); - //upgrade_2.resize(0, 0); - + f.close(); + return; } - } -} - -void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, - vector par_names, vector obs_names) -{ - //declare these helpers in here so they are thread safe... - class local_utils - { - public: static Eigen::DiagonalMatrix get_matrix_from_map(vector& names, unordered_map& dmap) { Eigen::VectorXd vec(names.size()); @@ -2531,7 +2501,7 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b if (verbose_level > 2) { ss.str(""); - ss << "thread_" << thread_id << "part_map.csv"; + ss << "thread_" << thread_id << ".part_map.csv"; f_thread.open(ss.str()); ss.str(""); } @@ -2612,16 +2582,20 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b } } - if (verbose_level > 2) - { - f_thread << t_count << "," << iter; - for (auto name : par_names) - f_thread << "," << name; - for (auto name : obs_names) - f_thread << "," << name; - f_thread << endl; - } - +// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_obs_names",obs_names); +// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_par_names",par_names); +// +// if (verbose_level > 2) +// { +// +// f_thread << t_count << "," << iter; +// for (auto name : par_names) +// f_thread << "," << name; +// for (auto name : obs_names) +// f_thread << "," << name; +// f_thread << endl; +// } + //reset all the solution parts par_resid.resize(0, 0); par_diff.resize(0, 0); @@ -2635,106 +2609,6 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b Am.resize(0, 0); - //now loop until this thread gets access to all the containers it needs to solve with -// while (true) -// { -// //if all the solution pieces are filled, break out and solve! -// if (((use_approx) || (par_resid.rows() > 0)) && -// (weights.size() > 0) && -// (parcov_inv.size() > 0) && -// (par_diff.rows() > 0) && -// (obs_resid.rows() > 0) && -// (obs_err.rows() > 0) && -// (obs_diff.rows() > 0) && -// ((!use_localizer) || (loc.rows() > 0)) && -// ((use_approx) || (Am.rows() > 0))) -// break; -// -// //get access to the localizer -// if ((use_localizer) && (loc.rows() == 0) && (loc_guard.try_lock())) -// { -// //get a matrix that is either the shape of par diff or obs diff -// if (loc_by_obs) -// { -// //loc = localizer.get_localizing_par_hadamard_matrix(num_reals, obs_names[0], par_names); -// loc = localizer.get_pardiff_hadamard_matrix(num_reals, key, par_names); -// } -// -// else -// { -// //loc = localizer.get_localizing_obs_hadamard_matrix(num_reals, par_names[0], obs_names); -// loc = localizer.get_obsdiff_hadamard_matrix(num_reals, key, obs_names); -// } -// loc_guard.unlock(); -// } -// -// //get access to the obs_diff container -// if ((obs_diff.rows() == 0) && (obs_diff_guard.try_lock())) -// { -// //piggy back here for thread safety -// //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_svd_pack() == PestppOptions::SVD_PACK::PROPACK) -// // use_propack = true; -// obs_diff = local_utils::get_matrix_from_map(num_reals, obs_names, obs_diff_map); -// obs_diff_guard.unlock(); -// } -// -// //get access to the residual container -// if ((obs_resid.rows() == 0) && (obs_resid_guard.try_lock())) -// { -// obs_resid = local_utils::get_matrix_from_map(num_reals, obs_names, obs_resid_map); -// obs_resid_guard.unlock(); -// } -// -// //get access to the obs noise container -// if ((obs_err.rows() == 0) && (obs_err_guard.try_lock())) -// { -// obs_err = local_utils::get_matrix_from_map(num_reals, obs_names, obs_err_map); -// obs_err_guard.unlock(); -// } -// -// //get access to the par diff container -// if ((par_diff.rows() == 0) && (par_diff_guard.try_lock())) -// { -// par_diff = local_utils::get_matrix_from_map(num_reals, par_names, par_diff_map); -// par_diff_guard.unlock(); -// } -// -// //get access to the par residual container -// if ((par_resid.rows() == 0) && (par_resid_guard.try_lock())) -// { -// par_resid = local_utils::get_matrix_from_map(num_reals, par_names, par_resid_map); -// par_resid_guard.unlock(); -// } -// -// //get access to the obs weights container -// if ((weights.rows() == 0) && (weight_guard.try_lock())) -// { -// weights = local_utils::get_matrix_from_map(obs_names, weight_map); -// weight_guard.unlock(); -// } -// -// //get access to the inverse prior parcov -// if ((parcov_inv.rows() == 0) && (parcov_guard.try_lock())) -// { -// parcov_inv = local_utils::get_matrix_from_map(par_names, parcov_inv_map); -// parcov_guard.unlock(); -// } -// -// //if needed, get access to the Am container - needed in the full glm solution -// if ((!use_approx) && (Am.rows() == 0) && (am_guard.try_lock())) -// { -// //Am = local_utils::get_matrix_from_map(num_reals, par_names, Am_map).transpose(); -// int am_cols = Am_map[par_names[0]].size(); -// Am.resize(par_names.size(), am_cols); -// Am.setZero(); -// -// for (int j = 0; j < par_names.size(); j++) -// { -// Am.row(j) = Am_map[par_names[j]]; -// } -// am_guard.unlock(); -// } -// } if ((use_localizer) && (loc.rows() == 0)) { while (true) { @@ -2776,216 +2650,242 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b } } - - - par_diff.transposeInPlace(); + Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, d_dash; + par_diff.transposeInPlace(); obs_diff.transposeInPlace(); obs_resid.transposeInPlace(); par_resid.transposeInPlace(); obs_err.transposeInPlace(); - //form the scaled obs resid matrix - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); - //Eigen::MatrixXd scaled_residual = weights * obs_resid; - - //form the (optionally) scaled par resid matrix - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); - - stringstream ss; - - double scale = (1.0 / (sqrt(double(num_reals - 1)))); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); - - //apply the localizer here... - if (use_localizer) - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "loc", loc); - if (use_localizer) - { - if (loc_by_obs) - par_diff = par_diff.cwiseProduct(loc); - else - obs_diff = obs_diff.cwiseProduct(loc); - } - - //Eigen::MatrixXd upgrade_1; -// ensemble_solution(iter,verbose_level,maxsing,thread_id,t_count, use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid,par_diff,Am,obs_resid, -// obs_diff,upgrade_1,obs_err,weights,parcov_inv); - - Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, d_dash; - - //---------------------------------- - //es-mda solution - //---------------------------------- - - if (!use_glm_form) - { - if (true) - { - // Low rank Cee. Section 14.3.2 Evenson Book - obs_err = obs_err.colwise() - obs_err.rowwise().mean(); - obs_err = sqrt(cur_lam) * obs_err; - Eigen::MatrixXd s0, V0, U0, s0_i; - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); - s0_i = s0.asDiagonal().inverse(); - Eigen::MatrixXd X0 = U0.transpose() * obs_err; - X0 = s0_i * X0; - Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; - rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); - - s1_2 = s1.cwiseProduct(s1); - s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); - Eigen::MatrixXd X1 = s0_i * U1; - X1 = U0 * X1; - - Eigen::MatrixXd X4 = s1_2i * X1.transpose(); - Eigen::MatrixXd X2 = X4 * obs_resid; - Eigen::MatrixXd X3 = X1 * X2; - - X3 = obs_diff.transpose() * X3; - upgrade_1 = -1 * par_diff * X3; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - upgrade_1.transposeInPlace(); - - } - else - { - // Use when ensemble size is larger than number of observation. - // This is a good option when number of observation is small - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); - obs_err = obs_err.colwise() - obs_err.rowwise().mean(); - - Eigen::MatrixXd s2_, s, V, U, cum_sum; - SVD_REDSVD rsvd; - Eigen::MatrixXd C; - C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); - Eigen::VectorXd s2; - - - rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - s2 = s.cwiseProduct(s); - s2_ = s2.asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); - - - Eigen::MatrixXd X1 = s2_ * U.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - X1 = X1 * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); - - X1 = U * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); - - X1 = obs_diff.transpose() * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); - - upgrade_1 = -1 * par_diff * X1; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - upgrade_1.transposeInPlace(); - } - - - } - - -// ---------------------------------- -// glm solution -// ---------------------------------- - else - { - - obs_resid = weights * obs_resid; - obs_diff = scale * (weights * obs_diff); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); - if (use_prior_scaling) - par_diff = scale * parcov_inv * par_diff; - else - par_diff = scale * par_diff; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); - SVD_REDSVD rsvd; - rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); - - Ut.transposeInPlace(); - obs_diff.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); - - Eigen::MatrixXd s2 = s.cwiseProduct(s); - - ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); - - Eigen::MatrixXd X1 = Ut * obs_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); - - Eigen::MatrixXd X2 = ivec * X1; - X1.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); - - Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; - X2.resize(0, 0); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); - upgrade_1 = -1.0 * par_diff * X3; + UpgradeThread::ensemble_solution(iter, verbose_level,maxsing, thread_id, + t_count, use_prior_scaling,use_approx, use_glm_form, + cur_lam,eigthresh, par_resid, par_diff, + Am, obs_resid,obs_diff, upgrade_1, + obs_err, weights, + parcov_inv, + obs_names,par_names); - if (use_prior_scaling) - { - //upgrade_1 = parcov_inv * upgrade_1; - } - - upgrade_1.transposeInPlace(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); - X3.resize(0, 0); - - Eigen::MatrixXd upgrade_2; - if ((!use_approx) && (iter > 1)) - { - if (use_prior_scaling) - { - par_resid = parcov_inv * par_resid; - } - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); - Eigen::MatrixXd x4 = Am.transpose() * par_resid; - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); - - par_resid.resize(0, 0); - - Eigen::MatrixXd x5 = Am * x4; - x4.resize(0, 0); - Am.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); - Eigen::MatrixXd x6 = par_diff.transpose() * x5; - x5.resize(0, 0); - - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); - Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; - x6.resize(0, 0); - - if (use_prior_scaling) - { - upgrade_2 = -1.0 * parcov_inv * par_diff * x7; - } - else - { - upgrade_2 = -1.0 * (par_diff * x7); - } - x7.resize(0, 0); - upgrade_1 = upgrade_1 + upgrade_2.transpose(); - local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); - upgrade_2.resize(0, 0); - - } - } +// par_diff.transposeInPlace(); +// obs_diff.transposeInPlace(); +// obs_resid.transposeInPlace(); +// par_resid.transposeInPlace(); +// obs_err.transposeInPlace(); +// +// //form the scaled obs resid matrix +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); +// //Eigen::MatrixXd scaled_residual = weights * obs_resid; +// +// //form the (optionally) scaled par resid matrix +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); +// +// stringstream ss; +// +// double scale = (1.0 / (sqrt(double(num_reals - 1)))); +// +// if (verbose_level > 1) { +// ss.str(""); +// ss << "solution scaling factor: " << scale; +// cout << ss.str() << endl; +// if (parcov_inv.size() < 10000) +// { +// Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); +// } +// } +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); +// //apply the localizer here... +// if (use_localizer) +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "loc", loc); +// if (use_localizer) +// { +// if (loc_by_obs) +// par_diff = par_diff.cwiseProduct(loc); +// else +// obs_diff = obs_diff.cwiseProduct(loc); +// } +// +// //Eigen::MatrixXd upgrade_1; +//// ensemble_solution(iter,verbose_level,maxsing,thread_id,t_count, use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid,par_diff,Am,obs_resid, +//// obs_diff,upgrade_1,obs_err,weights,parcov_inv); +// +// +// +// //---------------------------------- +// //es-mda solution +// //---------------------------------- +// +// if (!use_glm_form) +// { +// if (true) +// { +// // Low rank Cee. Section 14.3.2 Evenson Book +// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); +// obs_err = sqrt(cur_lam) * obs_err; +// Eigen::MatrixXd s0, V0, U0, s0_i; +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); +// s0_i = s0.asDiagonal().inverse(); +// Eigen::MatrixXd X0 = U0.transpose() * obs_err; +// X0 = s0_i * X0; +// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; +// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); +// +// s1_2 = s1.cwiseProduct(s1); +// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); +// Eigen::MatrixXd X1 = s0_i * U1; +// X1 = U0 * X1; +// +// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); +// Eigen::MatrixXd X2 = X4 * obs_resid; +// Eigen::MatrixXd X3 = X1 * X2; +// +// X3 = obs_diff.transpose() * X3; +// upgrade_1 = -1 * par_diff * X3; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// upgrade_1.transposeInPlace(); +// +// } +// else +// { +// // Use when ensemble size is larger than number of observation. +// // This is a good option when number of observation is small +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); +// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); +// +// Eigen::MatrixXd s2_, s, V, U, cum_sum; +// SVD_REDSVD rsvd; +// Eigen::MatrixXd C; +// C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); +// Eigen::VectorXd s2; +// +// +// rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); +// s2 = s.cwiseProduct(s); +// s2_ = s2.asDiagonal().inverse(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); +// +// +// Eigen::MatrixXd X1 = s2_ * U.transpose(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); +// +// X1 = X1 * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); +// +// X1 = U * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); +// +// X1 = obs_diff.transpose() * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); +// +// upgrade_1 = -1 * par_diff * X1; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// upgrade_1.transposeInPlace(); +// } +// +// +// } +// +// +//// ---------------------------------- +//// glm solution +//// ---------------------------------- +// else +// { +// +// obs_resid = weights * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_resid", par_diff); +// obs_diff = scale * (weights * obs_diff); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_diff", par_diff); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); +// if (use_prior_scaling) +// par_diff = scale * parcov_inv * par_diff; +// else +// par_diff = scale * par_diff; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); +// SVD_REDSVD rsvd; +// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); +// +// Ut.transposeInPlace(); +// obs_diff.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); +// +// Eigen::MatrixXd s2 = s.cwiseProduct(s); +// +// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); +// +// Eigen::MatrixXd X1 = Ut * obs_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); +// +// Eigen::MatrixXd X2 = ivec * X1; +// X1.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); +// +// Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; +// X2.resize(0, 0); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); +// upgrade_1 = -1.0 * par_diff * X3; +// +// if (use_prior_scaling) +// { +// //upgrade_1 = parcov_inv * upgrade_1; +// } +// +// upgrade_1.transposeInPlace(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); +// X3.resize(0, 0); +// +// Eigen::MatrixXd upgrade_2; +// if ((!use_approx) && (iter > 1)) +// { +// if (use_prior_scaling) +// { +// par_resid = parcov_inv * par_resid; +// } +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); +// Eigen::MatrixXd x4 = Am.transpose() * par_resid; +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); +// +// par_resid.resize(0, 0); +// +// Eigen::MatrixXd x5 = Am * x4; +// x4.resize(0, 0); +// Am.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); +// Eigen::MatrixXd x6 = par_diff.transpose() * x5; +// x5.resize(0, 0); +// +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); +// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; +// x6.resize(0, 0); +// +// if (use_prior_scaling) +// { +// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; +// } +// else +// { +// upgrade_2 = -1.0 * (par_diff * x7); +// } +// x7.resize(0, 0); +// +// upgrade_1 = upgrade_1 + upgrade_2.transpose(); +// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); +// upgrade_2.resize(0, 0); +// +// } +// } while (true) { @@ -4245,6 +4145,63 @@ map L2PhiHandler::calc_regul(ParameterEnsemble & pe) return phi_map; } +vector L2PhiHandler::get_violating_realizations(ObservationEnsemble& oe, const vector& viol_obs_names) +{ + stringstream ss; + vector viol_real_names; + if (viol_obs_names.size() == 0) + return viol_real_names; + + Eigen::MatrixXd resid = get_actual_obs_resid(oe); + map vmap; + vector nz_onames = oe_base->get_var_names(); + for (int i=0;i::iterator end = vmap.end(); + vector missing; + ObservationInfo* oi = pest_scenario->get_observation_info_ptr(); + vector nz_viol_obs_names; + map nz_viol_vmap; + for (auto& name : viol_obs_names) + { + if (vmap.find(name) == end) + missing.push_back(name); + else + { + if (oi->get_weight(name) > 0) + { + nz_viol_obs_names.push_back(name); + nz_viol_vmap[name] = vmap.at(name); + } + } + } +// if (missing.size() > 0) +// { +// ss.str(""); +// ss << "L2PhiHandler::get_violating_realizations(): the following violation obs names not found " << endl; +// for (auto& m : missing) +// ss << m << endl; +// throw runtime_error(ss.str()); +// } + Eigen::VectorXd real; + double sum; + vector rnames = oe.get_real_names(); + for (int i=0;i 1.0e-7) + { + viol_real_names.push_back(rnames[i]); + } + } + return viol_real_names; +} + void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector &names) { @@ -5606,6 +5563,8 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) message(1, "lambda decrease factor: ", dec_fac); message(1, "max run fail: ", ppo->get_max_run_fail()); + prep_drop_violations(); + //todo: add sanity checks for weights en //like conflict with da weight cycle table sanity_checks(); @@ -6134,28 +6093,52 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) //check that all obs are in conflict message(1, "dropping conflicted observations"); - if (in_conflict.size() == act_obs_names.size()) - { - if (cycle == NetPackage::NULL_DA_CYCLE) { + if (in_conflict.size() == act_obs_names.size()) { + if (cycle == NetPackage::NULL_DA_CYCLE) { throw_em_error("all non-zero weighted observations in conflict state, cannot continue"); + } else { + message(0, "all non-zero weighted observations in conflict state, continuing to next cycle"); + zero_weight_obs(in_conflict, false, false); + ph.update(oe, pe, weights); + return; + } + } + if (violation_obs.size() > 0) + { + set sviol_obs(violation_obs.begin(),violation_obs.end()); + set::iterator end = sviol_obs.end(); + vector temp,temp2; + for (auto& obs : in_conflict) + { + if (sviol_obs.find(obs) == end) + { + temp.push_back(obs); + } + else + { + temp2.push_back(obs); + } } - else + if (temp2.size() > 0) { - message(0,"all non-zero weighted observations in conflict state, continuing to next cycle"); - zero_weight_obs(in_conflict,false,false); - ph.update(oe,pe, weights); - return; + ss.str(""); + ss << "WARNING: the following observations are 'in conflict' but" << endl; + ss << " also contain 'drop_violation' conditions, " << endl; + ss << " so they are not being treated as conflicts...user beware..." << endl; + message(1,ss.str(),temp2); + in_conflict = temp; + } + } + if (in_conflict.size() > 0) { + zero_weight_obs(in_conflict); + if (ppo->get_ies_localizer().size() > 0) { + message(1, "updating localizer"); + if (localizer.get_use()) + localizer.get_orgmat_ptr()->clear_names(); + ofstream &frec = file_manager.rec_ofstream(); + use_localizer = localizer.initialize(performance_log, frec, true); } } - zero_weight_obs(in_conflict); - if (ppo->get_ies_localizer().size() > 0) - { - message(1, "updating localizer"); - if (localizer.get_use()) - localizer.get_orgmat_ptr()->clear_names(); - ofstream& frec = file_manager.rec_ofstream(); - use_localizer = localizer.initialize(performance_log, frec, true); - } } } @@ -6173,7 +6156,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) } - drop_bad_phi(pe, oe); + drop_bad_reals(pe, oe); if (oe.shape().first == 0) { throw_em_error(string("all realizations dropped as 'bad'")); @@ -6316,6 +6299,43 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) message(0, "initialization complete"); } +void EnsembleMethod::prep_drop_violations() +{ + violation_obs.clear(); + string section = "observation data external"; + string viol_col = "drop_violations"; + map viol_map = pest_scenario.get_ext_file_string_map(section,viol_col); + if (viol_map.size() == 0) + { + return; + } + + string sval; + for (auto& v : viol_map) + { + sval = pest_utils::strip_cp(pest_utils::lower_cp(v.second)); + if (sval == "true") + { + violation_obs.push_back(v.first); + } + } + stringstream ss; + ss << violation_obs.size() << " 'drop_violations' observations detected, see rec file for listing"; + message(1,ss.str()); + ofstream& f_rec = file_manager.rec_ofstream(); + f_rec << endl << "The following observations are being used as 'drop_violations' constraints:" << endl; + + for (auto& o : violation_obs) + { + f_rec << o << endl; + } + f_rec << "Note: any 'drop_violation' observations that are not inequalities will be treated as " << endl; + f_rec << " equality constraints such that any realization that does not match the observation" << endl; + f_rec << " value (numerically) exactly will be dropped. User beware..." << endl << endl; + +} + + void EnsembleMethod::check_and_fill_phi_factors(map>& group_to_obs_map,map>& group_map,map>& phi_fracs_by_real, vector& index, bool check_reals) { @@ -7499,7 +7519,7 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto frec << "lambda, scale value " << lam_vals[i] << ',' << scale_vals[i] << " obs ensemble saved to " << ss.str() << endl; } - drop_bad_phi(pe_lams[i], oe_lams[i], subset_idxs); + drop_bad_reals(pe_lams[i], oe_lams[i], subset_idxs); if (oe_lams[i].shape().first == 0) { message(1, "all realizations dropped as 'bad' for lambda, scale fac ", vals); @@ -7736,7 +7756,7 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto //append the remaining obs en oe_lam_best.append_other_rows(remaining_oe_lam); assert(pe_lams[best_idx].shape().first == oe_lam_best.shape().first); - drop_bad_phi(pe_lams[best_idx], oe_lam_best); + drop_bad_reals(pe_lams[best_idx], oe_lam_best); if (oe_lam_best.shape().first == 0) { throw_em_error(string("all realization dropped after finishing subset runs...something might be wrong...")); @@ -9171,19 +9191,40 @@ Eigen::MatrixXd EnsembleMethod::get_Am(const vector& real_names, const v return Am; } -void EnsembleMethod::drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector subset_idxs) +void EnsembleMethod::drop_bad_reals(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector subset_idxs) { + stringstream ss; //don't use this assert because _pe maybe full size, but _oe might be subset size bool is_subset = false; if (subset_idxs.size() > 0) is_subset = true; if (!is_subset) if (_pe.shape().first != _oe.shape().first) - throw_em_error("EnsembleMethod::drop_bad_phi() error: _pe != _oe and not subset"); + throw_em_error("EnsembleMethod::drop_bad_reals() error: _pe != _oe and not subset"); double bad_phi = pest_scenario.get_pestpp_options().get_ies_bad_phi(); double bad_phi_sigma = pest_scenario.get_pestpp_options().get_ies_bad_phi_sigma(); vector idxs = ph.get_idxs_greater_than(bad_phi, bad_phi_sigma, _oe, weights); + vector viol_reals = ph.get_violating_realizations(_oe,violation_obs); + map rmap = _oe.get_real_map(); + int idx; + if (viol_reals.size() > 0) { + ss.str(""); + ss << viol_reals.size() << " realizations meet 'drop_violations' conditions"; + message(2,ss.str()); + for (auto &v : viol_reals) { + if (v == BASE_REAL_NAME) { + message(3, "not dropping 'base' real even though it meets 'drop_violations' conditions"); + + } else + { + idx = rmap.at(v); + if (find(idxs.begin(), idxs.end(), idx) == idxs.end()) { + idxs.push_back(idx); + } + } + } + } if (pest_scenario.get_pestpp_options().get_ies_debug_bad_phi()) idxs.push_back(0); @@ -9194,7 +9235,7 @@ void EnsembleMethod::drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _ message(0, "dropping realizations as bad: ", idxs.size()); vector par_real_names = _pe.get_real_names(), obs_real_names = _oe.get_real_names(); - stringstream ss; + string pname; string oname; @@ -9222,7 +9263,7 @@ void EnsembleMethod::drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _ if (find(subset_idxs.begin(), subset_idxs.end(), pidx) == subset_idxs.end()) { ss.str(""); - ss << "drop_bad_phi() error: idx " << pidx << " not found in subset_idxs"; + ss << "drop_bad_reals() error: idx " << pidx << " not found in subset_idxs"; throw_em_error(ss.str()); } pname = full_pnames[pidx]; @@ -9247,12 +9288,12 @@ void EnsembleMethod::drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _ catch (const exception& e) { stringstream ss; - ss << "drop_bad_phi() error : " << e.what(); + ss << "drop_bad_reals() error : " << e.what(); throw_em_error(ss.str()); } catch (...) { - throw_em_error(string("drop_bad_phi() error")); + throw_em_error(string("drop_bad_reals() error")); } } } diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.h b/src/libs/pestpp_common/EnsembleMethodUtils.h index c9c9cfd7..b9c264f7 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.h +++ b/src/libs/pestpp_common/EnsembleMethodUtils.h @@ -111,6 +111,8 @@ class L2PhiHandler map get_actual_swr_map(ObservationEnsemble& oe, string real_name=""); map> get_meas_phi_weight_ensemble(ObservationEnsemble& oe, ObservationEnsemble& weights); + vector get_violating_realizations(ObservationEnsemble& oe, const vector& viol_obs_names); + private: string tag; map get_summary_stats(phiType pt); @@ -287,6 +289,14 @@ class UpgradeThread virtual void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, vector obs_names) { ; } + + static void ensemble_solution(const int iter, const int verbose_level,const int maxsing, const int thread_id, + const int t_count, const bool use_prior_scaling,const bool use_approx, const bool use_glm, + const double cur_lam,const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff, + const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1, + Eigen::MatrixXd& obs_err, const Eigen::DiagonalMatrix& weights, + const Eigen::DiagonalMatrix& parcov_inv, + const vector& act_obs_names,const vector& act_par_names); protected: PerformanceLog* performance_log; Localizer::How how; @@ -308,55 +318,35 @@ class UpgradeThread mutex par_diff_lock, am_lock, put_lock, obs_err_lock; mutex next_lock; -}; - -void ensemble_solution(const int iter, const int verbose_level,const int maxsing, const int thread_id, - const int t_count, const bool - use_prior_scaling,const bool use_approx, const bool use_glm, const double cur_lam, - const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff, - const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1, - Eigen::MatrixXd& obs_err, - const Eigen::DiagonalMatrix& weights, - const Eigen::DiagonalMatrix& parcov_inv); - -class CovLocalizationUpgradeThread : public UpgradeThread -{ -public: - using UpgradeThread::UpgradeThread; - void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, vector obs_names); }; +//void ensemble_solution(const int iter, const int verbose_level,const int maxsing, const int thread_id, +// const int t_count, const bool +// use_prior_scaling,const bool use_approx, const bool use_glm, const double cur_lam, +// const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff, +// const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1, +// Eigen::MatrixXd& obs_err, +// const Eigen::DiagonalMatrix& weights, +// const Eigen::DiagonalMatrix& parcov_inv, +// const vector& act_obs_names,const vector& act_par_names); + + +//class CovLocalizationUpgradeThread : public UpgradeThread +//{ +//public: +// using UpgradeThread::UpgradeThread; +// +// void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, vector obs_names); +//}; + class LocalAnalysisUpgradeThread: public UpgradeThread { public: using UpgradeThread::UpgradeThread; void work(int thread_id, int iter, double cur_lam,bool use_glm_form, vector par_names, vector obs_names); - - -//private: -// PerformanceLog* performance_log; -// Localizer::How how; -// vector keys; -// int count, total; -// -// unordered_map, vector>>& cases; -// -// ParameterEnsemble& pe_upgrade; -// Localizer& localizer; -// unordered_map& parcov_inv_map; -// unordered_map& weight_map; -// -// unordered_map& par_resid_map, & par_diff_map, & Am_map; -// unordered_map& obs_resid_map, & obs_diff_map, obs_err_map; -// -// mutex ctrl_lock, weight_lock, loc_lock, parcov_lock; -// mutex obs_resid_lock, obs_diff_lock, par_resid_lock; -// mutex par_diff_lock, am_lock, put_lock, obs_err_lock; -// mutex next_lock; - }; class EnsembleMethod @@ -456,11 +446,13 @@ class EnsembleMethod vector lam_mults; vector oe_org_real_names, pe_org_real_names; vector act_obs_names, act_par_names; + vector violation_obs; ParameterEnsemble pe, pe_base; ObservationEnsemble oe, oe_base, weights; Eigen::DiagonalMatrix obscov_inv_sqrt, parcov_inv_sqrt; bool oe_drawn, pe_drawn; + bool solve_glm(int cycle = NetPackage::NULL_DA_CYCLE); bool solve_mda(bool last_iter, int cycle = NetPackage::NULL_DA_CYCLE); @@ -477,7 +469,7 @@ class EnsembleMethod void initialize_restart(); - void drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector subset_idxs = vector()); + void drop_bad_reals(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector subset_idxs = vector()); void add_bases(); @@ -505,5 +497,7 @@ class EnsembleMethod map>& phi_fracs_by_real, vector& index, bool check_reals); + void prep_drop_violations(); + }; #endif diff --git a/src/libs/pestpp_common/Localizer.cpp b/src/libs/pestpp_common/Localizer.cpp index 72ebd800..d7b3e6a1 100644 --- a/src/libs/pestpp_common/Localizer.cpp +++ b/src/libs/pestpp_common/Localizer.cpp @@ -29,7 +29,12 @@ bool Localizer::initialize(PerformanceLog *performance_log, ofstream& frec, bool loc_typ = pest_scenario_ptr->get_pestpp_options().get_ies_loc_type(); if (loc_typ[0] == 'C') { - loctyp = LocTyp::COVARIANCE; + ss.str(""); + ss << "WARNING: covariance type localization deprecated, resetting to local-analysis"; + frec << ss.str() << endl; + cout << ss.str() << endl; + performance_log->log_event(ss.str()); + loctyp = LocTyp::LOCALANALYSIS; } else if (loc_typ[0] == 'L') { @@ -42,14 +47,28 @@ bool Localizer::initialize(PerformanceLog *performance_log, ofstream& frec, bool throw runtime_error("Localizer error: 'ies_localization_type' must start with 'C' (covariance) or 'L' (local analysis) not " + loc_typ); } - if (how_str[0] == 'P') - { - how = How::PARAMETERS; - } - else if (how_str[0] == 'O') - { - how = How::OBSERVATIONS; - } + how = How::PARAMETERS; + if (how_str[0] == 'P') + { + + } + else if (how_str[0] == 'O') + { + ss.str(""); + ss << "WARNING: localization 'how' by observations is deprecated, resetting to 'how' by parameters"; + frec << ss.str() << endl; + cout << ss.str() << endl; + performance_log->log_event(ss.str()); + how = How::PARAMETERS; + } +// if (how_str[0] == 'P') +// { +// how = How::PARAMETERS; +// } +// else if (how_str[0] == 'O') +// { +// how = How::OBSERVATIONS; +// } else { performance_log->log_event("Localizer error: 'ies_localize_how' or 'da_localize_how' must start with 'P' (pars) or 'O' (obs) not " + how_str); @@ -57,13 +76,13 @@ bool Localizer::initialize(PerformanceLog *performance_log, ofstream& frec, bool throw runtime_error("Localizer error: 'ies_localize_how' or 'da_localize_how' must start with 'P' (pars) or 'O' (obs) not " + how_str); } - if ((loctyp == LocTyp::COVARIANCE) && (how == How::OBSERVATIONS)) { - performance_log->log_event("Localizer error: covariance localization can only be used with localization by parameters"); - cout << "Localizer error: covariance localization can only be used with localization by parameters" << endl; - - throw runtime_error( - "Localizer error: covariance localization can only be used with localization by parameters"); - } +// if ((loctyp == LocTyp::COVARIANCE) && (how == How::OBSERVATIONS)) { +// performance_log->log_event("Localizer error: covariance localization can only be used with localization by parameters"); +// cout << "Localizer error: covariance localization can only be used with localization by parameters" << endl; +// +// throw runtime_error( +// "Localizer error: covariance localization can only be used with localization by parameters"); +// } string filename = pest_scenario_ptr->get_pestpp_options().get_ies_localizer(); autoadaloc = pest_scenario_ptr->get_pestpp_options().get_ies_autoadaloc(); sigma_dist = pest_scenario_ptr->get_pestpp_options().get_ies_autoadaloc_sigma_dist(); diff --git a/src/libs/pestpp_common/Pest.h b/src/libs/pestpp_common/Pest.h index c3131e6e..dfe23ad7 100644 --- a/src/libs/pestpp_common/Pest.h +++ b/src/libs/pestpp_common/Pest.h @@ -121,7 +121,7 @@ class Pest { protected: //this is the list of external file cols that have meaning... - set efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link" }; + set efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link","drop_violations" }; int n_adj_par = 0; string prior_info_string; ControlInfo control_info; diff --git a/src/libs/pestpp_common/SVDPackage.h b/src/libs/pestpp_common/SVDPackage.h index c490d0a0..899eed72 100644 --- a/src/libs/pestpp_common/SVDPackage.h +++ b/src/libs/pestpp_common/SVDPackage.h @@ -27,7 +27,7 @@ class SVDPackage { public: - SVDPackage(std::string _descritpion="undefined", int _n_max_sing=1000, double _eign_thres=1.0e-7); + SVDPackage(std::string _descritpion="undefined", int _n_max_sing=1000000, double _eign_thres=1.0e-7); SVDPackage(const SVDPackage &rhs) : description(rhs.description), n_max_sing(rhs.n_max_sing), eign_thres(rhs.eign_thres), performance_log(rhs.performance_log) {} virtual void solve_ip(Eigen::SparseMatrix& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix & U, Eigen::SparseMatrix& VT, Eigen::VectorXd &Sigma_trunc) = 0; virtual void solve_ip(Eigen::SparseMatrix& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix & U, Eigen::SparseMatrix& VT, Eigen::VectorXd &Sigma_trunc, double _eigen_thres) = 0; @@ -48,7 +48,7 @@ class SVDPackage class SVD_EIGEN : public SVDPackage { public: - SVD_EIGEN(int _n_max_sing = 1000, double _eign_thres = 1.0e-7) : SVDPackage("Eigen JacobiSVD", _n_max_sing, _eign_thres) {} + SVD_EIGEN(int _n_max_sing = 1000000, double _eign_thres = 1.0e-7) : SVDPackage("Eigen JacobiSVD", _n_max_sing, _eign_thres) {} SVD_EIGEN(const SVD_EIGEN &rhs) : SVDPackage(rhs) {} virtual void solve_ip(Eigen::SparseMatrix& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix& U, Eigen::SparseMatrix& VT, Eigen::VectorXd &Sigma_trunc); virtual void solve_ip(Eigen::SparseMatrix& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix& U, Eigen::SparseMatrix& VT, Eigen::VectorXd &Sigma_trunc, double _eigen_thres); @@ -62,7 +62,7 @@ class SVD_EIGEN : public SVDPackage class SVD_REDSVD : public SVDPackage { public: - SVD_REDSVD(int _n_max_sing = 1000, double _eign_thres = 1.0e-7) : SVDPackage("RedSVD", _n_max_sing, _eign_thres) {} + SVD_REDSVD(int _n_max_sing = 1000000, double _eign_thres = 1.0e-7) : SVDPackage("RedSVD", _n_max_sing, _eign_thres) {} SVD_REDSVD(const SVD_REDSVD &rhs) : SVDPackage(rhs) {} virtual void solve_ip(Eigen::SparseMatrix& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix& U, Eigen::SparseMatrix& VT, Eigen::VectorXd &Sigma_trunc); diff --git a/src/programs/pestpp-ies/pestpp-ies.cpp b/src/programs/pestpp-ies/pestpp-ies.cpp index ee554999..41257816 100644 --- a/src/programs/pestpp-ies/pestpp-ies.cpp +++ b/src/programs/pestpp-ies/pestpp-ies.cpp @@ -133,7 +133,7 @@ int main(int argc, char* argv[]) if (!restart_flag || save_restart_rec_header) { fout_rec << " pestpp-ies - a GLM iterative Ensemble Smoother" << endl - << "for PEST(++) datasets " << endl << endl; + << " for PEST(++) datasets " << endl << endl; fout_rec << " by the PEST++ developement team" << endl << endl << endl; fout_rec << endl; fout_rec << endl << endl << "version: " << version << endl;