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

\documentclass[twoside,10pt]{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

\newcommand{\rond}{$\bullet \ $}
\newcommand{\etoile}{$\star \ $}
\newcommand{\cercle}{$\circ \ $}
\newcommand{\carre}{$\Box \ $}

\newcommand{\bitem}[1]{\item[]{\bf #1} }

\begin{document}

\begin{titlepage}
%  The title page - top of the page with the title of the paper
\titrehp{SOPHYA user's guide }
%  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     \\
}
% \auteursall
%  The title page - bottom of the page with the paper number
\vspace{1cm}
\begin{center}
{\bf \Large SOPHYA Version: 2.2 (V\_Sep2010) }
\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 of 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 packages and libraries, 
encapsulating them in classes in order to facilitate their usage.
\par
Our main requirements in designing SOPHYA classes can be summarized as
follow:
\begin{itemize}
\item[\rond] Provide a comprehensive set of data containers, such as arrays and tables 
(tuple) covering the most common needs in scientific simulation and data analysis
softwares.
\item[\rond] Take advantage of the C++ language and define methods and operators 
for most basic operation, such as arithmetic operations, in a rather intuitive way, while
maintaining performances comparable to low level coding in other languages 
(C, Fortran, F90 \ldots) 
\item[\rond] Simplify memory management for programmers using the class library.
This has been a strong requirement for most SOPHYA classes. Automatic reference 
sharing and memory management is implemented in SOPHYA classes intended 
for large size objects. We recommend to allocate SOPHYA objects on the stack,
including when objects are returned by methods or functions. 
See section \ref{memgt} for more information.
\item[\rond] Archiving, importing (reading) and exporting (writing) data in a 
efficient and consistent way is a major concern in many scientific software 
and projects. SOPHYA provide a native data I/O or persistence system, 
(PPF, \ref{ppfdesc}) as well as import/export  services for ASCII and FITS formats.
\end{itemize}

% \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}. 
%%%
%%%
\subsection{Acknowlegments}
Many people have contributed to the development SOPHYA and/or the PI library 
and (s)piapp interactive analysis tool, in particular:

\begin{tabular}{lcl}
Reza Ansari & \hspace{5mm} & (LAL-Univ.Paris Sud, Orsay) \\
Eric Aubourg & & (DAPNIA-CEA/APC, Saclay) \\
Sophie Henrot-Versille & & (LAL-IN2P3/CNRS, Orsay) \\
Alex Kim & & (LBL, Berkeley) \\
Guy Le Meur & & (LAL-IN2P3/CNRS, Orsay) \\
Eric Lesquoy & & (DAPNIA-CEA, Saclay) \\
Christophe Magneville & & (DAPNIA-CEA, Saclay) \\
Bruno Mansoux & & (LAL-IN2P3/CNRS, Orsay) \\
Olivier Perdereau & & (LAL-IN2P3/CNRS, Orsay) \\
Nicolas Regnault & & (LPNHE-IN2P3/CNRS, Paris) \\
Benoit Revenu & & (APC/Univ.Paris 7, Paris) \\
Francois Touze & & (LAL-IN2P3/CNRS, Orsay) \\
\end{tabular}

We thank also the persons who have helped us by useful suggestions, among others : \\
S. Bargot, S. Plasczczynski, C. Renault and D. Yvon.

\subsection{SOPHYA modules}
\label{sopmodules}
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 
as 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). This module needs extensive 
checking and corrections.
\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. 
FITS is maintained by NASA and SAO and is available from: \\
\href{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html}
{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html}
\item[] {\bf LinAlg/} Interface with Lapack linear algebra package.
Lapack is a linear algebra package and can be downloaded from: \\ 
\href{http://www.netlib.org/lapack/}{http://www.netlib.org/lapack/}
\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a).
FFTW is a package for performing Fourier transforms, written in C.
Documentation and source code can be found at: \\
\href{http://www.fftw.org/}{http://www.fftw.org/}
\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}. 
Information on Xephem can be found at : \\
\href{http://www.clearskyinstitute.com/xephem/}{http://www.clearskyinstitute.com/xephem/}

\end{itemize}

The following modules contain each a set of related programs using the 
SOPHYA library.
\begin{itemize}
\item[] {\bf Tests/} Simple test programs. Many of these test programs can also be 
used as examples for using SOPHYA. 
\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
\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 and scripts:
\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.
\item[] {\bf DemoData/} data files (PPF, FITS \ldots) useful for the (s)piapp 
user's manual examples.
\end{itemize}   

The following modules contains additional programs or libraries, which rely on SOPHYA:
\begin{itemize}
\item[] {\bf AnaLC} contains program files extracted from the EROS/PEIDA software
and adapted to be compiled with SOPHYA. It can be used in particular to read and analyse
EROS light curve data files (fichiers de suivi). 
\item[] {\bf LUC}  ( {\bf L}ittle {\bf U}niverse {\bf C}alculator) is a module containing classes to 
perform basic computation related to the universe geometry (FRW metric).
\item[] {\bf SimLSS} Tools for the {\bf Sim}ulation and analysis of cosmological 
{\bf L}arge {\bf S}cale {\bf S}tructures
\item[] {\bf PIPhoto} is an interface module for (s)piapp and ImageMagick for 
import/export of graphics in image format (jpeg, gif \ldots) in piapp.
\item[] {\bf PMixer/} skymixer and related programs. This module is {\bf obsolete}.
\end{itemize}

\newpage 

\section{Using Sophya}
The organisation of SOPHYA directories and some of the associated 
utility programs are described in this section.
Basic usage of Sophya classes is described 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}
\label{directories}
The environment variable {\bf SOPHYABASE}  is used 
to define the path where the Sophya libraries and binaries 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 or .dylib on Darwin/MacOS)
\item \$SOPHYABASE/exe : Path for binary program 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, while some more specific flags,
are found in  { \tt \$SOPHYABASE/include/sspvflags.h }  

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. It also copy (or create symbolic link) for all SOPHYA include 
files as well as symbolic links for external libraries 
include files path in {\tt \$SOPHYABASE/include} (FitsIO, FFTW, XAstro \ldots). 

Object files for each module are grouped in a static archive library 
by the build procedure (libXXX.a for module
XXX, with XXX = BaseTools, TArray, HiStats,  FitsIOServer  \ldots).

When shared libraries are build, all stand alone SOPHYA modules
are grouped in  {\tt libsophya.so},  {\tt libextsophya.so} contains 
the interface modules with external libraries {\bf (FitsIOServer, LinAlg \ldots)},
while {\bf PI, PIext, PIGcont} modules  are grouped in  {\tt libPI.so}.
Alternatively, it is possible to group all modules in a single shared 
library {\tt libAsophyaextPI.so} (See \ref{build})

In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
should contain the Sophya shared library path 
({\tt \$SOPHYABASE/slb}). 
On Silicon Graphics machines with IRIX64 operating system, 
the default SOPHYA configuration correspond to the 64 bit architecture. 
The environment variable { \bf LD\_LIBRARY64\_PATH } replace in
this case the usual {\bf LD\_LIBRARY\_PATH} variable.
On IBM machines with AIX, the {\bf LIBPATH} environment variables 
contains the shared libraries search path.

When using the dynamic load services in SOPHYA ({\tt PDynLinkMgr}
class), in runcxx or (s)piapp applications for example, the shared 
library search path must contain the current working directory (
dot . in unix).

\subsection{Copy constructor and assignment operator}
\label{memgt}
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.

\subsection{Multi-thread programming with SOPHYA}
Multi-thread programming is usually safe as long as different threads DO NOT access
the same memory locations. SOPHYA is mainly organized as classes having only 
data members, with very few cases having static data members (memory locations
common to all instances of a class), or seldom, global variables. \\
Using different instances of a class in different threads is thus safe for most 
classes / methods in SOPHYA. \\
A major exception is the reference sharing mechanism, where different instances
of a class may  share memory locations. This reference sharing mechanism has been
made thread-safe in SOPHYA, from version V=2.1. \\
As a consequence, different execution threads can access non overlapping sub-arrays
and access (read/write) the corresponding elements without the need of 
mutex and thread synchonisation. Moreover, thread safe filling of 
NTuple and DataTable objects can be activated, on an object per object basis. \\
The {\tt ZThread, ZMutex \ldots} classes in the {\bf SysTools } module offer a relatively 
easy way of writing multi-threaded programs (See \ref{mtpsop}).

\subsection{the runcxx program}
\index{runcxx} \label{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 preprocessor macro {\tt KeepObj()} is defined by runcxx and let the user 
to save an objet to an PPF file, with the same name as the corresponding variable.
The example below shows how to compile, link and run a sample
code.
\begin{verbatim}
##  File example.icc :

Matrix mxa(3,3);
mxa = IdentityMatrix(1.);
cout << mxa ; 
// Save the object mxa in a PPF file named mxa
KeepObj(mxa);

## Executing this sample code
csh> runcxx -f example.icc
\end{verbatim}

\subsection{the scanppf program}
\index{scanppf} \label{scanppf}
{\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  2.0 Revision 0 (V_Jul2006) -- Jul 17 2006 14:13:27 cxx 
 Usage: scanppf [flags] filename 
 flags = -s -n -a0 -a1 -a2 -a3 -lh -lho -lmod 
   -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 
   -lmod : List initialized/registered modules 
\end{verbatim}

\subsection{the scanfits program}
\index{scanfits} \label{scanfits}
{\bf scanfits} is a SOPHYA program using the FitsIOServer 
\footnote{FitsIOServer module uses the cfitsio library. scanfits has to be linked with
with FitsIOServer module and cfitsio libraries, or libextsophya.so}
module which can be used 
to analyse the content of FITS files. It can list the FITS headers,  the appropriate 
SOPHYA-FITS handler (implementing {\tt FitsHandlerInterface}) class, and the list of
all registered  FITS handlers. 
\begin{verbatim}
csh> scanfits -h
 PIOPersist::Initialize() Starting Sophya Persistence management service 
SOPHYA Version  2.0 Revision 0 (V_Jul2006) -- Jul 17 2006 14:13:27 cxx 
 Usage: scanfits [flags] filename 
 flags = -V1 -lh -rd -header 
   -V1 : Scan using old (V1) code version 
   -lh : Print the list of registered handlers (FitsHandlerInterface) 
   -rd : try to read each HDU data using appropriate handler 
   -header : List header information 
\end{verbatim}

\subsection{the spiapp program}
\index{spiapp} \label{spiapp}
{\bf spiapp} is an interactive data analysis program, built on top of the SOPHYA 
library, , the PI GUI library. The interactive data analysis framework is defined 
by the classes in the {\bf PIext} module. spiapp has a c-shell like 
command interpreter and can be used to manipulate SOPHYA and display 
objects. Refer to the piapp user manual for more information.
\begin{verbatim}
csh> spiapp -h 
 SophyaInitiator::SophyaInitiator() BaseTools Init
 PIOPersist::Initialize() Starting Sophya Persistence management service 
SOPHYA Version  2.1 Revision 0 (V_Nov2007) -- Nov 24 2007 13:08:58 gcc 3.3 20030304 (Apple Computer, Inc. build 1495)

 piapp: Interactive data analysis and visualisation program 
 Usage: piapp [-nored] [-doublered] [-termread] [-term] 
              [-hidezswin] [-small] [-nosig] [-nosigfpe] [-nosigsegv] 
              [-tmpdir TmpDirectory] [-help2tex] [-exec file [args]] 
  -nored : Don't redirect stdout/stderr to piapp console
  -doublered : Redirect stdout/stderr to piapp console AND the terminal 
  -termread : Read commands on terminal (stdin)
  -term : equivalent to -nored -termread -small 
  -hidezswin : Hide Zoom/Stat/ColMap window 
  -small : Create small size main piapp window 
  -nosig : Don't catch SigFPE, SigSEGV 
  -nosigfpe -nosigsegv: Don t catch SigFPE / SigSEGV 
  -tmpdir TmpDirectory: defines TMDIR for temporary files 
  -help2tex: Create a LaTeX help file (piahelp.tex)
  -exec file [args] : Execute command file (last option)

\end{verbatim}


\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 (object persistence services).
\vspace*{5mm}

\subsection{Initialisation}
\index{SophyaInitiator}
A number of actions have to be taken before 
some of the services provided by SOPHYA become operational. This is the case
of SOPHYA persistence, as well as FITS I/O facilities. 
Initialisation of many SOPHYA modules is performed through an initialiser class,
which inherits from {\bf SophyaInitiator}. 
\par
Static instance of each initialiser class exist in the library and the various SOPHYA services
should be operational when the user code ({\tt main()}) starts, except for 
modules in the second or third shared libraries 
({\tt libextsophya.so libPI.so}). Indeed, a problem related 
to the initialisation of shared libraries arises on some systems
(Darwin/Mac OS X in particular) causing program crash at start-up, 
if static instance of initialiser class is present in the second shared library. 
The FitsIOServer module should thus be explicitly initialised in the user 
program. 
\par 
In cases where the run time loader does not perform correctly the static 
object initialisation, the initialiser class for the modules used in the 
program must be instanciated in the beginning of your main program: \\
{\tt TArrayInitiator , HiStatsInitiator , SkyMapInitiator , FitsIOServer \ldots}
%%%
\subsection{SOPHYA persistence}
\label{ppfdesc}
\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}
The most useful methods for using Sophya
persistence are listed below. A brief description of the PPF file format 
and some guidelines for writing writing delegate classes for handling 
object persistence can be found in the following paragraphs.
\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.

\subsubsection{PPF file format}
The PPF file consist of :
\begin{itemize}
\item[\rond] A file header (3x32=96 bytes) , containing  an identification string, 
PPF format version (currently V3), the file endiannes {\tt (BIG-ENDIAN , LITTLE-ENDIAN)}
and the file creation date and time. It should be noted that the present SOPHYA version 
(V=2.2) does not allow updating an existing PPF file.
\item[\rond] A collection of tagged data items and data structuring tags. 
the PPF tags are one byte (8 bits) long and may be followed by a length information
for arrays, or strings. The length information might be 4 bytes or 8 bytes long.
Characters, various length (1,2,4,8 bytes long) signed and unsigned integers, 
simple or double precision (4/8) bytes floating point and complex numbers are 
the basic data types handled by PPF streams.
\begin{verbatim}
tag=PPS_SIMPLE_Integer4  +  data (4 bytes)
tag=PPS_SIMPLE_Float4 + data (4 bytes)
tag=PPS_SIMPLE_Complex8 + data (2x8=16 bytes)
tag=PPS_STRING + Length (4 bytes) + data (Length bytes)
tag=PPS_SIMPLE_ARRAY4_UInt2 + Length (4 bytes) + data (2xLength bytes)
\end{verbatim}
The two tags  {\tt PPS\_OBJECT} and {\tt PPS\_ENDOBJECT} are used for identifying 
objects. 
\item[\rond] The file trailer contains some global statistics such as the total number 
of objects in file, the list of name tags, with their corresponding positions. 
These informations are grouped under the identifier tag {\tt  PPS\_NAMETAG\_TABLE }.
The last byte in file contains the value {\tt PPS\_EOF}
\end{itemize}
We show below the result of the analysis of a PPF file using 
 {\tt PPFBinaryInputStream::AnalyseTags()}. This can easily be done using the 
 {\bf runcxx} program.
The PPF file contains four objects: 
a {\tt TimeStamp} object, two {\tt NDataBlock} objects, the second one sharing its data with 
the first one. The second NDataBlock is only written as a reference to the first object. The fourth 
object is a {\tt RandomGenerator} which has an {\tt NDataBlock} object as data member. 
Notice the corresponding nested structure of object marker tags.
The first object and the last objects in the file are written associated with a {\bf NameTag}.
{\small 
\begin{verbatim}
 ---------------------------------------------------------- 
 PPFBinaryInputStream::AnalyseTags(Level= 49)
   FileName= toto.ppf
   Version= 3 FileSize= 8884 Creation Date= Fri Dec  7 15:30:58 2007

 NbPosTag=0 NbNameTag=2 NbObjs=4 NbTopLevObjs=3 NbRefs=1 MaxNest=2

<PPS_NAMETAG_MARK> tag at position 60
<PPS_OBJECT> tag at position 61 ClassId= c09dfe032b341cee  ObjectId= 10
  <PPS_SIMPLE> tag at position 72 DataType=INTEGER x4
  <PPS_SIMPLE> tag at position 77 DataType=INTEGER x8
  <PPS_SIMPLE> tag at position 80 DataType=FLOAT x8
<PPS_ENDOBJECT> tag at position 89  ObjectId= 10
<PPS_OBJECT> tag at position 92 ClassId= e670300b367585d1  ObjectId= 21
  <PPS_SIMPLE_ARRAY4> tag at position a3 DataType=UNSIGNED x8 NElts= 3
  <PPS_SIMPLE_ARRAY4> tag at position c0 DataType=INTEGER x4 NElts= 50
<PPS_ENDOBJECT> tag at position 18d  ObjectId= 21
<PPS_REFERENCE> tag at position 196  ObjectId= 21  OrigPos=92
<PPS_NAMETAG_MARK> tag at position 1a7
<PPS_OBJECT> tag at position 1a8 ClassId= eb6c427a5a30caed  ObjectId= 30
  <PPS_SIMPLE_ARRAY4> tag at position 1b9 DataType=UNSIGNED x4 NElts= 6
  <PPS_SIMPLE> tag at position 1d6 DataType=UNSIGNED x8
  <PPS_SIMPLE> tag at position 1df DataType=UNSIGNED x8
  <PPS_OBJECT> tag at position 1e8 ClassId= 6531a8f47336d4aa  ObjectId= 41
    <PPS_SIMPLE_ARRAY4> tag at position 1f9 DataType=UNSIGNED x8 NElts= 3
    <PPS_SIMPLE_ARRAY4> tag at position 216 DataType=FLOAT x8 NElts= 1024
  <PPS_ENDOBJECT> tag at position 221b  ObjectId= 41
<PPS_ENDOBJECT> tag at position 2224  ObjectId= 30
<PPS_NAMETAG_TABLE> tag at position 222d
<PPS_NAMETAG_ENTRY>  NameTag=rg-randgen NameTagMark Position=1a7
<PPS_NAMETAG_ENTRY>  NameTag=ts-timestamp NameTagMark Position=60
<PPS_EOF> tag at position 22b3 TagPos=222d
 PPFBinaryInputStream::AnalyseTags() - End - Total Number of Tags= 23
 ---------------------------------------------------------- 
\end{verbatim}
}  % fin de small 

\subsubsection{Writing PPF handlers}
Here are some guidelines for creating PPF handler classes to make 
a class persistent through the SOPHYA PPersist mechanism. The example
discussed here can be found in the {\bf Examples} module, in the directory 
{\bf MyPPF}. 
\begin{enumerate}
\item The class should inherit from the {\bf SOPHYA::AnyDataObj} class. 
   In the example here, we want to make the class {\bf Vfs} persistent
\begin{verbatim}   
class Vfs : public SOPHYA::AnyDataObj {
 ... 
}
\end{verbatim}
%%%%
\item The PPF handler class 
should inherit from {\bf SOPHYA::PPersist} . The pure virtual methods 
of PPersist class (DataObj() , SetDataObj() , ReadSelf(), WriteSelf()
must be implemented. 
\begin{verbatim}   
class PPFHandlerVfs : public SOPHYA::PPersist {
public:
  virtual SOPHYA::AnyDataObj* DataObj();
  virtual void       SetDataObj(SOPHYA::AnyDataObj &);
protected: 
  virtual void       ReadSelf(SOPHYA::PInPersist&);
  virtual void       WriteSelf(SOPHYA::POutPersist&) const;
}
\end{verbatim}
%%%
It is possible to use the template class {\tt SOPHYA::ObjFileIO<T>}, and 
specialize it for the tharget class (here {\tt SOPHYA::ObjFileIO<Vfs>}), 
by defining the two methods :  \\
\hspace*{5mm} {\tt  SOPHYA::ObjFileIO<Vfs>::ReadSelf(SOPHYA::PInPersist\&) } \\
\hspace*{5mm}  {\tt  SOPHYA::ObjFileIO<Vfs>::WriteSelf(SOPHYA::POutPersist\&) const }
%%%%
\item If it is NOT possible to have the target class inherit from {\bf AnyDataObj}, 
a wrapper class should be used. The same class can play the role of wrapper
AND the PPF handler. See the PPF handler / wrapper class for STL vectors : \\
\hspace*{5mm} {\tt   SOPHYA::PPFWrapperSTLVector<T>} in file BaseTools/ppfwrapstlv.h 
%%%
\item Implement the ReadSelf() and WriteSelf() methods of the PPF handler.
All the I/O services from the following classes can be used : 
\begin{itemize}
\item PPFBinaryIOStrem , PPFBinaryInputStream , PInPersist
\item PPFBinaryIOStrem , PPFBinaryOutputStream , POutPersist
\end{itemize}
Writing and reading of the embeded objects for which a handler has been register
ed can simply be performed by :  \\
\hspace*{5mm} {\tt   POutPersist::PutObject() } \\
\hspace*{5mm} {\tt   PInPersist::GetObject() }

{\bf Warning:} The services associated with nametags : \\
\hspace*{5mm}  {\tt PPFBinaryOutputStream::WriteNameTag() } \\
\hspace*{5mm}  {\tt PPFBinaryInputStream::GotoNameTag() } \\
are {\bf NOT} intented to be used in WriteSelf() , ReadSelf()
%%%%%
\item The new PPF handler, as well as the list of classes it can handle has to be
registered prior to use PPF read/write for the target classes. This must be 
performed during the initialization phase, for example at the beginning of the 
main() program. Another possibility is to use a module initializer 
(See {\bf SophyaInitiator} class in file BaseTools/sophyainit.h ) 
and declare a static instance of the class. Notice that this works only if 
the system loader handles correcly the call of constructor 
for the statically declared objects. 

The registration can be performed using the CPP macros defined in 
BaseTools/ppersist.h
\begin{verbatim}
  // First, register the PPF handler ObjFileIO<Vfs> 
  PPRegister(ObjFileIO<Vfs>);
  // Register the list of classes which can be handled by ObjFileIO<Vfs> 
  DObjRegister(ObjFileIO<Vfs>, Vfs);
\end{verbatim}
%%%%%%%%%%%%%%%%
\end{enumerate}

%%%%%%%%%%%%
\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{DVList, MuTyV and TimeStamp classes}
\index{DVList} \index{MuTyV} \index{TimeStamp}
\begin{figure}[hbt]
\dclsa{MuTyV}
\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.
{\bf MuTyV} objects can also hold {\bf TimeStamp } 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}

\begin{figure}[hbt]
\dclsbb{AnyDataObj}{TimeStamp}
\end{figure}
%
The {\bf TimeStamp} class represent date and time and provides
many standard operations, such as Initialisation from strings,
conversion to strings and time interval computations. \\
Usage example:
\begin{verbatim}
// Create a object with the current date and time and prints it to cout
TimeStamp ts; 
cout << ts << endl;
// Create an object with a specified date and time 
TimeStamp ts2("01/01/1905","00:00:00");
// Get the number of days since 0 Jan 1901
cout << ts2.ToDays() << endl;

// Combined use of TimeStamp and MuTyV 
string s;         
TimeStamp ts;      // Current date/time
MuTyV mvt = ts;
s = mvt;           // s contains the current date in string format
cout << s << endl;
\end{verbatim}

\subsection{\tcls{SegDataBlock} , \tcls{SwSegDataBlock}}
\index{\tcls{SegDataBlock}}  \index{\tcls{SwSegDataBlock}}
%%
\begin{figure}[hbt]
\dclsccc{AnyDataObj}{\tcls{SegDBInterface}}{ \tcls{SegDataBlock} }
\dclscc{\tcls{SegDBInterface}}{ \tcls{SwSegDataBlock} }
\end{figure}
\begin{itemize} 
\item[\rond] \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[\rond] \tcls{SwSegDataBlock} Implements the same \tcls{SegDBInterface}
using a data swapper object. Data swappers implement the interface defined in 
(\tcls{DataSwapperInterface} class.  \tcls{SwSegDataBlock} can 
thus be used to handle arrays with very large number of objects.
These classes handles reference sharing.
\end{itemize}
 
 \subsection{Pseudo-random number generators}
 \index{Random numbers}  \label{randgen}
\begin{figure}[hbt]
\dclsbb{ AnyDataObj }{ RandomGeneratorInterface } 
\vspace*{5mm}
\dclsccc{RandomGeneratorInterface}{DR48RandGen}{ ThSDR48RandGen }
\dclsbb{ RandomGeneratorInterface }{ FMTRandGen }
\end{figure}
The {\bf RandomGeneratorInterface } is a pure virtual class which 
defines the interface for random number generator classes and implements the 
generation algorithm for few probability distribution functions:
flat, Gaussian, Poisson, Exponential, 2D gaussian \ldots. The pure virtual method
{\tt Next() } should be implemented by inheriting  classes.
 \index{RandomGeneratorInterface}

Several implementations using different pseud-random generators are provided by SOPHYA:
\begin{enumerate}
\item {\bf DR48RandGen} uses the well known {\tt drand48() }  pseudo-random
sequence generator. It implements also the possibility of initializing the generator state, as well saving and 
restoring the generator state through the SOPHYA PPF mechanism. The  {\tt drand48() }  and thus 
{\bf DR48RandGen} objects are not NOT thread safe. Different instances of {\bf DR48RandGen}  share 
the same underlying 48 bit state.
\item {\bf ThSDR48RandGen} is a thread safe version of DR48RandGen. Using the constructor flag, in can 
be instanciated to reproduce the non thread-safe behaviour of DR48RandGen. 
Several instances of this class can be used in different threads without the risk of 
corrupting the internal state of the drand48() generator. However, in multi-thread applications,
 there is no guarantee to obtain the same sequence of numbers in each thread.  
\index{ ThSDR48RandGen} 
\item {\bf FMTRandGen } ids an implementation of RandomGeneratorInterface using 
the SIMD-oriented Fast Mersenne Twister (SFMT). It uses the code developed by 
Mutsuo Saito and Makoto Matsumoto from Hiroshima University. For more information on this 
method, refer to \\
\href{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html} 
{ttp://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html} \\
The instances of FMTRandGen are thread safe, each instance having its own internal state 
and produces an independent sequence of random numbers. 
\index{FMTRandGen }
\end{enumerate}

A number of c-like functions are defined in the file BaseTools/srandgen.h and can be used 
to generate random sequences: \\
{\tt drand01() , drandpm1() ,  GaussianRand() , PoissonRand \ldots} \\
These functions use a global instance  of RandomGeneratorInterface class which can be 
defined and accessed through the static methods:\\
\hspace*{5mm} { \tt RandomGeneratorInterface* RandomGeneratorInterface::GetGlobalRandGenP() } \\
\hspace*{5mm} { \tt void SetGlobalRandGenP(RandomGeneratorInterface* rgp) } \\

In multi-threaded programs, an instance of ThSDR48RandGen or FMTRandGen 
should be created in each thread running in parallel and needing to generate random sequences. 
All the three classes discussed here implement the automatic initialization method defined 
{\tt (AutoInit() ) } defined in the interface class and implement a number of other initialization
methods {\tt (SetSeed() )}.  SOPHYA provides also PPF handlers for these classes 
which can be used to save the complete  state of the object and the underlying generator to PPF files. 
The generator object state can be restored subsequently from the PPF file. This feature can be used 
to generate very long sequence of random numbers, distributed over several runs.
The example below illustrates the use of the generator classes and state save/restore 
through PPF files:
\begin{enumerate} 
\item We create and use a {\bf ThSDR48RandGen} object in a first program. Its sate is saved 
to the PPF file rg.ppf at some point in the program. 
\begin{verbatim}
// ------ program A 
// A.1- Create a thread safe generator based on drand48() 
ThSDR48RandGen rg;
// A.2- Auto initilize its state (using the system time) 
rg.AutoInit();
// A.3- compute and print a smal sequence of random numbers 
int N = 10;
for(int i=0; i<N; i++) 
  cout << " I=" << i << "  rand_flat01= " << rg.Flat01() 
       << " rand.gaussian= " << rg.Gaussian() << endl;
// A.4- Save the generator state for subsequent use
POutPersist po("rg.ppf");
po << rg;
// A.5- compute and print a second sequence of random numbers 
for(int i=0; i<N; i++) 
  cout << "++ I=" << i+N << "  rand_flat01= " << rg.Flat01() 
       << " rand.gaussian= " << rg.Gaussian() << endl;
\end{verbatim}

\item The pseudo-random generator object state is restored in a second program from the previously created 
PPF file, before being used :
 \begin{verbatim}  
// ------ program B
// B.1- Create and initialize the generator from the previously saved state
ThSDR48RandGen rgr;
PInPersist pin("rg.ppf");
pin >> rgr;
int N = 10;
// B.2- Compute and print a sequence of random number, 
//      should be compared to the sequence A.5
for(int i=0; i<N; i++) 
  cout << "-- I=" << i << "  rand_flat01= " << rgr.Flat01() 
       << " rand.gaussian= " << rgr.Gaussian() << endl;
\end{verbatim}
\end{enumerate}

%%%%%%%%%%%%
\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 from inline functions is
not in the 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{Arithmetic operations}
The four usual arithmetic operators ({\bf  + \, - \, * \, / }) are defined
to perform constant addition, subtraction, multiplication and division. 
The three operators ({\bf  + \, - \, / }) between two arrays of the same type
are defined to perform element by element addition, subtraction
and division. In order to avoid confusion with matrix multiplication, 
element by element multiplication is defined by overloading the
operator {\bf \, \&\& \, }, as shown in the example below:
\begin{verbatim}
TArray<int_4> a(4,3), b(4,3), c , d, e;
a = RegularSequence(1.,1.);
b = RegularSequence(10.,10.);
cout << a << b ;
c = a && b; 
d = c / a;
e = (c / b) - a;
cout << c << d << e;
\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.

\subsubsection{TVectors and std::vector}
SOPHYA \tcls{TVector} objects can be constructed or filled from \tcls{std::vector} objects.
TVector objects can also be converted to std::vector.
\begin{verbatim}
TVector(const vector<T>& v, short lcv=BaseArray::AutoVectorType);
sa_size_t TVector<T>::FillFr(const vector<T>& v,bool noresize=false);
sa_size_t TVector<T>::FillTo(vector<T>& v,bool addtoend=false);
vector<T> TVector<T>::ConvertTostdvec()

\end{verbatim}

In addition, \tcls{std::vector} objects can be written to, read from PPF streams
through the class \tcls{PPFWrapperSTLVector} (file ppfwrapstlv.h). The std::vector
PPF handler is currently registered  for the following data types: \\
\hspace*{5mm} {\tt string, int\_2, uint\_2, int\_4, uint\_4, int\_8, uint\_8 } \\
\hspace*{5mm} {\tt r\_4 , r\_8, complex$<$r\_4$>$ ,  complex$<$r\_8$>$ }

\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}
\hline \\
\multicolumn{2}{c}{ {\bf Range} {\tt (start, end, size, step) } } \\[2mm]
\hline \\
{\bf Range} {\tt r(7); } &  index range:  \hspace{2mm} 7 \\
{\bf Range} {\tt r(3,6); } &  index range:  \hspace{2mm} 3,4,5,6 \\
{\bf Range} {\tt r(3,7,2); } &  index range: \hspace{2mm} 3,5,7 \\
{\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 \\[2mm]
\hline 
\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{Persistence, IO}
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. Only 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} 

{\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 (dumped) again to the same  PPF stream. 
However, this behavior can be changed using the {\tt RenewObjId()} method,
as illustrated below.
\begin{verbatim}
{
POutPersist pos("mca.ppf");
TArray<int_4> ia(5,3);
ia = 8;
pos << ia;         // (1) 
ia = 16;
pos << ia;        // (2) Only a reference to the previously ia array is written
ia = 32;
ia.RenewObjId();    // We change the object Id
pos << ia;        //  (3) The complete array is dumped again 
}
\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, for the first two 
write operations. When we 
read the three array from the file mca.ppf, the same array elements
are obtained two times (all elements equal to 8), and a different array is obtained 
the third time 
\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{Cast without conversion}
Data conversion between arrays with different data type is handled transparently,
through the copy constructor or the assignment (=) operator . However, in rare cases, 
one wants to access the same memory locations, without data type conversion.
The template functions defined in {\tt arrctcast.h} can be used to access the same 
memory locations, by arrays with different data types. The SOPHYA/NDataBlock 
reference sharing mechanism is effective When using these functions.
Notice that the array size or stride may change during these cast operations. \\
 {\tt arrctcast.h} has been introduced in version V=2.1 (Nov 2007), and has not been 
 fully tested for non packed arrays.
\begin{verbatim}
// We define and initialize a  real array :
TArray<r_4> fa(5);
fa = RegularSequence(1.25,0.5);
cout << " fa= " << fa;
// We construct an integer array from fa, where the floating point values
// are converted to integer values
TArray<uint_2> ia(fa);
cout << "  ia= " << ia;
cout << "  ia (in hex)= " << hex << ia << dec;
// We can also access the fa memory locations interpreted as short integers
uint_2 ui2;
// Note that sfia size is double the fa size 
TArray<uint_2> sfia = ArrayCast(fa, ui2);
cout << "  sfia= " << sfia;
cout << "  sfia (in hex)= " << hex << sfia << dec;
\end{verbatim}
One of the most useful case of these array type cast without conversion 
correspond to accessing the real or imaginary part of a complex array. 
Two specific template functions {\tt SDRealPart()} and {\tt SDImagPart()}
are also defined in  {\tt arrctcast.h}. Two other functions {\tt ArrCastR2C()}
and {\tt ArrCastC2R() } are also defined for real to complex, and
complex to real cast. 
Their usage is shown in the next  paragraph on complex arrays.

\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 these functions create an array and copies
the data from the original complex array. They do not provide
shared memory access to real and imaginary parts. 
For shared memory access, use functions {\tt SDRealPart()} and
{\tt SDImagPart() } (see below). 
\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}

\noindent {\bf Shared data access :} It is possible to access a complex array 
elements (real and imaginary parts) through the template functions defined 
in {\tt arrctcast.h} and discussed above. Four specific template functions are 
defined for shared data access to complex arrays: 
\begin{itemize}
\item {\bf SDRealPart } return a float/double TArray, when applied to an
array with complex elements, corresponding the the {\bf real} part of the complex values.
\item {\bf SDImagPart } return a float/double TArray, when applied to an
array with complex elements, corresponding the the {\bf imaginary} part of the complex values.
\item {\bf ArrCastC2R } return a float/double TArray, when applied to an
array with complex elements. The data is shared between the two arrays, the real array 
containing real and imaginary parts as successive elements. 
\item {\bf ArrCastR2C } return a complex valued TArray, when applied to a float/double array. 
\end{itemize} 
The example below shows how to use these functions:

\begin{verbatim}
// We define a complex array 
TArray< complex<r_4> > za(5);
cout << " za= " << za;
// And two real arrays, corresponding to the real and imaginary parts 
TArray<r_4> rza = SDRealPart(za); 
TArray<r_4> iza = SDImagPart(za); 
// We initialize the real and imaginary parts of the complex array 
rza = RegularSequence(10.,2.);
iza = RegularSequence(5.,0.75);
cout << " rza=..., iza=... ----> za = " << za;
// The complex array seen as a real array (double size)
TArray<r_4> aza = ArrCastC2R(za);
cout << " za --> aza= " << aza;
TArray< complex<r_4> > zb;
zb = ArrCastR2C(aza);
KeepObj(zb);
\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)  (1,0)  (2,0)  (3,0) |
    | (0,1)  (1,1)  (2,1)  (3,1) |
    | (0,2)  (1,2)  (2,2)  (3,2) |
\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}
BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
TArray<r_4> X(4,2); 
X = RegularSequence(1,1);
cout << "Array X= " << X << endl;
TMatrix<r_4> X_C(X, true);
cout << "Matrix X_C (CMemoryMapping) = " << X_C << endl;
TMatrix<r_4> X_F(X_C, true);
X_F.TransposeSelf();    // we make it a FortranMemoryMapping matrix
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}

{\bf Note:} Only the {\tt TMatrix} class is sensitive to the memory organisation. It should also 
be noted that the first index for a matrix is the row index. For a {\tt TArray} object, the first 
index correspond always to the X axis. For a matrix in the {\tt CMemoryMapping}, the row
index is the Y index of the corresponding array, while for a matrix 
with {\tt FortranMemoryMapping}, the row index is the X index.
\begin{verbatim}
#  2D array  arr ,  matrix mx
TArray<r_4> arr; 
...
TMatrix<r_4> mx(arr);
...
###   CMemoryMapping  :
mx( row , col ) -> arr( ix=col, jy= row )
###   FortranMemoryMapping  :
mx( row , col ) -> arr( ix=row, jy= col )
\end{verbatim}
The following code fragment shows the difference between element access 
index order for the different memory organisation:
\begin{verbatim}
   TArray<int_4> a(5,3);
   a = RegularSequence(10,2);
   cout << "Array a = " << a << endl;
   TMatrix<r_4> mac(a);
   cout << "Matrix mac (CMemoryMapping) = " << mac << endl;
   TMatrix<r_4> maf(mac);
   maf.TransposeSelf();    // we make it a FortranMemoryMapping matrix
   cout << "Matrix maf (FortranMemoryMapping) = " << maf << endl;
   // ---- element access
   cout <<  " Array(ix,jy) : a(2,0)= " << a(2,0) << " a(1,2)= " << a(1,2) 
        <<   "  a(0,2)= " << a(0,2) << endl;
   cout <<  " mac(row,col) :  mac(2,0)= " << mac(2,0) << " mac(1,2)= " << mac(1,2) 
        <<   "  mac(0,2)= " << mac(0,2) << endl;
   cout <<  " maf(row,col) :  maf(2,0)= " << maf(2,0) << " maf(1,2)= " << maf(1,2) 
        <<   "  maf(0,2)= " << maf(0,2) << endl;
\end{verbatim}


\subsection{Special Square Matrices (SSQM) }
SOPHYA V2.2 introduces new classes for handling special square matrices (diagonal, symmetric, triangular \ldots).
The figure below shows the corresponding class hierarchy:
\begin{figure}[hbt]
\dclsccc{AnyDataObj}{\tcls{SpecialSquareMatrix}}{ \tcls{ DiagonalMatrix } }
\dclsc{ \tcls{ LowerTriangularMatrix } } 
\dclsc{ \tcls{ SymmetricMatrix } } 
\end{figure}

\index{ \tcls{DiagonalMatrix} }
\index{ \tcls{SymmetricMatrix} }
\index{ \tcls{TriangularMatrix} }

Theses classes provides methods similar to TArray/TMatrix objects for manipulating these special kind of matrices.
However, no sub-matrix extraction method is provided. The  following features are currently implemented:
\begin{itemize}
\item These classes implements reference sharing and behave in a way similar to TArray objects. The copy 
constructor shares the data, while the equal operator (=) performs an element by element copy. 
\item The {\bf FIO\_SpecialSquareMatrix$<$T$>$ } manages the persistence (I/O) in PPF format for these classes. 
\item Constructor from and conversion to general {\bf TMatrix$<$T$>$ } objects. 
 Non relevant elements are ignored when constructing from a general matrix and the special square matrix can be 
converted to the general matrix through the {\tt ConvertToStdMatrix () }.
\item Definition of element access operator {\tt mx(r,c) } as well as global assignment operator {\tt mx = a; mxb=mxa; }.
operator {\tt mx[k] } can be used to access the non zero element k. 
\item Addition, subtraction of a constant, multiplication or division by a constant, as well as  
\item Element by element addition, subtraction, multiplication and division methods between same type of SSQM 
matrices. As in the cases of general matrices, the four binary operators ( + , - , \&\& ,   /) are overloaded to perform 
these operations. 
\item Matrix multiplication operation between same type SSQM objects or between a general matrix and an SSQM
object is defined through the methods: \\
{\\ Multiply() , MultiplyXG(), MultiplyGX(),  X=D,S,T } \\
The binary multiplication operator ( * ) is then overloaded to perform matrix multiplication.   
\end{itemize}
The sample code below illustrates the use of these special matrices:
\begin{verbatim} 
// Create and initialize a diagonal matrix
DiagonalMatrix<double> diag(3);
diag(0,0)=1.;  diag(1,1)=2.;  diag(2,2)=3.;
// use of multiplication by a constant operator
diag *= 10.;
// check the diagonal matrix
cout << diag << diag.ConvertToStdMatrix() << endl;
// define and initialize a general matrix
TMatrix<double> mx(3,3);
mx=RegularSequence();
cout << mx;
// check the result of a matrix mutiplication
cout << mx*diag;
\end{verbatim} 

\newpage 

\section{Module HiStats}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{Histo}
\dclscc{Histo}{HProf}
\dclsbb{AnyDataObj}{Histo2D}
\dclsbb{AnyDataObj}{HistoErr}
\dclsbb{AnyDataObj}{Histo2DErr}
\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 DataTable} which are
more or less the same as the binary or ascii FITS tables.

\subsection{Histograms}
\subsubsection{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}

\subsubsection{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}

\subsubsection{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.

\subsubsection{Histograms HistoErr and Histo2DErr}
\index{HistoErr}
The {\bf HistoErr} are basic histograms where the number of entries for each bin is kept.
Methods to compute of the mean and the variance in each bin are provided.
The {\bf Histo2DErr} is the same for $2$ dimensions.

\subsection{Data tables (tuples)}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{NTuple}
\dclsccc{AnyDataObj}{BaseDataTable}{DataTable}
\dclscc{BaseDataTable}{SwPPFDataTable}
\end{figure}

\subsubsection{NTuple}
\index{NTuple}
{\bf 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 interactive analysis tool can be used to create 1D, 2D and 3D plots using the date from 
NTuples and DataTables. \\
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}

{\bf XNTuple} provide additional functionalities, compared to NTuple. However,
this class is deprecated and superseded by classes inheriting from  BaseDataTable.
It is 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. 

\subsubsection{DataTables}
\label{datatables}
\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}


\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 this class. A similar class
SwFitsDataTable (\ref{SwFitsDataTable}), using 
FITS files as swap space is also provided in the FitsIOServer module.

\index{DataTableRow}
The class {\bf DataTableRow } is an auxiliary class which simplifies the manipulation
of BaseDataTable object rows. 
The example below show how to create and filling a table, using a PPF stream as 
swap space. In addition, we have used a {\tt DataTableRow} to prepare data 
for  each table line.
\begin{verbatim}
   #include "swppfdtable.h"
   // ...
   {
   // ---------- Create an output PPF stream (file)
   POutPersist po("swdtable.ppf");
   // ------------------
   // Create a table with 3 columns, using the above stream as swap space
   SwPPFDataTable dtrow(po, 64);
   dtrow.AddStringColumn("sline");
   dtrow.AddIntegerColumn("line");
   dtrow.AddDateTimeColumn("datime");   
   //
   TimeStamp ts, ts2;  // Initialize current date and time
   string sline;
   //---- Create a table row with the required structure
   DataTableRow row = dtrow.EmptyRow();
   // ----- Fill the table
   for(int k = 0; k<2500; k++) {
     sline = "L-";
     sline += (string)MuTyV(k);
     row["sline"] = sline;
     row[1] = k;
     ts2.Set(ts.ToDays()+(double)k);
     row["datime"] = ts2;
     dtrow.AddRow(row);
  }
  //------ Write the table itself to the stream, before closing the file
  po << PPFNameTag("SwTable") << dtrow;
  }
\end{verbatim}
%%
The  previously created table can easily be read in, as shown below:
%%
\begin{verbatim}
   #include "swppfdtable.h"
   // ...
   {
   // ------ Create the input  PPF stream (file)
   PInPersist pin("swdtable.ppf");
   // ------ Read in the SwPPFDataTable object 
   SwPPFDataTable dtr;
   pin >> PPFNameTag("SwTable") >> dtr;
   // ---- Create a table row with the required structure
   DataTableRow row = dtr.EmptyRow();
   // ---- Acces and print two of the table rows :
   cout << dtr.GetRow(6, row) << endl;
   cout << dtr.GetRow(17, row) << endl;   
  }
\end{verbatim}
   
\subsection{Persistence and visualisation }
%%
Histogram and NTuple/DataTable objects have PPF handlers and 
can be written to or read from PPF files. Sophya interactive analysis tool 
(spiapp) can automatically display and operate on all of these objects.
The following example shows how to write the previously created objects
into such a file~:
\begin{verbatim}
{
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}

\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{Optimisation and parameter fitting}
\subsubsection{GeneralFit: non linear 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.

\subsubsection{Simplex method}
The class {\bf MinZSimplex} implements the simplex method for non linear 
optimization / minimization. See the SOPHYA reference manual and the 
Tests/tsimplex.cc for more information.

\subsection{Interpolation}
\index{Interpolation} \index{SLinInterp1D}
The class {\bf SLinInterp1D} (file slininterp.h) can be used for linear 1D interpolation. 
\begin{verbatim}
#include "slininterp.h"
#include "srandgen.h"
  ...
vector<double> ys;
double xmin = 0.5;
double xmax = 0.;
for(int i=0; i<=12; i++) {
  xmax = xmin+i*0.1;
  yreg.push_back(sin(xmax)*cos(2.2*xmax));
}
SLinInterp1D interpYR(xmin, xmax, yreg);
cout << interpYR
for(int i=0; i<=12; i++) {
  double x = drand01()*2.;
  cout << " Interpol result for X=" << x << " -> " << interpYR(x) 
       << " ?= " << sin(x)*cos(2.2*x) << endl;
}
\end{verbatim}
The {\bf CSpline} and {\bf CSpline2} classes perform one dimensional and two dimensional spline interpolation.

\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, such as thread creation and management, dynamic loading and 
resource usage information. For example, yhe class {\bf Periodic} provides the 
necessary services needed to implement the execution of a periodic action.

\subsection{Resource usage (CPU, memory \ldots) }
 The class {\bf ResourceUsage} \index{ResourceUsage}
and {\bf Timer} \index{Timer} provides access to information 
about various resource usage (memory, CPU, ...).
The class {\bf Timer} \index{time (CPU, elapsed)} and c-functions
{\tt InitTim() , PrtTim(const char * Comm) } can be used to print 
the amount of CPU and elapsed time in programs.

The following sample code illustrates the use of {\bf ResourceUsage} :
\begin{verbatim}
  // How to check resource usage for a given part of the program
  ResourceUsage res;
  // --- Part of the program to be checked : Start
  // ...
  res.Update();
  cout << " Memory size increase (KB):" << res.getDeltaMemorySize() << endl;
  cout << " Resource usage info : \n" << res << endl;
\end{verbatim}

\subsection{Thread management classes}
\index{ZThread} \index{ZMutex}  \label{thrcls}
A basic interface to POSIX threads is also provided
through the \index{threads} {\bf ZThread}, {\bf ZMutex} and {\bf ZSync}
classes. The best way to use thread management classes is by inheriting
from {\bf ZThread} and redefining the {\tt run() } method. 
It is also possible to use the default {\tt run() } implementation and associate 
a function to perform the action, as in the example below :
\begin{verbatim}
  // The functions to perform computing
  void fun1(void *arg) { }
  void fun2(void *arg) { }
  // ...
  ZThread zt1;
  zt1.setAction(fun1, arg[1]);
  ZThread zt2;
  zt2.setAction(fun2, arg[2]);
  cout << " Starting threads ... " << endl;
  zt1.start();
  zt2.start();
  cout << " Waiting for threads to end ... " << endl;  
  zt1.join();
  zt2.join();
\end{verbatim}
The classes {\bf ZMutex} \index{mutex} and {\bf ZSync} can be used 
to perform synchronisation and signaling between threads. Each {\bf ZMutex} object
contains a pair of mutex and condition variable.
The {\bf ParallelExecutor } class can be used for parallel simultaneous execution 
of several function of objects inheriting from {\bf  ParallelTaskInterface}. 
Some hints about parallel programming and usage of these classes can be found in 
\index{ParallelExecutor} \index{ ParallelTaskInterface }
\ref{mtpsop}. 

\begin{figure}[hbt]
\dclsa{ParallelExecutor}
\dclsbb{ParallelTaskInterface}{ParallelTaskFunction}
\dclsbb{ZThread}{ParalExThread}
\dclsa{ZMutex}
\end{figure}


\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}
\index{Commander}
\index{Expression evaluation}
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. The latter can be used to perform evaluation
in reverse polish notation (old HP calculators), 
while the former uses the usual algebraic (c-language like) expressions.
The command language provides variable manipulation through the usual 
{\tt \$varname} vector variable and arithmetic expression extensions, as well 
as the control and test blocs. A description of this command interpreter syntax 
and possibilities can be found in the (s)piapp user's guide.
\begin{verbatim}
#include "commander.h"
...
Commander cmd;
char* ss[3] = {"foreach f ( AA bbb CCCC ddddd )", "echo $f" , "end"};
for(int k=0; k<3; k++) {
  string line = ss[k];
  cmd.Interpret(line);
}
\end{verbatim}
 
\newpage
\section{Module SkyMap}
\begin{figure}[hbt]
\dclsbb{AnyDataObj}{PixelMap}
\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
\dclsc{SphereThetaPhi}
\dclsc{SphereECP}
\dclsb{LocalMap}
\caption{partial class diagram for spherical map 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{3D geometry}
Some of the classes in this module simplify the angle manipulation, or the 3D geometry
and rotation calculations, such as the {\bf Vector3d} and {\bf UnitVector} classes.
\index{Vector3d}
The classes {\bf Vector3d} and {\bf UnitVector}  are useful for angle conversions.
\index{Angle}
{\bf LongLat} class and {\bf Angle} class (defined in file vector3d.h) can be used for
angle conversions :
\begin{verbatim}
   // Example to convert 0.035 radians to arcsec
   double vr = 0.035;
   cout << "Angle rad= " << vr << " ->arcsec= " << Angle(vr).ToArcSec() << endl;
   // Example to convert 2.3 arcmin to radian - we use the conversion operator
   double vam = 2.3;
   cout << "Angle arcmin= " << vam << " ->rad= " 
        << (double)Angle(vam, Angle::ArcMin) << endl;
\end{verbatim}
%%%
\subsection {Spherical maps}
SkyMap module provides three classes for representing data with pixels distributed over complete ($4 \pi$ steradians) spheres. These classes implements the common interface defined 
in the base class \tcls{SphericalMap} with three different algorithms or 
pixelization scheme. 
The spherical maps can be instanciated for the followind data types: \\
\hspace*{5mm} int\_4 , r\_4 (float) , r\_8 (double) , complex$<$r\_4$>$ , complex$<$r\_8$>$. \\
The SphereHEALPix can in addition be instanciated for T=uint\_2.

\begin{enumerate}
\item \index{\tcls{SphereHEALPix}} 
{\bf \tcls{SphereHEALPix}} implements the HEALPix 
({\bf H}ierarchical {\bf E}qual {\bf A}rea iso{\bf L}atitude {\bf Pix}elization) scheme, 
developped originaly by K. Gorski \& E. Hivon. Refer to 
\href{http://healpix.jpl.nasa.gov/}{HEALPix home page} for 
detailed information about this pixelisation scheme and related software
\footnote{HEALPix home page: http://healpix.jpl.nasa.gov/ }.
FITS read/write for SphereHEALPix objects is handled by the \tcls{FITS\_SphereHEALPix}
class in module FitsIOServer.
\item \index{\tcls{SphereThetaPhi}} 
{\bf \tcls{SphereThetaPhi}} represents spheres pixelized following an algorithm 
developed at LAL-ORSAY, for SOPHYA. 
The sphere is divided into a number of rings or slices
along the parallels, corresponding to different values of the angle $\theta$.
Each slice is then divided into a number of pixels, with an aspect ratio close 
to one (square pixels). Pixels are exactly iso-latitude and very uniform in surface,
over the sphere. 

\item \index{\tcls{SphereECP}}

The {\bf \tcls{SphereECP}} class correspond to the cylindrical projection.
Like SphereThetaPhi class, the sphere is divided into a number of rings
or slices, and each ring is divided into a constant number of pixels
along the $\varphi$ direction. Although the \tcls{SphereECP} does not have 
equal area pixels when used for complete spheres, it can be used for 
representing partial or full spherical maps. 
\end{enumerate}

The example below shows creating and filling of a 
SphereHEALPix with nside = 8 
(The sphere will have $12 \times 8 \times 8= 768$ pixels) :

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

Pixels at an angular posistion can be directly accessed through the operator \\
\hspace*{15mm} T \tcls{SphericalMap}::operator()($\theta, \varphi$) :
\begin{verbatim}
#include "vector3d.h"
#include "spherethetaphi.h"
...
// Create a sphere with 40 arcmin resolution
int M = SphereThetaPhi<r_4>::ResolToSizeIndex(
        Angle(40., Angle::ArcMin).ToRadian() );
SphereThetaPhi<r_4> sphtp(M);
double tet, phi;
for (int k=0; k< sphtp.NbPixels(); k++) {
  sphtp.PixThetaPhi(k, tet, phi);
  sphtp(k) = cos(5.*tet)*sin(7.*phi);
}
cout << sphtp;
// To save sphtp to file sphtp (if executed through runcxx)
KeepObj(sphtp);
\end{verbatim}

\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{Persistence and visualisation }
The different SkyMap objects have PPF handlers and written to or read from PPF files.
The following example shows how to save 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 interactive analysis tool (spiapp) can automatically display and operate
on SkyMap objects. 

\newpage
\section{Samba and SkyT}
\subsection{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}

\subsection{Module SkyT}
\index{RadSpectra} \index{SpectralResponse} 
\begin{center}
{\bf \large THIS MODULE NEEDS EXTENSIVE CHECKING/DEBUGGING AND REDESIGN } 
\end{center}
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}
This module provides classes for handling file input-output in FITS 
\footnote{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html} 
format using the cfitsio library. Its
design is similar to  the SOPHYA persistence (see Module BaseTools). 
Delegate classes or handlers perform the actual read/write from/to fits files. 
\par
Compared to the SOPHYA native persistence (PPF format),
FITS format has the advantage of being used extensively, and handled 
by a many different software tools. It is a de facto standard in 
astronomy and astrophysics. 
However, FITS lacks some of the features present in SOPHYA PPF, and although 
many SOPHYA objects can be saved in FITS files, FITS persistence has 
some limitations. For example, FITS does not currently handle complex arrays.
\subsection{FITS streams}
\index{FITS} \index{FitsInOutFile} 
%%
The class {\bf FitsInOutFile} can be seen as a wrapper class for the cfitsio library functions.  
This class has been introduced in 2005 (V=1.9), when the module has been 
extensively  changed. In order to keep backward compatibility, the old fits wrapper 
classes ({\bf  FitsFile, FitsInFile, FitsOutFile}) has been changed to inherit from
 {\bf FitsInOutFile}.  The use of class FitsFile and specific services of these old classes 
 should  be avoided, but FitsInFile, FitsOutFile can be safely considered a specialisation
 of FitsInOutFile for read/input and write/output operations respectively.
 Most c-fitsio errors are converted to an exception: {\bf FitsIOException}.
 \par
 File names are passed to cfitsio library. It is thus possible to use cfitsio file name conventions,
 such as {\bf ! } as the first character, for overwriting existing files (when creating files).
 The diagram below shows the class hierarchy for cfitsio wrapper classes.
\begin{figure}[hbt]
\dclsa{FitsInOutFile}
\dclscc{FitsFile}{FitsInFile}
\dclscc{FitsFile}{FitsOutFile}
\end{figure}
%%%%
\subsection{FITS handlers and I/O operators}
\index{FitsManager}
Handlers classes inheriting from {\bf FitsHandlerInterface} perform write/read operations
for AnyDataObj objects to/from FitsInOutFile streams. The {\bf FitsManager} class provides
top level services for object read/write in FITS files. 
\par  In most cases, 
\hspace{5mm} {\tt FitsInOutFile\& $<<$ } \, and \,  {\tt FitsInOutFile\& $>>$ } \hspace{5mm}
operators can be used to write and read objects.  
When reading objects from a fits file using the {\tt FitsInOutFile\& $>>$ } operator, 
the fits file is positioned on the next HDU, after reading. Also, if the {\bf FitsInOutFile}
object is positioned on a first empty HDU (without data, naxis=0), reading in objects
corresponding to a binary or ascii table using the operator $>>$  will skip automatically 
the empty HDU and position the fits file on the second HDU, before trying to read in
the object.
\par
The two main types of fits data structures, images and tables 
{\tt (IMAGE\_HDU , BINARY\_TBL , ASCII\_TBL)} are handled by the generic handlers: \\
{\bf \tcls{FitsArrayHandler}}  and {\bf FitsHandler$<$BaseDataTable$>$}.  
\par
A number of more specific handlers are also available, in particular for NTuple,
\tcls{SphereHealPix} and \tcls{SphereThetaPhi}.  \\[2mm]
{\bf Warning:} Some handlers were written with the old FitsIOServer classes.
They  inherit from the intermediate class {\bf FitsIOHandler} and
have been adapted to the new scheme. \\[2mm] 
%%%
The examples below illustrates the usage of FitsIOServer classes. They can be compiled
and executed using runcxx, without the {\tt include} lines: \\[1mm]
\hspace*{5mm} {\tt csh> runcxx -import SkyMap -import FitsIOServer -inc fiosinit.h } 
%%%
\begin{enumerate}
\item Saving an array and a HealPix map to a Fits file
\index{HEALPix-FITS}
\begin{verbatim}
#include "fitsioserver.h"
#include "fiosinit.h"
// ....
{
// Make sure FitsIOServer module is initialised :
FitsIOServerInit();
//  Create and open a fits file named myfile.fits
FitsInOutFile fos("myfile.fits", FitsInOutFile ::Fits_Create);
// Create and save a 15x11 matrix of integers
TMatrix<int_4> mxi(15, 11);
mxi = RegularSequence(10.,10.);
fos << mxi;
// Save a HEALPix spherical map using FitsManager services 
SphereHEALPix<r_8> sph(16);
sph = 48.3;
FitsManager::Write(fos, sph);
// --- The << operator could have been used instead :    fos << sph;
}
\end{verbatim}
%%%%
%%%%
\item Reading objects and the header from the previously created fits file:
\begin{verbatim}
{
FitsIOServerInit();   //  Initialisation 
// ----  Open the fits file named myfile.fits
FitsInFile fis("myfile.fits");
//---- print file information on cout
cout << fis << endl;
//--- Read in the array 
TArray<int_4> arr; 
fis >> arr;
arr.Show();
//--- Position on second HDU 
fis.MoveAbsToHDU(2);
//--- read and display header information
DVList hdu2;
fis.GetHeaderRecords(hdu2, true, true);
cout << hdu2;
//--- read in the HEALPix map 
SphereHEALPix<r_8> sph;
FitsManager::Read(fis, sph);
// --- The >> operator could have been used instead :    fis >> sph;
sph.Show();
}
\end{verbatim}
%%%%%%%
%%% 
\item DataTable objects can be read from and written to FITS files as ASCII or
binary tables. The example below show reading the DataTable created in the example
in section \ref{datatables} from a PPF file and saving it to a fits file.
\index{DataTable-FITS}
\begin{verbatim}
#include "swfitsdtable.h"
// ....
{
FitsIOServerInit();   //  FitsIOServer Initialisation 
//--- Reading in DataTable object from PPF file
PInPersist pin("dtable.ppf");
DataTable dt;
pin >> dt;
dt.Show(); 
//--- Saving table to FITS
FitsInOutFile fos("!dtable.fits", FitsInOutFile ::Fits_Create);
fos << dt;
}
\end{verbatim}
%%%%
\end{enumerate}
%%%
A partial class diagram of FITS persistence handling classes is shown below. The 
class {\tt FitsIOhandler} conforms to the old FitsIOServer module design and 
should not be used anymore. 
\begin{figure}[hbt]
\dclsbb{FitsHandlerInterface}{FitsArrayHandler$<$T$>$}
\dclsb{\tcls{FitsHandler}}
\dclscc{FitsIOhandler}{FITS\_NTuple}
\dclsc{FITS\_SphereHEALPix}
% \dclsb{FITS\_LocalMap}
\end{figure}

\subsection{Fits Array I/O in chunk }
Array read/write from/to FITS files (Image HDU) is rather straightforward, in particular thanks to 
the overloaded stream operators, as in the examples given above and here:
\index{Array-FITS}
\begin{verbatim}
//  ---- Writing an array to fits file
//  Create and open a fits file named myfile.fits
FitsInOutFile fos("!myarr.fits", FitsInOutFile ::Fits_Create);
TArray<r_8> a(50,30);
// ... fill the array
// Write the array to file
fos << a;

// ---- Reading an array from file myarr.fits
// Open the fits file for reading
FitsInOutFile fis("myarr.fits", FitsInOutFile::Fits_RO);
TArray<r_8> a;
// Read in the array 
fis >> a;
\end{verbatim}

It is also possible to write arrays FITS Image HDU in chunks  or read arrays in chunk, 
using the following methods : \\
\begin{verbatim}
void FitsArrayHandler<T>::ReadAtOffset(FitsInOutFile& is, sa_size_t* offset);
void FitsArrayHandler<T>::WriteAtOffset(FitsInOutFile& os, sa_size_t* offset);
\end{verbatim} 
The sample code below show how a FITS file corresponding to a 3D data cube of $5\times4\times3$  elements 
can be created and written as a series of 2D slices:
\begin{verbatim}
// Make sure FitsIOServer module is initialised :
FitsIOServerInit();
// Open/create the output file, named mycube.fits (overwrite it)
FitsInOutFile fos("!mycube.fits", FitsInOutFile ::Fits_Create);
// Define the cube size 
LONGLONG size[3] = {5,4,3};
// Create a 2D array corresponding to an x-y cube slice 
TArray<int_4> slice(size[0], size[1]);
// Create an IMAGE HDU: cube size and slice data type
fos.CreateImageHDU(FitsTypes::ImageType(slice(0,0)), 3, size);
// initialize the array content 
slice = RegularSequence();
// Create a Fits handler for the array 
FitsArrayHandler<int_4> fh(slice);
// Write the 3 slices to the file 
sa_size_t offset[3];
offset[0]=0; offset[1]=0; 
for(int k=0;k<size[2]; k++) {
  slice += (int_4)(k*100);
  offset[2]=k;
  fh.WriteAtOffset(fos,offset);
}
\end{verbatim}
The example below show how to read arrays from Image HDU in chunk:
\begin{verbatim}
// Make sure FitsIOServer module is initialised :
FitsIOServerInit();
// Open the existing file mycube.fits for reading
FitsInOutFile fis("mycube.fits", FitsInOutFile::Fits_RO);
// Get the array size
int naxis=BASEARRAY_MAXNDIMS;
LONGLONG axes[BASEARRAY_MAXNDIMS];
fis.GetImageHDUInfo(naxis, axes);
// Read in the data as 2D slices 
TArray<int_4> slice(axes[0], axes[1]);
FitsArrayHandler<int_4> fh(slice);
sa_size_t offset[3];
offset[0]=0; offset[1]=0; 
for(int k=0;k<axes[2]; k++) {
  offset[2]=k;
  fh.ReadAtOffset(fis,offset);
  cout << " Slice k=" << k << " Sx=" << axes[0] << " Sy=" << axes[1] << endl;
  slice.Print(cout,slice.Size());  cout << endl;
}
\end{verbatim}
 
\subsection{SwFitsDataTable and other classes}
\label{SwFitsDataTable}
\index{SwFitsDataTable} 
The {\bf SwFitsDataTable} class implements the BaseDataTable interface
using a FITS file as swap space. Compared to SwPPFDataTable, they can be 
used in R/W mode (reading from the table, when it is being created / filled).
They can be used in a way similar to DataTable and SwPPFDataTable. 
When creating the table, a {\tt FitsInOutFile } stream, opened for writing has
to be passed to the creator. No further operation is needed. 
\begin{verbatim}
// ....
FitsInOutFile so("!myswtable.fits", FitsInOutFile::Fits_Create);
SwFitsDataTable dt(so, 16);
// define table columns
dt.AddIntegerColumn("X0_i");
dt.AddFloatColumn("X1_f");
// ... Fill the table
r_8 x[5];
for(int i=0; i<63; i++) {
    x[0] = (i%9)-4.;  x[1] = (i/9)-3.;  
    dt.AddLine(x);
}
\end{verbatim}
The class {\bf FitsBTNtuIntf } provide an alternative tool to read FITS tables. 
{\bf FitsABTColRd} , {\bf FitsABTWriter } and {\bf FitsImg2DWriter } can also 
be used to manipulate FITS files.
\par 
The {\bf scanfits} program can be used to check FITS files and analyse their
content (See \ref{scanfits}).

%%%%   
\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  using the FFTW package (available from {\tt http://www.fftw.org}).
Depending on the configuration options during library build, only 
transforms on double precision data might be available.
Usage example:
\begin{verbatim}
int nx=500, ny=500;
TArray<r_8> img(nx, ny);
// fill in the image 
....
// Compute Fourier transform of the image
TArray< complex<r_8> > fourimg;
FFTWServer ffts;                     
ffts.setNormalize(true);  
ffts.FFTForward(img, fourimg);
// perform some filtering in Fourier domain
// modifying the fourier coefficients fourimg
...
// Compute back the new (filtered) image
TArray<r_8> img2(nx, ny);
ffts.FFTBackward(fourimg, img2);
\end{verbatim}

\section{Multi-thread and parallel programming with SOPHYA} 
 \label{mtpsop}
A general introduction to parallel programming model in the Shared Memory multi-Processor
systems and the POSIX thread model is well beyond the scope of this 
document. Interested readers can find some introductory texts on the following 
web pages:
\begin{enumerate}
\item Symmetric multiprocessing, \\
\href{http://en.wikipedia.org/wiki/Symmetric\_multiprocessing } {http://en.wikipedia.org/wiki/Symmetric\_multiprocessing }
\item POSIX Threads Tutorial {\it  Mark Hays, Software Interest Group }\\
\href{ http://math.arizona.edu/\~swig/documentation/pthreads/ } {http://math.arizona.edu/\~swig/documentation/pthreads/} 
\item POSIX Threads Programming {\it Blaise Barney, Lawrence Livermore National Laboratory } \\
\href{ https://computing.llnl.gov/tutorials/pthreads/ } {https://computing.llnl.gov/tutorials/pthreads/  } 
\end{enumerate} 
In this section, we present some tips and example, focusing on the use of SOPHYA thread management classes. 
Few examples of multithread programs using these classes can be found in the {\bf Tests} module : \\
\hspace{10mm} {\tt zthr.cc  , tmtdt.cc , tmtrnd.cc }

\subsection{General guidelines} 
\begin{enumerate}
\item In general, most classes in SOPHYA are inherently thread-safe, as they do not have global (static) data class 
members and methods operate on instance data members. This is in particular the case of data objects 
in SOPHYA: \\
\tcls{NDataBlock}. \tcls{TArray}, Histo, DataTable, NTuple, \tcls{SphericalMap} \ldots \\
This means that different instances of these classes can be used within threads without any precaution. 
The SOPHYA reference sharing mechanism being itself thread safe. 
\item Once the handlers have been registered, the SOPHYA persistence mechanism is in principle thread safe, 
as long as I/O is performed on different files in different threads. Synchronisation through mutex should
be implemented if the same PPF file (PInPersist/POutPersist object) is being used in different threads. 
\item If you want to use the same NTuple or DataTable or SwPPFDataTable object in different threads, you have 
to call the method {\tt SetThreadSafe()} on the corresponding object:
\begin{verbatim}
---> Set fg=true for thread safe operations 
void NTuple::SetThreadSafe(bool fg);
void BaseDataTable::SetThreadSafe(bool fg);  
\end{verbatim}
\item Some computation functions or classes in SOPHYA use currently global data and are not thread-safe.
FFTPackServer, FFTWServer and SphericalTransform servers are not thread safe.  FFTWServer is not currently 
thread safe as the creation of fftw\_plan is not thread safe in fftw library. 
\item the cfitsio library is not thread safe, making the services of the FitsIOServer NON thread-safe 
\item Many functions in the standard libraries are NOT thread-safe. This is the case of the usual random generators.
Refer to section \ref{randgen} for using random generators in multi-thread applications.
\end{enumerate}

\subsection{Using ZThread and ZMutex }
\begin{itemize}
\item Define and implement a class inheriting from the {\bf ZThread} class (zthread.h). 
This class should overwrite  the {\tt virtual void ZThread::run() } method of the base class 
which has to perform the actual computation. 
At the end of the computation in the {\tt run()} method, the {\tt setRC(int rc) } can be called to set 
the return code of the thread object. 
\item If the actual computation is performed by a function having a single {\tt void *} argument, 
the ZThread base class can directly be used.  See the example in paragraph \ref{thrcls}. 
\item The actual  thread is created by the {\tt start()} method. This method has to be called 
on each  ZThread object to perform the computation. The {\tt start() } method return immediately 
after the creation of the underlying thread object. 
\item The method {\tt join() } can be used to wait for the thread termination.
\item {\bf ZMutex} objects should be used to control to data (memory location) shared 
by several threads.  These ZMutex objects would usually be attributes (data members) of the 
ZThread objects or classes used by the ZThread objects.  Consider a variable {\tt int count }
which is accessed by two different threads. Thread A changes (increment or decrement for example) 
the value of {\tt count}, while thread B waits until {\tt count} reaches a given threshold, before 
performing some actions and resetting {\tt count}. Access to {\tt count} should be protected by 
a ZMutex object; The actual code would then look like:
\begin{verbatim}
// count and its associated ZMutex declaration 
int count = 0;
int threshold=10;
ZMutex mex;
.......
// ----- Code executed by thread A 
mex.lock();
count++;
mex.unlock();
mex.signal();
// Alternatively, mex.broadcast() can be called to wake up 
// several threads waiting on the ZMutex condition variable
.......
// ----- Code executed by thread B
mex.lock();
while (count<threshold) mex.wait();
// here, count has reached or exceeded the threshold
// Do some action, and reset count
count=0;
mex.unlock();
\end{verbatim} 
\end{itemize}

\subsection{Using ParallelExecutor }
The {\bf ParallelExecutor } class (parlex.h) objects can be used to perform a computation divided into 
parts which can be performed independently and in parallel in different threads. A typical 
would be the application of a complex function to all elements in an array. 
\begin{itemize}
\item Define and implement a class inheriting from the {\bf ParallelTaskInterface} class. This class should 
implement the pure virtual method {\tt int  execute(int tid)}. {\tt tid} is the thread identifier, or rank in the 
set of parallel threads assigned to perform the computation. The first thread has a rank=0 and 
{\tt ParallelTaskInterface::getNbParallelThreads() } method can be used to find the total 
number of threads assigned to the task in the parallel execution context. 
\item Create an {\bf ParallelExecutor} object passing to it the ParallelTask object and the number of
threads, or alternatively, a vector of ParallelTask object pointers. 
\item The set of threads are started by calling the {\tt ParallelExecutor::start() } method. The created threads
will be in a waiting state, until the {\tt ParallelExecutor::execute() } is called.
\item The call to {\tt execute() } triggers the execution of the {\tt ParallelTask::execute(tid) } methods
for all of the parallel threads in the context. The {\tt execute() }  method returns when all the threads
have finished executing the {\tt ParallelTask::execute(tid) } .
\item The threads will then be again in a waiting state, ready for a new call to {\tt execute() }. 
The threads will terminate when the  ParallelExecutor object is destroyed.
\end{itemize}

\begin{verbatim}
// Class MyTask implements ParallelTaskInterface
MyTask ptask(...);
// Create a parallel execution context with 4 threads
ParallelExecutor parex(ptask,4);
// Start the threads 
parex.start();
// Loop over several data sets that 
int NDS=5;
for(int i=0; i<NDS; i++ {
  // Prepare the data 
  ...
  // Do the parallel computation 
  parex.execute();
} 
\end{verbatim}

\newpage
\section{Building and installing Sophya}
\subsection{Supported platforms}
Presently, the SOPHYA and PI libraries and tools have been compiled and tested with 
the following compiler/platform pairs:

\begin{center}
\begin{tabular}{|l|l|}
\hline
OS & compiler \\
\hline
Linux  (SL5, 64)          &     g++  4.1  \\
Linux  (SL5, 64)          &     icc - Intel compiler (10.1)   \\
Linux Ubuntu             &    g++ 4.4 \\
MacOSX/Darwin 10.3 (PowerPC) \, 10.4  (Intel) &   Apple g++ 3.3  \, 4.0 \\
MacOSX/Darwin 10.5 (Intel)  &  Apple   g++ 4.0 \\
MacOSX/Darwin 10.6 (Intel)   Snow Leopard &  Apple  g++ 4.2 \\
IBM AIX 5.3      &     xlC   8.0   \\
HP/Compaq/DEC Tru64 ( OSF1)  &     cxx  6.1 \, 6.3  \\
\hline
\end{tabular}
\end{center}

The SOPHYA and PI library and associated tools (spiapp, runcxx \ldots) were ported and tested 
on the Silicon Graphics (SGI)  IRIX system and the corresponding native c++ compiler ( IRIX64, CC  7.3 ),
up to SOPHYA V2.1 (V\_Nov2007). 

\subsection{Library and makefile structure}
% 
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}).
Alternatively, it is possible to group all modules in a single shared 
library {\tt libAsophyaextPI.so}.
\par 
Each library module has a {\tt Makefile} which compiles the source files
and build the correspond static (archive) library using the compilation
rules and flags defined in \\ 
\hspace*{5mm} {\tt \$SOPHYABASE/include/sophyamake.inc}. \\
Each program module has a {\tt Makefile} which compiles and link the 
corresponding programs using the compilation rules and libraries
defined in {\$SOPHYABASE/include/sophyamake.inc}. 
The top level Makefile in BuildMgr/ compiles each library modules
and builds shared libraries. 
\par
Some of the modules in the Sophya package uses external libraries. The 
{\bf FitsIOServer} is an example of such a module, where the {\tt libcfitsio.a}
is used. The list of all Sophya modules using external libraries is
presented in section \ref{sopmodules}.
The external libraries should be installed before the configure step
(see below) and the compilation of the corresponding Sophya modules. 
\par 
The series of Makefiles use the link to {\tt sophyamake.inc} in BuildMgr.
There are also the {\tt smakefile} series which uses the explicit path, using
{\tt \$SOPHYABASE} environment variable.

\subsection{Build and install procedure} 
\label{build}
The build procedure has two main steps: 
\begin{enumerate}
\item The configure step (BuildMgr/configure) setup the directory structure and
the necessary configuration file. Refer to section \ref{directories} for 
the description of SOPHYA directory tree and files.
\item The make step compiles the different sources files, create the library and optionaly
builds all or some of the associated executables. 
\end{enumerate}

{\tt BuildMgr/configure } is a c-shell script with a number of arguments. Please note
 that this is NOT an autoconf generated script, but a custom c-shell script and it does 
not create the makefiles.  The c-shell (or tcsh) must be accessible to run this script.  
\begin{verbatim}
csh> ./configure -h
configure [-sbase SOPHYABASE] [-scxx SOPHYACXX] [-incln] 
  [-minc mymake.inc]  [-compopt 'cc/cxxOptions'] 
  [-arch64] [-arch32] [-sasz64] [-ldble128] 
  [-nofpic] [-nothsafe] [-boundcheck] [-sodebug]
  [-extp dir1 -extp dir2 ...] [-extip dir1 -extip dir2 ... ]
  [-extlp dir1 -extlp dir2 ... ] [-followlink] 
  [-noextlib] [-noext fits] [-noext fftw] [-noext lapack] [-noext astro]
  [-alsofftwfloat] [-usefftw2] [-uselapack2] 
  [-wgrdl] [-epip mdir1 -epip mdir2 ...]  [-noPI] [-slballinone]
  (See SOPHYA manual/web pages for a detailed description of configure options)
\end{verbatim}
\begin{itemize}
\item General parameters and options: \\[1mm]
{\bf -sbase {\it path} }  define SOPHYA installation base directory. \$SOPHYABASE is used 
if not specified \\
{\bf -scxx {\it cmd} } selects the C++ compiler. \$SOPHYACXX  is used if not specified. \\
{\bf -incln }  creates symbolic link for include files, instead of copying them. \\
{\bf -minc {\it filename} } give an explicit name for the file used as the seed to create  the file \\
{\tt \$SOPHYABASE/include/sophyamake.inc}. \hspace{2mm}
If  -minc is not specified, one of the files in BuildMgr/ directory is selected, based on the 
system name and the compiler: \\ 
{\tt  Linux\_g++\_make.inc , OSF1\_cxx\_make.inc , AIX\_xlC\_make.inc \ldots}  \\
{\bf -compopt {\it opt} } can be used to specify additional compiler options 
\item Conditional compilation flags related to SOPHYA configuration (see {\tt machdefs.h} )\\[1mm]
{\bf -arch64 }  select/force 64 bits compilation mode. This is the default on 64 bits Linux systems and AIX. \\
{\bf -arch32  } select/force 32 bits compilation mode. useful on 64 bits Linux systems or AIX. \\
{\bf -sasz64 } select 8 byte (64 bits) size for array indices, useful if you need to allocate an manipulate 
large arrays, with more than $2^32$ elements along a single array dimension. \\
{\bf -ldble128 } set the flags activating {\tt long double} (128 bits) arrays.  \\
{\bf -nofpic } disable -fPIC (Position Independent Code) generation flag \\
{\bf -nothsafe } disable thread-safety protection procedures in SOPHYA : reference sharing \ldots. 
 ( activate the conditional compilation flag {\tt SO\_NOTHSAFE } ) \\
{\bf -boundcheck } compile SOPHYA with bound checking activated when accessing array elements
( conditional compilation flag {\tt SO\_BOUNDCHECKING } ) \\
{\bf -sodebug } activate conditional compilation flag {\tt SOPHYA\_DEBUG } for debugging the library. 
\item  External libraries \\[1mm]
{\bf -extp {\it path} } Adds the specified path to the search path of the external libraries 
include files and archive library. \\
{\bf -extip {\it path} } Adds the specified path to the search path of the external libraries 
include files. \\
{\bf -extlp {\it path} } Adds the specified path to the search path of the external libraries 
archive (libxxx.a). \\
{\bf -followlink}  follow symbolic links when searching for include files and libraries
(-L option of the find command). default: don't follow symbolic links. \\
{\bf -noextlib } Disable compilation of all modules referencing external libraries. \\
{\bf -noext {\it modname} }  Disable compiling of the specified module (with reference to external
library : {\tt -noext fits , -noext fftw  -noext lapack -noext astro } \\
{\bf -usefftw2 } Use FFTW V2 instead of the default FFTW V3 - A preprocessor
flag will be defined in sspvflags.h \\
{\bf -uselapack2 } Use Lapack V2 istead of the default V3 - A preprocessor
flag will be defined in sspvflags.h \\
{\bf -alsofftwfloat } compile single precision (float) version  of the Fourier 
transform methods (module IFFTW, class FFTWServer). A preprocessor
flag will be defined in sspvflags.h and float version of the FFTW library 
(libfftw3f.a) will be linked with SOPHYA, in addition to the default double
precision library (libfftw3.a). \\
{\bf -slballinone}  Use of  a single shared library for all SOPHYA, PI and 
external library interface modules. A compilation flag 
will be defined in sspvflags.h . See also target {\tt slballinone} below. 
\item PI and piapp specific options  \\[1mm]
{\bf -wgrdl } compile and link piapp with GNU readline library \\
{\bf -epip {\it path} }  Adds the specified path to the search path for include 
files and libraries for Motif, and GNU readline, if applicable. \\
{\bf  -noPI } ignore PI modules 
\end{itemize}

\subsubsection{Configure steps }

The configure script performs the following actions :
\begin{enumerate}
\item Creates directory tree under {\tt \$SOPHYABASE }
\item Copy include files to {\tt \$SOPHYABASE/include/ }  (or creates symbolic link) 
\item Search for external libraries include files and create the necessary links 
in {\tt \$SOPHYABASE/include}
\item Search for external libraries (-lfits \ldots) and add the corresponding 
directories to the library search path, in {\tt sophyamake.inc}
\item Creates the file { \tt \$SOPHYABASE/include/sophyamake.inc }
\item Creates the file {\tt  machdefs.h} from { BaseTools/machdefs\_mkmf.h } and 
{\tt sspvflags.h }
\item Creates the object list files for shared library creation
\end{enumerate}

\subsubsection{Compile and install, an example } 

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 XAstroPack package.

\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
# Step 1.b : Check the generated file  $SOPHYABASE/include/sophyamake.inc
csh> ls -lt *.inc 
csh> more sophyamake.inc 
\end{verbatim}
If necessary, edit the generated file {\tt sophyamake.inc } in order to modify 
compilation flags, library list. The file is rather short and self documented.
\begin{verbatim}
# 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 perform the steps 2.a to 2.f using the targets {\tt all} and {\tt slball} 
defined  in the Makefile
\begin{verbatim}
# Step 2.a ... 2.f
csh> make all slball
\end{verbatim}

It is also possible to group all modules in a single shared library using 
the target {\tt slballinone}. 
\begin{verbatim}
# Step 2.a ... 2.f
csh> make all slballinone
\end{verbatim}

At this step, all libraries should have been made. Programs using
Sophya libraries can now be built:
\begin{verbatim}
# To compile some of the test programs
csh> make basetests
# To compile runcxx , scanppf , scanfits 
csh> make prgutil
# To build (s)piapp (libPI.so is needed)
csh> make piapp 
\end{verbatim}

If no further modification or update of source files is foreseen, it is possible 
to remove all .o files:
\begin{verbatim}
#  To clean $SOPHYABASE/obj directory :
csh> make cleanobj
\end{verbatim}


\subsection{Notes}
\begin{itemize}
\item[{\bf Makefile}] List of top level Makefile build targets 
\begin{verbatim}
> libs extlibs PI = all 
> slb slbext slbpi = slball    (OR = slballinone)
> clean  cleanobj
> tests prgutil prgmap progpi = prgall 
> basetests piapp (ou progpi)  pmixer 
\end{verbatim}
\item[{\bf MacOS X}] A high performance mathematic and signal processing
library, including LAPACK and BLAS is packaged in Darwin/MacOS X (10.3, 10.4) : \\
\hspace*{5mm} {\bf -framework Accelerate}
\item[{\bf Tru64/OSF}] An optimised math library with LAPACK and BLAS might 
optionaly be installed {\bf (-lcxlm) }. On our system, this libray contained Lapack V2.
So we used the LAPACK, as compiled from the public sources, and linked with 
the Tru64 native BLAS.
\item[{\bf  IRIX64}] We used the math library with LAPACK V2 and BLAS 
from SGI : {\bf -lcomplib.sgimath}
\item[{\bf AIX}] There seem to be a problem on AIX when several shared 
libraries are used. We have been able to run SOPHYA programs either 
using static libraries, or a single shared library (libAsophyaextPI.so) 
if extlibs and PI are needed, in addition to stand alone SOPHYA modules.
It has not been possible to link SOPHYA with fortran libraries
\item[{\bf Mgr}] This module contains makefiles and build scripts
that were used in SOPHYA up to version 1.7 (2004) : OBSOLETE.
\end{itemize}

\subsection{Files and scripts in BuildMgr/ }
\begin{itemize}
\item[] {\bf Makefile:} Top level Makefile for building SOPHYA.
{\tt smakefile} is similar to Makefile, except that it uses 
the smakefiles in each module.
\item[] {\bf mkmflib:}  c-shell script for creation of library module 
Makefile / smakefile. \\
\hspace*{5mm} {\tt  ./mkmflib -sbase /tmp/sbase SUtils } 
\item[] {\b mkmfprog:} 
c-shell script for creation of programs module Makefile / smakefile \\
\hspace*{5mm}  {\tt ./mkmfprog -sbase /tmp/sbase ProgPI }
\item[] {\bf domkmf:} c-shell script - calls mkmflib for all modules \\
\hspace*{5mm} {\tt ./domkmf -sbase /tmp/sbase}
\item[] {\bf xxx\_make.inc:} Configuration files for different compilers and OS
{\tt ( Linux\_g++\_make.inc , OSF1\_cxx\_make.inc \ldots )}. 
These files are used to generate {\tt sophyamake.inc} 
\end{itemize}

\subsection{Test programs}
The module {\bf Tests} contains a number of test programs:
\begin{enumerate}
\item  {\bf tobjio } tests part of SOPHYA persistence (PPF files) 
\begin{verbatim}
csh> tobjio A
 ===> check the created file tobjio.ppf using scanppf 
csh> scanppf tobjio.ppf
csh> tobjio B
 ===> OK if RC=0 
\end{verbatim}
\item {\bf arrt } and {\bf carrt} perform some checks on TArray modules, including I/O, 
sub-array extraction, array conversion and memory organisation, as well as array operations.
\begin{verbatim}
csh> arrt check 
csh> carrt ac
csh> carrt oso 
csh> carrt odo  
 ===> OK if RC=0 
csh> carrt ace
 ===> OK if RC=9  (throw an exception)
\end{verbatim}
\item {\bf spar } is a TArray performance measurement tool 
\begin{verbatim}
csh> time spar 2 10 800 1000 
csh> time spar 2 1000 80 100 
\end{verbatim}
\item {\bf tssqmx } checks DiagonalMatrix, SymmetricMatrix and LowerTriangularMatrix classes
\begin{verbatim}
csh> tssqmx d 
csh> tssqmx s
csh> tssqmx t
 ===> OK if RC=0 
\end{verbatim}
\item  {\bf ttimestamp } and {\bf tnt } perform some checks on various SOPHYA objects ( DVList, TimeStamp, NTuple \ldots)
\begin{verbatim}
csh> ttimestamp 
 ===> OK if RC=0 
csh> tnt d 
 ===> DVList test, check the output
csh> tnt n 
 ===> NTuple test, check the output
csh> tnt DT 
 ===> DataTable test, check the output
csh> tnt DT 
 ===> SwPPFDataTable test, check the output
\end{verbatim}
\item {\bf tfft } checks the operation of FFTPackServer and FFTWServer in 1D transforms
\begin{verbatim}
# FFTPackServer tests in double 
csh> tfft 1531 P D 0 50 0.0002
csh> tfft 1220 P D 0 50 0.0002
 ===> OK if RC=0 
# FFTPackServer tests in float 
csh> tfft 1531 P F 0 50 0.0002
csh> tfft 1220 P F 0 50 0.0002
 ===> OK if RC=0 
# FFTWServer tests in double 
csh> tfft 1531 W D 0 50 0.0002
csh> tfft 1220 W D 0 50 0.0002
 ===> OK if RC=0 
\end{verbatim}
\item {\bf lpk} checks the LinAlg interface module
\begin{verbatim}
csh> lpk inverse 100,100 0
 ===> OK if RC=0 
csh> lpk linsolve 100,100 0
 ===> OK if RC=0 
csh> lpk lss 100,50 0
 ===> OK if RC=0 
csh> lpk svd 100,50 0
 ===> OK if RC=0 
## These checks can be combined in a single call to lpk 
csh> lpk all 100,50 0
\end{verbatim}
\item {\bf tmtdt } tests several aspects of SOPHYA : threads (ZThread, ZMutex \ldots), DataTables \ldots 
\begin{verbatim}
csh> tmtdt NT 50000 2
 ===> OK if RC=0 
csh> tmtdt DT 50000 2
 ===> OK if RC=0 
csh> tmtdt SWPPF 50000 2
 ===> OK if RC=0 
# cfitsio is NOT thread safe, test SwFitsDataTable with a single thread
csh> tmtdt SWFITS 50000 1
\end{verbatim}

\end{enumerate}

\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. From version V2.2 , all SOPHYA exceptions inherit from 
the standard C++ exception class ( {\bf std::std::exception}). The method {\tt const char* what()} is 
thus implemented and returns the corresponding error message.
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}