annintro.tex
% TEX SET-UP
% TO BE PUBLISHED IN THE PROCEEDINGS OF WORKSHOP ON AUTOMATED
% DATA ANALYSIS IN ASTRONOMY
% Please do not change the page layout, except \hoffset and/or
% \voffset, where changes may be required depending on default printer
% settings. Leave a margin of 3.8 cm on the left and 5 cm at the top

\magnification=\magstep0
\hsize=13.5 cm               %  horizontal size of printed page
\vsize=19.0 cm               %  vertical size of printed page
\baselineskip=12 pt plus 1 pt minus 1 pt  % The line spacing
\parindent=0.5 cm  % The paragraph indent
\hoffset=1.3 cm      % The horizontal offset (may need to be changed)
\voffset=2.5 cm      % The vertical offset (may need to be changed)
%
\font\twelvebf=cmbx10 at 12truept % Set bold font for Title
\font\twelverm=cmr10 at 12truept % Set large font for Name
\overfullrule=0pt
%\nopagenumbers    %  Actual page nos. will be inserted by the Editor
\input psfig.tex
%
% The headlines
% The changes in the headlines should be made just before the Abstract
\newtoks\leftheadline \leftheadline={\hfill {\eightit Authors' name}
\hfill}
\newtoks\rightheadline \rightheadline={\hfill {\eightit the running title}
 \hfill}
% Do not change the headline on the first page of paper.
\newtoks\firstheadline \firstheadline={
%{\eightrm To appear in  {\eightit Automated Data Analysis in Astronomy}, 
%R.\ Gupta, H.P.\ Singh, C.A.L.\ Bailer-Jones (eds.),
%Narosa Publishing House, New Delhi, India, 2001} 
\hfill}
\def\makeheadline{\vbox to 0pt{\vskip -22.5pt
\line{\vbox to 8.5 pt{}\ifnum\pageno=1\the\firstheadline\else%
\ifodd\pageno\the\rightheadline\else%
\the\leftheadline\fi\fi}\vss}\nointerlineskip}
%
% Defining 8-pt fonts for figure captions and references
\font\eightrm=cmr8  \font\eighti=cmmi8  \font\eightsy=cmsy8
\font\eightbf=cmbx8 \font\eighttt=cmtt8 \font\eightit=cmti8
\font\eightsl=cmsl8
\font\sixrm=cmr6    \font\sixi=cmmi6    \font\sixsy=cmsy6
\font\sixbf=cmbx6
%
%for switching to eight point type \eightpoint
\def\eightpoint{\def\rm{\fam0\eightrm}
\textfont0=\eightrm \scriptfont0=\sixrm \scriptscriptfont0=\fiverm
\textfont1=\eighti  \scriptfont1=\sixi  \scriptscriptfont1=\fivei
\textfont2=\eightsy \scriptfont2=\sixsy \scriptscriptfont2=\fivesy
\textfont3=\tenex   \scriptfont3=\tenex \scriptscriptfont3=\tenex
\textfont\itfam=\eightit  \def\it{\fam\itfam\eightit}%
\textfont\slfam=\eightsl  \def\sl{\fam\slfam\eightsl}%
\textfont\ttfam=\eighttt  \def\tt{\fam\ttfam\eighttt}%
\textfont\bffam=\eightbf  \scriptfont\bffam=\sixbf
\scriptscriptfont\bffam=\fivebf \def\bf{\fam\bffam\eightbf}%
\normalbaselineskip=10pt plus 0.1 pt minus 0.1 pt
\normalbaselines
\abovedisplayskip=10pt plus 2.4pt minus 7pt
\belowdisplayskip=10pt plus 2.4pt minus 7pt
\belowdisplayshortskip=5.6pt plus 2.4pt minus 3.2pt \rm}
%
% define the displayed equations to be indented 1.5 cm from left
% as required by the Proceedings. Hopefully this will work for all
% equations. With this definition using the normal $$...$$ should
% produce the equations with correct indentation, but it will be necessary
% to use \eqno to put equation numbers (though it will be possible to put
% blank eq. nos.). Further, eq. nos using \eqalignno will not work
%
\def\leftdisplay#1\eqno#2$${\line{\indent\indent\indent%
$\displaystyle{#1}$\hfil #2}$$}
\everydisplay{\leftdisplay}
%
% Some useful definitions
% less than or order of \la
\def\frac#1#2{{#1\over#2}}
\def\la{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr<\cr%
\sim\cr}}}}}
% greater than or order of \ga
\def\ga{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr>\cr\sim%
\cr}}}}}
%
\input psfig.tex
%
%to generate boldface characters
\def\pmb#1{\setbox0=\hbox{$#1$}\kern-0.015em\copy0\kern-\wd0%
\kern0.03em\copy0\kern-\wd0\kern-0.015em\raise0.03em\box0}
%
% AUTHOR DEFINITIONS
\def\teff{T$_{\rm eff}$}
\def\logg{$\log g$}
\def\met{[M/H]}

\def\yv{{\bf y}}
\def\vv{{\bf v}}
\def\xv{{\bf x}}
\def\pv{{\bf p}}
\def\qv{{\bf q}}
\def\wv{{\bf w}}

\def\sig68{$\sigma_{68}$}
\def\sigrms{$\sigma_{RMS}$}
\def\sigint{$\sigma_{int}$}
%
%Beginning of Document%
%
%
\pageno=1
%
\vglue -28 pt
%
\leftline{\eightrm To appear in  {\eightit Automated Data Analysis in Astronomy},
	  R.\ Gupta, H.P.\ Singh, C.A.L.\ Bailer-Jones (eds.),} 
\leftline{\eightrm Narosa Publishing House, New Delhi, India, 2001}
%
\vglue 68 pt  %Leave some space on page 1 before the title
% The title
%
\leftline{\twelvebf  An introduction to artificial neural networks}
% if more than one line is required for the title, then use next two lines ...
%
\smallskip
% end of title
\vskip 46 pt  % Space between title and author(s) name(s).
%
\leftline{\twelverm Coryn A.L.\ Bailer-Jones} % Name of Authors
\vskip 4 pt
\leftline{\eightit Max-Planck-Institut f\"ur Astronomie, K\"onigstuhl 17, 69117 Heidelberg, Germany
}
\leftline{\eightit email: calj@mpia-hd.mpg.de}
%\leftline{\eightit  to reduce the number of lines}
%
% If authors are from different institutes, repeat the above lines
% for each institution. For authors from same institution write the
% names in one line.
%
\vskip 0.5 cm
\leftline{\twelverm Ranjan Gupta} % Name of Authors
\vskip 4 pt
\leftline{\eightit IUCAA, Post Bag 4, Ganeshkhind, Pune-411007, India}
\leftline{\eightit email: rag@iucaa.ernet.in}
\vskip 0.5 cm
\leftline{\twelverm Harinder P.\ Singh}
\vskip 4 pt
\leftline{\eightit Department of Physics, Sri Venkateswara College, Benito Juarez Road, Dhaula Kuan}
\leftline{\eightit New Delhi-110021, India}
\leftline{\eightit email: hps@ttdsvc.ernet.in}

\vskip 20 pt % leave some space between author(s) names(s) and abstract
%
%
% The leftheadline should include the Authors' name, for two authors use
% \&  (e.g. I. M. Author \& I. M. Co-author) for three or more authors
% use et al.,
\leftheadline={\hfill {\eightit Bailer-Jones, Gupta \& Singh} \hfill}
% Use a short running title as the rightheadline
\rightheadline={\hfill {\eightit An introduction to artificial neural networks
}  \hfill}

% Abstract begins
%
{\parindent=0cm\leftskip=1.5 cm

{\bf Abstract.}
\noindent
Artificial neural networks are algorithms which have been developed to
tackle a range of computational problems. These range from modelling
brain function to making predictions of time-dependent phenomena to
solving hard (NP-complete) problems.  In this introduction we describe
a single, yet very important, type of network known as a feedforward
network. This network is a mathematical model which can be trained to
learn an arbitrarily complex relationship between a data and a
parameter domain, so can be used to solve interpolation and
classification problems. We discuss the structure, training and
interpretation of these networks, and their implementation, taking the
classification of stellar spectra as an example.
\smallskip
%
\vskip 0.5 cm  %  Space between Abstract and Key words
%
{\it Key words:} neural networks, nonlinear models, classification,
interpolation, data modelling

}                                 %  End of abstract

% Beginning of document
%
%
% Beginning of a section heading
%
% for the first section leave 20 pt space, for subsequent sections just
% leave bigskip (i.e. 12 pt)
\vskip 20 pt
\centerline{\bf 1.\ Biological neural networks}
\bigskip
\noindent
Artificial neural networks were originally introduced as very
simplified models of brain function, so it is initially instructive to
consider the broad analogies between artificial neural networks and
biological ones. Of course, the latter are considerably more complex,
and many artificial neural network models today bear very little
resemblance to biological ones. For a history of the development of
artificial neural networks, see, for example, the introduction of
Hertz, Krogh and Palmer (1991).

The human brain consists of billions of interconnected neurons. These
are cells which have specialized membranes which allow the
transmission of signals to neighbouring neurons. Fig.~1 shows the
structure of a neuron.  A single {\it axon} extends from the basic
cell body.  This is the output from the neuron, and typically divides
into many sub-branches before terminating at a {\it
synapse}. Electrical pulses are transmitted along the axon to the
synapses by the transfer of Na$^{+}$ ions. The arrival of a voltage
pulse stimulates the release of neurotransmitting chemicals across the
synaptic cleft towards the postsynaptic cell, which is the receiving
part of the next neuron.  This postsynaptic cell passes the signal via
the {\it dendrite} to the main part of the neuron body.  The inputs
from the different dendrites are then combined to produce an output
signal which is passed along the axon, and so the process
continues. However, a signal is only produced in the axon if there are
enough inputs of sufficient strength to overcome some threshold value,
and then the output is some nonlinear function of the input stimuli.

\noindent
\midinsert
\vskip 8.3cm
\special{psfile=neuron.ps voffset=270mm hoffset=-5mm
vscale=50mm hscale=50mm angle=270}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 1.} Structure of a single neuron in the brain.
Information is passed into the neuron from other neurons via the
dendrites to the left. Given appropriate input stimuli, the cell sends
an output signal along the axon to the synapses, where the signal is
transmitted to other neurons.  }
\endinsert

The human brain consists of about $10^{11}$ neurons, and each neuron
has between a few and a few thousand synapses on its dendrites, giving
a total of about $10^{14}$ synapses (connections) in the brain. The
`strength' of the synaptic connection between neurons can be
chemically altered by the brain in response to favourable and
unfavourable stimuli, in such a way as to adapt the organism to
function optimally within its environment. The synapses are therefore
believed to be the key to learning in biological systems.



\vskip 12 pt
\centerline{\bf 2.\ Artificial neural networks}
\bigskip
\noindent
``Artificial neural networks'' is a relatively loose term referring to
mathematical models which have some kind of distributed architecture,
that is, consist of processing nodes (analogous to neurons) with
multiple connections (analogous to dendrites and axons). These
connections generally have adaptable parameters which modify the
signals which pass along them. There are numerous types of artificial
neural networks for addressing many different types of problems, such
as modelling memory, performing pattern recognition, and predicting
the evolution of dynamical systems. Most networks therefore perform
some kind of data modelling, and they may be split into two broad
classes: {\it supervised} and {\it unsupervised}. The former refers to
networks which attempt to learn the relationship between a data and a
parameter domain, while the latter refers to networks used to find
``natural'' groupings with a data set independently of external
constraints. What they have in common is the idea of learning about a
problem through relationships intrinsically present in data, rather
that through a set of predetermined rules.  An introduction to
several types of neural networks (but by no means all) is given by
Hertz, Krogh and Palmer (1991). A less mathematical overview is
provided by Beale \& Jackson (1990).

In this article we introduce what is one of the most important types
of supervised neural networks, called a {\it feedforward multilayer
perceptron}. The term {\it perceptron} is historical, and refers to
the function performed by the nodes. {\it Feedforward} means that
there is a definite input and output, and a flow of data in one
direction. This is in contrast to {\it recurrent neural networks} in
which data flows in a loop: these are important when time plays a
relevant role in the problem (see, e.g., Bailer-Jones 1998a). A
comprehensive treatment of feedforward networks is provided by Bishop
(1996), where more details on many of the themes discussed below can
be found.

\noindent
\midinsert
\vskip 9.5cm
\special{psfile=ann_2layer_2.eps voffset=15mm hoffset=0mm vscale=43mm hscale=43mm angle=0}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 2.} Feedforward artificial neural network architecture.
This particular example shows two hidden layers. Not explicitly
shown are the bias nodes: each layer has an extra node that holds a
constant value (and has no inputs) and provides an offset to the
subsequent layers.  These are necessary for the network to model
functions properly.}
\endinsert

Fig.~2 shows a feedforward neural network with several inputs, two
hidden layers, and several outputs.  Each node in the input layer
holds a value, $x_i$. In our example application, the input vector,
$(x_1,x_2,\ldots,x_i,\ldots)$, is a stellar spectrum, and the output
vector, $(y_1,y_2,\ldots,y_l,\ldots)$, is a set of stellar parameters,
e.g.\ \teff, \logg\ and
\met. Each of the input nodes connects to every node in the next layer
of nodes, the first `hidden' layer, and each of these connections has
a weight, $w_{i,j}$, associated with it.  A node in the hidden layer
forms a weighted sum of its inputs, and passes this through a {\it
nonlinear transfer function}, such that the output from the $j^{th}$
hidden node is

$$ p_j = \tanh \left( \sum_i w_{i,j} x_i \right) \ \ . \eqno(1) $$

\noindent
These values are passed to a second hidden layer which performs
a similar processing, the output from that layer being the vector $\qv$

$$ q_k = \tanh \left( \sum_j w_{j,k} p_j \right) \ \ . \eqno(2) $$

\noindent
The output layer then performs a simple sum of its inputs, so that the
network output, $\yv$, is

$$ y_l =  \sum_k w_{k,l} q_k \ \ . \eqno(3) $$

\noindent
The tanh function in the hidden layers provides the nonlinear
capability of the network.  Other nonlinear functions are possible;
the sigmoidal function ($1/(1-\exp[-\sum wx])$) is often used. Both
functions map an infinite possible input range onto a finite output
range, $-1$ to $+1$ in the case of tanh. This imitates the transfer
function of neurons.

% force a page break
\vskip 30 pt

%\vskip 12 pt
\centerline{\bf 3.\ Network training}
\bigskip
\noindent
The weights, $\wv$, appearing in equations 1--3 are the free parameters of the
network. Clearly, in order for the network to yield appropriate
outputs for given inputs, the weights must be set to suitable
values. This is done by {\it training} the network on a set of input
vectors for which the ideal outputs (or {\it targets}) are already
known. This is {\it supervised learning}. The idea is that the network
will encapsulate the relationship between, in our example, a wide
range of types of stellar spectra and their associated physical
parameters. The network does not come up with a new classification
system, but rather will {\it interpolate} the training data to give a
{\it generalization} of the relationship between spectra
and their classifications. We can then use this trained network to
obtain reliable classifications for unclassified spectra.

Supervised learning proceeds by minimizing an error function with
respect to all of the network weights. The error function typically
used is the sum-of-squares error, which for a single input vector,
$n$, is

$$ e^{\{n\}} = \frac{1}{2}\sum_l \beta_l(y_l^{\{n\}} - T_l^{\{n\}})^2 \ \ , \eqno(4) $$

\noindent where $T_l$ is the target output value for the $l^{th}$ output node.
The $\beta_l$ terms allow us to assign different weights to different
outputs, and thereby give more priority to determining certain outputs
correctly.

The most popular method for training the network is with the {\it
backpropagation algorithm} (Rumelhart, Hinton and Williams 1986a,b),
in which we determine the gradient of the error function with respect
to each of the weights in the network.  The network is first
initialized by setting the weights to small random values. We then
pass one of the training vectors (spectra) through the network and
evaluate the outputs, $y_l$, as well as all the intermediate node
values in the network (the vectors {\bf p} and {\bf q}).
Differentiating equation 4 with respect to a weight, $w_{k,l}$, in the
final weights layer and substituting for $y_l$ from equation 3 gives

$$ \frac{\partial e}{\partial w_{k,l}} =
\beta_l (y_l - T_l) q_k \ \ , \eqno(5) $$

\noindent 
where we have used the fact that the weights are causally independent
of one another and of the node values. (We have dropped the $n$
superscript to make things less crowded.) We can determine the
gradient of $e$ with respect to a weight, $w_{j,k}$, in the previous
weights layer in the same way, but now $q_k$ is a function of
$w_{j,k}$,

$$ \frac{\partial e}{\partial w_{j,k}} =
 \sum_l \beta_l (y_l - T_l) \frac{\partial y_l}{\partial w_{j,k}} = 
 \frac{\partial q_k}{\partial w_{j,k}} \sum_l \beta_l (y_l - T_l) w_{k,l}
\ \ . \eqno(6) $$

\noindent 
The remaining partial derivative can then be evaluated from equation
2.  Thus we see that the error gradient with respect to any weight in the
network can be evaluated by {\it propagating} the error dependency back
through the network to that weight. This can be continued for any
number of layers, up to the input layer, giving us the complete error
gradient vector, $\partial e/ \partial {\bf w}$, where ${\bf w}$ is
the set of all network weights.
 
Having evaluated the error gradient vector, there are a number of ways in
which it can be used to train the network. The most common is the {\it
gradient descent} process, in which we adjust the weight
vector in the direction of the negative of the gradient vector, i.e.\

$$ \Delta {\bf w} = - \mu \frac{\partial e}{\partial {\bf w}} \ \ . \eqno(7) $$

\noindent The factor $\mu$ determines how large a step is made, and typically
has a value between 0.1 and 1.0. The gradient is then recalculated using the new
values of the weights, and over many iterations the weights should
move towards a value which represents a small value of the error, $e$.
For training to strictly converge (i.e.\ reach an exact minimum),
$\mu$ would have to decrease to zero during training, although in
practice a sufficiently small error can often be achieved without this.

Of course, we are generally interested in getting the network to
generalize its input--output mapping for a range of types of objects
(e.g.\ different classes of stars), so we must correspondingly use a
range of objects which are representative of the problem. We therefore
apply the above training algorithm successively to each vector in a
training sample. We may either update the weights after each vector is
passed, or we can keep the weights fixed and update only after the
errors have been calculated for all vectors in the training set. In
this latter case, known as {\it batch} training, we update using the
the average error,

$$ E = \frac{1}{N} \sum_{n=1}^{n=N} e^{\{n\}} \ \ , \eqno(8) $$

\noindent
and the corresponding average error gradient vector.

The network training is a nonlinear minimization process in $W$
dimensions, where $W$ is the number of weights in the network.  As $W$
is typically large, this can lead to various complications.  One of
the most important is the problem of local minima. Clearly, training
using the gradient descent technique will stop upon reaching the
bottom of a minimum, but this may only be a {\it local} minimum, not
the {\it global} minimum, and may correspond to a much larger
error. To help avoid local minima, a {\it momentum term} is sometimes
added to the weight update equation 7, in which a fraction, $\nu$, of
the previous weight update is added to the current one. This provides
some ability to ``escape'' from small local minima: the larger is
$\nu$ the larger are the minima which can be ignored, although
possible at the expense of slower convergence.  A brute force approach
to avoiding local minima is to retrain the network several times from
different initial random weights, and see whether each case converges
to the same solution. If the majority do, this could be taken as the
best solution. Alternatively, we could use all of these networks
together in a {\it committee}, in which the outputs from the networks
are averaged.

A common difficulty is knowing when to stop training, as true
convergence in all but very simple (or fortuitous) cases is rarely
achieved. Without searching the entire weights space, it is generally
impossible to know whether the minimum reached is the global one or
just a local one, and if local, how much better the global minimum
is. Furthermore, the minimum may have a long, almost horizontal bottom
(thinking analogously in three dimensions), so it will often be
necessary to terminate training when the gradient falls below some
small value. Setting the convergence gradient too small will slow
convergence; too large and the search may terminate prematurely.

Fig.~3 shows several examples of how the error in equation 8 changes
during training. The dotted line shows the error on the actual
training data. However, we really want to test whether the network has
generalized its mapping, rather than ``memorized'' the specific
training data.  The solid line therefore shows the error evaluated on
a separate test data set, which is not used to modify the
weights. Under certain circumstances it is possible for the network to
fit the noise in the training set, in which case it is said to have
{\it overtrained}. This can be seen in Fig.~3c, where the error on the
test data starts to rise at some point, although the error on the
training data continues to fall. Frequently, both errors continue to
drop monotonically until levelling out (Fig.~3a), which is
desirable. (A difference in the errors finally achieved may reflect
differences in the training and test data sets.)  Fig.~3b shows an
example of where the error on the test data rises briefly, and then
drops: we should not, therefore, cease training at the first sign of a
rising test error. In certain cases the error can be very noisy
(Fig.~3d), oscillating around a minimum value, probably due to $\mu$
in equation 7 being held constant. Anecdotal evidence suggests that
this is more likely to occur when multiple outputs are used, although
it depends on the minimization algorithm used.

\noindent
\midinsert
\vskip 8.5cm
\special{psfile=trainerr.ps voffset=-50mm hoffset=50mm vscale=40mm hscale=40mm angle=0}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 3.} Examples of the evolution of the network error during
training. The dotted (solid) line is the error evaluated on the
training (test) data (except in case D where the lower error is for
the training data).  (A),(B) and (C) are for networks with a single
continuous output. (D) is for a 50:5:5:3 network with probabilistic
outputs. All cases used gradient descent with constant $\mu$ and $\nu$
terms. See Bailer-Jones (1996) for more details.  }
\endinsert

Gradient descent only uses knowledge of the local gradient, which will
not, in general, point towards the global minimum, thus requiring many
steps to reach this minimum. Numerous other optimization techniques
exist to overcome this problem. One is the method of {\it conjugate
gradients}, in which a more efficient gradient is defined using the
second order derivatives (albeit implicitly) of the error with respect
to the weights. Several other methods make direct use of the second
derivatives (or {\it Hessian} matrix) to obtain more rapid
convergence.  We refer the interested reader to any number of books on
optimization techniques for further details.

The mapping produced by a neural network is an interpolation of the
training data. Whenever data are interpolated, there needs to be some
kind of complexity control to prevent overfitting the data. For
example, if we imagine fitting a curve through five noisy data points
in a two-dimensional space, then -- unless the data are co-linear -- a
fourth order polynomial model will give a better fit (a perfect one in
fact) than a straight line. This is unlikely to be appropriate,
though, as the data contain noise. Hence there is a trade off between
obtaining a good fit and not fitting the noise, and a compromise can
be achieved through the use of a {\it regularization} technique.  A
very basic approach is to split the training sample into two parts,
and during training on one half, monitor the total error on the other
half. As seen in Fig.~3, the test error {\it may} start to rise after
some number of iterations, indicating that the network has begun to
fit the noise (or {\it overfit} the data). Training could then be
stopped just before the error starts to rise. However, experience
shows that in some cases the training error never rises, even after a
very large number of iterations have occurred and a small error has
been achieved.

Another approach to regularization is {\it weight decay}. A high order
fit to data is characterized by large curvature of the mapping
function, which in turn corresponds to large weights. We can penalise
large weights by adding the following term to the training error in
equation 8

$$ \alpha \frac{1}{2} \sum_i w_i^2 \ \ , \eqno(9) $$

\noindent
where the sum is over all weights in the network and $\alpha$ is some
constant.  

Many neural network applications in astronomy have used gradient
descent for training, to the extent that some people think gradient
descent and neural networks to be synonymous.  This is not the case:
the feedforward neural network is a general means of representing a
nonlinear relationship between two domains, with the feature that the
error gradient is easily calculable in terms of the free parameters
(the weights) of the model. What optimization technique is then used
to determine the best solution for the weights is really a separate
issue. Gradient descent is often used because it is easy to code, but
it may not be the best choice.

\vskip 12 pt
\centerline{\bf 4.\ A Bayesian perspective}
\bigskip
\noindent
A rather different approach to interpolating data with a feedforward
network is provided within a Bayesian probabilistic framework.  If we
go back to first principles, then we see that the problem we are
trying to answer with the neural network is fundamentally
probabilistic: Given the set of training data, $D$, and some
background assumptions, $H$, what is the most probable output vector,
$\yv^*$, for some given input vector, $\xv$ (not in $D$)? When phrased
this way, we see that we are not really interested in the network
weights, and, moreover, an ``optimal'' set of network weights is not
actually required.  For example, if we have trained a network and
found a good mapping, we may retrain it and find a different set of
weights which give a mapping which is almost as good. Thus, while a given
input will yield much the same output as before, the weights may be
very different. In the optimization approach we ignore all but one
solution (except when using committees), whereas in a Bayesian
approach we combine all such solutions by the process of {\it
marginalization}. Specifically, if the output of the network with
weights $\wv$ and input $\xv$ is written $\yv(\xv, \wv)$, then

$$ \yv^* = \int \yv(\xv, \wv) P(\wv|D,H) d \wv  \ \ . \eqno(10) $$

\noindent
In other words, we weight each possible output by the probability,
$P(\wv|D,H)$, that the training data would produce those weights.
This {\it posterior} probability for the weights is given by Bayes'
Theorem

$$ P(\wv|D,H) \propto P(D|\wv,H) P(\wv|H) \ \ . \eqno(11) $$

\noindent
The first term on the right hand side of this equation is the {\it
likelihood}. It is the probability that a network with a given set of
weights produces the targets in the training data, D, when the inputs
from these data are fed in. If we assume a Gaussian error model for
the discrepancies between outputs and targets, this is just

$$ P(D|\wv,H) \propto \exp{-E} \ \ , \eqno(12) $$

\noindent
where $E$ is the sum of squares error from equations 8 and~4.  The
other term in equation 11 is the {\it prior probability} over the
weights, and represents our prior belief of what the weights should
be. In this Bayesian framework, the weight decay regularizer mentioned
earlier (equation 9) is naturally incorporated into this prior.

From equations 11 and 12 we can evaluate all of the terms under the
integral in equation 10. However, this integral must usually be
evaluated numerically with Monte Carlo methods. This involves
replacing the integral with a summation over $\yv(\xv, \wv)$ for a
large number of values of $\wv$ drawn from the distribution
$P(\wv|D,H)$. In this Bayesian approach, the optimization of the
weights is replaced with the integration over all possible weight
values, so no ``optimum'' weight vector is (or has to be)
calculated. While this approach is theoretically more satisfying than
an optimization, it is time consuming, as the marginalization must be
carried out for each new input vector we wish to classify. Partly for
this reason there exists an approximation to this approach which,
although still within the Bayesian framework, returns the problem to
one of optimization. There are several advantages to the Bayesian
approach, such as automatic complexity control (avoidance of
overfitting) and a natural way for determining the uncertainty in each
network output. More on the Bayesian approach for feedforward neural
networks can be found in (MacKay 1995).  For a general introduction to
Bayesian methods, Sivia (1996) is recommended.


\vskip 12 pt
\centerline{\bf 5.\ Implementation issues}
\bigskip
\noindent
Now that the basic principles of neural networks have been laid out,
it is necessary to ask how they are used in practice to solve a
problem.  We illustrate this with the problem of stellar spectral
classification.

\vskip 12 pt
\centerline{\bf 5.1.\ Training and testing data}
\bigskip
\noindent
The first thing which must be specified is the training data set,
consisting of spectra preclassified on the classification system of
interest. It must be emphasized that the network cannot be more
accurate than these preclassifications, so considerable effort must
often be spent in assembling the training data.  We must also ensure
that the examples in the training data set cover the full range of
types to which we will want to apply the trained network:  We cannot
expect a network trained only on B stars to classify M stars well, as
the characteristic features of M stars are not present in B star
spectra. Conversely, we do expect a network trained only on early and
late B stars to classify mid-B stars well, as B stars all have
qualitatively similar spectra.

Once a network has been trained, we want to evaluate its performance,
and this must be done using a separate data set. As mentioned earlier,
it is possible for a network to fit the noise and thus memorize the
specific training data set rather then generalize from it.  This can
be a particular problem if the data are noisy or if the network is
very complex (i.e.\ has many hidden nodes so is flexible to model many
nonlinearities). We can test for this by training the network on only
a subset of the available data and then evaluating it on the rest, ensuring
that the two data subsets have similar distributions of objects.

\vskip 12 pt
\centerline{\bf 5.2.\  Form of the inputs}
\bigskip
\noindent
The input to the network is some measured feature vector. In the case
of stellar classification this will be a spectrum or a list of fluxes
through various filters, or both. We could also input other
measurements we believe to be relevant, such as the parallax or
interstellar extinction. The inputs need not be contiguous or in any
particular order. A stellar spectrum typically contains a large amount
of redundant information (from the point of view of general
classification), so using the entire spectrum as the input may not be
necessary. If we can reduce the number of input nodes, we reduce the
number of weights and hence the dimensionality of the mapping
function. This results in faster training, and may reduce the chance
of getting stuck in local minima. It also increases the data/weights
ratio (see section 5.5) and density of the data, thus providing a more
reliable mapping. There are a number of approaches to {\it
dimensionality reduction}. The simplest is to remove any inputs which
we know, or suspect, to have little or no influence on the problem in
hand. A more sophisticated approach is with Principal Component
Analysis (PCA), which finds the most discriminatory linear
combinations of the inputs. By removing the less significant principal
components, we can compress the data with negligible loss of
discriminatory information, and filter out some noise too. PCA is
discussed in the article by Singh, Bailer-Jones \& Gupta in these
proceedings.

\vskip 12 pt
\centerline{\bf 5.3.\ Output modes}
\bigskip
\noindent
The output vector is the set of classifications.  For example, we may
have one output for MK spectral type (SpT) and another for luminosity
class (LC). We could alternatively use separate networks for each
parameter.  As the network outputs are real numbers, these classes are
represented on a continuous scale. Thus the network will be able to
produce intermediate classes. This is useful for the spectral type
parameter, as spectral types are simply points on what is really a
continuous scale (closely related to effective temperature). We refer
to this way of representing the outputs as {\it continuous
mode}. Luminosity classes, on the other hand, have more claim to being
discrete, for which we may want to use the network in {\it
probabilistic mode}, whereby each class is represented by a separate
output.  By confining each output to the range 0--1 (e.g.\ by replacing
equation 3 with the function $0.5 \times [1 + \tanh]$), we can
interpret the output as the probability that the input is of that
class.  For a five class problem, the target vector for an object in
class 2 would be $(0,1,0,0,0)$. When we apply the trained network to
new data, we would assign the class with the largest probability,
although if that probability is below some threshold, or not much
larger than the probability of some other class, then we may prefer to
label this object as uncertain.

Clearly, the probabilistic mode offers more information than the
continuous one. For example, if the input is a composite of two
classes (e.g.\ two stars of different brightnesses), we can represent
both in the target vector (e.g.\ 0,0.25,0,0.75,0). Thus if the network
thinks an input is composite, it can (in principle at least) tell us
this.  The only way we could model binary star spectra in
continuous mode is to have two spectral type outputs, although then
we also need some means of knowing when the input spectrum is just a
single star.  See Weaver (2000) for an attempt to classify binary
spectra with neural networks. A disadvantage of the probabilistic mode
is that it may result in many more network weights that a network in
continuous mode.

\vskip 12 pt
\centerline{\bf 5.4.\ How many hidden nodes and layers are required?}
\bigskip
\noindent
Deciding upon the appropriate number of hidden nodes and layers is
largely a matter of experience. With many problems, sufficient
accuracy can be obtained with one or two hidden layers and 5--10
hidden nodes in those layers. There is a theorem which states that a
network with one hidden layer can approximate any continuous function
to arbitrary accuracy, provided it has a sufficient number of hidden
nodes (Hornick, Stinchcombe and White 1990). In practice,
such a large number of nodes may be required that it is more efficient
to go to a second hidden layer. This can be seen in the work of
Bailer-Jones et al.\ (1998b). Fig.~4 shows how the error in
determining the spectral type (in continuous mode) varies as a
function of the number of hidden nodes in a network with just one
hidden layer.  The error is already reasonably small for a network
with just one hidden node, which is equivalent to a network with no
hidden layer. Adding a few hidden nodes improves the error, but once
there are 5--10 nodes, the error does not decrease further. However,
if we add a second hidden layer, then a network with 5 nodes in each
of these layers produces a $\sigma_{\rm 68}$ error of 0.85 SpT, a
statistically significant improvement. A similar improvement using a
second hidden layer was seen again in a related (but more complex)
problem (Bailer-Jones 2000).

\noindent
\midinsert
\vskip 6.5cm
\special{psfile=err03.ps voffset=185mm hoffset=60mm vscale=30mm hscale=30mm angle=270}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 4.} Variation of the network classification error as a
function of the number of hidden nodes, $q$. The solid line is for
\sig68\ and the dotted line for \sigrms; they represent the core and
outliers of the error distribution respectively (see Bailer-Jones et
al.\ 1998b).  The network architecture is 50:$q$:1 with a continuous
spectral type output.  The error bars are $3 \times \delta$, where
$\delta$ is the standard deviation in \sig68\ due to a finite sample
size. The mean RMS error in the target classifications themselves is
0.65 SpT.  }
\endinsert

\vskip 12 pt
\centerline{\bf 5.5.\ Ensuring a well-determined network solution}
\bigskip
\noindent
An important consideration in setting up the network is the ratio of
data to weights. If we train a network with one output using $N$
training vectors, we have a total of $N$ error measurements with which
to determine the $W$ weights. Thus we have a system of $N$ equations
with $W$ unknowns, and, if the weights and training vectors are
independent of each other, we require $N > W$ to find a solution.
With fewer data the solution for the weights will be
underdetermined. Thus a 10:5:5:1 network (indicating 10 inputs, 5
nodes in each of two hidden layers, and one output) has 91 weights
(including the bias nodes), so will require at least as many training
examples. The data/weight ratio estimation is more complex with
multiple outputs (especially if they are correlated, as will be the
case in probabilistic mode), but if there are two independent outputs,
then the number of error measures is $2N$.

In practice, this ``overdetermination'' requirement may not be so
strict.  For example, Gulati et al.\ (1994) used only 55 training
examples to train a 161:64:64:55 probabilistic network. This network
has about 18,000 free parameters (each output is parametrized by about
15,000 weights), yet only of order 55 error measures were available to
determine them. In spite of this, a reasonably small error on a test
data set was obtained. We suspect that many weights were essentially
not trained, and, by staying at their small initial values, played no
role in the network. Additionally, correlations between the inputs
(flux bins of stellar spectra) could result in correlated weights in
the input--hidden layer. Both of these reduce the effective number of
free parameters. Other cases of apparently good results with formally
underdetermined weights appear in the literature, implying that
correlated or unused weights are not uncommon. Nonetheless, in many of
these cases simpler networks could probably have been used without
loss of accuracy.


\vskip 12 pt
\centerline{\bf 6.\ Example Applications}
\bigskip
\noindent
The paper by Gupta et al.\ in these proceedings gives an application
of a 93:16:16:17 network to classify 2000 bright sources from IRAS LRS
spectra.

Fig.~5 shows the performance of a committee of ten 50:5:1 networks
applied to the spectral classification problem, trained using gradient
descent. It shows the correlation between the network classifications
and the ``true'' classifications (those from the training catalogue),
and the histogram of the error residuals, for a test data set. There
is an even scatter about the ideal classifications in the left hand
plot, except for the hottest (early type) stars, where the network
assigns cooler classes.

\noindent
\midinsert
\vskip 5.8cm
\special{psfile=exampleres.eps voffset=0mm hoffset=22mm vscale=70mm hscale=70mm angle=0}
\vskip 0.3cm
{\eightpoint % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 5.} Results from a committee of ten 50:5:1 networks with a
single continuous spectral type output. The left plot shows the
correlation between network and catalogue classifications (perfect
classifications would fall on the diagonal line), and the residuals
from this are plotted on the right.  The spectral types O3 to M9 have
been coded by the integers 1 to 57. The summary results are
$\sigma_{68} = 1.03$ SpT and $\sigma_{RMS} = 1.35$ SpT. The internal
error of the networks is 0.41 SpT. This is a measure of the average
differences in classifications which members of the committee assign
to a given spectrum.  }
\endinsert

A similar plot to that in Fig.~5 is not very useful for a network used
in probabilistic mode when there are only a handful of classes. In
such cases it is more useful to look at the {\it confusion matrix}.
Examples are shown in Fig.~6 for networks with three outputs,
corresponding to the three luminosity classes III, IV and V. In
probabilistic mode, we can also examine the confidence with which the
network classifies objects. An example is shown in Fig.~7, where we
plot the cumulative fraction of objects classified with a probability
greater than some amount.  This plot is for the 50:5:5:3 network with
the line+continuum spectra from Fig.~6.  Thus while 98.1\% of the
class V objects were correctly classified (i.e.\ a fraction 0.981 were
classified with a probability, $p>1/3$), only 77.5\% have a
probability of over 0.90. Part of the reason for this is that a
network with sigmoidal outputs can never produce an output of exactly
0.0 or 1.0, as that requires very large weights, something which
regularization tries to avoid. (See Bailer-Jones 1996 and Bailer-Jones
et al.\ 1998b for more details on these examples.)


\noindent
\midinsert
\vskip 8.5cm
\special{psfile=lumres.ps voffset=268mm hoffset=-5mm vscale=50mm hscale=50mm angle=270}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 6.}  Confusion matrices for four different network models
differing in the number of PCA inputs and whether or not the continuum
has been subtracted from the spectra.  Each matrix lists the
percentage of spectra which have been classified correctly and
incorrectly: in the top-left matrix the network classifies 1.5\% of
spectra which are class V in the catalogue as class III, whereas it
correctly classifies 97.7\% of class Vs. In all four cases the
networks' do not believe in the existence of luminosity class IV
(random assignments would lead to about 33\% in the middle row).  }
\endinsert

\noindent
\midinsert
\vskip 13.5cm
\special{psfile=conf08.ps voffset=-10mm hoffset=30mm vscale=50mm hscale=50mm angle=0}
\vskip 0.3cm
{\eightpoint   % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 7.}  Confidence of correct network classifications based
on network output probability.  Each figure is a cumulative plot of
the fraction of spectra which have been classified with a probability
$>p$, for the luminosity classes III, IV and V.}
\endinsert


\vskip 12 pt
\centerline{\bf 7.\ Understanding the network}
\bigskip
\noindent
The underlying principle of a neural network is that it learns its
input--output mapping from the relationship intrinsically present in a
preclassified data set. Unlike rule-based classification algorithms,
the networks are not explicitly told anything about the relevance of
certain inputs in determining the output. This relevance information
is expressed in the optimized weight vectors, and in principle we can
look at the values of the weights to find out which inputs are
significant in determining the outputs.  In practice, however, this
may be rather difficult on account of there being more than one hidden
node (which is necessary to get a nonlinear mapping). For example,
weights leading from a relevant input to different nodes in the hidden
layer could be both small and large in magnitude. There is no
constraint on the network training to provide easily-interpretable
weights, although such a one may be possible. Furthermore, there may
be many minima in the weights space which have approximately
equal errors, and the corresponding ``optimal'' weight vectors may
give conflicting opinions about the relevance of the inputs. A
degeneracy in the weights is made more likely if there are redundant,
or correlated, inputs, which is often the case with stellar
spectra. Ultimately, obtaining both an accurate model of a complex
problem {\it and} a good understanding are conflicting demands, as
understanding often requires simplification.


\noindent
\midinsert
\vskip 17.0cm
\special{psfile=node05.ps voffset=-30mm hoffset=-10mm vscale=65mm hscale=65mm angle=0}
\vskip 0.3cm
{\eightpoint % Switch to 8 pt fonts for figure captions
\noindent
{\bf Figure 8.}  Correlation between the node values and outputs for a
network trained to give spectral types in continuous mode.  The column
on the left (right) represents the first (second) hidden layer of
nodes in a 50:5:5:1 network.  Each point is the node value produced by
one input vector in an test data set. The transfer function in each
layer (including the subsequent output layer) is a sigmoid function,
which restricts the node values to the range 0--1.  }
\endinsert

We can nonetheless get some idea of {\it how} a network produces its
classifications. Fig.~8 shows the correlations between the values on
the hidden nodes of a network and the outputs from that
network. Looking at the first hidden layer, we see that no single node
shows a low-scatter, monotonic relationship with spectral type. If one
did, then that alone could be used for spectral type classification
and we could dispense with the network altogether. Generally, the
correlations are nonlinear, and with large scatter.  However, node 1
does show linearity over a limited part of the spectral range (SpT =
28--42, or F0--K0). This in interesting, as it means that for classes
F and G this node could serve as the output node and we would get
reasonable classifications. Looking back at equation 1, we see that if
we were to re-parametrize the spectral type as $tanh^{-1}$ SpT, we
could classify F and G stars using a simple linear regression on the
spectral indices (as $tanh^{-1}$ is monotonic over the range of
interest).  Nodes 3--5 also show linear correlations over small
portions of the spectral type range. However, all nodes show either
some insensitivity (due to saturation of the node values) or large
scatter in the correlation, for at least some part of the spectral
type range.

When these values are passed to the second hidden layer, the situation
is somewhat different. The correlations are still nonlinear, but there
is much less scatter. Furthermore, node 5 has a low scatter, {\it
linear} correlation for spectral types earlier than about F3 (SpT =
30), and node 4 shows a linear correlation for types later than about
K0 (SpT = 42). Node 2 (or 3) alone could reasonably classify all intermediate
types. Thus we see that there is some localization in the network,
with certain nodes taking on the job of identifying certain ranges of
the target variable. There is, of course, no significance in the
numbering of the nodes within a layer, as the nodes could be shuffled
with no change in the final output.


\vskip 12 pt
\centerline{\bf 8.\ Summary}
\bigskip
\noindent
We have introduced the basic concepts of feedforward neural networks.
Although the type introduced, the multilayer perceptron, is only one
of many types, it is very important as judged by its wide and numerous
applications.  We have described the backpropagation algorithm for
training the network, as well as some of the ways in which it can be
implemented to optimize the weights. We also briefly looked at a
Bayesian approach to modelling with a network, which marginalizes over
the weights rather than optimizing them. A number of implementation
issues were discussed and example applications were presented.

Although it is interesting to see how these models were originally
inspired by attempts to model brain behaviour, these artificial neural
networks really have very little in common with biological neural
networks, and, in our opinion, are better understood in purely
mathematical terms. Whilst they may be powerful models, there is
nothing ``mysterious'' about them. They are simply parametrized
nonlinear models: the input--output function of a two-hidden layer
model is described entirely by equations 1--3. They are no more a
``black box'' than is a straight line fit to two-dimensional data. The
nonlinearity and multidimensionality of neural networks are present
because real-world problems are often complex. If we have difficulty
understanding how they work, then it is because an {\it accurate}
description of real-word problems often requires that we avoid the
simplifying assumptions which typically make easily-understood
descriptions only approximate.

We finish by stressing that the purpose of these models is to
generalize the relationship between a data domain and a parameter
domain which is intrinsically present in a preclassified data set.
They are useful because they can model complex relationships in a
large number of parameters, and generalize this to areas of the input
space where training data are sparse. This has potential applications
to large, multiparameter astronomical surveys (see Bailer-Jones, these
proceedings). But the models can only be as good as the data on which
they are trained (although they can be a lot worse if trained badly),
so construction of an accurate and relevant training data set is of
fundamental importance.

% Sample references
\bigskip
\centerline{\bf References}
\bigskip
{\eightpoint\parindent=0pt\everypar={\hangindent=0.5 cm}
% References in the format of the Proceedings
% using 8 pt fonts
% leave one line blank between two references to force a paragraph break
%scussions.

Bailer-Jones, C.A.L., 1996, PhD thesis, Univ.\ Cambridge

Bailer-Jones, C.A.L., 2000, A\&A, 357, 197

Bailer-Jones, C.A.L., MacKay, D.J.C., Withers, P.J., 1998a, Network:
Computation in Neural Systems, 9, 531

Bailer-Jones, C.A.L., Irwin, M., von Hippel T., 1998b, MNRAS, 298, 361

Beale, R., Jackson, T., {\it Neural Computing: An Introduction}, 1990,
Institute of Physics: Bristol

Bishop, C.M., {\it Neural Networks for Pattern Recognition}, 1995,
Oxford University Press: Oxford

Gulati, R.K., Gupta, R., Gothoskar, P., Khobragade, S., 1994, ApJ, 426, 340

Hertz, J., Krogh, A., Palmer, R.G., 1991, {\it Introduction to the
Theory of Neural Computation}, Addison-Wesley: Wokingham, England

Hornick, K., Stinchcombe, M., White, H., 1990, Neural Networks, 3, 551

MacKay, D.J.C., 1995, Network: Computation in Neural Systems, 6, 469

Rumelhart, D.E., Hinton, G.E., Williams, R.J., 1986a, in 
D.E.\ Rumelhart, J.L.\ McClelland J.L, the PDP Research Group, eds, 
Parallel Distributed Processing: Explorations in the
  Microstructure of Cognition, MIT Press: Boston, p.\ 318

Rumelhart, D.E., Hinton, G.E., Williams, R.J., 1986b, Nature, 323, 533

Sivia, D.S., 1996, {\it Data Analysis. A Bayesian Tutorial}, Oxford
University Press: Oxford

Weaver, W.B., 2000, ApJ, 541, 298

}                                         % End of references
% leave one line blank before the closing braces.
\end

psfig.tex
% Psfig/TeX Release 1.8
% dvips version
%
% All psfig/tex software, documentation, and related files
% in this distribution of psfig/tex are 
% Copyright 1987, 1988, 1991 Trevor J. Darrell
%
% Permission is granted for use and non-profit distribution of psfig/tex 
% providing that this notice is clearly maintained. The right to
% distribute any portion of psfig/tex for profit or as part of any commercial
% product is specifically reserved for the author(s) of that portion.
%
% *** Feel free to make local modifications of psfig as you wish,
% *** but DO NOT post any changed or modified versions of ``psfig''
% *** directly to the net. Send them to me and I'll try to incorporate
% *** them into future versions. If you want to take the psfig code 
% *** and make a new program (subject to the copyright above), distribute it, 
% *** (and maintain it) that's fine, just don't call it psfig.
%
% Bugs and improvements to trevor@media.mit.edu.
%
% Thanks to Greg Hager (GDH) and Ned Batchelder for their contributions
% to the original version of this project.
%
% Modified by J. Daniel Smith on 9 October 1990 to accept the
% %%BoundingBox: comment with or without a space after the colon.  Stole
% file reading code from Tom Rokicki's EPSF.TEX file (see below).
%
% More modifications by J. Daniel Smith on 29 March 1991 to allow the
% the included PostScript figure to be rotated.  The amount of
% rotation is specified by the "angle=" parameter of the \psfig command.
%
% Modified by Robert Russell on June 25, 1991 to allow users to specify
% .ps filenames which don't yet exist, provided they explicitly provide
% boundingbox information via the \psfig command. Note: This will only work
% if the "file=" parameter follows all four "bb???=" parameters in the
% command. This is due to the order in which psfig interprets these params.
%
%  3 Jul 1991	JDS	check if file already read in once
%  4 Sep 1991	JDS	fixed incorrect computation of rotated
%			bounding box
% 25 Sep 1991	GVR	expanded synopsis of \psfig
% 14 Oct 1991	JDS	\fbox code from LaTeX so \psdraft works with TeX
%			changed \typeout to \ps@typeout
% 17 Oct 1991	JDS	added \psscalefirst and \psrotatefirst
%

% From: gvr@cs.brown.edu (George V. Reilly)
%
% \psdraft	draws an outline box, but doesn't include the figure
%		in the DVI file.  Useful for previewing.
%
% \psfull	includes the figure in the DVI file (default).
%
% \psscalefirst width= or height= specifies the size of the figure
% 		before rotation.
% \psrotatefirst (default) width= or height= specifies the size of the
% 		 figure after rotation.  Asymetric figures will
% 		 appear to shrink.
%
% \psfigurepath#1	sets the path to search for the figure
%
% \psfig
% usage: \psfig{file=, figure=, height=, width=,
%			bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=, angle=, silent=}
%
%	"file" is the filename.  If no path name is specified and the
%		file is not found in the current directory,
%		it will be looked for in directory \psfigurepath.
%	"figure" is a synonym for "file".
%	By default, the width and height of the figure are taken from
%		the BoundingBox of the figure.
%	If "width" is specified, the figure is scaled so that it has
%		the specified width.  Its height changes proportionately.
%	If "height" is specified, the figure is scaled so that it has
%		the specified height.  Its width changes proportionately.
%	If both "width" and "height" are specified, the figure is scaled
%		anamorphically.
%	"bbllx", "bblly", "bburx", and "bbury" control the PostScript
%		BoundingBox.  If these four values are specified
%               *before* the "file" option, the PSFIG will not try to
%               open the PostScript file.
%	"rheight" and "rwidth" are the reserved height and width
%		of the figure, i.e., how big TeX actually thinks
%		the figure is.  They default to "width" and "height".
%	The "clip" option ensures that no portion of the figure will
%		appear outside its BoundingBox.  "clip=" is a switch and
%		takes no value, but the `=' must be present.
%	The "angle" option specifies the angle of rotation (degrees, ccw).
%	The "silent" option makes \psfig work silently.
%

% check to see if macros already loaded in (maybe some other file says
% "\input psfig") ...
\ifx\undefined\psfig\else\endinput\fi

%
% from a suggestion by eijkhout@csrd.uiuc.edu to allow
% loading as a style file:
\edef\psfigRestoreAt{\catcode`@=\number\catcode`@\relax}
\catcode`\@=11\relax
\newwrite\@unused
\def\ps@typeout#1{{\let\protect\string\immediate\write\@unused{#1}}}
\ps@typeout{psfig/tex 1.8}

%% Here's how you define your figure path.  Should be set up with null
%% default and a user useable definition.

\def\figurepath{./}
\def\psfigurepath#1{\edef\figurepath{#1}}

%
% @psdo control structure -- similar to Latex @for.
% I redefined these with different names so that psfig can
% be used with TeX as well as LaTeX, and so that it will not 
% be vunerable to future changes in LaTeX's internal
% control structure,
%
\def\@nnil{\@nil}
\def\@empty{}
\def\@psdonoop#1\@@#2#3{}
\def\@psdo#1:=#2\do#3{\edef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \expandafter\@psdoloop#2,\@nil,\@nil\@@#1{#3}\fi}
\def\@psdoloop#1,#2,#3\@@#4#5{\def#4{#1}\ifx #4\@nnil \else
       #5\def#4{#2}\ifx #4\@nnil \else#5\@ipsdoloop #3\@@#4{#5}\fi\fi}
\def\@ipsdoloop#1,#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@ipsdoloop\fi\@nextwhile#2\@@#3{#4}}
\def\@tpsdo#1:=#2\do#3{\xdef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \@tpsdoloop#2\@nil\@nil\@@#1{#3}\fi}
\def\@tpsdoloop#1#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@tpsdoloop\fi\@nextwhile#2\@@#3{#4}}
% 
% \fbox is defined in latex.tex; so if \fbox is undefined, assume that
% we are not in LaTeX.
% Perhaps this could be done better???
\ifx\undefined\fbox
% \fbox code from modified slightly from LaTeX
\newdimen\fboxrule
\newdimen\fboxsep
\newdimen\ps@tempdima
\newbox\ps@tempboxa
\fboxsep = 3pt
\fboxrule = .4pt
\long\def\fbox#1{\leavevmode\setbox\ps@tempboxa\hbox{#1}\ps@tempdima\fboxrule
    \advance\ps@tempdima \fboxsep \advance\ps@tempdima \dp\ps@tempboxa
   \hbox{\lower \ps@tempdima\hbox
  {\vbox{\hrule height \fboxrule
          \hbox{\vrule width \fboxrule \hskip\fboxsep
          \vbox{\vskip\fboxsep \box\ps@tempboxa\vskip\fboxsep}\hskip 
                 \fboxsep\vrule width \fboxrule}
                 \hrule height \fboxrule}}}}
\fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file reading stuff from epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%   Revised by Tomas Rokicki to accept bounding boxes with no
%      space after the colon, 18 Jul 1990.
%   Portions modified/removed for use in PSFIG package by
%      J. Daniel Smith, 9 October 1990.
%
\newread\ps@stream
\newif\ifnot@eof       % continue looking for the bounding box?
\newif\if@noisy        % report what you're making?
\newif\if@atend        % %%BoundingBox: has (at end) specification
\newif\if@psfile       % does this look like a PostScript file?
%
% PostScript files should start with `%!'
%
{\catcode`\%=12\global\gdef\epsf@start{%!}}
\def\epsf@PS{PS}
%
\def\epsf@getbb#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\ps@stream=#1
\ifeof\ps@stream\ps@typeout{Error, File #1 not found}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\not@eoftrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \if@psfile
	  \read\ps@stream to \epsf@fileline
       \else{
	  \obeyspaces
          \read\ps@stream to \epsf@tmp\global\let\epsf@fileline\epsf@tmp}
       \fi
       \ifeof\ps@stream\not@eoffalse\else
%
%   Check the first line for `%!'.  Issue a warning message if its not
%   there, since the file might not be a PostScript file.
%
       \if@psfile\else
       \expandafter\epsf@test\epsf@fileline:. \\%
       \fi
%
%   We check to see if the first character is a % sign;
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox:' and the `(atend)' specification was not found.
%   That is, the only way to stop is when the end of file is reached,
%   or a `%%BoundingBox: llx lly urx ury' line is found.
%
          \expandafter\epsf@aux\epsf@fileline:. \\%
       \fi
   \ifnot@eof\repeat
   }\closein\ps@stream\fi}%
%
% This tests if the file we are reading looks like a PostScript file.
%
\long\def\epsf@test#1#2#3:#4\\{\def\epsf@testit{#1#2}
			\ifx\epsf@testit\epsf@start\else
\ps@typeout{Warning! File does not start with `\epsf@start'.  It may not be a PostScript file.}
			\fi
			\@psfiletrue} % don't test after 1st line
%
%   We still need to define the tricky \epsf@aux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12\global\let\epsf@percent=%\global\def\epsf@bblit{%BoundingBox}}
%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.  We continue searching if `(at end)'
%   was found after the `%BoundingBox:'.
%
\long\def\epsf@aux#1#2:#3\\{\ifx#1\epsf@percent
   \def\epsf@testit{#2}\ifx\epsf@testit\epsf@bblit
	\@atendfalse
        \epsf@atend #3 . \\%
	\if@atend	
	   \if@verbose{
		\ps@typeout{psfig: found `(atend)'; continuing search}
	   }\fi
        \else
        \epsf@grab #3 . . . \\%
        \not@eoffalse
        \global\no@bbfalse
        \fi
   \fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsf@grab #1 #2 #3 #4 #5\\{%
   \global\def\epsf@llx{#1}\ifx\epsf@llx\empty
      \epsf@grab #2 #3 #4 #5 .\\\else
   \global\def\epsf@lly{#2}%
   \global\def\epsf@urx{#3}\global\def\epsf@ury{#4}\fi}%
%
% Determine if the stuff following the %%BoundingBox is `(atend)'
% J. Daniel Smith.  Copied from \epsf@grab above.
%
\def\epsf@atendlit{(atend)} 
\def\epsf@atend #1 #2 #3\\{%
   \def\epsf@tmp{#1}\ifx\epsf@tmp\empty
      \epsf@atend #2 #3 .\\\else
   \ifx\epsf@tmp\epsf@atendlit\@atendtrue\fi\fi}


% End of file reading stuff from epsf.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trigonometry stuff from "trig.tex"
\chardef\letter = 11
\chardef\other = 12

\newif \ifdebug %%% turn me on to see TeX hard at work ...
\newif\ifc@mpute %%% don't need to compute some values
\c@mputetrue % but assume that we do

\let\then = \relax
\def\r@dian{pt }
\let\r@dians = \r@dian
\let\dimensionless@nit = \r@dian
\let\dimensionless@nits = \dimensionless@nit
\def\internal@nit{sp }
\let\internal@nits = \internal@nit
\newif\ifstillc@nverging
\def \Mess@ge #1{\ifdebug \then \message {#1} \fi}

{ %%% Things that need abnormal catcodes %%%
	\catcode `\@ = \letter
	\gdef \nodimen {\expandafter \n@dimen \the \dimen}
	\gdef \term #1 #2 #3%
	       {\edef \t@ {\the #1}%%% freeze parameter 1 (count, by value)
		\edef \t@@ {\expandafter \n@dimen \the #2\r@dian}%
				   %%% freeze parameter 2 (dimen, by value)
		\t@rm {\t@} {\t@@} {#3}%
	       }
	\gdef \t@rm #1 #2 #3%
	       {{%
		\count 0 = 0
		\dimen 0 = 1 \dimensionless@nit
		\dimen 2 = #2\relax
		\Mess@ge {Calculating term #1 of \nodimen 2}%
		\loop
		\ifnum	\count 0 < #1
		\then	\advance \count 0 by 1
			\Mess@ge {Iteration \the \count 0 \space}%
			\Multiply \dimen 0 by {\dimen 2}%
			\Mess@ge {After multiplication, term = \nodimen 0}%
			\Divide \dimen 0 by {\count 0}%
			\Mess@ge {After division, term = \nodimen 0}%
		\repeat
		\Mess@ge {Final value for term #1 of 
				\nodimen 2 \space is \nodimen 0}%
		\xdef \Term {#3 = \nodimen 0 \r@dians}%
		\aftergroup \Term
	       }}
	\catcode `\p = \other
	\catcode `\t = \other
	\gdef \n@dimen #1pt{#1} %%% throw away the ``pt''
}

\def \Divide #1by #2{\divide #1 by #2} %%% just a synonym

\def \Multiply #1by #2%%% allows division of a dimen by a dimen
       {{%%% should really freeze parameter 2 (dimen, passed by value)
	\count 0 = #1\relax
	\count 2 = #2\relax
	\count 4 = 65536
	\Mess@ge {Before scaling, count 0 = \the \count 0 \space and
			count 2 = \the \count 2}%
	\ifnum	\count 0 > 32767 %%% do our best to avoid overflow
	\then	\divide \count 0 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 0 < -32767
		\then	\divide \count 0 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\ifnum	\count 2 > 32767 %%% while retaining reasonable accuracy
	\then	\divide \count 2 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 2 < -32767
		\then	\divide \count 2 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\multiply \count 0 by \count 2
	\divide \count 0 by \count 4
	\xdef \product {#1 = \the \count 0 \internal@nits}%
	\aftergroup \product
       }}

\def\r@duce{\ifdim\dimen0 > 90\r@dian \then   % sin(x+90) = sin(180-x)
		\multiply\dimen0 by -1
		\advance\dimen0 by 180\r@dian
		\r@duce
	    \else \ifdim\dimen0 < -90\r@dian \then  % sin(-x) = sin(360+x)
		\advance\dimen0 by 360\r@dian
		\r@duce
		\fi
	    \fi}

\def\Sine#1%
       {{%
	\dimen 0 = #1 \r@dian
	\r@duce
	\ifdim\dimen0 = -90\r@dian \then
	   \dimen4 = -1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 90\r@dian \then
	   \dimen4 = 1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 0\r@dian \then
	   \dimen4 = 0\r@dian
	   \c@mputefalse
	\fi
%
	\ifc@mpute \then
        	% convert degrees to radians
		\divide\dimen0 by 180
		\dimen0=3.141592654\dimen0
%
		\dimen 2 = 3.1415926535897963\r@dian %%% a well-known constant
		\divide\dimen 2 by 2 %%% we only deal with -pi/2 : pi/2
		\Mess@ge {Sin: calculating Sin of \nodimen 0}%
		\count 0 = 1 %%% see power-series expansion for sine
		\dimen 2 = 1 \r@dian %%% ditto
		\dimen 4 = 0 \r@dian %%% ditto
		\loop
			\ifnum	\dimen 2 = 0 %%% then we've done
			\then	\stillc@nvergingfalse 
			\else	\stillc@nvergingtrue
			\fi
			\ifstillc@nverging %%% then calculate next term
			\then	\term {\count 0} {\dimen 0} {\dimen 2}%
				\advance \count 0 by 2
				\count 2 = \count 0
				\divide \count 2 by 2
				\ifodd	\count 2 %%% signs alternate
				\then	\advance \dimen 4 by \dimen 2
				\else	\advance \dimen 4 by -\dimen 2
				\fi
		\repeat
	\fi		
			\xdef \sine {\nodimen 4}%
       }}

% Now the Cosine can be calculated easily by calling \Sine
\def\Cosine#1{\ifx\sine\UnDefined\edef\Savesine{\relax}\else
		             \edef\Savesine{\sine}\fi
	{\dimen0=#1\r@dian\advance\dimen0 by 90\r@dian
	 \Sine{\nodimen 0}
	 \xdef\cosine{\sine}
	 \xdef\sine{\Savesine}}}	      
% end of trig stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\psdraft{
	\def\@psdraft{0}
	%\ps@typeout{draft level now is \@psdraft \space . }
}
\def\psfull{
	\def\@psdraft{100}
	%\ps@typeout{draft level now is \@psdraft \space . }
}

\psfull

\newif\if@scalefirst
\def\psscalefirst{\@scalefirsttrue}
\def\psrotatefirst{\@scalefirstfalse}
\psrotatefirst

\newif\if@draftbox
\def\psnodraftbox{
	\@draftboxfalse
}
\def\psdraftbox{
	\@draftboxtrue
}
\@draftboxtrue

\newif\if@prologfile
\newif\if@postlogfile
\def\pssilent{
	\@noisyfalse
}
\def\psnoisy{
	\@noisytrue
}
\psnoisy
%%% These are for the option list.
%%% A specification of the form a = b maps to calling \@p@@sa{b}
\newif\if@bbllx
\newif\if@bblly
\newif\if@bburx
\newif\if@bbury
\newif\if@height
\newif\if@width
\newif\if@rheight
\newif\if@rwidth
\newif\if@angle
\newif\if@clip
\newif\if@verbose
\def\@p@@sclip#1{\@cliptrue}


\newif\if@decmpr

%%% GDH 7/26/87 -- changed so that it first looks in the local directory,
%%% then in a specified global directory for the ps file.
%%% RPR 6/25/91 -- changed so that it defaults to user-supplied name if
%%% boundingbox info is specified, assuming graphic will be created by
%%% print time.
%%% TJD 10/19/91 -- added bbfile vs. file distinction, and @decmpr flag

\def\@p@@sfigure#1{\def\@p@sfile{null}\def\@p@sbbfile{null}
	        \openin1=#1.bb
		\ifeof1\closein1
	        	\openin1=\figurepath#1.bb
			\ifeof1\closein1
			        \openin1=#1
				\ifeof1\closein1%
				       \openin1=\figurepath#1
					\ifeof1
					   \ps@typeout{Error, File #1 not found}
						\if@bbllx\if@bblly
				   		\if@bburx\if@bbury
			      				\def\@p@sfile{#1}%
			      				\def\@p@sbbfile{#1}%
							\@decmprfalse
				  	   	\fi\fi\fi\fi
					\else\closein1
				    		\def\@p@sfile{\figurepath#1}%
				    		\def\@p@sbbfile{\figurepath#1}%
						\@decmprfalse
	                       		\fi%
			 	\else\closein1%
					\def\@p@sfile{#1}
					\def\@p@sbbfile{#1}
					\@decmprfalse
			 	\fi
			\else
				\def\@p@sfile{\figurepath#1}
				\def\@p@sbbfile{\figurepath#1.bb}
				\@decmprtrue
			\fi
		\else
			\def\@p@sfile{#1}
			\def\@p@sbbfile{#1.bb}
			\@decmprtrue
		\fi}

\def\@p@@sfile#1{\@p@@sfigure{#1}}

\def\@p@@sbbllx#1{
		%\ps@typeout{bbllx is #1}
		\@bbllxtrue
		\dimen100=#1
		\edef\@p@sbbllx{\number\dimen100}
}
\def\@p@@sbblly#1{
		%\ps@typeout{bblly is #1}
		\@bbllytrue
		\dimen100=#1
		\edef\@p@sbblly{\number\dimen100}
}
\def\@p@@sbburx#1{
		%\ps@typeout{bburx is #1}
		\@bburxtrue
		\dimen100=#1
		\edef\@p@sbburx{\number\dimen100}
}
\def\@p@@sbbury#1{
		%\ps@typeout{bbury is #1}
		\@bburytrue
		\dimen100=#1
		\edef\@p@sbbury{\number\dimen100}
}
\def\@p@@sheight#1{
		\@heighttrue
		\dimen100=#1
   		\edef\@p@sheight{\number\dimen100}
		%\ps@typeout{Height is \@p@sheight}
}
\def\@p@@swidth#1{
		%\ps@typeout{Width is #1}
		\@widthtrue
		\dimen100=#1
		\edef\@p@swidth{\number\dimen100}
}
\def\@p@@srheight#1{
		%\ps@typeout{Reserved height is #1}
		\@rheighttrue
		\dimen100=#1
		\edef\@p@srheight{\number\dimen100}
}
\def\@p@@srwidth#1{
		%\ps@typeout{Reserved width is #1}
		\@rwidthtrue
		\dimen100=#1
		\edef\@p@srwidth{\number\dimen100}
}
\def\@p@@sangle#1{
		%\ps@typeout{Rotation is #1}
		\@angletrue
%		\dimen100=#1
		\edef\@p@sangle{#1} %\number\dimen100}
}
\def\@p@@ssilent#1{ 
		\@verbosefalse
}
\def\@p@@sprolog#1{\@prologfiletrue\def\@prologfileval{#1}}
\def\@p@@spostlog#1{\@postlogfiletrue\def\@postlogfileval{#1}}
\def\@cs@name#1{\csname #1\endcsname}
\def\@setparms#1=#2,{\@cs@name{@p@@s#1}{#2}}
%
% initialize the defaults (size the size of the figure)
%
\def\ps@init@parms{
		\@bbllxfalse \@bbllyfalse
		\@bburxfalse \@bburyfalse
		\@heightfalse \@widthfalse
		\@rheightfalse \@rwidthfalse
		\def\@p@sbbllx{}\def\@p@sbblly{}
		\def\@p@sbburx{}\def\@p@sbbury{}
		\def\@p@sheight{}\def\@p@swidth{}
		\def\@p@srheight{}\def\@p@srwidth{}
		\def\@p@sangle{0}
		\def\@p@sfile{} \def\@p@sbbfile{}
		\def\@p@scost{10}
		\def\@sc{}
		\@prologfilefalse
		\@postlogfilefalse
		\@clipfalse
		\if@noisy
			\@verbosetrue
		\else
			\@verbosefalse
		\fi
}
%
% Go through the options setting things up.
%
\def\parse@ps@parms#1{
	 	\@psdo\@psfiga:=#1\do
		   {\expandafter\@setparms\@psfiga,}}
%
% Compute bb height and width
%
\newif\ifno@bb
\def\bb@missing{
	\if@verbose{
		\ps@typeout{psfig: searching \@p@sbbfile \space  for bounding box}
	}\fi
	\no@bbtrue
	\epsf@getbb{\@p@sbbfile}
        \ifno@bb \else \bb@cull\epsf@llx\epsf@lly\epsf@urx\epsf@ury\fi
}	
\def\bb@cull#1#2#3#4{
	\dimen100=#1 bp\edef\@p@sbbllx{\number\dimen100}
	\dimen100=#2 bp\edef\@p@sbblly{\number\dimen100}
	\dimen100=#3 bp\edef\@p@sbburx{\number\dimen100}
	\dimen100=#4 bp\edef\@p@sbbury{\number\dimen100}
	\no@bbfalse
}
% rotate point (#1,#2) about (0,0).
% The sine and cosine of the angle are already stored in \sine and
% \cosine.  The result is placed in (\p@intvaluex, \p@intvaluey).
\newdimen\p@intvaluex
\newdimen\p@intvaluey
\def\rotate@#1#2{{\dimen0=#1 sp\dimen1=#2 sp
%            	calculate x' = x \cos\theta - y \sin\theta
		  \global\p@intvaluex=\cosine\dimen0
		  \dimen3=\sine\dimen1
		  \global\advance\p@intvaluex by -\dimen3
% 		calculate y' = x \sin\theta + y \cos\theta
		  \global\p@intvaluey=\sine\dimen0
		  \dimen3=\cosine\dimen1
		  \global\advance\p@intvaluey by \dimen3
		  }}
\def\compute@bb{
		\no@bbfalse
		\if@bbllx \else \no@bbtrue \fi
		\if@bblly \else \no@bbtrue \fi
		\if@bburx \else \no@bbtrue \fi
		\if@bbury \else \no@bbtrue \fi
		\ifno@bb \bb@missing \fi
		\ifno@bb \ps@typeout{FATAL ERROR: no bb supplied or found}
			\no-bb-error
		\fi
		%
%\ps@typeout{BB: \@p@sbbllx, \@p@sbblly, \@p@sbburx, \@p@sbbury} 
%
% store height/width of original (unrotated) bounding box
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\ps@bbw{\number\count203}
		\edef\ps@bbh{\number\count204}
		%\ps@typeout{ psbbh = \ps@bbh, psbbw = \ps@bbw }
		\if@angle 
			\Sine{\@p@sangle}\Cosine{\@p@sangle}
	        	{\dimen100=\maxdimen\xdef\r@p@sbbllx{\number\dimen100}
					    \xdef\r@p@sbblly{\number\dimen100}
			                    \xdef\r@p@sbburx{-\number\dimen100}
					    \xdef\r@p@sbbury{-\number\dimen100}}
%
% Need to rotate all four points and take the X-Y extremes of the new
% points as the new bounding box.
                        \def\minmaxtest{
			   \ifnum\number\p@intvaluex<\r@p@sbbllx
			      \xdef\r@p@sbbllx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluex>\r@p@sbburx
			      \xdef\r@p@sbburx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluey<\r@p@sbblly
			      \xdef\r@p@sbblly{\number\p@intvaluey}\fi
			   \ifnum\number\p@intvaluey>\r@p@sbbury
			      \xdef\r@p@sbbury{\number\p@intvaluey}\fi
			   }
%			lower left
			\rotate@{\@p@sbbllx}{\@p@sbblly}
			\minmaxtest
%			upper left
			\rotate@{\@p@sbbllx}{\@p@sbbury}
			\minmaxtest
%			lower right
			\rotate@{\@p@sbburx}{\@p@sbblly}
			\minmaxtest
%			upper right
			\rotate@{\@p@sbburx}{\@p@sbbury}
			\minmaxtest
			\edef\@p@sbbllx{\r@p@sbbllx}\edef\@p@sbblly{\r@p@sbblly}
			\edef\@p@sbburx{\r@p@sbburx}\edef\@p@sbbury{\r@p@sbbury}
%\ps@typeout{rotated BB: \r@p@sbbllx, \r@p@sbblly, \r@p@sbburx, \r@p@sbbury}
		\fi
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\@bbw{\number\count203}
		\edef\@bbh{\number\count204}
		%\ps@typeout{ bbh = \@bbh, bbw = \@bbw }
}
%
% \in@hundreds performs #1 * (#2 / #3) correct to the hundreds,
%	then leaves the result in @result
%
\def\in@hundreds#1#2#3{\count240=#2 \count241=#3
		     \count100=\count240	% 100 is first digit #2/#3
		     \divide\count100 by \count241
		     \count101=\count100
		     \multiply\count101 by \count241
		     \advance\count240 by -\count101
		     \multiply\count240 by 10
		     \count101=\count240	%101 is second digit of #2/#3
		     \divide\count101 by \count241
		     \count102=\count101
		     \multiply\count102 by \count241
		     \advance\count240 by -\count102
		     \multiply\count240 by 10
		     \count102=\count240	% 102 is the third digit
		     \divide\count102 by \count241
		     \count200=#1\count205=0
		     \count201=\count200
			\multiply\count201 by \count100
		 	\advance\count205 by \count201
		     \count201=\count200
			\divide\count201 by 10
			\multiply\count201 by \count101
			\advance\count205 by \count201
			%
		     \count201=\count200
			\divide\count201 by 100
			\multiply\count201 by \count102
			\advance\count205 by \count201
			%
		     \edef\@result{\number\count205}
}
\def\compute@wfromh{
		% computing : width = height * (bbw / bbh)
		\in@hundreds{\@p@sheight}{\@bbw}{\@bbh}
		%\ps@typeout{ \@p@sheight * \@bbw / \@bbh, = \@result }
		\edef\@p@swidth{\@result}
		%\ps@typeout{w from h: width is \@p@swidth}
}
\def\compute@hfromw{
		% computing : height = width * (bbh / bbw)
	        \in@hundreds{\@p@swidth}{\@bbh}{\@bbw}
		%\ps@typeout{ \@p@swidth * \@bbh / \@bbw = \@result }
		\edef\@p@sheight{\@result}
		%\ps@typeout{h from w : height is \@p@sheight}
}
\def\compute@handw{
		\if@height 
			\if@width
			\else
				\compute@wfromh
			\fi
		\else 
			\if@width
				\compute@hfromw
			\else
				\edef\@p@sheight{\@bbh}
				\edef\@p@swidth{\@bbw}
			\fi
		\fi
}
\def\compute@resv{
		\if@rheight \else \edef\@p@srheight{\@p@sheight} \fi
		\if@rwidth \else \edef\@p@srwidth{\@p@swidth} \fi
		%\ps@typeout{rheight = \@p@srheight, rwidth = \@p@srwidth}
}
%		
% Compute any missing values
\def\compute@sizes{
	\compute@bb
	\if@scalefirst\if@angle
% at this point the bounding box has been adjsuted correctly for
% rotation.  PSFIG does all of its scaling using \@bbh and \@bbw.  If
% a width= or height= was specified along with \psscalefirst, then the
% width=/height= value needs to be adjusted to match the new (rotated)
% bounding box size (specifed in \@bbw and \@bbh).
%    \ps@bbw       width=
%    -------  =  ---------- 
%    \@bbw       new width=
% so `new width=' = (width= * \@bbw) / \ps@bbw; where \ps@bbw is the
% width of the original (unrotated) bounding box.
	\if@width
	   \in@hundreds{\@p@swidth}{\@bbw}{\ps@bbw}
	   \edef\@p@swidth{\@result}
	\fi
	\if@height
	   \in@hundreds{\@p@sheight}{\@bbh}{\ps@bbh}
	   \edef\@p@sheight{\@result}
	\fi
	\fi\fi
	\compute@handw
	\compute@resv}

%
% \psfig
% usage : \psfig{file=, height=, width=, bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=}
%
% "clip=" is a switch and takes no value, but the `=' must be present.
\def\psfig#1{\vbox {
	% do a zero width hard space so that a single
	% \psfig in a centering enviornment will behave nicely
	%{\setbox0=\hbox{\ }\ \hskip-\wd0}
	%
	\ps@init@parms
	\parse@ps@parms{#1}
	\compute@sizes
	%
	\ifnum\@p@scost<\@psdraft{
		%
		\special{ps::[begin] 	\@p@swidth \space \@p@sheight \space
				\@p@sbbllx \space \@p@sbblly \space
				\@p@sbburx \space \@p@sbbury \space
				startTexFig \space }
		\if@angle
			\special {ps:: \@p@sangle \space rotate \space} 
		\fi
		\if@clip{
			\if@verbose{
				\ps@typeout{(clip)}
			}\fi
			\special{ps:: doclip \space }
		}\fi
		\if@prologfile
		    \special{ps: plotfile \@prologfileval \space } \fi
		\if@decmpr{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile.Z \space }
			}\fi
			\special{ps: plotfile "`zcat \@p@sfile.Z" \space }
		}\else{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile \space }
			}\fi
			\special{ps: plotfile \@p@sfile \space }
		}\fi
		\if@postlogfile
		    \special{ps: plotfile \@postlogfileval \space } \fi
		\special{ps::[end] endTexFig \space }
		% Create the vbox to reserve the space for the figure
		\vbox to \@p@srheight true sp{
			\hbox to \@p@srwidth true sp{
				\hss
			}
		\vss
		}
	}\else{
		% draft figure, just reserve the space and print the
		% path name.
		\if@draftbox{		
			% Verbose draft: print file name in box
			\hbox{\frame{\vbox to \@p@srheight true sp{
			\vss
			\hbox to \@p@srwidth true sp{ \hss \@p@sfile \hss }
			\vss
			}}}
		}\else{
			% Non-verbose draft
			\vbox to \@p@srheight true sp{
			\vss
			\hbox to \@p@srwidth true sp{\hss}
			\vss
			}
		}\fi	



	}\fi
}}
\psfigRestoreAt



