\begin{center}
{\bf Due Date: 17:00, Wednesday, October 13, 2010} \\
{\bf TA: Xukai Shen (xukai@sep.stanford.edu)}
\end{center}
\title{Lab 3 - Model fitting by least squares}

\def\figdir{./Fig}
\author{Shunzhi}
\email{shunzhi@Qing}
\righthead{Iterative Methods}
\lefthead{GEE - Lab 3}


\begin{abstract}
In this computer exercise you are given a module which performs a
  simple velocity transform.  First you will improve the linear operator
	by implementing linear interpolation instead of nearest-neighbor 
	interpolation.  In the second part of the assignment, you
  will apply an iterative least-square optimization for velocity
  transform inversion.
\end{abstract}

\section{Introduction}
Chapter 4 of Basic Earth Imaging 
\footnote{http://sepwww.stanford.edu/sep/prof/bei/toc\_html/}
introduces the concept of velocity spectra. Performing NMO and stack
with different stacking velocities transforms a CMP gather into velocity space. Primary and multiple reflection events usually have
different apparent velocities. In this sense, they can be separated 
from the primaries in the velocity space. If we could forward model 
the CMP gather from its velocity-space representation, it would allow 
us to model multiples only and to subtract them from the original gather.

\plot{velvel}{width=6.0in,height=4.0in}{        
  Iteration between spaces.
  Left are model spaces.
  Right are data spaces.
  Right derived from left.
  Lower model space derived from upper data space.
}

Figure \ref{fig:velvel}, borrowed (with some modifications) from BEI,
illustrates this idea. An idealized velocity spectrum (the top left plot)
consists of four high-velocity spikes and a low-velocity spike. In
this figure, the horizontal axis in the velocity transform domain is
slowness, and the vertical axis is time. The top right plot shows a
model CMP gather, derived from that model. The bottom left figure is
the result of a velocity transform, and the bottom right figure is a CMP
gather derived from it. 

Both the forward modeling and the velocity transform were performed with
module \texttt{velsimp}. To better understand
this program, recall the following pseudo-code
\par\noindent
\vbox{\begin{tabbing}
indent \= char \= char \= char \kill
\>do $v$ \{\\
\>do $\tau$ \{\\
\>do $x$  \{\\
\>      \>$t = \sqrt{ \tau^2 + x^2/v^2 }$ \\
\>      \>if hyperbola superposition \\
\>      \>      \> data$(t,x)$ = data$(t,x)$ + vspace$(\tau,v)$ \\
\>      \>else if velocity analysis \\
\>      \>      \> vspace$(\tau,v)$ = vspace$(\tau,v)$ + data$(t,x)$\\
\>      \>      \}\}\}
\end{tabbing}}

We can see from Figure \ref{fig:velvel} that the adjoint process does
take us back to the velocity space, but the result is far from
perfect. The underlying velocity-space model is not reconstructed, and
we are unable to predict the CMP gather from it. In this assignment,
you will improve this situation first by enhancing the velocity transform 
operator, and then by applying an iterative least-square optimization.

\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 Lab3.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 Lab3.tar}
to create your {\tt Lab3} directory.

In your { \tt Lab3 } directory you will find three main programs
({\tt Veltran.f90, Velmod.f90,} and {\tt DPtest.f90}), and two subroutines
({\tt velsimp.lop} and {\tt halfdif.lop}).  You will also find a Makefile, a 
SConstruct file ({\tt SConstruct}) and a \LaTeX file ({\tt paper.tex})

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 Lab3} directory
and then type {\tt tar -xf Lab3.tar}.

{\bf Answer the following questions. }
  \begin{itemize}

	\item Begin by typing {\tt make velvel.view}.  A figure with four panels
		should appear.  The upper-left-hand panel is the ``known'' slowness,
		spikes corresponding to five hyperbolic events, four of which have
		the same slowness.  For the remainder of the assignment, we will refer
		to the velocity transform operator as $\bf H$.  In other words,
		\[ \bf Hm = d \]
		creates $(t,x)$ data, $\bf d$, by applying the operator $\bf H$ to
		a model, $\bf m$, in the space of zero-offset traveltime and slowness
		$(\tau,s)$. 
		$\bf H$ ``pushes'', or ``spreads'' each point in model space $(\tau,s)$	
		to the corresponding hyperbola (apex at $\tau$) in data space $(t,x)$. 
		Conversely, the adjoint operation,
		\[ \bold H^T \bf d = m, \]
		maps data, hopefully consisting of a series of events with hyperbolic
		moveout, to model space, which hopefully consists of a series of well
		resolved points, each corresponding to a single event in data space. \\

		In short, we hope the operator $\bf H$ is unitary, or orthogonal:
		\[ \bold H^T \bf H = I, \]
		but as you will see, it is not.  This fact is the basis for the lab.

	\item Type {\tt make velvel.view} again.  First look at the upper-right-hand
		plot, or the data modeled from the initial model vector:
		\[ \bf d_{\mbox{mod}} = Hm_0.\]
		Does the data look as you would expect?  If not, explain the source of the
		trouble. \\
	\boxit{YOUR ANSWER.}
		
	\item Now look at the lower-left-hand plot, or the velocity scan derived
		by applying the adjoint of $\bf H$ to the modeled data:
		\[ \bold m = \bold H^T \bold d_{\mbox{mod}} = \bold H^T \bf Hm_0.\]
		Compare this panel to the starting velocity model.  Obviously, there
		are considerable differences.  Notice that the original spikes in
		velocity turn into laterally-smeared spikes.  What causes this 
		smearing? \\
	\boxit{YOUR ANSWER.}


	\item Focus still on the lower-left-hand plot.  Can you think of any 
		other artifacts that might appear?  Do you see them in the lower-left-hand
		panel?  {\bf Hint:} Type {\tt Grey < mod2.H | Tube } to view this panel
		in grayscale. \\
	\boxit{YOUR ANSWER.}


	\item Last, focus on the lower-right-hand plot, which is the ($t,x$) data
		modeled from the ugly velocity panel in the lower left.  These data look
		quite different from the ``theoretical'' data, or the data in the upper-right-hand
		corner, which was modeled directly from spikes in velocity space.
		What do you think are the two most important differences between this
		data and the theoretical data?  Explain the cause of both differences. \\
	\boxit{YOUR ANSWER.}


        \item The modeled CMP gather in the top right plot of Figure \ref{fig:velvel} 
		appears to have a step-like moveout. This is clearly an artifact of the 
		nearest-neighbor interpolation, used in module \texttt{velsimp}. 
		Put simply, the operator loops over slowness, 
		offset, and zero-offset traveltime.  Nonzero-offset traveltimes are computed 
		via the hyperbola equation: 
		\[t = \sqrt{ \tau^2 + (xs)^2. } \]
		Rarely does this value fall on a traveltime ``bin''. \\
	
		Your next task is to modify the \texttt{velsimp}  module so that linear 
		interpolation between two adjacent samples in time is used instead of 
		nearest-neighbor interpolation.  When you are done, Figure \ref{fig:velvel-lint}
		will reflect your changes.  \\

		{\bf Hints:} 
		Notice that the {\tt velsimp\_init} interface has a logical argument,
		{\tt lint}.  If {\tt lint = .false.}, nearest-neighbor interpolation
		is used.  If {\tt lint = .true.}, linear interpolation is used. 
		For an example of the algorithm, review operator \texttt{lint1} from 
		Chapter 1.

	\item Once you are finished, check to see if you have coded the adjoint
		correctly by typing {\tt make dottest}.  Modify program \texttt{DPtest.f90} 
		so that \texttt{velsimp} is called using linear 
		interpolation.  You need not report your results, but if you have coded
		the adjoint incorrectly, your later results will have problems.
		{\bf Hint: } To get an intuition for calling the {\tt dot\_test}
		function, look at the source code: {\tt \~prof/gee/Lib/dottest.f90}.
	\boxit{YOUR ANSWER.}

	\item Now type {\tt make velvel-lint.view} to view your results.  A 
		two-frame movie will come onto the screen.  The first frame shows
		your nearest-neighbor results, while the second frame shows your
		improved, linear interpolation results.  Compare the CMP data 
		modeled directly from the slowness spikes.  Has linear interpolation
		improved the hyperbola modeling?  \\
	\boxit{YOUR ANSWER.}


\plot{velvel-lint}{width=6.0in,height=4.0in}{        
	Same as Figure \ref{fig:velvel}, only we use linear interpolation in the
	velocity transform operator.
  Iteration between spaces.
  Left are model spaces.
  Right are data spaces.
  Right derived from left.
  Lower model space derived from upper data space.
}

\item At this point, you should be wondering why the addition of linear
	interpolation failed to make better-looking lower panels on Figure
	\ref{fig:velvel-lint}.  Certainly, smoother hyperbolas were generated,
	but say you were an interpreter, picking velocities from real data.
	The upper-right-hand panel of Figure \ref{fig:velvel-lint} is made up
	of perfect hyperbolas, yet the velocity scan is nearly useless for 
	the real goal, which is picking velocities! \\

	Least squares estimation (inversion) comes to the rescue.  Rather than
	simply applying the adjoint and stopping, we will search for the model
	(velocity panel) which optimally explains the data (CMP gather), given 
	a physical connection ($\bf H$) between the data and model.

\item Modify program \texttt{Veltran} so that
  in addition to applying the adjoint operator, it is enabled to call 
	subroutine \texttt{solver} to run an iterative least-square optimization 
	by the conjugate-gradient method.  You can find examples of using the
  \texttt{solver\_smp} subroutine in Chapter 2 of the textbook. 

\item After you have modified the program, run {\tt make velinv-lint.view}
  to recreate Figure \ref{fig:velinv-lint}.  You will see a two frame movie.
	The first frame is simply Figure \ref{fig:velvel-lint}, while the second
	is the result using optimization.
  
  \plot{velinv-lint}{width=6.0in,height=4.0in}{ Iterative approach
		to solving problem, using linear interpolation along the time axis.  
		Left are model spaces.  Right are data spaces.  Right derived from 
		left.  Lower model space derived from upper data space after 20 
		conjugate-gradient iterations.} 

	\item Has optimization improved the quality of the velocity scan
		(lower-left-hand panel)?  Do you see model space artifacts in the
		inversion result?  Explain differences with the adjoint processing
		result. Explain why artifacts are still or are no longer present. \\
	\boxit{YOUR ANSWER.}

		
	\item In practice, the data modeled from the inversion result provides
		us with an excellent diagnostic tool.  In short, if application of the 
		modeling operator to the model we've just estimated doesn't give 
		something that looks like the input data, something has gone wrong.
		Type {\tt make diff.view} and examine the top row of the plot. From
		left to right, we have the input CMP data, the modeled CMP data, and
		the difference.  Describe the difference mathematically.  \\
	\boxit{YOUR ANSWER.}


		\plot{diff}{width=6.0in,height=4.0in}{  Top, from left to right:
			Input CMP gather;  CMP gather modeled from estimated velocity;
			Difference.  Bottom, from left to right: Input slowness model;
			Estimated slowness model; Difference.
		}
	
 	 \item Type {\tt make diff.view} and examine the bottom row of the
		plot.  From left to right, we have ``known'' slowness model, estimated
		slowness model, and difference.  From this information and your answer
		to the previous question, does the inverse problem appear to 
		have a null space?  If so, describe the components of the null space 
		that you directly observe. \\
	\boxit{YOUR ANSWER.}


	\item In the method of conjugate directions, is the residual power
   	guaranteed to decrease at each iteration? \\
	\boxit{YOUR ANSWER.}

		
   \item Is the gradient power guaranteed to decrease? \\
	\boxit{YOUR ANSWER.}

		
    \item How would you define the rate of convergence? \\
	\boxit{YOUR ANSWER.}

    \end{itemize}

  {\bf EXTRA CREDIT }:
  \begin{itemize}
	\item Who was Shunzhi. Describe his relationship with Qing Dynasty. Is he the first Emperor of Qing Dynasty, if not, who is(are) before him.	
	\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 Formulate at least one (research) question/idea on the topic of this lab. You are not expected to have the answer. This can be on the method or on the problem. If you can't think up a research question, forumate a question of something that is unclear to you about what you have learned so far in GEE.
	\boxit{YOUR ANSWER.}

   \end{itemize}

\section{DONE}
When you are all done, print your paper and hand in to your TA.
Clean up your directory by typing \texttt{make clean}. 



