Running Pest++ with WRF-HYDRO

[jbensabat]

Hello I am running pest++ to calibrate a WRF-HYDRO mode. Prior to running the calibration phase I ran the sensitivity analysis (Morris) to identify the most sensitive parameters. Looking at the pest++ iterative process I see that it almost does change the parameter values (really insignificant changes) while there are gaps between measured and simulated values. Has anybody encountered such kind of behavior ? is it due to specific settings (line search and or partial derivative calculation?)
any help would really be appreciated.
best
jac

[jtwhite79]

Are you saying that pestpp-glm isnt able to reduce the objective function across several iterations? Is anything reported in the rec file regarding upgrade length and if upgrade vectors are being truncated to enforce parameter bounds and/or parameter change limits? What version of pestpp-glm are you using? @mmorphew recently made some changes to pestpp-glm in the way it handles parameter bound enforcement and that might help in this situation...

[jbensabat]

hi and thanks for your reply. Yes the changes in the parameters is almost insignificant with regard to the parameter values. I have attached the rec file. I am running version 5.2.18. I have checked the parameters are not close to their bounds (at least the vast majority).
best
jac
HMU17_GLM_03_01-08_01_2020.zip

[mmorphew]

I've looked through the .rec file. I don't think universal-shrinkage is the issue here. I do think it's related to the other issue I fixed regarding lambda search adjustment. The first iteration has significant phi reduction, but then there's no phi improvement onward. This tells me that it may be lambda-related.

@jbensabat , I would try building the development branch version of pestpp-glm if you can and running it instead. The development branch improves how lambda values are tested. You may also want to implement a wider lambda search in your pest control file using the ++ lambdas() option as well as the ++ lambda_scale_vec() option.

[jtwhite79]

Actually I think your changes are rolled into 5.2.19 if you wanna just download precompiled bins

[mmorphew]

Oh, even better! Thank you for letting me know.

[jbensabat]

Hi thanks
Just for my curiosity the lambda parameter is picked from a set of predefined values or is it determined via a process of secondary iterations
best
jac

[mmorphew]

If ++lambdas() is not supplied in the control file, it will use a predefined set of values (I believe it's 0.1,1,10,100, or something similar). If the maximum or minimum of those values provides the best phi reduction, the 5.2.19 version of pestpp-glm will add a new higher/lower value to the next iteration (by a power of 10) so that the algorithm doesn't get stuck.

I've gotten sidetracked on some other projects, but I plan on expanding lambda functionality to include pest_hp behavior where lambdas are selected and scaled based on the phi value of the initial forward. But the current behavior should be enough to allow pestpp-glm to find the correct values of lambda to give you phi reduction.

[jtwhite79]

This is what pestpp-ies does. @mmorphew maybe we could sync up how lambda is handled between the glm and ies?

[mmorphew]

Absolutely, would be great to have it be identical between the two. I will try to claw out some free time to look through the source code again - it's been a couple of months so I'm rusty.

[jbensabat]

hi
thanks
could you please provide an example on how to set ++lambdas and ++ lambda_scale_vec in the control file ?
best
jac

[mmorphew]

PEST++ parameters can be included on any line in the control file. You can paste and edit the following lines anywhere in the control file. I typically put my pest++ control options before the parameter group section of the control file.

++ lambdas(94023,100000,202566,1000000,10000000)
++ lambda_scale_fac(1.25,1.5,1.75,2)

The values above are fairly arbitrary and linked to some specific testing I was doing. Depending on how much compute you have, I would just do a wide range across several orders of magnitude to try to capture where the algorithm is happy with lambda. You can also include lambda_scale_vec values between 0 and 1, I was just testing higher values in that particular run.

[jbensabat]

Hi and thanks a lot for your help. If I understand correctly there are three sets of lambda related parameters
lambdas -> the actual value of the parameter in the LM iteration
lambda_scale_fac
lambda_scale_vec
could you please provide an explanation on how these parameters are used ? what is the difference between  lambda_scale_fac and lambda_scale_vec ? Are they used to set the value of lambda for the next iteration or for the current iteration ?
thanks
jac

[mmorphew]

Hello, sincere apologies for the confusion. Lambda_scale_vec was a typo on my part and does not exist as a parameter. The parameter is lambda_scale_fac, short for lambda scale factor. Lambda values are multiplied by this factor, so the total number of lambdas tested in an iteration is effectively the length of lambdas times the length of lambda_scale_fac (a scale of 1 is always included regardless of user input).

Conceptually, I think of these scaling factors as a way to catch a lambda value that is "good," but simply at the wrong scale for our current position in the objective function space. These scaling factors effectively increase/decrease the magnitude of the upgrade vector, which may lead to a better phi reduction than the original lambda. Practically, these scaling factors can be an easy way to generate more lambda values to test. If we set our lambdas by powers of 10, then common lambda scale factors like [0.5,0.75,1.25,1.5], etc. allow us to catch some lambda values in between those powers of 10 rather easily.

Re: current or next iteration? Lambda_scale_fac will not change between iterations. If we reach a point where the best lambda in the lambdas array is the highest or lowest value, then the lambdas array is expanded in the relevant direction for the next iteration.

All pestpp options can be found in the user manual:
https://github.com/usgs/pestpp/blob/master/documentation/pestpp_users_manual.md#643-pest-control-variables

[jbensabat]

thanks a lot
now it is clear
thanks
jac

[jbensabat]

Hi
I have built version 5.2.19 and updated the pst file to include the lambda definitions
I get the following runtime error:
processing control file HMU17_GLM_03_01-08_01_2020.pst
control file parsing error: '* control data keyword' cant be used with '++' args
control file parsing error: '* control data keyword' cant be used with '++' args
Error processing control file: HMU17_GLM_03_01-08_01_2020.pst

control file parsing error: '* control data keyword' cant be used with '++' args

Error condition prevents further execution:
control file parsing error: '* control data keyword' cant be used with '++' args

I am attaching the pst file
how should i resolve this incompatibility ?
best
jac

HMU17_GLM_03_01-08_01_2020.zip

[jbensabat]

I resolved it just by replacing ++ by *
thanks
jac

[mmorphew]

Hi Jac,

Looks like you're using an enhanced pest file. Its documentation can be found here: https://github.com/usgs/pestpp/blob/master/documentation/pestpp_users_manual.md#418-keyword-and-external-file-control-file-format

In an enhanced pest file, the preceding ++ is dropped entirely from any pestpp-specific settings. For instance, here is a modified enhanced pest control file from the pastas python package that I created when testing my lambda changes to GLM. You can see that lambdas and lambda_scale_fac is provided without any prior symbols. (Not saying you need to use this many lambdas and this many scaling factors. This is probably a bit overkill.)

pcf version=2
* control data keyword
pestmode                                 estimation
noptmax                               100
svdmode                                 1
maxsing                          10000000
eigthresh                           1e-06
eigwrite                                1
lambdas                        1e-09,1e-08,1e-07,1e-06,1e-05,0.0001,0.001,0.01,0.1,1,10,100,1000,10000,100000,1000000,10000000
lambda_scale_fac               0.25,0.5,0.75,1,1.25,1.5,1.75,2
* parameter groups external
pest.pargp_data.csv
* parameter data external
pest.par_data.csv
* observation data external
pest.obs_data.csv
* model command line
python forward_run.py
* model input external
pest.tplfile_data.csv
* model output external
pest.insfile_data.csv

[jbensabat]

Hi and thanks for your comments.
I do not know if you have tested that, but the determination of the lambda parameter is basically a one-dimensional optimization problem (like a line-search). The problem could be formulated as lambda = argmin ||(F(x-(J^TJ+lambdaI)^(-1) J^Tx)||^2. Using a derivative free approach such as Brent algorithm (working on log (lambda) instead of lambda) could provide a faster and more accurate estimate of this parameter.
best
jac

[mmorphew]

Hi Jac,

There are some other methods of determining lambda, as Jeremy and I were discussing above. Getting a more sophisticated set of lambdas to test has been on my list for a while and I should hopefully get to it around July. Pull requests welcome if you are interested in implementing and testing the Brent algorithm for GLM!