Some thoughts on PEST integration

[jtwhite79]

The primary “forward run” process we need for PEST integration involves:

• load an existing PASTAS model instance
  • “current” PASTAS model structure (ie stress models, noise models, burn-in/spin-up settings, etc),
  • “current” forcing time series
  • “current “parameter” values
• Run the PASTAS model forward (including any burn-in/spin up)
• Save the time series contributions from each stressmodel, noise model, as well as the aggregate result (ie the standard output) and other derived outputs like signatures, etc. This would be all the simulated results, not just the results that align with actual observations – in the pest world, we commonly carry a lot more simulated results than we have observations for, which is nice for diagnostics, etc.

This is what Ive tried to do when I was hooking up pastas to pest before. With this in mind, it would sweet if we could pickle a pastas model instance (or save all the settings/definitions in a config, which I think is in the works) so that we don’t need to have all the py code to re-create the pastas model instance each time we run the model. I had some hackery to force pastas to run over the same time period every time - It would awesome if the pastas model would just run from a start to end datetime without supplying any observed series...

Also, I felt like I was living dangerously trying to hot-wire pastas to simply use parameter values from a file (a file that pest fills with the current values at runtime), rather than using the internal PASTAS calibration routine(s). So it would be nice if there was a clean and standard way to tell a pastas model to read parameter values from a plain-ole CSV file (with name-value pairs). For this to work, I had to hope that pastas was using the parameter info in the order I supplied it (and the “A” and “a” parameter naming is a problem in the pest world bc pest is case insensitive).

I think these are the primary obstacles to easily making a big ‘ole bowl of PEST-o-PASTA haha!

Maybe some of these things are already dealt with and/or I was just doing it wrong (highly likely) so please let me know if Ive missed somethings along the way...

[dbrakenhoff]

Thanks for starting this discussion @jtwhite79! This is something we're very interested in, and I'm sure my colleague @martinvonk has more thoughts and ideas on this topic. But I had some thoughts and questions I wanted to ask/share, so here goes:

I'm assuming the PEST run is performed by some Python script? If so, the following code snippets might be useful, if you weren't already doing it that way...

Concerning the storing/loading of existing Pastas models:

Pastas supports saving/loading JSON-style files. It is a bit of a pain to modify stuff in this file directly, so it is recommended to do that stuff via the Pastas Model interface.

ml.to_file("model.pas")  # store as JSON file
ml = ps.io.load("model.pas") # load JSON file

There is something of an alternative, through the pastastore package, which allows storing Pastas Models and time series separately in some sort of database-like structure (which can be some folder on your PC). This package allows you to reading and writing models as dictionaries (which is very fast), instead of building the full Pastas Model instances each time. That would allow modifying parameters/settings quickly. However, since you need to run the Pastas Model at some point anyway, it might be more efficient to just use the Pastas Model class directly? The pastastore package also has a YAML interface for building Pastas models, explained here, which is a mini-language that allows you to build Pastas Models and define settings in a human-readable text file. I'm not sure if these methods are fully what you are looking for, but they may already contain some of the functionality we would need to better support PEST runs?

With respect to the parameters, setting parameters is done using:

ml  # some pastas model

for param, value in csv_value_pairs.items():
    ml.set_parameter(param, optimal=value)

# or maybe since PEST is case-insensitive
for i, (param, value) in enumerate(csv_value_pairs.items()):
    ipar = ml.parameters.index[i]
    ml.set_parameter(ipar, optimal=value)

And pastas does have a strict parameter order, based on the order in which StressModels were added to the original model, so you can count on that being the case.

Running a model without using the internal calibration routine can also be done in multiple ways:

1. Supply ml.simulate(p=p) with an array p containing parameter values (in the order Pastas expects).
2. Set the optimal parameters using ml.set_parameter(<name>, optimal=<value>). Then call ml.simulate().

Many methods in pastas support tmin/tmax arguments so that should force Pastas to compute the results for the simulation period you're interested in.

Once we've settled on an approach to better supporting PEST runs, it would definitely be a good idea to document a clear and intuitive way to set parameters or tweak a model when using external calibration routines. Or maybe we can create some kind of PestSolver class, that basically just reads/writes external csv files...? Just a random thought.

Anyway, I'm curious if you were already using the logic above? If not, would that make it running PEST-o-PASTA (great name by the way :p ) easier? Let us know!

[jtwhite79]

That is more or less how I was getting at the problem, just wanted to make sure that assuming that order of parameters was safe - I think saw a warning somewhere about messing with the parameter dataframe that a pastas model has attached to it...but maybe I was also using an older version of pastas.

I like the idea of a helper class to as an entry point where we could automate and hide all the gory details of preparing for and running pest/pest++ would be great. There is a monster class in pyemu called PstFrom (which @martinvonk has been experimenting with I think) that is supposed to help setup for pest by building all the components one needs to run pest/pest++, so I think if we could link from what a pastas model knows about the parameters (their initial values and ranges) and the observations (values and noise expectations), then we could make things pretty easy...altho there might be some trickery around how do deal with forcing uncertainty, but we can cross that bridge when we come to it.

(I also cant take credit for that name - I stole it from @raoulcollenteur !)

[dbrakenhoff]

That is more or less how I was getting at the problem, just wanted to make sure that assuming that order of parameters was safe - I think saw a warning somewhere about messing with the parameter dataframe that a pastas model has attached to it...but maybe I was also using an older version of pastas.

Yea, you're right about the recommendation about messing with the parameters dataframe directly. That is usually not recommended, because certain parameters (initial values and bounds) are also stored in the StressModels they're associated with, and are read from there when (re-)initializing the model. Also setting the parameters directly means the bounds might be violated.

For the PEST case however, the optimal parameters are always stored only in the Model.parameters dataframe, so it's okay to modify the dataframe directly (and if the bounds are used, to respect those somehow as well).

altho there might be some trickery around how do deal with forcing uncertainty, but we can cross that bridge when we come to it.

This is an interesting step we need to think about. Maybe some of the pastastore logic could come in handy here, since it stores time series and models separately that allows you to modify the forcing time series, which are pulled back into the model when you load it from the database.

Anyway, plenty to think about! I'll get @martinvonk and the rest of the Pastas developers to look at this, and start thinking of a way to implement this helper class. We'll post updates here, and feel free to share any other ideas/implementations/code if you think of anything.

[martinvonk]

Alright, I've got a working example on the pastas-plugins repository: https://github.com/pastas/pastas-plugins/blob/pestas/docs/examples/pest_pastas.ipynb. All logic is in the solver.py file which contains the following:

Pastas has a BaseSolver class from which each other solver (LeastSquares, LmFitSolve etc.) inherits. This BaseSolver contains some method that are essential. I created a PestSolver class that inherits from the BaseSolver. The PestSolver contains all the logic for setting up a pastas model and the pest file structure via pyemu.

• PestSolver.setup_model does three things:

1. Save the Pastas model to a pas file. This pas file can be loaded again into a model during a pest optimization run to simulate with updated parameters.
2. Save the (to be) optimized parameters to a CSV (parameters_sel.csv). These can be loaded into the model during a pest optimization run. Note that pest parameters are case insensitive so _A is renamed to _g (for gain). Otherwise _A and _a parameters are confused.
3. Save the observations to a CSV (simulation.csv). This file is needed to load the original observations in pest, and write the newly simulated values of each optimization run to.

Note: I know that the parameters are also available in the pas file but I figured this was faster and pyemu nice support for reading CSV files.

• PestSolver.setup_files sets up the pest structure via the pyemu PstFrom class:

1. PstFrom class and required directories are initialized
2. Parameters are added to the class
3. Observations are added to the class
4. A python script file _run_pastas_model.py (now still separate from script that's run pyemu issue 527) is added that loads the Pastas model from the pas file, updates the parameters from pest optimization run and simulates again with the new parameters.
5. Write a pest control file

• PestSolver.run starts the pest optimization with the executable.

There are many different versions of pest(++); with different output files and required handling. For each pest version we need another class with a solve method. On the pastas-plugins repository I now have two classes; PestGlmSolver and PestIemSolver from which only the PestGlmSolver currently works. All those classes inherit from the PestSolver class, because they only differ in the solving and interpretation of the results. The solve method is called with pastas.Model.solve(). The PestGlmSolver solve method calls PestSolver.setup_model, PestSolver.setup_files, and PestSolver.run. If the optimization is successful, the optimal parameters are read from a pest output file. Handling of all the other important output files from pest (covariance matrices etc.) is yet to be implemented for this class. That's also the difficult part for me because I don't know how to interpret them, yet 🙈. It also differs per pest version so that's where we can learn from your expertise @jtwhite79 👼.

I hope this makes things a bit more clear for everyone, and we can work from here. Next step is (I think?) to compare the result from the notebook to a simple pastas LeastSquares model result (including the standard deviation etc.). That should iron out some bugs that are still present in my code. I don't know if pestpp-glm is the best/most simple version for such a comparison?

P.S. For me, what helped was a small example I created to compare SciPy least_squares to the pyemu pestpp-glm setup. This showed me more clearly how to set up a pest optimization with pyemu. Notebook can be found here: https://github.com/pastas/pastas-plugins/blob/pestas/docs/examples/pest_linear_regression.ipynb

[martinvonk]

Pest also does not like negative parameters, since it does log transformations by default. We don't have those negative parameters a lot, But one example: _f in the linear recharge model is negative... We have to think about how we want to handle that.

[cnicol-gwlogic]

Pest also does not like negative parameters, since it does log transformations by default. We don't have those negative parameters a lot, But one example: _f in the linear recharge model is negative... We have to think about how we want to handle that.

FYI all Pest versions offer a user input parameter "offset" (specified per paremeter), so for example, you can specify an offset such that your _f parameter never becomes negative. Eg if the parameter minimum bound is -3, make that offset value -4. Then pest deals in offset space (param_value+4 - your bounds/initial values are all specified as +4), but pest internally adds your -4 to the parameter when writing to model input files.

Anyway, I am keen to try this pest interface you've made @martinvonk. I see you have a branch where you are working on this. I'll give it a shot in the coming month or two and contribute where I can offer anything useful. Appreciate the work on it.

[martinvonk]

Hi @cnicol-gwlogic,

Awesome, you're willing to test it and contribute :) Please note that the branch is still very much work in progress. Also, we're still testing and comparing the outcomes of pestpp and default pastas, and we are getting some results we can't explain yet.

Good suggestion for the offset. What I'm now thinking of is just multiplying the parameter with both a negative pmin and pmax with -1 when giving the parameter to pest. This way the standard deviations are also more easily interpretable. Additionally, this helps with our other negative parameters; _A if up=False.