Intro to pyEMU¶

This notebook provides a quick run through some of the capabilities of pyemu for working with PEST(++). This run through is very "high level". We will not go into detail, merely show some of the core functionality as a primer for subsequent tutorials. We assume the reader has at least some understanding of PEST(++), common file structures and workflows.

We will make use of an existing PEST(++) interface. You do not need to be familiar with the details of the setup for the purposes of the current tutorial.

Throughout the notebook we will:

• introduce how to access and edit an existing PEST control file and setup using the Pst class.
• explore some of the built-in methods for customing PEST setups and post-processing outcomes.
• introduce geostats in pyemu.
• introduce methods for handling matricies with the Matrix class, as special instances with the Cov and Jco classes.
• introduce classes that facilitate generating and handling parameter and obsveration Ensembles.

Here we do not demonstrate how to setup a PEST interface from scratch. See the "part2_pstfrom_pest_setup" tutorial for a demonstration on how to use the PstFrom class to do so.

1. Admin¶

In the tutorial directory there is a folder named handling_files. This folder contains the PEST(++) dataset which we will access. Feel free to explore it if you wish.

Amongst other things, it includes:

• PEST(++) control file: freyberg_pp.pst
• A Jacobian matrix produced by PEST: freyberg_pp.jcb
• Parameters assigned using pilot points. These are recorded in the pilot point file: hkpp.dat.tpl

2. The Pst class¶

pyEMU encapsulates the PEST control file in the Pst class. Let's read the freyberg_pp.pst control file in the f_d folder.

From the pst instance, we can access all the "*" sections of the control file. Each of these is stored as an attribute. PEST variable names are used for consistency.

For example, the * parameter data section can be accessed by calling pst.parameter_data. This returns the parameter data as a Pandas DataFrame, making it easy to access and edit.

The same for * observation data:

You can access, slice and edit pst.parameter_data or pst.observation_data just as you would a DataFrame. You can add columns, these will not be writen to the control file. pyemu is clever like that.

It's best not to mess around with parameter names (parnme). If you do, you will need to make sure that the corresponding names appear in a .tpl file. (The same applies for observation names and .ins files.)

You can edit values like so:

2.1. Control Data Section¶

The * control data section is handled by a special class that tries to prevent mistakes. pyemu will not allow you to assign an illogical value to a control data variable. This avoids silly mistakes. But it also requires that the user must know what type(s) can be passed to each data variable.

For example, NOPTMAX can only be assigned an integer. If you try to assign text, you will get an error. Try replacing the value 0 in the next cell with the string "zero":

We can inspect all control data values using the pst.control_data.formatted_values attribute:

2.2. PEST++ Options¶

PEST++ options are stored in a dictionary in which the keys are the PEST++ Control Variable name (see the PEST++ user manual for names of these variables and their descriptions). Values must be asigned according to what PEST++ expects as input.

You can access the existing PEST++ options like so:

You can change the values of existing control variables like so:

Or add new PEST++ variables like so:

2.3. Writing the .pst control file¶

All of these edits are kept in memory untill explicitly written to a .pst file. This is accomplished with the Pst.write() method.

The control file can be written as version 1 or 2 (see the PEST++ user manual for descriptions of versions). Unlike the original PEST version 1, version 2 control files have each of the "*" sections stored in external csv files. This makes them easier to access and manipulate, either programatically or using common spreadsheet software. PEST and PEST_HP only accept version 1. Only PEST++ accepts version 2. If there are more than 10k parameters, version 2 is written by default.

You can specify the version by passing the relevant argument. Run the cells below then inspect the folder to see the differences.

2.4. Adding Parameters/Observations from .tpl/.ins files¶

In another tutorial we demonstrate the use of the PstFrom to automate construction of a complete PEST interface, starting from scratch.

Alternatively, it is possible to construct a PEST control file from template (.tpl) and instruction (.ins) files. The Pst class includes methods to read .tpl/.ins files and add parameters/observations to the control file. This enables construction of a PEST dataset, or simply adding new ones to an existing control file. It is particularily usefull for pesky model input/output files with inconvenient file structures (e.g. that are not array or tabular formats).

The cell bellow writes a "dummy" template file to the tutorial folder for demonstration purposes.

Adding the parameters from "special_pars.dat.tpl" to the control file is a simple matter of calling the .add_parameters() method. This method adds the parameters to the pst.parameter_data section, and updates other relevant sections of the pst control file; it also returns a DataFrame. Note that parameter names come from the .tpl file. However, intial values, bounds and group name are assigned default values; you will need to specify correct values before writing the control file to disk.

You can then adjust parameter data details:

2.5. Tying Parameters¶

We may on ocasion need to tie parameters in the control file. In the pyemu world, tied parametes are specified in the Pst.parameter_data dataframe. Start by adding a partied column and, for parameters you want to tie, changing "partrans" to "tied" and adding the name of the parameter to tie to in the "partied" column.

We will demonstrate step-by-step by tying "special_par2" to "special_par1" (the parameters we just added from the .tpl file):

3. Utilities¶

pyemu has several built-in methods to make your PEST-life easier. Several of these handle similar tasks as utilities from the PEST-suite, such as adjusting observation weights and assigning prior information equations. Others provide usefull tables or plots that summarize details of the PEST setup and/or outcomes.

3.1. Par and Obs Summaries¶

You can access quick summaries of observation and paramaeter names, group names, etc thorugh the respective pst attributes:

You can write a parameter or observation summary table wth the Pst.write_par_symmary_table() and Pst.write_obs_symmary_table() methods, respectively. Quite usefull when preparing those pesky reports.

These methods return a Pandas DataFrame and (by default) write the table to an external file. Parameters and observations are summarized by group name.

3.2. Phi and residuals¶

The Pst class tries to load a residuals file iduring construction. It looks for a file in the same folder as the control file and with the same base name, but with the extension ".rei". Alterantaively, you can specify the name of the residual file when constructing the Pst. (e.g. pyemu.Pst("controlfile.pst", resfile="residualfile.rei"))

If that file is found, you can access some pretty cool stuff. The Pst.res attribute is stored as a Pandas DataFrame.

Of course, all of this relies on PEST(++) having been run at least once before hand to record the residuals file. For the purposes of this tutorial, we have already done so. When we constructed pst at the beggining of this notebook, pyemu also loaded the residuals file.

Inspect it by running the cell bellow. As you can see, the DataFrame lists all observations and group names, their modelled and measured values, weights and of course the residual:

A somewhat clunky (and meaningless) look at everyones favourite "good fit" plot:

Or a clunky look at the residuals for selected observations:

There are built in routines for some common plots. These can be called with the Pst.plot() method and specifying the kind argument. For example, a 1to1 plot for each observation group:

The weighted sum of square residuals (Phi) is also stored in the respective Pst attribute.

We can access the components of the objective function as a dictionary. These allow us to breakdown the contributions to Phi from each observation group:

They can also be displayed with a plot, like so:

These values can be recalculated for different observation weights by simply changing the weights in the pst.observation_data. No need to re-run PEST!

3.3. Adjusting Weights for "Visibility"¶

Prior to estimating parameters using PEST(++), a user must decide how to weight observations. In some cases, it is wise to weight observations strictly according to the inverse of the standard deviation of measurement noise. Certainly, observations of greater credibility should be given greater weights than those of lower credibility.

However, when history-matching, this approach can result in the loss of information contained in some observations, due to their contribution to the objective function being over-shadowed by the contribution from other obsevrations. An alternative approach is to weight observations (or observation groups) such that, at the start of the history-matching process, they each contribute the same amount to the objective function. The information content of each of these groups is thereby given equal right of entry to the parameter estimation process. (This matter is extensively discussed in the PEST Book.)

As stated above, a practical means of accommodating this situation is to weight all observation groups such that they contribute an equal amount to the starting measurement objective function. In this manner, no single group dominates the objective function, or is dominated by others; the information contained in each group is therefore equally “visible” to PEST(++).

The Pst.adjust_weights() method provides a mechanism to fine tune observation weights according to their contribution to the objective function. (Side note: the PWTADJ1 utility from the PEST-suite automates this same process of "weighting for visibility".)

Some caution is required here. Observation weights and how these pertain to history-matching versus how they pertain to generating an observation ensemble for use with pestpp-ies is a frequent source of confusion.

• when history-matching, observation weights listed in the control file determine their contribution to the objective function, and therefore to the parameter estiamtion process. Here, observation weights may be assigned to reflect observation uncertainty, the balance required for equal "visibility", or other modeller-defined (and perhaps subjective...) measures of observation worth.
• when undertaking generating an observation ensemble, weights should reflect the inverse of the standard deviation of measurement noise. Unless instructed otherwise, pestpp-ies will generate the observation ensemble using observation weights in the PEST control file. Therefore, when history-matching with pestpp-ies and using weights that do not reflect observation uncertainty, it is important to provide pestpp-ies with a previously prepared observation ensemble (we will demonstrate this further on).

OK, so let's adjust some observation weights using Pst.adjust_weights(). This method allows us to adjust weights for individual observations, or for entire observation groups. We do so by passing a dictionary with observation names (or group names) as keys, and the correspding value they contribute to the objective function as values.

Let's get started.

Now, let's create the dictionary of non-zero weighted observation groups. We will specify that we want each group to contribute a value of 100 to the objective function. (Why 100? No particular reason. Could just as easily be 1000. Or 578. Doesn't really matter. 100 is a nice round number though.)

Now we can simply pass this dictionary as an argument:

And voila. Run the cell below to see that phi components form each group are (roughly) 100. The same approach can be implemented for individual observations (see the obs_dict argument in pst.adjust_weights()).

3.4. Discrepancy based weight adjustment¶

In a perfect (model and algorithm) world, we would acheive a final objective function that is equal to the number of (non-zero weighted) observations. But because of model error and simplifying assumptions in the algorithms we use for history matching, this is rarely the case. More often, the final objective function is much larger than the number of observations. This implies that we were not able to "fit" as well as we thought we could (where "thought" is incapsulated in the observations weights in the control file, representing the inverse of measurment noise). This really matters when we do posterior uncertainty analyses following a PEST run (this will be discussed further in the FOSM and data-worth notebooks).

The simpliest way to try to rectify this situation is to adjust the weights in the control file so that the resulting contribution to the objective function from each observation (or optional observation group) is equal to 1 (or the number of members of the group). This is related to Morozov's discrepancy principal (google it!). pyEMU has a built in routine to help with this: Pst.adjust_weights_discrepancy() - great name!

Note 1: dont make this adjustment until after you are through with history matching! The point is for weights to represent the inverse of observation uncertainty, which includes both measurement error and model error.

Note 2: the PWTADJ2 utility from the PEST-suite acomplishes a similar task.

So we were expecting Phi to be equal to nnz_obs (number of non zero observations). This did not happen due to the original_ceiling argument being set to True (which is the default value).

What this means is that, for some observations, weights would have to be increased to achieve a contribution to Phi of 1.0. Which is illogical if the original weight is assumed to be the inverse of the measurement noise. In some cases, this requirement may not apply; such as when observatons are weighted for visibility.

4. Geostatistics in pyEMU¶

The pyemu.geostats module provides tools for implementing geostatistics in the pyemu world. These have similar functions to the PPCOV* utilities from the PEST suite.

In the PEST world, geostatistics are used (1) to describe how parameters vary in space and/or time and (2) when assigning parameter values from pilot points to model input files. Front-end users will mostly be interested in the former. A separate tutorial notebook delves into geostatistics in greater detail (see the "intro to geostats" notebook). The "part1 pest setup" and "part2 pstfrom pest setup" tutorials provide examples of their use in the wild.

At the heart of geostatistics is some kind of model expressing the variability of properties in a field. This is a "variogram". pyemu supports spherical, exponential and gaussian variograms.

A pyemu geostatistical structure object that mimics the behavior of a PEST geostatistical structure (see the PEST manual for details). The object contains variogram instances, information on parameter transforms and (optionally) nugget information. Multiple nested variograms are supported.

Using the geostatistical structure (or a variogram), a covariance matrix describing the statistical relationship between spatialy distributed parameters can be constructed. These covariance matrices are used (1) to describe prior parameter probability distributions, and (2) specifying the weights of prior information equations as a regularisation device.

For example, let's create a covariance matrix for a set of pilot point parameters from the gs structure. This requires a list of X and Y coordinates. Let's make use of the coordinates from a pilot points file in the f_d folder.

First, read the pilot points file to get X and Y values:

Now simply pass the respective values from the pilot points file as arguments to gs.covariance_matrix(). Conveniently, we can also take the parameter names from the pilot point file. This returns a pyemu.Cov object (which we will discuss later). The same can be complished with a single variogram (e.g. v.covariance_matrix())

We can do the same thing for a 1-D sequence (think time-series), to get a covariance matrix for parameters distributed in time:

5. Prior Information Equations¶

The mathematical term for the process through which a unique solution is sought for a nonunique inverse problem is “regularisation”. The goal of regularised inversion is to seek a unique parameter field that results in a suitable fit between model outputs and field measurements, whilst minimizing the potential for wrongness in model predictions. That is, out of all the ways to fit a calibration dataset, regularized inversion seeks the parameter set of minimum error variance. (Regularisation is discussed in greater detail in a separate notebook.)

One way to seek a parameter field of minimum error variance is to seek a parameter field that allows the model to fit the calibration dataset, but whose values are also as close as possible to a set of “preferred parameter values”. Ideally, preferred parameter values should also be initial parameter values as listed in the “parameter data” section of the PEST control file. These preferred parameter values are normally close to the centre of the prior parameter probability distribution. At the same time, scales and patterns of departures from the preferred parameter set that are required for model outputs to fit a calibration dataset should be achieved in ways that make “geological sense”.

PEST provides a user with a great deal of flexibility in how Tikhonov constraints can be introduced to an inversion process. The easiest way is to do this is through the use of prior information equations. When prior information equations are employed, Tikhonov constraints are expressed through preferred values that are assigned to linear relationships between parameters. (Equality is the simplest type of linear relationship.) Weights must be assigned to these equations. As is described in PEST documentation, when PEST is run in “regularisation” mode, it makes internal adjustments to the weights that are assigned to any observations or prior information equations that belong to special observation groups that are referred to as “regularisation groups”.

Note: in a similar manner, the PEST-utilities ADDREG1 and ADDREG2 also automate the addition of prior information equations to a PEST-control file.

5.1. Preffered value or Zero Order Tikhonov¶

pyemu provides utilities to apply preferred value prior equations to a PEST control file. Note though, pyemu doesn't call it "preferred value"! Rather, it uses the mathematical term "Zero Order Tikhonov".

Before we do so, we can inspect the control file * prior information section. As you may have guessed, this is accessed using the pst.prior_information attribute, which returns a DataFrame:

Now, as you can see below, prior information equations have been added for all adjustable parameters. The "preferred value" is obtained from the parameter initial value in the * parameter_data section. The weight is calculated from the parameter bounds (this behaviour can be changed with arguments in pyemu.helpers.zero_order_tikhonov()).

5.2. Preferred difference or First Order Pearson Tikhonov¶

We may wish (almost certaintly) to express a preference for similarity between parameters. For example, hydraulic properties of two points close together are more likley to be similar to each other, than two points which are far apart. We describe this relationship using geostatistics, encapsulated in a covariance matrix.

As previously described, pyemu.geostats module provides tools for generating such matrices. The PEST suite also includes many utilities for this purpose (see the PPCOV* set of utilities.) Let's use the cov covariance matrix we constructed earlier for the set of pilot points.

Now, we can assign prior information equations for preferred difference. Note that the preferred difference = 0, which means our preferred difference regularization is really a preferred homogeneity condition! If observation data doesn't say otherwise, parameters which are close together should be similar to each other.

The weights on the prior information equations are the Pearson correlation coefficients implied by the covariance matrix.

5.3. Custom Prior Information Equations¶

Pst.add_pi_equation() is a helper to construct prior information equations. We demonstrate below using the "special_par"s we added earlier.

6. Matrices¶

The pyemu.Matrix class is the backbone of pyemu. The Matrix class does some fancy things in the background to make manipulating matrices and linear algebra easier. The class overloads all common mathematical operators and also uses an "auto-align" functionality to line up matrices for multiplication, addition, etc. The pyemu.Jco class provides a wrapper to deal with PEST Jacobian matrices. It functions the same as the Matrix class.

The pyemu.Cov class is a special class designed specifically to handle covariance matrices. (We have already seen examples of the pyemu.Cov earlier.) The Cov class makes some assumptions, such as symmetry (and accordingly that matrix rows and columns names are equal) and provides additional functionality.

All classes provide functionality to record matrices to, as well as read from, external files in PEST-compatible formats (e.g. ASCII or binary). Cov can also write PEST uncertainty (*.unc) files. The Cov class has additional functionality that allows covariance matrices to be constructed from information contained in PEST control files (e.g. from parameter or observation bounds).

As you can see below, the Cov object retains the parameter names as row and column names.

We can check if the matrix is diagonal with the .isdiagonal attribute.

How about the cov covariance matrix we generated earlier for spatialy correlated parameters?

Recall that, a matrix is "diagonal" when all the entries off the diagonal are zero. This means that elements of the matrix are uncorrelated. By generating a covariance matrix from parameter data, the only information we have is the uncertainty of each individual parameter (expressed by the parameter bounds). There is no information on correlation between parameters. Therefore cov_pb.isdiagonal is True. When we generated cov, we specified correlation between parameters; so, off-diagonal entries in cov are non-zero. Therefore, cov.isdiagonal is False.

The values of the matrix are accessed in the .x attribute:

Note that parcov.x is 1-dimensional. Again, this is because parcov is "diagonal". So .x only returns the diagonal entries. On the other hand, cov.x is 2-dimensional:

You can access the full 2-dimensional matrix with the .as_2d attribute (see all the off-diagonals are zero):

It may be easier to visualize this in a plot. Because parcov has many parameters it may be hard to see. So let's first use the get() method to get a submatrix. We can use this method to extract specific row and column names. We will get the same pilot point parameters that are in the cov matrix to compare.

The singular values are the diagonal entries of the matrix, arranged in descending order. These can be accessed with the .s attribute. This still returns a Matrix object (technically a vector). Values in the matrix can be accessed in the same fashion as described above.

"Right singular vectors" are obtained with .v attribute:

Matrices can be written to, and read from, external files. These can be ASCII or binary type files. As you saw above, we can access numpy arrays with the .as_2d attribute. They can also be converted to a Pandas Dataframe.

Jacobian matrices, recorded by PEST(++) during parameter estimation, can be read and manipulated with the Jco class. This calss has the same functioanlity as the Matrix class.

7. Linear Analysis or FOSM¶

FOSM stands for "First Order, Second Moment". You may also see this referred to as "linear analysis" (e.g. in PEST documentation). We will delve into FOSM in more detail in another tutorial. Here we merely provide a brief introduction.

The pyemu.Schur object is one of the primary object for FOSM in pyemu. Instantiating a Schur object requires, at minimum, a Jacobian matrix and a PEST(++) control file. From these, pyemu builds the prior parameter covariance matrix (from parameter bounds) and the observation noise covariance matrix (from observation weights). Alternatively, the parameter and observation covariance matrices can be provided explicitly.

Optionaly, observation names can be specified as forecasts. The Schur object extracts the corresponding rows from the Jacobian matrix to serve as forecast sensitivity vectors.

As we saw earlier, there is a Jacobian matrix in the tutorial folder, recorded in the file named freyberg_pp.jcb.

The prior parameter covariance matrix is stored in the .parcov attribute:

The same for the observation noise covariance matrix:

The posterior parameter covariance matrix is calculated and stored in the .posterior_parameter attribute:

Let's record the prior and posterior coavariance matrix to external files (we will use these later):

The Schur object found the "++forecasts()" optional pestpp argument in the control, found the associated rows in the Jacobian matrix file and extracted those rows to serve as forecast sensitivity vectors. lterantively, we can also pass a list of observation names to use as forecasts when instantiating the Schur object. (e.g.pyemu.Schur("jacobian.jcb", forecasts=[obsname1, obsname2, etc...]).)

pyemu makes FOSM easy to undertake. It has lot's of usefull functionality. We will not go into further detail here. See the "intro to fosm" tutorial for a deeper dive.

8. Ensembles¶

The pyemu.ParameterEnsemble and pyemu.ObservationEnsemble ensemble classes store parameter or observation ensembles, respectively. These classes are DataFrames under the hood, allowing us to use all the baked in Pandas conveniences.

We will focus on the pyemu.ParameterEnsemble class. Similar concepts apply to pyemu.ObservationEnsemble.

Ensembles can be read from (and written to) external files, or generated by pyemu using one of several methods. These methods draw stochastic values from (multivariate) (log) gaussian, uniform and triangular distributions for parameters in a Pst control file. Much of what we do is predicated on the gaussian distribution, so let's use that here by generating an ensemble with the .from_gaussian_draw() method.

Note that these draw methods use initial parameter values in the control file (the Pst.parameter_data.parval1 attribute) the $\boldsymbol{\mu}$ (mean) prior parameter vector. Unless otherwise specified, the parameter bounds are assumed to represent six standard deviations of the parameter probability distribution (i.e. ~99% of parameter values lie within the lower and upper bounds).

We can express parameter correlation by passing a covariance matrix to the cov argument. In doing so, the covariance matrix describes the second moement (the standard deviation) of the gaussian distribution. This allows us to draw parameter ensembles respecting parameter covariance.

Because the Ensemble is stored as a DataFrame, we can easily access it. Each row is an individual realisation. Each column is a parameter.

Let's plot a histogram of values generated for one of the parameters:

As you can see, parameter bounds may be violated when drawing the ensemble. The .enforce() method goes through it and makes sure that bounds are respected. It does this by assigning the bound value to any parameter values which exceed the bound.

As you can see in the subsequent plot, bounds are now respected:

There are also built in functions to automate plotting:

You can record an ensemble to external files (.csv of binary):

And of course, read from external files (csv, binary) or even from pandas DataFrames:

We can even form an empirical covariance matrix from an ensemble!

8.1 Bayes Linear Monte Carlo¶

We can use the bayes linear posterior parameter covariance matrix (aka Schur compliment) to "precondition" the realizations using linear algebra so that they hopefully yield a lower phi. The trick is we just need to pass this posterior covariance matrix to the draw method. Note this covariance matrix is the second moment of the posterior (under the FOSM assumptions) and the final parameter values is the first moment.

In other GMDSI educational material, and in the PEST Roadmaps, this approach is sometimes referred to as "using the linearized posterior parameter distribution".

Applying Bayes linear Monte Carlo requires that we have previously calibrated a model and calculated a post-calibration Jacobian matrix. From the Jacobian, we can obtain the post-calibration parameter covariance. The "calibrated" parameter values represent the mean of the posterior parameter probability distribution. By centering the distribution on values that already provide a good fit with measurment data, we are increasing the likelihood that the realisatiosn that we draw will also fit measured data well. This can be usefull in reducing subsequent history-matching computation time.

Or for comparing histograms from several ensembles. (Can also be used for observation ensembles; see additional method arguments):