\begin{center}
{\bf Due Date: 17:00, Monday, October 25, 2010} \\
{\bf TA: Xukai Shen (xukai@sep.stanford.edu)}
\end{center}
\title{Lab 4 - Missing data estimation and multiple realizations}

\def\figdir{./Fig}
\author{Hongwu Emperor}
\email{hongwu@ming.com}
\righthead{Missing data \& multiple realizations}
\lefthead{GEE - Lab 4}

\begin{abstract} 
	Geostatistics has gained much recent attention for its ability to solve
	missing data problems not with a single ``answer'', but with a suite of 
	so-called equiprobable models, or models that fit the data equally well.
	The geostatistical approach gives better results in fluid flow simulation
	(for instance), and also provides a useful tool for multiple scenario
	analysis. \\

	\cite{Clapp.sep.105.bob3} showed that multiple realizations can easily
	be incorporated into the existing SEP framework to solve missing data problems. \\

  This lab has two parts.  In the first, you will modify a conjugate gradient
	algorithm to allow multiple realizations, and then test it on the 1-D inverse
  linear interpolation problem given in Chapter 3 of the book.  In the second
  part, you will use your new solver to tackle the ``Wells not matching the
  seismic map'' problem of Chapter 3.  You will improve the ``extension'' of
	the seismic data into empty regions of the model space, first by writing a 
	program to de-trend the seismic data with a least squares best-fitting plane,  
	and second, by using your new solver to come up with multiple realizations of 
	the extended seismic data.
\end{abstract} 

\section{Introduction}
	In this course, we rely heavily on least squares, and sometimes too blindly
	believe that least squares provides {\em the} answer.
	Given data, a model, and a linear operator connecting the two -- $\bf Lm = d$
	-- the least squares solution to the problem takes the familiar form:
	\[ \bold m_{LS} = \left(\bold L^T \bold L \right)^{-1} \bold L^T \bf d.\]
%	However, many assumptions underlie this ``optimal'' result.  First, the
%	operator $\bf L$ must be overdetermined.  Second, and more importantly,
%	the result is unique only if $\bold L$ has no nullspace.  

	Assuming that we have computed an orthonormal basis for the nullspace of 
	$\bf L$, we can write the least squares solution in a more general form:
	\[ \bold m_{LS} = \left(\bold L^T \bold L \right)^{-1} \bold L^T \bold d
	                + \sum_{i=1}^p a_i \bold m^n_i.\]
	The $\bold m^n_i$ are the basis vectors for the nullspace of $\bf L$, which we 
	assume has dimension $p$.  The $a_i$ are arbitrary scalar coefficients.

	In other words, the space of all possible models consists of the least
	squares solution plus an arbitrary linear combination of the nullspace
	vectors. 
	
	In the first chapter of the book, you were introduced to simple binning
	operators.  It is simplest to start with the adjoint of binning: it maps
	a vector of irregularly-spaced data (with fixed coordinates) to a regular
	grid.  We can easily conceptualize the nullspace of the binning operator:
	any point in model space to which no data point maps (empty bins)
	is in the nullspace of the binning operator.  

	Here is the crux of the problem: we can set the model points in the 
	empty bins to any value we like, and still fit the data equally well.

	We use regularization to constrain the nullspace to something ``reasonable''.
	For instance, using $\nabla$ as a regularization forces the empty bins
	to obey Poisson's equation ($\nabla^T \nabla m = 0$).  As you may recall from
	college physics, solutions to Poisson's equation have minimum energy.
	In one sense, the minimum energy solution is the safest solution.  If a 
	given empty bin is far from any data, the safe estimate of the value in the
	empty bin is zero, assuming the data are zero mean.

	Unfortunately, after solving a missing data problem, the result often doesn't 
	look realistic.  Acquisition footprint still betrays the known data locations.
	
	For some time, Geostatistics has provided an interesting alternative to the
	over-smooth model maps: multiple realizations.  Instead of an ``answer'', the
	geostatistician returns a suite of equiprobable (i.e., which fit the data 
	equally well) answers, all of which have a realistic texture in the empty 
	regions.  The realistic texture was imposed by computing the variance of the
	known data and drawing the empty bin values from that distribution.  In fluid
	flow problems, the models with texture provided more accurate results than
	the smooth, minimum energy models described above.

	In short, geostatistics provides models where the empty bins have the same
	covariance {\em and} variance as the known data.  Model covariance is 
	essentially the way in which neighboring model points depend on one another.
	We specify the model covariance through the regularization operator.
	By minimizing $\| \nabla^T \nabla \bold m \|^2$, we impose the assumption 
	that neighboring model points are similar in magnitude, leading to a smooth
	model.  Model variance is more difficult to conceptualize, because we must
	treat the model points as random variables.  The minimum energy model is
	essentially the average of many different realizations.

	\cite{Clapp.sep.105.bob3} showed how our existing estimation methodology
	can be extended to create multiple realizations.  Read his paper.  Morgan Brown has 
	provided a slightly different explanation to achieve the same goal.  
	Start with the fitting 
	goals for the regularized missing data problem:
	\begin{eqnarray}
		\bf Lm &\approx& \bf d \\
		\bf \epsilon Am &\approx& \bf 0 
	\end{eqnarray}
	Recall from above that we can replace $\bf m$ with 
	$\bf m + \sum_{i=1}^p a_i \bold m^n_i.$ 
	\begin{eqnarray}
		\bf Lm &\approx& \bf d \\
 		{\bf \epsilon Am + \epsilon A} \sum_{i=1}^p a_i \bold m^n_i &\approx& \bf 0
	\end{eqnarray}
	The nullspace vectors are annihilated by the modeling operator, $\bf L$.
	Now expand the fitting goals into the objective function that we minimize
	with least squares.
	\begin{equation}	
		Q( \bold m) = \| {\bf Lm - d} \|^2
		            + \epsilon^2 \| {\bf Am} \|^2
		            + \epsilon^2 \left|\left| \bold A \sum_{i=1}^p a_i \bold m^n_i \right|\right|^2 \label{eq:objfunc}
	\end{equation}
	The third term of equation (\ref{eq:objfunc}) is confusing, but recalling some
	basic least squares theory leads to simplification.  After minimization,
	both model and data residuals should be white.  Rather than worrying about 
	choosing the $a_i$ and $\bold m^n_i$ explicitly, we can simply replace the
	third term with a vector of random numbers multiplied by a scalar, $\sigma \bf v$.
	\begin{equation}	
		Q( \bold m) = \| {\bf Lm - d} \|^2
		            + \epsilon^2 \| {\bf Am} \|^2
		            + \sigma^2 \| \bold v \|^2 \label{eq:objfunc2}
	\end{equation}	
	This is Clapp's \cite{Clapp.sep.105.bob3} strategy.  To implement it,
	we only need to initialize the model residual vector (in our conjugate gradient
	solver) with $\sigma \bf v$.
	Running the optimization with many different choices of $\bf v$
	gives a suite of equiprobable realizations -- exactly what we were looking
	to achieve.

\section{Your Assignment}

	Visit the GEE class webpage ({\tt http://sep/doku.php?id=sep:courses:gp211}),
	find and download the source code for this lab, {\tt Lab4.tar}
	If you have not already, log into {\tt glad} using your class account.
	Save this file in your home directory, then type {\tt tar -xf Lab4.tar}
	to create your {\tt Lab4} directory.

	Begin the exercise by editing the {\tt paper.tex} file
	to put your name on the paper and to answer the questions below.
	If you aren't familiar with {\texttt emacs} or {\texttt vi}, there is
	an intuitive editor called {\texttt gedit} which you can use by typing
	{\tt gedit paper.tex}.

	 To see the results of your editing, type {\tt SCons}.
	 Hopefully the document will rebuild again with your changes. View it with your favourite pdf viewer.
	 If you really mess things up and the document does not rebuild, you can
	 always restore the original document by removing your {\tt Lab4} directory
	 and then type {\tt tar -xf Lab4.tar}.

	This lab has two parts.  In the first, you will use module {\tt solver\_reg\_mod}
	to perform multiple realizations, and then test it on the 1-D inverse
	linear interpolation problem given in Chapter 3 of the book.  In the second
	part, you will use the same solver to tackle the ``Wells not matching the
	seismic map'' problem of Chapter 3.  First you will attempt to improve the
	first step of the process, which is to extend the seismic map into all
	regions of the model plane, by writing a program to de-trend the seismic
	data with a least squares best-fitting plane, and then use the regularized
	solver to come up with multiple realizations of the extended seismic data.

	\subsection{Inverse Linear Interpolation}
	In this section, we use the simple 1-D example in section {\em Inverse
	Linear Interpolation}.  The data is an irregularly sampled sine wave.  
	The model is an estimate of the function on a regular, and much finer,
	grid.  The problem is highly underdetermined, so we regularize with the
	gradient.  In fitting goals, we have:
	\begin{eqnarray*}
		\bf Lm &\approx& \bf d \\
		\bf \nabla m &\approx& \bf 0
	\end{eqnarray*}
	$\bf L$ is linear interpolation; a mapping of the data between the 
	irregular and regular grids.

	\begin{enumerate}
		\item Type {\tt make multreal.view}.  The frames of the movie show three
			things. In yellow, we have the known data points. In blue is the ``normal'' 
			inversion result; the one we get without adding random numbers to the
			model residual.  In green are eight successive equiprobable realizations.
			At this point, the realizations all look like the normal result.
		\item Open the file {\tt invint.f90} for editing. 
			Modify the call to {\tt solver\_reg} to obtain multiple equiprobable
			realizations. Hint: mind the parameter {\tt xnull}.
		\item Once you are finished, type {\tt make multreal.view} to view your
			result.  Figure \ref{fig:multreal} displays the results in a different
			format.
		\item Do the realizations look ``realistic''?  In other words, if you 
			didn't know where the known data was, could you be fooled into thinking
			that some of the regions had real data? \\
			\boxit{YOUR ANSWER:}
	\end{enumerate}

	\plot{multreal}{width=6.0in,height=4.0in}{ Eight equiprobable realizations
		which fit the known data equally well. }

	\subsection{ Wells not matching the seismic map }

	This section follows directly from the section in Chapter 3 of the same
	name.  Review that section to refresh your memory on the main ideas.

	The most important (and problematic) step in the processing flow in the
	book was the infill and extension of the irregularly sampled seismic
	data to the model grid.  In this section of the lab, you will estimate 
	multiple realizations for the extended seismic data using what you have
	learned in the previous section.

	\begin{enumerate}

	\item 
		To see what the gridded data looks like, type { \tt make wellseisfold.view } on 
		the command line (also in Figure \ref{fig:wellseisfold} of this paper).  The
		left hand panels show the gridded seismic and well data, while the right hand 
		panels shows the {\em fold}, or the number of data points which fell into
		each grid cell. 

	\plot{wellseisfold}{width=6.0in,height=4.5in}{Top row, left to right: Binned well data;
		Fold of well data.  Bottom row, left to right: Binned seismic data; Fold of seismic
		data. }

	\item
		Now that we have the seismic data on a regular grid, we need to extend it to fill
		the entire grid.  Recall that we need to know the {\em difference} between the 
		binned seismic data and the binned well data at each well.  Unfortunately, after
		binning the seismic data, we have empty bins, or places with zero fold.  If we have 
		no data, we can't very well compute the difference!

		Recall that the book advocated first imposing a complicated boundary condition 
		before estimating the missing seismic data.  Clearly, this complication is 
		undesirable.  However, there was little way around this problem, since the solution 
		to missing data problems regularized with $\nabla$ or $\nabla^2$ tend to zero, and 
		the data is not zero mean. 

		If we could somehow estimate the ``trend'' of the seismic data, and then subtract 
		the trend from the seismic data, we could solve the missing data problem with
		no boundary condition, because the de-trended data would then be zero mean.  

		As a zeroth order solution, we could simply subtract the mean of the seismic data.
		In this section of the lab, you will go a step beyond the trivial constant-valued
		trend and write a least squares program to compute the plane that best fits the 
		seismic data.

	\item
		The least squares problem we're solving to estimate the plane coefficients is this 
		one:
		\begin{equation}	
			\left[ \begin{array}{ccc}
			        x_1 & y_1 & 1 \\
			        x_2 & y_2 & 1 \\
			        x_3 & y_3 & 1 \\
			            & \vdots& \\
			        x_n & y_n & 1 
			       \end{array}
			\right]
			\left[ \begin{array}{c}
			        A \\
			        B \\
			        C
			       \end{array}
			\right] \approx
			\left[ \begin{array}{c}
			        z_1 \\
			        z_2 \\
			        z_3 \\ 
			        \vdots \\ 
			        z_n \\ 
			       \end{array}
			\right] \label{eq:planefit}
		\end{equation}
		The goal is to minimize the sum of squared differences between the seismic data 
		($z_i$) and the plane $z = Ax + By + C$ at all $x_i$ and $y_i$.  As expected, 
		you will need to write the operator suggested by equation (\ref{eq:planefit}),
		and the adjoint.  To do this, open the module {\tt planefit\_mod.f90} in a file
		editor.  You need only write the operator and adjoint where prompted.  

	\item
		Next, in a file editor, open the program {\tt FitPlane.f90}, which drives the
		plane estimation process.  Read over the documentation to get an idea of what 
		the program is doing.  You are prompted to do two things.  First, you must
		make a call to {\tt solver\_smp()} to estimate the coefficients of the plane.
		For this you need to code the operator {\tt planefit\_lop}.
		Second, you must use the plane coefficients to compute the value of the plane
		on a regular grid, which is then output.
		
	\item
		When you are finished, type {\tt make seisplane.H}.  The {\tt FitPlane.x} 
		program runs the dot product test and displays the plane coefficients to the
		screen.  Make sure your operator passes the dot product test before going on!

	\item
		Next, put the correct solver call in {\tt lapfill\_mod.f90}.  Once you're done, 
		type {\tt make seisfill.view}.  You'll see the
		original binned seismic data, the plane you've estimated, and the results
		of solving a missing data problem to ``extend'' the seismic data - where
		the data is de-trended with the mean, and de-trended by your estimated
		seismic plane.  (Note: To solve the missing data problem, we use a version
		of the {\tt lapfill} routine given in the book, only modified to use
		{\tt solver\_reg\_mod}. )

		Moreover, you'll see that the movie has five frames.  The first frame
		does not add random numbers to the model residual.  The next three frames
		are equiprobable realizations -- computed using the {\tt solver\_reg\_mod}
		routine.  The fifth frame is the average of six separate realizations.

		The same result is shown in your paper copy as Figure \ref{fig:seisfill}.

		\plot{seisfill}{width=6.5in,height=8.5in}{ Comparison of subtracting
			mean and subtracting plane for seismic de-trending and extending 
			into empty regions.  Also, 3 equiprobable realizations for each
			case. }

	\item
		Do the realizations have the same texture as the known data in regions 
		of missing data?  If not, adjust the {\tt scale} parameter in file
		{\tt seisfill.P} until it does.  {\tt scale} corresponds to $\sigma$ in
		equation (\ref{eq:objfunc2}). \\
		\boxit{YOUR ANSWER:}
	\item
		Clearly propose a viable strategy for estimating $\sigma$ without using trial
		and error.  Get creative! If your strategy is
		viable, implement it and test it on this problem. \\
		\boxit{YOUR ANSWER:}
	\item
		We are now able to compute multiple, equiprobable estimates of the 
		extended seismic data.  As you recall from the book, the final step 
		of the calibration problem is to tie the extended seismic data to
		the well data by computing a smooth, additive correction function 
		(depth).  We estimate this correction function by solving a second 
		least squares problem.  The data is the depth error (extended seismic
		minus well depth) at the well locations. \\

		Clearly explain how you might use an ensemble of equiprobable extended
		seismic estimates to improve the estimation of the correction function. \\
		\boxit{YOUR ANSWER:}
	\item
		Does your seismic plane seem to capture the general trend that you 
		observe in the known seismic data?  Note: All figures in the plot
		are clipped to the same level.  In your opinion, does de-trending with 
		the seismic plane lead to a better extended seismic panel than 
		de-trending with the mean? \\
		\boxit{YOUR ANSWER:}
	\item
		We could easily fit the known seismic data with a more complicated 
		conic surface, i.e., $z = Ax^2 + Bx + Cy^2 + Dy + E$.  While we could 
		capture more complex trends by adding more parameters, such a
		strategy is dangerous.  Why?  Explain why simply de-trending with
		the mean might even be a safer approach than de-trending with the
		seismic plane. \\
		\boxit{YOUR ANSWER:}
	\item
		Examine equation (\ref{eq:planefit}).  Recall that in least squares,
		if we seek to ``invert'' $\bf Lm=d$, we really are computing the 
		inverse of $\bold L^T \bf L$.  Under what circumstances will 
		$\bold L^T \bf L$ be poorly conditioned?  What effects do you expect
		that poor conditioning would cause?  How could you overcome poor
		conditioning? \\
		\boxit{YOUR ANSWER:}

        {\bf EXTRA CREDIT }:
	\item Who was Hongwu Emperor? What's his given name and temple name? Hong long does Ming Dynasty last?
		\boxit{YOUR ANSWER:}
	\item With a green pen, correct all grammar and spelling mistakes
         on this Lab. Also mark any out of date instructions or otherwise
         confusing elements in the lab.

	\item Comment on the 6 realizations and their average. Would it be sufficient to describe the solution of the inverse problem as a mean and a standard deviation ?
		\boxit{YOUR ANSWER:}

	\end{enumerate}

\bibliographystyle{seg}
\bibliography{lab4}

\section{DONE}
When you are finished modifying the latex file. Compile the latex file into a pdf using the {\tt scons} command. Print your paper. Hand your copy in to your TA. Clean up your directory by typing \texttt{make clean}, note that \texttt{scons} should still compile your paper.


