<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Model background
Here is an example based on the model of Freyberg, 1988.  The synthetic model is a 2-dimensional MODFLOW model with 1 layer,  40 rows, and 20 columns.  The model has 2 stress periods
Step1:  The plot shows the Freyberg (1988) model domain.  The colorflood is the hydraulic conductivity ($\frac{m}{d}$).  Red and green cells coorespond to well-type and river-type boundary conditions. Blue dots indicate the locations of water levels used for calibration.
Using pyemu
Step2:  Drawing from the prior
Now we need a prior-realized ParameterEnsemble, which stores a pandas.DataFrame under the hood.
draw
The ParameterEnsemble class has several draw type methods to generate stochastic values from (multivariate) (log) gaussian, uniform and triangular distributions.  Much of what we do is predicated on the gaussian distribution, so let's use that here.
The gaussian draw accepts a cov arg which can be a pyemu.Cov instance.  If this isn't passed, then the draw method constructs a diagonal covariance matrix from the parameter bounds (assuming a certain number of standard deviations represented by the distance between the bounds - the sigma_range argument)
Step3:  draw also accepts a num_reals argument to specify the number of draws to make
Step4:  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.  To change that, we need to update the parameter values in the control file
Step5:  plotting
Since ParameterEnsemble stores a pandas.DataFrame, it has all the cool methods and attributes we all love.  Let's compare the results of drawing from a uniform vs a gaussian distribution.  The actual dataframe is stored under the private attribute ParameterEnsemble._df
Step6:  The gaussian histo go beyond the parameter bound - bad times.  Luckily, ParameterEnsemble includes an enforce method to apply parameter bounds
Step7:  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 (which we parrep'ed into the Pst earlier)
<ASSISTANT_TASK:>
Python Code:
import flopy
# load the model
model_ws = os.path.join("Freyberg","extra_crispy")
ml = flopy.modflow.Modflow.load("freyberg.nam",model_ws=model_ws)
# Because this model is old -- it predates flopy's modelgrid implementation. 
# And because modelgrid has been implemented without backward compatability 
# the modelgrid object is not constructed properly. 
# - We will use some sneaky pyemu to get things to how they should be 
import pyemu
sr = pyemu.helpers.SpatialReference.from_namfile(
    os.path.join(model_ws, ml.namefile), 
    delc=ml.dis.delc, 
    delr=ml.dis.delr
)
ml.modelgrid.set_coord_info(
    xoff=sr.xll,
    yoff=sr.yll,
    angrot=sr.rotation,
    proj4=sr.proj4_str,
    merge_coord_info=True,
)
# plot some model attributes
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(111,aspect="equal")
ml.upw.hk.plot(axes=ax,colorbar="K m/d",alpha=0.3)
ml.wel.plot(axes=ax)  # flopy possibly now only plots BCs in black
ml.riv.plot(axes=ax)
# plot obs locations
obs = pd.read_csv(os.path.join("Freyberg","misc","obs_rowcol.dat"),delim_whitespace=True)
obs_x = [ml.modelgrid.xcellcenters[r-1,c-1] for r,c in obs.loc[:,["row","col"]].values]
obs_y = [ml.modelgrid.ycellcenters[r-1,c-1] for r,c in obs.loc[:,["row","col"]].values]
ax.scatter(obs_x,obs_y,marker='.',label="obs")
#plot names on the pumping well locations
wel_data = ml.wel.stress_period_data[0]
wel_x = ml.modelgrid.xcellcenters[wel_data["i"],wel_data["j"]]
wel_y = ml.modelgrid.ycellcenters[wel_data["i"],wel_data["j"]]
for i,(x,y) in enumerate(zip(wel_x,wel_y)):
    ax.text(x,y,"{0}  ".format(i+1),ha="right",va="center", font=dict(size=15), color='r')
ax.set_ylabel("y")
ax.set_xlabel("x")
ax.add_patch(rect((0,0),0,0,label="well",ec="none",fc="r"))
ax.add_patch(rect((0,0),0,0,label="river",ec="none",fc="g"))
ax.legend(bbox_to_anchor=(1.5,1.0),frameon=False)
plt.savefig("domain.pdf")
Explanation: Model background
Here is an example based on the model of Freyberg, 1988.  The synthetic model is a 2-dimensional MODFLOW model with 1 layer,  40 rows, and 20 columns.  The model has 2 stress periods: an initial steady-state stress period used for calibration, and a 5-year transient stress period.  The calibration period uses the recharge and well flux of Freyberg, 1988; the last stress period use 25% less recharge and 25% more pumping.
The inverse problem has 761 parameters: hydraulic conductivity of each active model cell, calibration and forecast period recharge multipliers, storage and specific yield, calibration and forecast well flux for each of the six wells, and river bed conductance for each 40 cells with river-type boundary conditions.  The inverse problem has 12 head obseravtions, measured at the end of the steady-state calibration period.  The forecasts of interest include the sw-gw exchange flux during both stress periods (observations named sw_gw_0 and sw_gw_1), and the water level in well cell 6 located in at row 28 column 5 at the end of the stress periods (observations named or28c05_0 and or28c05_1).  The forecasts are included in the Jacobian matrix as zero-weight observations. The model files, pest control file and previously-calculated jacobian matrix are in the freyberg/ folder
Freyberg, David L. "AN EXERCISE IN GROUND‐WATER MODEL CALIBRATION AND PREDICTION." Groundwater 26.3 (1988): 350-360.
End of explanation
import pyemu
pst = pyemu.Pst(os.path.join("Freyberg","freyberg.pst"))
Explanation: The plot shows the Freyberg (1988) model domain.  The colorflood is the hydraulic conductivity ($\frac{m}{d}$).  Red and green cells coorespond to well-type and river-type boundary conditions. Blue dots indicate the locations of water levels used for calibration.
Using pyemu
End of explanation
pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=200)
Explanation: Drawing from the prior
Now we need a prior-realized ParameterEnsemble, which stores a pandas.DataFrame under the hood.
draw
The ParameterEnsemble class has several draw type methods to generate stochastic values from (multivariate) (log) gaussian, uniform and triangular distributions.  Much of what we do is predicated on the gaussian distribution, so let's use that here.
The gaussian draw accepts a cov arg which can be a pyemu.Cov instance.  If this isn't passed, then the draw method constructs a diagonal covariance matrix from the parameter bounds (assuming a certain number of standard deviations represented by the distance between the bounds - the sigma_range argument)
End of explanation
pe.head()
Explanation: draw also accepts a num_reals argument to specify the number of draws to make:
End of explanation
pst.parrep(pst.filename.replace(".pst",".par"))
pst.parameter_data.parval1
pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=200)
pe.head()
Explanation: 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.  To change that, we need to update the parameter values in the control file:
End of explanation
pe = pyemu.ParameterEnsemble.from_uniform_draw(pst=pst,num_reals=1000)
ax = plt.subplot(111)
pe._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=1000)
pe._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
Explanation: plotting
Since ParameterEnsemble stores a pandas.DataFrame, it has all the cool methods and attributes we all love.  Let's compare the results of drawing from a uniform vs a gaussian distribution.  The actual dataframe is stored under the private attribute ParameterEnsemble._df:
End of explanation
pe = pyemu.ParameterEnsemble.from_uniform_draw(pst=pst,num_reals=1000)
ax = plt.subplot(111)
pe._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=1000)
pe.enforce(how="reset")
pe._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
Explanation: The gaussian histo go beyond the parameter bound - bad times.  Luckily, ParameterEnsemble includes an enforce method to apply parameter bounds:
End of explanation
# get the list of forecast names from the pest++ argument in the pest control file
jco = os.path.join("Freyberg","freyberg.jcb")
sc = pyemu.Schur(jco=jco)
pe_post = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,cov=sc.posterior_parameter, num_reals=1000)
pe_post.enforce()
ax = plt.subplot(111)
pe._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
pe_post._df.loc[:,"rch_1"].plot(kind="hist",bins=20,ax=ax,alpha=0.5)
Explanation: 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 (which we parrep'ed into the Pst earlier)
End of explanation
<END_TASK>