source: Sophya/trunk/SophyaLib/Manual/sophya.tex@ 2919

\documentclass[twoside,11pt]{article}
%  Package standard : Utilisation de caracteres accentues, mode francais et graphique
\usepackage{url}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage{graphicx}
% package a mettre pour faire du pdf
\usepackage{palatino}

%  Extension de symboles mathematiques
\usepackage{amssymb}

%  Definition pour Docs Sophya
\usepackage{defsophya}

%  Constitution d'index
\usepackage{makeidx}
    
\usepackage[ps2pdf,bookmarks,bookmarksnumbered,%
              urlcolor=blue,citecolor=blue,linkcolor=blue,%
                pagecolor=blue,%hyperindex,%
                 colorlinks=true,hyperfigures=true,hyperindex=true
           ]{hyperref}

\makeindex     %  Constitution d'index

\begin{document}

\begin{titlepage}
%  The title page - top of the page with the title of the paper
\titrehp{Sophya  \\ An overview }
%  Authors list
\auteurs{
R. Ansari            &  ansari@lal.in2p3.fr       \\
E. Aubourg           &  aubourg@hep.saclay.cea.fr \\
G. Le Meur           &  lemeur@lal.in2p3.fr       \\
C. Magneville        &  cmv@hep.saclay.cea.fr     \\
S. Henrot-Versille   &  versille@in2p3.fr     
}
% \auteursall
%  The title page - bottom of the page with the paper number
\vspace{1cm}
\begin{center}
{\bf \Large Sophya Version: 1.929 (V\_Fev2006) }
% Document revision 1.0 }
\end{center}
\titrebp{1}
\end{titlepage}

\tableofcontents

\newpage

\section{Introduction}

{\bf SOPHYA} ({\bf SO}ftware for {\bf PHY}sics {\bf A}nalysis) 
is a collection of C++ classes designed for numerical and
physics analysis software development. Our goal is to provide
easy to use, yet powerful classes which can be used by scientists.
Although some of the SOPHYA modules (SkyMap, Samba, SkyT) 
have been designed with the specific goal CMB data analysis, most
modules presented here have a much broader scope and can be
used in scientific data analysis and modeling/simulation. 
Whenever possible, we use existing numerical package and libraries, 
encapsulating them whenever
possible. 
\par
\vspace*{2mm}
This documents
presents only a brief overview of the class library, 
mainly from the user's point of view. A more complete description
can be found in the reference manual, available from the SOPHYA
web site: % {\bf http://www.sophya.org}.
\href{http://www.sophya.org}{http://www.sophya.org}. 
\par
\vspace*{2mm}
The source directory tree 
\footnote{ CVS: cvsserver.lal.in2p3.fr:/exp/eros/CVSSophya} 
is organised into a number of modules. 

\begin{itemize}
\item[] {\bf BuildMgr/}  Scripts for code management, 
makefile generation and software installation
\item[] {\bf BaseTools/} General architecture support classes such 
as {\tt PPersist, NDataBlock<T>}, and few utility classes 
such as the dynamic variable list manager ({\tt DVList}) as well
as the basic set of exception classes used in SOPHYA.
\item[] {\bf TArray/} template numerical arrays, vectors and matrices \\
({\tt TArray<T> TMatrix<T> TVector<T> } \ldots) 
\item[] {\bf NTools/} Some standard numerical analysis tools 
(linear, and non linear parameter fitting, FFT, \ldots)
\item[] {\bf HiStats/}  Histogram-ming and data set handling classes (tuples) \\
({\tt Histo Histo2D NTuple DataTable} \ldots)
\item[] {\bf SkyMap/} Local and full sky maps, and some 3D geometry 
handling utility classes. \\
({\tt PixelMap<T>, LocalMap<T>, SphericalMap<T>, \ldots})
\item[] {\bf SUtils/} This module contains few utility classes, such as the
{\tt DataCard} class, as well as string manipulation functions in C and C++.
\item[] {\bf SysTools/} This module contains classes implementing 
an interface to various OS specific services, such 
threads and dynamic link/shared library handling.

\end{itemize}

The modules listed below are more tightly related to the 
CMB (Cosmic Microwave Background) data analysis problem:
\begin{itemize}
\item[] {\bf SkyT/} 
classes for spectral emission and detector frequency response modelling \\
({\tt SpectralResponse, RadSpectra, BlackBody} \ldots)
\item[] {\bf Samba/} Spherical harmonic analysis, noise generators \ldots 
\end{itemize}

The following modules contain the interface classes with 
external libraries:
\begin{itemize}
\item[] {\bf FitsIOServer/} Classes for handling file input-output 
in FITS format using the cfitsio library.
\item[] {\bf LinAlg/} Interface with Lapack linear algebra package
\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a)
\item[] {\bf XAstroPack/} Interface to some common astronomical 
computation libraries. Presently, this module uses an external library
extracted from the {\bf Xephem } source code. The corresponding source 
code is also available from SOPHYA cvs repository, module {\bf XephemAstroLib}. 
\item[] {\bf MinuitAdapt/} Wrapper classes to CERN minimization routines (Minuit). 

\end{itemize}

The following modules contain each a set of related programs using the 
SOPHYA library.
\begin{itemize}
\item[] {\bf Tests/} Simple test programs
\item[] {\bf PrgUtil/} Various utility programs (runcxx, scanppf, scanfits, \ldots)
\item[] {\bf PrgMap/} Programs performing operations on skymaps: projections, 
power spectrum in harmonic space, \ldots
\item[] {\bf PMixer/} skymixer and related programs
\end{itemize}

As a companion to SOPHYA, the {\bf (s)piapp} interactive data analysis
program is built on top of SOPHYA and the {\bf PI} GUI class library 
and application framework. The {\bf PI} ({\bf P}eida {\bf Interactive})
development started in 1995, in the EROS \footnote{EROS: {\bf E}xp\'erience 
de {\bf R}echerche d'{\bf O}bjets {\bf S}ombres - http://eros.in2p3.fr}
microlensing search collaboration, with PEIDA++ \footnote {PEIDA++: 
The EROS data analysis class library - 
http://www.lal.in2p3.fr/recherche/eros/PeidaDoc/}.
The {\bf PI} documentation and the {\bf piapp} user's guide are available
from \href{http://www.sophya.org}{http://www.sophya.org}.
%\href{http://www.sophya.org}{http://www.sophya.org}. 
The {\bf PI} is organized as the following modules:
\begin{itemize}
\item[] {\bf PI/} Portable GUI class library and application development 
framework kernel.
\item[] {\bf PIGcont/} Contour-plot drawing classes.
\item[] {\bf PIext/} Specific drawers and adapters for SOPHYA objects,
and the {\bf piapp} interactive data analysis framework.
\item[] {\bf ProgPI/} interactive analysis tool main program and pre-loaded 
modules.
\end{itemize}   

Modules containing examples and demo programs:
\begin{itemize}
\item[] {\bf Examples/} Sample SOPHYA codes and example programs and
makefiles.
\item[] {\bf DemoPIApp/} Sample scripts and programs for (s)piapp 
interactive analysis tools.
\end{itemize}   
\newpage 

\section{Using Sophya}
Basic usage of Sophya classes are described in in the following sections.
Complete Sophya documentation can be found at our web site
{\bf http://www.sophya.org}. 

\subsection{Directories, environment variables, configuration files}
The environment variable {\bf SOPHYABASE}  is used 
to define the path where the Sophya libraries and executable are installed.
\begin{itemize}
\item \$SOPHYABASE/include : Include (.h) files
\item \$SOPHYABASE/lib : Path for the archive libraries (.a)
\item \$SOPHYABASE/slb: Shared library path (.so)
\item \$SOPHYABASE/exe : Executable file path
\end{itemize}

In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
should contain the Sophya shared library path 
({\tt \$SOPHYABASE/slb }). \\
On Silocon Graphics machines (IRIX64) , the default SOPHYA configuration
correspond to the 64 bit architecture. 
The environment variable { \bf LD\_LIBRARY64\_PATH } defines
the shared library path in this case and should contain ({\tt \$SOPHYABASE/slb }.

The configure script creates links for external libraries include files in : 
\begin{itemize}
\item \$SOPHYABASE/include/FitsIO : c-fitsio library include files
\item \$SOPHYABASE/include/FFTW : FFTW library include files
\item \$SOPHYABASE/include/XAstro : XEphem library include files
\end{itemize}

The directory { \tt \$SOPHYABASE/include/SophyaConfInfo/ } contains files
describing the installed configuration of SOPHYA software.

The file { \tt \$SOPHYABASE/include/machdefs.h } contains definitions
(flags, typedef) used in SOPHYA. 

The file { \tt \$SOPHYABASE/include/sophyamake.inc } contains the 
compilation commands and flags used for building the software. 
Users can use most of compilation and link commands defined in this file:
 {\tt \$CCOMPILE , \$CXXCOMPILE . \$CXXLINK \ldots}. 
 (See module Example).

The configure script (BuildMgr/configure) creates the directory tree and the
above files.

\subsection{the runcxx program}
\index{runcxx}
{\bf runcxx} is a simple program which can be used to compile, link
and run simple C++ programs. It handles the creation of a
complete program file, containing the basic set C++ include files,
the necessary include files for SOPHYA SysTools, TArray, HiStats 
and NTools modules, and the main program with exception handling.
Other Sophya modules can be included using the {\tt -import} flag.
Use of additional include files can be specified using the 
{\tt -inc} flag.
\begin{verbatim}
csh> runcxx -h
 PIOPersist::Initialize() Starting Sophya Persistence management service 
SOPHYA Version  1.9 Revision 0 (V_Mai2005) -- May 31 2005 15:11:32 cxx 
 runcxx : compiling and running of a piece of C++ code 
   Usage: runcxx [-compopt CompileOptions] [-linkopt LinkOptions] 
                 [-tmpdir TmpDirectory] [-f C++CodeFileName] 
                 [-inc includefile] [-inc includefile ...] 
                 [-import modulename] [-import modulename ...] 
                 [-uarg UserArg1 UserArg2 ...] 
   if no file name is specified, read from standard input 
   modulenames: SkyMap, Samba, SkyT, FitsIOServer, 
                LinAlg, IFFTW, XAstroPack 
\end{verbatim}
Most examples in this manual can be tested using runcxx. The 
example below shows how to compile, link and run a sample
code.
\begin{verbatim}
//  File example.icc
Matrix a(3,3);
a = IdentityMatrix(1.);
cout << a ; 
// Executing this sample code
csh> runcxx -f example.icc
\end{verbatim}

\subsection{the scanppf program}
{\bf scanppf} is a simple SOPHYA application which can be used to check 
PPF files and list their contents. It can also provide the list of all registered
PPF handlers.
\begin{verbatim}
csh> scanppf -h
 PIOPersist::Initialize() Starting Sophya Persistence management service 
SOPHYA Version  1.9 Revision 0 (V_Mai2005) -- May 31 2005 15:11:32 cxx 
 Usage: scanppf [flags] filename 
 flags = -s -n -a0 -a1 -a2 -a3 -lh -lho 
   -s[=default} : Sequential reading of objects 
   -n : Object reading at NameTags 
   -a0...a3 : Tag List with PInPersist.AnalyseTags(0...3) 
   -lh : List PPersist handler classes 
   -lho : List PPersist handler and dataobj classes 
\end{verbatim}


\newpage

\section{Copy constructor and assignment operator}
In C++, objects can be copied by assignment or by initialisation.
Copying by initialisation corresponds to creating an object and
initialising its value through the copy constructor.
The copy constructor has its first argument as a reference, or
const reference to the object's class type. It can have  
more arguments, if default values are provided.
Copying by assignment applies to an existing object and
is performed through the assignment operator (=). 
The copy constructor implements this for identical type objects:
\begin{verbatim}
class MyObject {
public:
  MyObject();      // Default constructor
  MyObject(MyObject const & a);  // Copy constructor
  MyObject & operator = (MyObject const & a) // Assignment operator
}
\end{verbatim}
The copy constructors play an important role, as they are 
called when class objects are passed by value, 
returned by value, or thrown as an exception. 
\begin{verbatim}
// A function declaration with an argument of type MyObject,
// passed by value, and returning a MyObject
MyObject f(MyObject x) 
{
  MyObject r;
  ...
  return(r);  // Copy constructor is called here
}
// Calling the function :
MyObject a;
f(a);        // Copy constructor called for a
\end{verbatim}
It should be noted that the C++ syntax is ambiguous for the 
assignment operator. {\tt MyObject x; x=y; } and 
{\tt MyObject x=y;} have different meaning.
\begin{verbatim}
MyObject a;        // default constructor call
MyObject b(a);     // copy constructor call
MyObject bb = a;   // identical to bb(a) : copy constructor call
MyObject c;        // default constructor call
c = a;             // assignment operator call
\end{verbatim}

As a general rule in SOPHYA, objects which implements 
reference sharing on their data members have a copy constructor
which shares the data, while the assignment operator copies or
duplicate the data.

\newpage
\section{Module BaseTools}

{\bf BaseTools} contains utility classes such as 
{\tt DVlist}, an hierarchy of exception classes for Sophya, a template
class {\tcls{NDataBlock}} for handling reference counting on numerical
arrays, as well as classes providing the services for implementing simple
serialisation. 
\vspace*{5mm}

\subsection{SOPHYA persistence}
\index{PPersist} \index{PInPersist} \index{POutPersist}
\begin{figure}[hbt]
\dclsa{PPersist}
\dclsccc{PPFBinarIOStream}{PPFBinaryInputStream}{PInPersist}
\dclscc{PPFBinaryOutputStream}{POutPersist}
\caption{partial class diagram for classes handling persistence in Sophya}
\end{figure}
A simple persistence mechanism is defined in SOPHYA. Its main 
features are:
\begin{itemize}
\item[] Portable file format, containing the description of the data structures
and object hierarchy. \\
{\bf PPF} {\bf P}ortable {\bf P}ersistence file {\bf F}ormat.
\index{PPF}
\item[] Handling of read/write for multiply referenced objects.
\item[] All write operations are carried using sequential access only. This
holds also for read operations, unless positional tags are used. 
SOPHYA persistence services can thus be used to transfer objects
through network links.
\item[] The serialisation (reading/writing) for objects for a given class
is implemented through a handler object. The handler class inherits
from {\tt PPersist} class.
\item[] A run time registration mechanism is used in conjunction with 
RTTI (Run Time Type Identification) for identifying handler classes
when reading {\bf PInPersist} streams, or for associating handlers
with data objects {\bf AnyDataObject} for write operations.
\end{itemize}
A complete description of SOPHYA persistence mechanism and guidelines
for writing delegate classes for handling object persistence is beyond
the scope of this document. The most useful methods for using Sophya
persistence are listed below:
\begin{itemize}
\item[] {\tt POutPersist::PutObject(AnyDataObj \& o)} \\ 
Writes the data object {\bf o} to the output stream.
\item[] {\tt POutPersist::PutObject(AnyDataObj \& o, string tagname)} \\
Writes the data object {\bf o} to the output stream, associated with an
identification tag {\bf tagname}.
\item[] {\tt PInPersist::GetObject(AnyDataObj \& o)} \\ 
Reads the next object in stream into {\bf o}. An exception is 
generated for incompatible object types.
\item[] {\tt PInPersist::GetObject(AnyDataObj \& o, string  tagname)} \\
Reads the object associated with the tag {\bf tagname} into {\bf o}. 
An exception is generated for incompatible object types.
\end{itemize}
The operators {\tt operator << (POutPersist ...) } and 
{\tt operator >> (PInPersist ...) } are often overloaded 
to perform {\tt PutObject()} and {\tt GetObject()} operations.
the {\bf PPFNameTag} (ppfnametag.h) class can be used in conjunction with 
{\tt << >> } operators to write objects with a name tag or to retrieve
an object identified with a name tag. The example below shows the
usage of these operators:
\begin{verbatim}
// Creating and filling a histogram
Histo hw(0.,10.,100);
...
// Writing histogram to a PPF stream
POutPersist os("hw.ppf");
os << PPFNameTag("myhisto") << hw;

// Reading a histogram from a PPF stream
PInPersist is("hr.ppf");
is >> PPFNameTag("myhisto")  >> hr;
\end{verbatim}

The {\bf scanppf} program can be used to list the content of a PPF file.

\subsection{\tcls{NDataBlock}}
\index{\tcls{NDataBlock}}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{\tcls{NDataBlock}}
\dclsbb{PPersist}{\tcls{FIO\_NDataBlock}}
\end{figure}
The {\bf \tcls{NDataBlock}} is designed to handle reference counting 
and sharing of memory blocs (contiguous arrays) for numerical data 
types. Initialisation, resizing, basic arithmetic operations, as
well as persistence handling services are provided. 
The persistence handler class ({\tt \tcls{FIO\_NDataBlock}}) insures
that a single copy of data is written for multiply referenced objects,
and the data is shared among objects when reading. 
\par
The example below shows writing of NDataBlock objects through the 
use of overloaded operator $ << $ :
\begin{verbatim}
#include "fiondblock.h"
// ...
POutPersist pos("aa.ppf");
NDataBlock<r_4> rdb(40);
rdb = 567.89;
pos << rdb;
// We can also use the PutObject method
NDataBlock<int_4> idb(20);
idb = 123;
pos.PutObject(idb);
\end{verbatim}
The following sample programs show the reading of the created PPF file : 
\begin{verbatim}
PInPersist pis("aa.ppf");
NDataBlock<r_4> rdb;
pis >> rdb;
cout << rdb;
NDataBlock<int_4> idb;
cout << idb;
\end{verbatim}

\subsection{Using DVList}
\index{DVList} \index{MuTyV}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{DVList}
\dclsbb{PPersist}{\tclsc{ObjFileIO}{DVList}}
\end{figure}
The {\bf DVList} class objects can be used to create and manage list
of values, associated with names. A list of pairs of (MuTyV, name(string))
is maintained by DVList objects. {\bf MuTyV} is a simple class
capable of holding string, integer, float or complex values, 
providing easy conversion methods between these objects.
\begin{verbatim}
// Using MuTyV objects
MuTyV s("hello");  // string type value
MuTyV x;
x = "3.14159626";    // string type value, ASCII representation for Pi
double d = x;        // x converted to double = 3.141596
x = 314;             // x contains the integer value = 314
// Using DVList
DVList  dvl;
dvl("Pi") = 3.14159626;            // float value, named Pi
dvl("Log2") = 0.30102999;          // float value, named Log2
dvl("FileName") = "myfile.fits";   // string value, named myfile.fits
// Printing DVList object
cout << dvl;
\end{verbatim}

\subsection{\tcls{SegDataBlock} , \tcls{SwSegDataBlock}}
\begin{figure}[hbt]
\dclsccc{AnyDataObj}{\tcls{SegDBInterface}}{ \tcls{SegDataBlock} }
\dclscc{\tcls{SegDBInterface}}{ \tcls{SwSegDataBlock} }
\end{figure}
\begin{itemize} 
\item[] \tcls{SegDataBlock}  handles arrays of object of
type {\bf T} with reference sharing in memory. The array can be extended
(increase in array size) with fixed segment size. It implements the interface
defined by tcls{SegDBInterface}.
\item[] \tcls{SwSegDataBlock} Implements the same tcls{SegDBInterface}
using a data swapper (\tcls{DataSwapperInterface} . 
Can thus be used for very large objects.
\end{itemize}
 
\newpage
\section{Module TArray}
\index{\tcls{TArray}}
{\bf TArray} module contains template classes for handling standard 
operations on numerical arrays. Using the class {\tt \tcls{TArray} },
it is possible to create and manipulate up to 5-dimension numerical
arrays {\tt (int, float, double, complex, \ldots)}. The include
file {\tt array.h} declares all the classes and definitions 
in module TArray. {\bf Array} is a typedef for arrays 
with double precision floating value elements. \\
{\tt typedef TArray$<$r\_8$>$ Array ; }

\begin{figure}[hbt]
\dclsccc{AnyDataObj}{BaseArray}{\tcls{TArray}}
\dclsbb{PPersist}{\tcls{FIO\_TArray}}
\end{figure}

The development of this module started around 1999-2000,
after evaluation of a number of publicly  available 
C++ array hadling packages, including TNT, Lapack++, Blitz++, 
as well as commercial packages from RogueWave (math.h++ \ldots). 
Most of these packages provide interesting functionalities, however, 
not any one package seemed to fulfill most of our requirements.
\begin{itemize}
\item Capability to handle {\bf large - multidimensional - dense} 
arrays, for numerical data types. Although we have used templates, for
data type specialisation, the actual code, apart inline functions is
not in header files. Instead, we use explicit instanciation, and the 
compiled code for the various numerical types of arrays is the
library .
\item The shape and size of the arrays can be defined and changed
at run time. The classes ensure the memory management of the 
created objets, with reference sharing for the array data. 
The default behaviour of the copy constructor is to share the data,
avoiding expensive memory copies.
\item The package provides transparent management of sub-arrays
and slices, in an intuitive way, somehow similar to what is
available in Mathlab or Scilab.
\item The memory organisation for arrays, specially matrices 
(row-major or column major) can be 
controled. This provide compatibility when using existing C or
Fortran coded numerical libraries.
\item The classes provide efficient methods to perform basic arithmetic
and mathematical operations on arrays. In addition, operator overload 
provides intuitive programming for element acces and most basic
arithmetic operations.
\item Conversion can be performed between arrays with different
data types. Copy and arithmetic operations can be done transparently 
between arrays with different memory organisation patterns.   
\item This module does not provide more complex operations 
such as FFT or linear algebra. Additional libraries are used, with interface
classes for these operations.
\item ASCII formatted I/O, for printing and read/write operations to/from text files.
\item Efficient binary I/O for object persistence (PPF format), or import/export 
to other data formats, such as FITS are provided by helper or handler classes. 
\end{itemize}

\subsection{Using arrays}
\index{Sequence} \index{RandomSequence} \index{RegularSequence} 
\index{EnumeratedSequence}
The example below shows basic usage of arrays, creation, initialisation
and arithmetic operations. Different kind of {\bf Sequence} objects
can be used for initialising arrays. 

\begin{figure}[hbt]
\dclsbb{Sequence}{RandomSequence}
\dclsb{RegularSequence}
\dclsb{EnumeratedSequence}
\end{figure}

The example below shows basic usage of arrays: 
\index{\tcls{TArray}} 
\begin{verbatim}
// Creating and initialising a 1-D array of integers
TArray<int> ia(5);
EnumeratedSequence es;
es = 24, 35, 46, 57, 68;
ia = es;
cout << "Array<int> ia = " << ia;
// 2-D array of floats 
TArray<r_4> b(6,4), c(6,4);
// Initializing b with a constant
b = 2.71828;
// Filling c with random numbers  
c = RandomSequence(); 
// Arithmetic operations
TArray<r_4> d = b+0.3f*c;
cout << "Array<float> d = " << d;
\end{verbatim}

The copy constructor shares the array data, while the assignment operator
copies the array elements, as illustrated in the following example:
\begin{verbatim}
TArray<int> a1(4,3);
a1 = RegularSequence(0,2);
// Array a2 and a1 shares their data 
TArray<int> a2(a1);
// a3 and a1 have the same size and identical elements
TArray<int> a3;
a3 = a1;
// Changing one of the a2 elements
a2(1,1,0) = 555;
// a1(1,1) is also changed to 555, but not a3(1,1)
cout << "Array<int> a1 = " << a1;
cout << "Array<int> a3 = " << a3;
\end{verbatim} 

\subsection{Matrices and vectors}
\index{\tcls{TMatrix}}  \index{\tcls{TVector}}
\begin{figure}[hbt]
\dclsccc{\tcls{TArray}}{\tcls{TMatrix}}{\tcls{TVector}}
\end{figure}
Vectors and matrices are 2 dimensional arrays. The array size
along one dimension is equal 1 for vectors. Column vectors
have {\tt NCols() = 1} and row vectors have {\tt NRows() = 1}.
Mathematical expressions involving matrices and vectors can easily 
be translated into C++ code using {\tt TMatrix} and 
{\tt TVector} objects. {\bf Matrix} and {\bf Vector} are 
typedefs for double precision float matrices and vectors.
The operator {\bf *} beteween matrices is redefined to 
perform matrix multiplication. One can then write: \\
\begin{verbatim}
  // We create a row vector
  Vector v(1000, BaseArray::RowVector);
  // Initialize values with a random sequence
  v = RandomSequence();
  // Compute the vector length (norm)
  double norm = (v*v.Transpose()).toScalar();
  cout << "Norm(v) = " << norm << endl;
\end{verbatim}

This module contains basic array and matrix operations 
such as the Gauss matrix inversion algorithm
which can be used to solve linear systems, as illustrated by the 
example below:
\begin{verbatim}
#include "sopemtx.h"
// ...
// Creation of a random 5x5 matrix
Matrix A(5,5);
A = RandomSequence(RandomSequence::Flat);
Vector X0(5);
X0 = RandomSequence(RandomSequence::Gaussian);
// Computing B = A*X0
Vector B = A*X0;
// Solving the system A*X = B 
Vector X;
LinSolve(A, B, X);
// Checking the result
Vector diff = X-X0;
cout << "X-X0= " << diff ;
double min,max;
diff.MinMax(min, max);
cout << " Min(X-X0) = " << min << " Max(X-X0) = " << max << endl; 
\end{verbatim}

{\bf Warning: } The operations defined in {\tt sopemtx.h}, such as 
matrix inversion and linear system solver use a basic Gauss pivot
algorithm which are not adapted for large matrices ($>\sim 100x100$).
The services provided in other modules, such as {\bf LinAlg} should
be preferred in such cases.

\subsection{Working with sub-arrays and Ranges}
\index{Range}
A powerful mechanism is included in array classes for working with 
sub-arrays. The class {\bf Range} can be used to specify range of array 
indexes in any of the array dimensions. Any regularly spaced index
range can be specified, using the {\tt start} and {\tt end} index
and an optional step (or stride). It is also possible to specify
the {\tt start} index and the number of elements.
\begin{itemize}
\item {\bf Range::all()} {\tt = Range(Range::firstIndex(), Range::lastIndex())}  \\
return a Range objects representing all  valid indexes along the 
corresponding axe.
\item {\bf Range::first()} {\tt = Range(Range::firstIndex())}  \\
return a Range object representing the first valid index
\item {\bf Range::last()} {\tt = Range(Range::lastIndex())} 
return a Range object representing the last valid index
\item {\bf Range(idx)} represents a single index ({\bf = idx})
\item {\bf Range(first, last)} represents the range of indices
{\bf first} $\leq$ index $\leq$ {\bf last}.
The static method {\tt Range::lastIndex()} can be used 
to specify the last valid index.
\item {\bf Range(first, last, step)} represents the range of index
which is equivalent to \\  {\tt for(index=first; index <= last; index += step) }
\item { \bf Range (first, last, size, step) } the general form can be used 
to specify an index range, using the number of elements. 
It is possible to specify a range of index, ending with the last valid index. 
For example \\
\hspace*{5mm} 
{\tt Range(Range::lastIndex(), Range::lastIndex(), 3, 2) }  \\
defines  the index range: \hspace*{5mm}  last-4, last-2, last.

\begin{center}
\begin{tabular}{ll}
\multicolumn{2}{c}{ {\bf Range} {\tt (start, end, size, step) } } \\[2mm]
\hline \\
{\bf Range} {\tt r(3,6); } &  index range:  \hspace{2mm} 3,4,5,6 \\
{\bf Range} {\tt r(3,6,0,1); } &  index range: \hspace{2mm} 3,4,5,6 \\
{\bf Range} {\tt r(7,0,3,1); } &  index range: \hspace{2mm} 7,8,9 \\
{\bf Range} {\tt r(10,0,5,2); } &  index range: \hspace{2mm} 10,12,14,16,18 \\
\end{tabular}
\end{center}
\end{itemize}

The method {\tt TArray<T>SubArray(Range ...)} can be used 
to extract subarrays and slices. The operator {\tt operator() (Range rx, Range ry ...)}
is also overloaded for sub-array extraction.
For matrices, {\tt TMatrix<T>::Row()} and {\tt TMatrix<T>::Column()} 
extract selected matrix rows and columns.

The example illustrates the sub-array extraction using Range objects:
\begin{verbatim}
  // Creating and initialising a 2-D (6 x 4) array of integers
  TArray<int> iaa(6, 4);
  iaa = RegularSequence(1,2);
  cout << "Array<int> iaa = \n" << iaa;
  // We extract a sub-array - data is shared with iaa
  TArray<int> iae = iaa(Range(1, Range::lastIndex(), 3) ,
                        Range::all(), Range::first() );
  cout << "Array<int> iae=subarray(iaa) = \n" << iae;
  // Changing iae elements changes corresponding iaa elements
  iae = 0;
  cout << "Array<int> iae=0 --> iaa = \n" << iaa;

\end{verbatim}

In the following example, a simple low-pass filter, on a one 
dimensional stream (Vector) has been written using sub-arrays:  

\begin{verbatim}
// Input Vector containing a noisy periodic signal
  Vector in(1024), out(1024);
  in = RandomSequence(RandomSequence::Gaussian, 0., 1.);
  for(int kk=0; kk<in.Size(); kk++)
    in(kk) += 2*sin(kk*0.05);
// Compute the output vector by a simple low pass filter
  Vector out(1024);
  int w = 2;
  for(int k=w; k<in.Size()-w; k++)
    out(k) = in(Range(k-w, k+w).Sum()/(2.*w+1.);
\end{verbatim}

\subsection{Input, Output}
Arrays can easily be saved to, or restored from files in different formats.
SOPHYA library can handle array I/O to ASCII formatted files, to PPF streams, 
as well as to files in FITS format.
FITS format input/output is provided through the classes in 
{\bf FitsIOServer} module. Onnly arrays with data types 
supported by the FITS standard can be handled during 
I/O operations to and from FITS streams (See the FitsIOServer section
for additional details).

\subsubsection{PPF streams} 

SOPHYA persistence (PPF streams) handles reference sharing, and multiply 
referenced objects are only written once. A hierarchy of arrays and sub-arrays
written to a PPF stream is thus completely recovered, when the stream is read.
The following example illustrates this point:
\begin{verbatim}
{
// Saving an array with a sub-array into a POutPersist file
Matrix A(3,4);
A = RegularSequence(10,5);
// Create a sub-array of A 
Matrix AS = A(Range(1,2), Range(2,3));
// Save the two arrays to a PPF stream
POutPersist pos("aas.ppf");
pos << A << AS; 
}
{
// Reading arrays from the previously created PPF file aas.ppf
PInPersist pis("aas.ppf");
Matrix B,BS;
pis >> B >> BS;
// BS is a sub-array of B, modifying BS changes also B
BS(1,1) = 98765.;
cout << " B , BS after BS(1,1) = 98765. " 
     << B << BS << endl; 
}
\end{verbatim} 
The execution of this sample code creates the file {\tt aas.ppf} and 
its output is reproduced here. Notice that the array hierarchy is 
recovered. BS is a sub-array of B, and modifying BS changes also
the corresponding element in B.
\begin{verbatim}
 B , BS after BS(1,1) = 98765.

--- TMatrix<double>(NRows=3, NCols=4) ND=2 SizeX*Y*...= 4x3 ---
10 15 20 25
30 35 40 45
50 55 60 98765

--- TMatrix<double>(NRows=2, NCols=2) ND=2 SizeX*Y*...= 2x2 ---
40 45
60 98765
\end{verbatim} 

\centerline{\bf Warning: }

There is a drawback in this behaviour: only a single 
copy of an array is written to a file, even if the array is modified, 
without being resized and written to a PPF stream. 
\begin{verbatim}
{
POutPersist pos("mca.ppf");
TArray<int_4> ia(5,3);
ia = 8;
pos << ia;
ia = 16;
pos << ia;
ia = 32;
pos << ia;
}
\end{verbatim}

Only a single copy of the data is effectively written to the output 
PPF file, corresponding to the value 8 for array elements. When we 
read the three array from the file mca.ppf, the same array elements
are obtained three times (all elements equal to 8):
\begin{verbatim}
{
PInPersist pis("mca.ppf");
TArray<int_4> ib;
pis >> ib;
cout << " First array read from mca.ppf : " << ib;
pis >> ib;
cout << " Second array read from mca.ppf : " << ib;
pis >> ib;
cout << " Third array read from mca.ppf : " << ib;
}
\end{verbatim}

\subsubsection{ASCII streams} 

The {\bf WriteASCII} method can be used to dump an array to an ASCII 
formatted file, while the {\bf ReadASCII} method can be used to decode 
ASCII formatted files. Space or tabs are the possible separators. 
Complex numbers should be specified as a pair of comma separated 
real and imaginary parts, enclosed in parenthesis. 

\begin{verbatim}
{
// Creating array A and writing it to an ASCII file (aaa.txt)
Array A(4,6);
A = RegularSequence(0.5, 0.2);
ofstream ofs("aaa.txt");
A.WriteASCII(ofs);
}
{
// Decoding the ASCII file aaa.txt
ifstream ifs("aaa.txt");
Array B;
sa_size_t nr, nc;
B.ReadASCII(ifs,nr,nc);
cout << " Array B; B.ReadASCII() from file " << endl;
cout << B ;
}
\end{verbatim}


\subsection{Complex arrays}
The {\bf TArray} module provides few functions for manipulating
arrays of complex numbers (single and double precision).
These functions are declared in {\tt matharr.h}.
\begin{itemize}
\item[\bul] Creating a complex array through the specification of the 
real and imaginary parts. 
\item[\bul] Functions returning arrays corresponding to real and imaginary 
parts of a complex array: {\tt real(za) , imag(za) } 
({\bf Warning:} Note that the present implementation does not provide
shared memory access to real and imaginary parts.)
\item[\bul] Functions returning arrays corresponding to the module, 
phase, and module squared of a complex array:
 {\tt phase(za) , module(za) , module2(za) } 
\end{itemize}

\begin{verbatim}
  TVector<r_4> p_real(10, BaseArray::RowVector);
  TVector<r_4> p_imag(10, BaseArray::RowVector);
  p_real = RegularSequence(0., 0.5);
  p_imag = RegularSequence(0., 0.25);
  TVector< complex<r_4> > zvec = ComplexArray(p_real, p_imag);
  cout << " :: zvec= " << zvec;
  cout << " :: real(zvec) = " << real(zvec) ;
  cout << " :::: imag(zvec) = " << imag(zvec) ;
  cout << " :::: module2(zvec) = " << module2(zvec) ;
  cout << " :::: module(zvec) = " << module(zvec) ;
  cout << " :::: phase(zvec) = " << phase(zvec) ;
\end{verbatim}

The decoding of complex numbers from an ASCII formatted stream 
is illustrated by the next example. As mentionned already, 
complex numbers should be specified as a pair of comma separated 
real and imaginary parts, enclosed in parenthesis. 

\begin{verbatim}
csh> cat zzz.txt
(1.,-1) (2., 2.5) -3. 12.
-24. (-6.,7.) 14.2 (8.,64.)

//  Decoding of complex numbers from an ASCII file
//  Notice that the << operator can be used instead of ReadASCII
TArray< complex<r_4> > Z;
ifstream ifs("zzz.txt");
ifs >> Z;
cout << " TArray< complex<r_4> > Z from file zzz.txt " << Z ;
\end{verbatim}


\subsection{Memory organisation}
{\tt \tcls{TArray} } can handle numerical arrays with various memory 
organisation, as long as the spacing (steps) along each axis is 
regular. The five axis are labeled X,Y,Z,T,U. The examples below
illustrates the memory location for a 2-dimensional, $N_x=4 \times N_y=3$.
The first index is along the X axis and the second index along the Y axis.
\begin{verbatim}
    | (0,0)  (0,1)  (0,2)  (0,3) |
    | (1,0)  (1,1)  (1,2)  (1,3) |
    | (2,0)  (2,1)  (2,2)  (2,3) |
\end{verbatim}
In the first case, the array is completely packed 
($Step_X=1, Step_Y=N_X=4$), with zero offset,
while in the second case, $Step_X=2, Step_Y=10, Offset=10$:
\begin{verbatim}
    | 0  1  2  3  |         | 10 12 14 16 |
Ex1 | 4  5  6  7  |     Ex2 | 20 22 24 26 |
    | 8  9  10 11 |         | 30 32 34 36 | 
\end{verbatim} 

For matrices and vectors, an optional argument ({\tt MemoryMapping}) 
can be used to select the memory mapping, where two basic schemes 
are available: \\
{\tt CMemoryMapping} and {\tt FortranMemoryMapping}. \\
In the case where {\tt CMemoryMapping} is used, a given matrix line
is packed in memory, while the columns are packed when 
{\tt FortranMemoryMapping} is used. The first index when addressing 
the matrix elements (line number index) runs along
the Y-axis if  {\tt CMemoryMapping} is used, and along the X-axis
in the case of {\tt FortranMemoryMapping}. 
Arithmetic operations between matrices
with different memory organisation is allowed as long as 
the two matrices have the same sizes (Number of rows and columns). 
The following code example and the corresponding output illustrates 
these two memory mappings. The {\tt \tcls{TMatrix}::TransposeSelf() }
method changes effectively the matrix memory mapping, which is also 
the case of {\tt \tcls{TMatrix}::Transpose() } method without argument.

\begin{verbatim}
TArray<r_4> X(4,2); 
X = RegularSequence(1,1);
cout << "Array X= " << X << endl;
TMatrix<r_4> X_C(X, true, BaseArray::CMemoryMapping);
cout << "Matrix X_C (CMemoryMapping) = " << X_C << endl;
TMatrix<r_4> X_F(X, true, BaseArray::FortranMemoryMapping);
cout << "Matrix X_F (FortranMemoryMapping) = " << X_F << endl;
\end{verbatim}
This code would produce the following output (X\_F = Transpose(X\_C)) :
\begin{verbatim}
Array X= 
--- TArray<f>  ND=2 SizeX*Y*...= 4x2 ---
1, 2, 3, 4
5, 6, 7, 8

Matrix X_C (CMemoryMapping) = 
--- TMatrix<f>(NRows=2, NCols=4) ND=2 SizeX*Y*...= 4x2 ---
1, 2, 3, 4
5, 6, 7, 8

Matrix X_F (FortranMemoryMapping) = 
--- TMatrix<f>(NRows=4, NCols=2) ND=2 SizeX*Y*...= 4x2 ---
1, 5
2, 6
3, 7
4, 8
\end{verbatim}

\newpage 

\section{Module HiStats}
\begin{figure}[hbt]
\dclsccc{AnyDataObj}{Histo}{HProf}
\dclsbb{AnyDataObj}{Histo2D}
\dclsbb{AnyDataObj}{Ntuple}
\caption{partial class diagram for histograms and ntuples}
\end{figure}

{\bf HiStats} contains classes for creating, filling, printing and
doing various operations on one or two dimensional histograms
{\tt Histo} and {\tt Histo2D} as well as profile histograms {\tt HProf}. \\
This module also contains {\tt NTuple} and {\tt XNTuple} which are
more or less the same that the binary FITS tables.

\subsection{1D Histograms}
\index{Histo}
For 1D histograms, various numerical methods are provided such as
computing means and sigmas, finding maxima, fitting, rebinning,
integrating \dots \\
The example below shows creating and filling a one dimensional histogram
of 100 bins from $-5.$ to $+5.$ to create a Gaussian normal distribution
with errors~:
\begin{verbatim}
#include "histos.h"
// ...
Histo H(-0.5,0.5,100);
H.Errors();
for(int i=0;i<25000;i++) {
  double x = NorRand();
  H.Add(x);
}
H.Print(80);
\end{verbatim}

\subsection{2D Histograms}
\index{Histo2D}
Much of these operations are also valid for 2D histograms. 1D projection
or slices can be set~:
\begin{verbatim}
#include "histos2.h"
// ...
Histo2D H2(-1.,1.,100,0.,60.,50);
H2.SetProjX(); // create the 1D histo for X projection
H2.SetBandX(25.,35.); // create 1D histo  projection for 25.<y<35.
H2.SetBandX(35.,45.); // create 1D histo  projection for 35.<y<45.
H2.SetBandX(40.,55.); // create 1D histo  projection for 40.<y<55.
//... fill H2 with what ever you want
H2.Print();
Histo *hx = H2.HProjX();
  hx->Print(80);
Histo *hbx2 = HBandX(1); // Get the second X band (35.<y<45.)
  hbx2->Print(80);
\end{verbatim}

\subsection{Profile Histograms}
\index{HProf}
Profiles histograms {\bf HProf} contains the mean and the 
sigma of the distribution
of the values filled in each bin. The sigma can be changed to
the error on the mean. When filled, the profile histogram looks
like a 1D histogram and much of the operations that can be done on 1D histo
may be applied onto profile histograms.

\subsection{Data tables (tuples)}
\index{NTuple}
NTuple are memory resident tables of 32 or 64 bits floating values 
(float/double).They are arranged in columns. Each line is often called an event.
These objects are frequently used to analyze data.
The piapp graphicals tools can plot a column against an other one
with respect to various selection cuts. \\
Here is an example of creation and filling~:
\begin{verbatim}
#include "ntuple.h"
#include "srandgen.h"
// ...
char* nament[4] = {"i","x","y","ey"};
r_4 xnt[4];
NTuple NT(4,nament);
for(i=0;i<5000;i++) {
  xnt[0] = i+1;
  xnt[1] = 5.*drandpm1();       // a random value between -5 and +5
  xnt[2] = 100.*exp(-0.5*xnt[1]*xnt[1]) + 1.;
  xnt[3] = sqrt(xnt[2]);
  xnt[2] += xnt[3] * NorRand(); // add a random gaussian error
  NT.Fill(xnt);
}
\end{verbatim}

XNTuple provide additional functionalities, compared to NTuple. 
They are deprecated and are only kept  for backward compatibility
and should not be used anymore. Use DataTable and 
SwPPFDataTable instead. 
Object of type XNTuple handle various types
of column values (double,float,int,string,...) and can handle
very large data sets, through swap space on disk. 

\index{DataTable} 

The class {\bf DataTable} extends significantly the functionalities provided by 
NTuple. DataTable is a memory resident implementation of the interface 
{\bf BaseDataTable } which organizes the data as a 2-D table. User can define
the name and data type of each column. Data is added to the table as rows.
The table is extended as necessary when adding rows. 
The sample code below shows an example of DataTable usage :
\begin{verbatim}
   #include "datatable.h"
   // ...
   DataTable dt(64);
   dt.AddFloatColumn("X0_f");
   dt.AddFloatColumn("X1_f");
   dt.AddDoubleColumn("X0X0pX1X1_d");
   double x[5];
   for(int i=0; i<63; i++) {
     x[0] = (i/9)-4.;  x[1] = (i/9)-3.;  x[2] = x[0]*x[0]+x[1]*x[1];
     dt.AddLine(x);
   }
   // Printing table info
   cout << dt ;
   // Saving object into a PPF file
   POutPersist po("dtable.ppf");
   po << dt ;

\end{verbatim}

\begin{figure}[hbt]
\dclsccc{AnyDataObj}{BaseDataTable}{DataTable}
\dclscc{BaseDataTable}{SwPPFDataTable}
\end{figure}

\index{SwPPFDataTable}
The class {\bf SwPPFDataTable} implements the BaseDataTable interface
using segmented data blocks with swap on PPF streams. Very large data sets
can be created and manipulated through tis class

\subsection{Writing, viewing \dots }

All these objects have been design to be written to or read from a persistent file.
The following example shows how to write the previously created objects
into such a file~:
\begin{verbatim}
//-- Writing
{
char *fileout = "myfile.ppf";
string tag;
POutPersist outppf(fileout);
tag = "H";   outppf.PutObject(H,tag);
tag = "H2";   outppf.PutObject(H2,tag);
tag = "NT";   outppf.PutObject(NT,tag);
}  // closing ``}'' destroy ``outppf'' and automatically close the file !
\end{verbatim}

Sophya graphical tools (spiapp) can automatically display and operate
all these objects.

\newpage
\section{Module NTools}

This module provides elementary numerical tools for numerical integration,
fitting, sorting and ODE solving. FFTs are also provided (Mayer,FFTPack).

\subsection{Fitting}
\index{Fitting} \index{Minimisation}
Fitting is done with two classes {\tt GeneralFit} and {\tt GeneralFitData}
and is based on the Levenberg-Marquardt method.
\index{GeneralFit} \index{GeneralFitData}
GeneralFitData is a class which provide a description of the data
to be fitted. GeneralFit is the fitter class. Parametrized functions
can be given as classes which inherit {\tt GeneralFunction}
or as simple C functions. Classes of pre-defined functions are provided
(see files fct1dfit.h and fct2dfit.h). The user interface is very close
from that of the CERN {\tt Minuit} fitter.
Number of objects (Histo, HProf \dots ) are interfaced with GeneralFit
and can be easily fitted. \\
Here is a very simple example for fitting the previously created NTuple
with a Gaussian~:
\begin{verbatim}
#include "fct1dfit.h"
// ...

// Read from ppf file
NTuple nt;
{
PInPersist pis("myfile.ppf");
string tag = "NT";  pis.GetObject(nt,tag);
}

// Fill GeneralData
GeneralData mGdata(nt.NEntry());
for(int i=0; i<nt.NEntry(); i++)
  mGdata.AddData1(xnt[1],xnt[2],xnt[3]); // Fill x, y and error on y
mGData.PrintStatus();

// Function for fitting : y = f(x) + noise
Gauss1DPol mFunction; // gaussian + constant

// Prepare for fit
GeneralFit mFit(&mFunction);  // create a fitter for the choosen function
mFit.SetData(&mGData);        // connect data to the fitter

// Set and initialise the parameters (that's non-linear fitting!)
//           (num par, name, guess start, step, [limits min and max])
mFit.SetParam(0,"high",90.,1..);
mFit.SetParam(1,"xcenter",0.05,0.01);
mFit.SetParam(2,"sigma",sig,0.05,0.01,10.);
                              // Give limits to avoid division by zero
mFit.SetParam(3,"constant",0.,1.);

// Fit and print result
int rcfit = mFit.Fit();
mFit.PrintFit();
if(rcfit>0) {)
  cout<<"Reduce_Chisquare = "<<mFit.GetChi2Red()
      <<" nstep="<<mFit.GetNStep()<<" rc="<<rcfit<<endl;
} else {
  cout<<"Fit_Error, rc = "<<rcfit<<"  nstep="<<mFit.GetNStep()<<endl;
  mFit.PrintFitErr(rcfit);
}

// Get the result for further use
TVector<r_8> ParResult = mFit.GetParm();
cout<<ParResult;
\end{verbatim}

Much more usefull possibilities and detailed informations might be found
in the HTML pages of the Sophya manual.

\subsection{Polynomial}
\index{Polynomial} \index{Poly} \index{Poly2}
Polynomials of 1 or 2 variables are supported ({\tt Poly} and {\tt Poly2}).
Various operations are supported~:
\begin{itemize}
\item elementary operations between polynomials $(+,-,*,/) $
\item setting or getting coefficients
\item computing the value of the polynomial for a given value
      of the variable(s),
\item derivating
\item computing roots (degre 1 or 2)
\item fitting the polynomial to vectors of data.
\end{itemize}
Here is an example of polynomial fitting~:
\begin{verbatim}
#include "poly.h"
// ...
Poly pol(2);
pol[0] = 100.; pol[1] = 0.; pol[2] = 0.01; // Setting coefficients
TVector<r_8> x(100);
TVector<r_8> y(100);
TVector<r_8> ey(100);
for(int i=0;i<100;i++) {
  x(i)  = i;
  ey(i) = 10.;
  y(i)  = pol((double) i) + ey(i)*NorRand();
  ey(i) *= ey(i)
}

TVector<r_8> errcoef;
Poly polfit;
polfit.Fit(x,y,ey,2,errcoef);

cout<<"Fit Result"<<polfit<<endl;
cout<<"Errors :"<<errcoef;
\end{verbatim}

Similar operations can be done on polynomials with 2 variables.

\subsection{Integration, Differential equations}
\index{Integration}
The NTools module provide also simple classes for numerical integration 
of functions and differential equations.
\begin{figure}[hbt]
\dclsbb{Integrator}{GLInteg}
\dclsb{TrpzInteg}
\end{figure}

\index{GLInteg} \index{TrpzInteg}
{\bf GLInteg} implements the integration through Gauss-Legendre method
and {\bf TrpzInteg} implements trapeze integration. For {\bf TrpzInteg},
number of steps specify the number of trapeze, and integration step, 
their width.
The sample code below illustrates the use of TrpzInteg class:
\begin{verbatim}
#include "integ.h"
// ......................................................
// Function to be integrated
double myf(double x)
{
// Simple a x + b x^2  (a=2 b=3)
return (x*(2.+3.*x));
}
// ......................................................

// Compute Integral(myf, 2., 5.) between xmin=2., xmax=5.
TrpzInteg trpz(myf, 2., 5.);
// We specify an integration step
trpz.DX(0.01);
// The integral can be computed as trpz.Value() 
double myf_integral = trpz.Value();
// We could have used the cast operator :
cout << "Integral[myf, 2., 5.]= " << (double)trpz << endl;
// Limits can be specified through ValueBetween() method 
cout << "Integral[myf, 0., 4.]= " << trpz.ValueBetween(0.,4.) << endl;
\end{verbatim}

\subsection{Fourier transform (FFT)}
\index{FFT} \index{FFTPackServer}
An abstract interface for performing FFT operations is defined by the 
{\bf FFTServerInterface} class. The {\bf FFTPackSever} class implements
one dimensional FFT, on real and complex data. FFTPackServer uses an
adapted and extended version of FFTPack (available from netlib), 
translated in C, and can operate on single and double precision
({\tt float, double}) data.

The sample code below illustrates the use of FFTServers:
\begin{verbatim}
#include "fftpserver.h"
   // ...
TVector<r_8> in(32);
TVector< complex<r_8> > out;
in = RandomSequence();
FFTPackServer ffts;
ffts.setNormalize(true);  // To have normalized transforms
cout << " FFTServer info string= " << ffts.getInfo() << endl;
cout << "in= " << in << endl;
cout << " Calling ffts.FFTForward(in, out) : " << endl;
ffts.FFTForward(in, out);
cout << "out= " << out << endl;
\end{verbatim}

% \newpage
\section{Module SUtils}
Some utility classes and C/C++ string manipulation functions are gathered
in {\bf SUtils} module.
\subsection{Using DataCards}
\index{DataCards}
The {\bf DataCards} class can be used to read parameters from a file.
Each line in the file starting with \@ defines a set of values
associated with a keyword. In the example below, we read the 
parameters corresponding with the keyword {\tt SIZE} from the 
file {\tt ex.d}. We suppose that {\tt ex.d} contains the line: \\
{\tt @SIZE 400 250} \\
\begin{verbatim}
#include "datacards.h"
// ...
// Initialising DataCards object dc from file ex.d
DataCards dc( "ex.d" );
// Getting the first and second parameters for keyword size
// We define a default value 100
int size_x = dc.IParam("SIZE", 0, 100);
int size_y = dc.IParam("SIZE", 1, 100);
cout << " size_x= " << size_x << " size_y= " << size_y << endl;
\end{verbatim}
 
\section{Module SysTools}
The {\bf SysTools} module contains classes implementing interface to some 
OS specific services. The class {\bf ResourceUsage} \index{ResourceUsage}
and {\bf Timer} {\index{Timer} provides access to information 
about various resource usage (memory, CPU, ...).
The class {\bf Periodic} provides the necessary services needed to 
implement the execution of a periodic action.

\subsection{Thread management classes}
A basic interface to POSIX threads \index{thread} is also provided
through the \index{ZThread} {\bf ZThread}, {\bf ZMutex} and {\bf ZSync}
classes.
 
\subsection{Dynamic linker and C++ compiler classes}
\index{PDynLinkMgr}
The class {\bf PDynLinkMgr} can be used for managing shared libraries
at run time. The example below shows the run time linking of a function:\\
{\tt extern "C" { void myfunc(); } } \\
\begin{verbatim}
#include "pdlmgr.h"
// ...
string soname = "mylib.so";
string funcname = "myfunc";
PDynLinkMgr dyl(soname);
DlFunction f = dyl.GetFunction(funcname); 
if (f != NULL) { 
// Calling the function 
  f();    
}
\end{verbatim}

\index{CxxCompilerLinker}
The {\bf CxxCompilerLinker} class provides the services to compile C++ code and building 
shared libraries, using the same compiler and options which have 
been used to create the SOPHYA shared library. 
The sample program below illustrates using this class to build
the shared library (myfunc.so) from the source file myfunc.cc :
\begin{verbatim}
#include "cxxcmplnk.h"
// ...
string flnm = "myfunc.cc";
string oname, soname;
int rc;
CxxCompilerLinker cxx;
// The Compile method provides a default object file name
rc = cxx.Compile(flnm, oname);
if (rc != 0 ) { // Error when compiling ... }
// The BuildSO method provides a default shared object file name
rc = cxx.BuildSO(oname, soname);
if (rc != 0 ) { // Error when creating shared object ... }
\end{verbatim}

\subsection{Command interpreter}
The class {\bf Commander} can be used in interactive programs to provide
c-shell like command interpreter and scripting capabilties. 
Arithmetic expression evaluation is implemented through the {\bf CExpressionEvaluator}
and {\bf RPNExpressionEvaluator} classes. 
 
\newpage
\section{Module SkyMap}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{PixelMap}
\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
\dclsc{SphereThetaPhi}
\dclsb{LocalMap}
\caption{partial class diagram for pixelization classes  in Sophya}
\end{figure}
The {\bf SkyMap} module provides classes for creating, filling, reading pixelized spherical  and 2D-maps. The types of values stored in pixels can be int, float, double , complex etc. according to the specialization of the template type. 
\subsection {Spherical maps}
There are two kinds of spherical maps according pixelization algorithms. SphereHEALPix represents spheres pixelized following the HEALPIix algorithm (E. Hivon, K. Gorski)
\footnote{see the HEALPix Homepage: http://www.eso.org/kgorski/healpix/ }
, SphereThetaPhi represents spheres pixelized following an algorithm developed at LAL-ORSAY. The example below shows creating and filling of a SphereHEALPix with nside = 8 (it will be 12*8*8= 768 pixels) :
\index{\tcls{SphereHEALPix}} 
\index{\tcls{SphereThetaPhi}} 

\begin{verbatim}
#include "spherehealpix.h"
// ...
SphereHEALPix<double> sph(8);
for (int k=0; k< sph.NbPixels(); k++) sph(k) = (double)(10*k);
\end{verbatim}

SphereThetaPhi is used in a similar way with an argument representing number of slices in theta (Euler angle) for an hemisphere.
\index{\tcls{SphereThetaPhi}}

\subsection {Local maps}
\index{\tcls{LocalMap}}
A local map is a 2 dimensional array, with i as column index and j as row index. The map is supposed to lie on a plan tangent to the celestial sphere in a point whose coordinates are (x0,y0) on the local map and (theta0, phi0) on the sphere. The range of the map is defined by two values of angles covered respectively by all the pixels in x direction and all the pixels in y direction (SetSize()). Default value of (x0, y0) is middle of the map, center of pixel(nx/2, ny/2).

Internally, a map is first defined within this reference plane and tranported until the point (theta0, phi0) in such a way that both axes are kept parallel to meridian and parallel lines of the sphere. The user can define its own map with axes rotated with respect to reference axes (this rotation is characterized by angle between the local parallel line and the wanted x-axis-- method SetOrigin(...))

The example below shows creating and filling of a LocalMap with 4 columns and 5 rows. The origin is set to default. The map covers a sphere portion defined by two angles of 30. degrees (methods \textit{SetOrigin()} and \textit{SetSize()} must be called in order to completely define the map).
\begin{verbatim}
#include "localmap.h"  
//..............       
  LocalMap<r_4> locmap(4,5);
  for (int k=0; k<locmap.NbPixels();k++) locmap(k)=10.*k;
  locmap.SetOrigin();
  locmap.SetSize(30.,30.);
\end{verbatim}

\subsection{Writing, viewing \dots }

All these objects have been design to be written to or read from a persistant file.
The following example shows how to write the previously created objects
into such a file~:
\begin{verbatim}
//-- Writing

#include "fiospherehealpix.h" 
//................

char *fileout = "myfile.ppf";
POutPersist outppf(fileout);
FIO_SphereHEALPix<r_8> outsph(sph);
outsph.Write(outppf);
FIO_LocalMap<r_8> outloc(locmap);
outloc.Write(outppf);
// It is also possible to use the << operator
POutPersist os("sph.ppf");
os << outsph;
os << outloc;
\end{verbatim}

Sophya graphical tools (spiapp) can automatically display and operate
all these objects.

\newpage
\section{Module Samba}
\index{Spherical Harmonics}
\index{SphericalTransformServer}
The module provides several classes for spherical harmonic analysis. The main class is \textit{SphericalTranformServer}. It contains methods for analysis and synthesis of spherical maps. The following example fills a vector of Cl's, generate a spherical map from these Cl's. This map is analysed back to Cl's...
\begin{verbatim}
#include "skymap.h"
#include "samba.h"
....................

// Generate input spectra  a + b* l + c * gaussienne(l, 50, 20)
int lmax = 92;
Vector clin(lmax);
for(int l=0; l<lmax; l++) {
  double xx = (l-50.)/10.;
  clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
}

// Compute map from spectra
SphericalTransformServer<r_8> ylmserver;
int m = 128;  // HealPix pixelisation parameter
SphereHEALPix<r_8> map(m);
ylmserver.GenerateFromCl(map, m,  clin, 0.);
// Compute power spectrum from map
Vector clout = ylmserver.DecomposeToCl(map, lmax,  0.);
\end{verbatim}

\newpage
\section{Module SkyT}
\index{RadSpectra} \index{SpectralResponse}
The SkyT module is composed of two types of classes:
\begin{itemize}
\item{} one which corresponds to an emission spectrum of
radiation, which is called RadSpectra
\item{} one which corresponds to the spectral response
of a given detector (i.e. corresponding to a detector
filter in a given frequency domain), which is called
SpectralResponse.
\end{itemize}
\begin{figure}[hbt]
\dclsbb{RadSpectra}{RadSpectraVec}
\dclsb{BlackBody}
\dclsccc{AnyDataObj}{SpectralResponse}{SpecRespVec}
\dclsc{GaussianFilter}
\caption{partial class for SkyT module}
\end{figure}

\begin{verbatim}
#include "skyt.h"
// ....
// Compute the flux from a blackbody at 2.73 K through a square filter
BlackBody myBB(2.73);
// We define a square filter from 100 - 200 GHz
SquareFilter mySF(100,200);
// Compute the filtered integrated flux :
double flux = myBB.filteredIntegratedFlux(mySF);
\end{verbatim}

A more detailed description of SkyT module can be found in: 
{\it The SkyMixer (SkyT and PMixer modules) - Sophya Note No 2. }
available also from Sophya Web site.

\newpage
\section{Module FitsIOServer}
\begin{figure}[hbt]
\dclsbb{FitsFile}{FitsInFile}
\dclsb{FitsOutFile}
\end{figure}
\index{FITS} \index{FitsInFile} \index{FitsOutFile}
This module provides classes for handling file input-output in FITS format using the cfitsio library. It works like the SOPHYA persistence (see Module SysTools), using delegate objects, but its design is simpler. The following example  writes a matrix (see module TArray) and a spherical map (see module SkyMap)  on a FITS file and reads back from FITS file and creates new objects :
\begin{verbatim}
#include "spherehealpix.h"   
#include "fitsspherehealpix.h" 
#include "fitstarray.h" 
#include "tmatrix.h" 
//...........................

int m=...;
SphereHEALPix<r_8> sph(m);
................
int dim1=...;
int dim2=...;
TMatrix<r_8> mat(dim1,dim2);
............

FITS_SphereHEALPix<r_8> sph_temp(sph);
FITS_TArray<r_8> mat_temp(mat);
// writing

FitsOutFile os("myfile.fits");
sph_temp.Write(os);
mat_temp.Write(os);

// reading
FitsInFile is("myfile.fits");
sph_temp.Read(is);
mat_temp.Read(is);
SphereHEALPix<r_8> new_sph=(SphereHEALPix<r_8>)sph_temp;
TMatrix<r_8> new_mat=(TMatrix<r_8>)mat_temp;
................ 
      
\end{verbatim}

The operators {\tt operator << (FitsOutFile ...)} and 
{\tt operator >> (FitsInFile ...)} are defined in order 
to facilitate the FITS file operations:
\begin{verbatim}
// Writing an array object to a FITS file
#include "fitstarray.h" 
FitsOutFile fio("arr.fits");
Matrix m(20,30);
m = 12345.;
fio << m;
// .....
// Reading a binary table to a XNTuple
#include "fitsxntuple.h" 
XNTuple xn;
FitsInFile fii("table.fits");
fii >> xn;
\end{verbatim}

The class {\bf FITS\_AutoReader} provides a limited FITS files reading
and decoding capabilities. A partial class diagram of FITS persistence 
handling classes is shown below:
\begin{figure}[hbt]
\dclsbb{FitsIOhandler}{FITS\_TArray}
\dclsb{FITS\_NTuple}
% \dclsb{FITS\_XNTuple}
\dclsb{FITS\_SphereHEALPix}
% \dclsb{FITS\_LocalMap}
\end{figure}

\newpage
\section{LinAlg and IFFTW modules}
An interface to use LAPACK library (available from {\tt http://www.netlib.org})
is implemented by the {\bf LapackServer} class, in module LinAlg.
\index{LapackServer}.
The sample code below shows how to use SVD (Singular Value Decomposition)
through LapackServer:
\begin{verbatim}
#include "intflapack.h"
// ...
// Use FortranMemoryMapping as default
BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
// Create an fill the arrays A and its copy AA
int n = 20;
Matrix A(n , n), AA;
A = RandomSequence(RandomSequence::Gaussian, 0., 4.);  
AA = A;  // AA is a copy of A
// Compute the SVD decomposition
Vector S;  // Vector of singular values
Matrix U, VT;
LapackServer<r_8> lpks;
lpks.SVD(AA, S, U, VT);
// We create a diagonal matrix using S
Matrix SM(n, n); 
for(int k=0; k<n; k++) SM(k,k) = S(k);
// Check the result : A = U*SM*VT
Matrix diff = U*(SM*VT) - A;
double min, max;
diff.MinMax(min, max);
cout << " Min/Max difference Matrix (?=0) , Min= " << min 
     << " Max= " << max << endl;
\end{verbatim}

\index{FFTWServer}
The {\bf FFTWServer} class (in module FFTW) implements FFTServerInterface class 
methods, for one dimensional and multi-dimensional Fourier 
transforms on double precision data using the FFTW package
(available from {\tt http://www.fftw.org}).

\newpage
\section{Building and installing Sophya}
\subsection{supported platforms}
Presently, the Sophya library has been tested with the following
compiler/platform pairs:

\begin{center}
\begin{tabular}{|l|l|}
\hline
OS & compiler \\
\hline
HP/Compaq/DEC Tru64 ( OSF1)  &     cxx  (6.1 , 6.3)  \\
Linux (RH)          &     g++  (3.2)  \\
Linux  (SCL)          &     icc  (8.1)   (Intel compiler) \\
SGI IRIX64       &     CC   (7.3)   \\
MacOSX/Darwin 10.3  &         g++ 3.3  \\
\hline
\end{tabular}
\end{center}

Some of the modules in the Sophya package uses external libraries. The 
{\bf FitsIOServer} is the example of such a module, where the {\tt libcfitsio.a}
is used.
par
The object files from a given Sophya module are grouped in an archive library 
with the module's name ({\tt libmodulename.a}). All Sophya modules 
 are grouped in a single shared library ({\tt libsophya.so}), while the 
modules with reference to external libraries are grouped in 
({\tt libextsophya.so}). The {\bf PI} and {\bf PIext} modules are 
grouped in ({\tt libPI.so}).

\subsection{Installation} 

The build procedure has two main steps: \\
- The configure step (BuildMgr/configure) setup the directory structure and
the necessary configuration file. \\
- The make step compiles the different sources files, create the library and optionaly
builds all or some of the associated executables. 

\par
{\tt BuildMgr/configure } is a c-shell script with a number of arguments:
\begin{verbatim}
csh> ./configure -h
configure [-sbase SOPHYABASE] [-scxx SOPHYACXX] [-incln] 
  [-minc mymake.inc] 
  [-extp dir1 -extp dir2 ...] [-extip dir1 -extip dir2 ... ] 
  [-extlp dir1 -extlp dir2 ... ]
  [-noextlib -noext fits -noext fftw -noext lapack ]
  [-noext astro -noext minuit]
\end{verbatim}
\begin{itemize}
\item[] -sbase : define SOPHYA installation base directory. \$SOPHYABASE is used 
if not specified.
\item[] -scxx : selects the C++ compiler. \$SOPHYACXX  s used 
if not specified.
\item[] -incln : creates symbolic link for include files, instead of copying them. 
\item[] -minc : give an explicit name for the file used to generate 
\$SOPHYABASE/include/sophyamake.inc.
\item[] -extp : Adds the specied path to the search path of the external libraries 
include files and archive library.
\item[] -extip : Adds the specied path to the search path of the external libraries 
include files.
\item[] -extp : Adds the specied path to the search path of the external libraries 
archive (libxxx.a).
\item[] -noextlib : Disable compiling of modules referencing external libraries.
\item[] -noext : Disable compiling of the specified module (with reference to external
library. 
\end{itemize}

In the example below, we assume that we want to install Sophya from a 
released (tagged) version in the source directory {\tt \$SRC} in the
{\tt /usr/local/Sophya} directory, using {\tt g++}. We assume that 
the external libraries can be found in {\tt /usr/local/ExtLibs/}. 
We disable the compilation of the MinuitAdapt and XAstrPack packages.

\vspace*{3mm}
\begin{verbatim}
# Create the top level  directory
csh> mkdir /usr/local/Sophya/
csh> cd $SRC/BuildMgr
# Step 1.a : Run the configuration script 
csh> ./configure -sbase /usr/local/Sophya -scxx g++  -extp /usr/local/ExtLibs/ \
-noext astro -noext minuit 
# Step 1.b : Check the generated file  $SOPHYABASE/include/
# Step 2.a: Compile the modules without external library reference
csh> make libs
# Step 2.b: Compile the modules WITH  external library reference (optional) 
csh> make extlibs
# Step 2.c: Build libsophya.so
csh> make slb
# Step 2.d: Build libextsophya.so (optional)
csh> make slbext
# Step 2.e: Compile the PI and PIext modules (optional)
csh> make PI
# Step 2.f: Build the corresponding shared library libPI.so (optional)
csh> make slbpi
\end{verbatim}

To compile all modules and build the shared libraries, it is possible 
to use:
\begin{verbatim}
# Step 2.a ... 2.f
csh> make all slball 
\end{verbatim}

At this step, all libraries should have been made. Programs using
Sophya libraries can now be built:
\begin{verbatim}
# To compile test programs
csh> cd ../PrgUtil
csh> make 
# To build (s)piapp (libPI.so is needed)
csh> cd ../ProgPI
csh> make 
csh> cd ..
\end{verbatim}

\subsection{Mgr module}
This module contains scripts which can be used for generating the 
makefiles for each module.
\begin{itemize}
\item {\bf Makefile} Top level Makefile for building the libraries.
\item {\bf Makefile.h} contains the definition of compilation flags for the 
different compilers and systems. This file is used for building the 
library and generating {\bf MakefileUser.h} (to be included in makefiles).
\item {\bf Makefile.slb} contains the rules for building shared libraries
for the different compilers and systems. (to be included in makefiles)
\item {\bf crerep\_sophya} c-shell script for creating the directory tree
under {\tt \$SOPHYABASEREP} and {\tt \$SOPHYADEVREP}
\item {\bf install\_sophya} c-shell script for installing the Sophya package.
Usually from {\tt \$SOPHYADEVREP} to {\tt \$SOPHYABASEREP}
\item {\bf mkmflien} c-shell script for making symbolic links or copying
include files to {\tt \$SOPHYADEVREP/Include} or {\tt \$SOPHYABASEREP/Include}
\item {\bf mkmf} c-shell script for generating module makefiles and the 
top level makefile (named GNUmakefile)
\item {\bf mkmflib} c-shell script for generating each library module 
makefile (named GNUmakefile)
\item {\bf mkmfprog} c-shell script for generating makefile for a module
containing the source for executable programs (named GNUmakefile).
\item {\bf mkmfPI} c-shell script for generating makefile for PI and PIext
modules (named GNUmakefile)
\item {\bf libdirs} List of Sophya modules without reference to external
libraries.
\item {\bf extlibdirs} List of Sophya modules with reference to external
libraries.

\end{itemize}
 
\newpage
\appendix
\section{SOPHYA Exceptions}
\index{Exception classes} \index{PThrowable} \index{PError} \index{PException}
SOPHYA library defines a set of exceptions which are used 
for signalling error conditions. The figure below shows a partial 
class diagram for exception classes in SOPHYA.
\begin{figure}[hbt]
\dclsbb{PThrowable}{PError}
\dclscc{PError}{AllocationError}
\dclscc{PError}{NullPtrError}
\dclscc{PError}{ForbiddenError}
\dclscc{PError}{AssertionFailedError}
\dclsbb{PThrowable}{PException}
\dclscc{PException}{IOExc}
\dclscc{PException}{SzMismatchError}
\dclscc{PException}{RangeCheckError}
\dclscc{PException}{ParmError}
\dclscc{PException}{TypeMismatchExc}
\dclscc{PException}{MathExc}
\dclscc{PException}{CaughtSignalExc}
\caption{partial class diagram for exception handling in Sophya}
\end{figure}

For simple programs, it is a good practice to handle 
the exceptions at least at high level, in the {\tt main()} function. 
The example below shows the exception handling and the usage 
of Sophya persistence.
 
\input{ex1.inc}


\newpage
\addcontentsline{toc}{section}{Index}
\printindex 
\end{document}