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

Last change on this file since 3851 was 3839, checked in by ansari, 15 years ago

MAJ sophya.tex pour V2.2 en cours, Reza 09/08/2010

File size: 116.6 KB
Line 
1\documentclass[twoside,10pt]{article}
2% Package standard : Utilisation de caracteres accentues, mode francais et graphique
3\usepackage{url}
4\usepackage[latin1]{inputenc}
5\usepackage[T1]{fontenc}
6\usepackage[english]{babel}
7\usepackage{graphicx}
8% package a mettre pour faire du pdf
9\usepackage{palatino}
10
11% Extension de symboles mathematiques
12\usepackage{amssymb}
13
14% Definition pour Docs Sophya
15\usepackage{defsophya}
16
17% Constitution d'index
18\usepackage{makeidx}
19
20\usepackage[ps2pdf,bookmarks,bookmarksnumbered,%
21 urlcolor=blue,citecolor=blue,linkcolor=blue,%
22 pagecolor=blue,%hyperindex,%
23 colorlinks=true,hyperfigures=true,hyperindex=true
24 ]{hyperref}
25
26\makeindex % Constitution d'index
27
28\newcommand{\rond}{$\bullet \ $}
29\newcommand{\etoile}{$\star \ $}
30\newcommand{\cercle}{$\circ \ $}
31\newcommand{\carre}{$\Box \ $}
32
33
34\begin{document}
35
36\begin{titlepage}
37% The title page - top of the page with the title of the paper
38\titrehp{Sophya \\ An overview }
39% Authors list
40\auteurs{
41R. Ansari & ansari@lal.in2p3.fr \\
42E. Aubourg & aubourg@hep.saclay.cea.fr \\
43G. Le Meur & lemeur@lal.in2p3.fr \\
44C. Magneville & cmv@hep.saclay.cea.fr \\
45}
46% \auteursall
47% The title page - bottom of the page with the paper number
48\vspace{1cm}
49\begin{center}
50{\bf \Large Sophya Version: 2.2 (V\_Sep2010) }
51\end{center}
52\titrebp{1}
53\end{titlepage}
54
55\tableofcontents
56
57\newpage
58
59\section{Introduction}
60
61{\bf SOPHYA} ({\bf SO}ftware for {\bf PHY}sics {\bf A}nalysis)
62is a collection of C++ classes designed for numerical and
63physics analysis software development. Our goal is to provide
64easy to use, yet powerful classes which can be used by scientists.
65Although some of the SOPHYA modules (SkyMap, Samba, SkyT)
66have been designed with the specific goal of CMB data analysis, most
67modules presented here have a much broader scope and can be
68used in scientific data analysis and modeling/simulation.
69Whenever possible, we use existing numerical packages and libraries,
70encapsulating them in classes in order to facilitate their usage.
71\par
72Our main requirements in designing SOPHYA classes can be summarized as
73follow:
74\begin{itemize}
75\item[\rond] Provide a comprehensive set of data containers, such as arrays and tables
76(tuple) covering the most common needs in scientific simulation and data analysis
77softwares.
78\item[\rond] Take advantage of the C++ language and define methods and operators
79for most basic operation, such as arithmetic operations, in a rather intuitive way, while
80maintaining performances comparable to low level coding in other languages
81(C, Fortran, F90 \ldots)
82\item[\rond] Simplify memory management for programmers using the class library.
83This has been a strong requirement for most SOPHYA classes. Automatic reference
84sharing and memory management is implemented in SOPHYA classes intended
85for large size objects. We recommend to allocate SOPHYA objects on the stack,
86including when objects are returned by methods or functions.
87See section \ref{memgt} for more information.
88\item[\rond] Archiving, importing (reading) and exporting (writing) data in a
89efficient and consistent way is a major concern in many scientific software
90and projects. SOPHYA provide a native data I/O or persistence system,
91(PPF, \ref{ppfdesc}) as well as import/export services for ASCII and FITS formats.
92\end{itemize}
93
94% \vspace*{2mm}
95This documents
96presents only a brief overview of the class library,
97mainly from the user's point of view. A more complete description
98can be found in the reference manual, available from the SOPHYA
99web site: % {\bf http://www.sophya.org}.
100\href{http://www.sophya.org}{http://www.sophya.org}.
101%%%
102%%%
103\subsection{Acknowlegments}
104Many people have contributed to the development SOPHYA and/or the PI library
105and (s)piapp interactive analysis tool, in particular:
106
107\begin{tabular}{lcl}
108Reza Ansari & \hspace{5mm} & (LAL-Univ.Paris Sud, Orsay) \\
109Eric Aubourg & & (DAPNIA-CEA/APC, Saclay) \\
110Sophie Henrot-Versille & & (LAL-IN2P3/CNRS, Orsay) \\
111Alex Kim & & (LBL, Berkeley) \\
112Guy Le Meur & & (LAL-IN2P3/CNRS, Orsay) \\
113Eric Lesquoy & & (DAPNIA-CEA, Saclay) \\
114Christophe Magneville & & (DAPNIA-CEA, Saclay) \\
115Bruno Mansoux & & (LAL-IN2P3/CNRS, Orsay) \\
116Olivier Perdereau & & (LAL-IN2P3/CNRS, Orsay) \\
117Nicolas Regnault & & (LPNHE-IN2P3/CNRS, Paris) \\
118Benoit Revenu & & (APC/Univ.Paris 7, Paris) \\
119Francois Touze & & (LAL-IN2P3/CNRS, Orsay) \\
120\end{tabular}
121
122We thank also the persons who have helped us by useful suggestions, among others : \\
123S. Bargot, S. Plasczczynski, C. Renault and D. Yvon.
124
125\subsection{SOPHYA modules}
126\label{sopmodules}
127The source directory tree
128\footnote{ CVS: cvsserver.lal.in2p3.fr:/exp/eros/CVSSophya}
129is organised into a number of modules.
130
131\begin{itemize}
132\item[] {\bf BuildMgr/} Scripts for code management,
133makefile generation and software installation
134\item[] {\bf BaseTools/} General architecture support classes such
135as {\tt PPersist, NDataBlock<T>}, and few utility classes
136such as the dynamic variable list manager ({\tt DVList}) as well
137as the basic set of exception classes used in SOPHYA.
138\item[] {\bf TArray/} template numerical arrays, vectors and matrices \\
139({\tt TArray<T> TMatrix<T> TVector<T> } \ldots)
140\item[] {\bf NTools/} Some standard numerical analysis tools
141(linear, and non linear parameter fitting, FFT, \ldots)
142\item[] {\bf HiStats/} Histogram-ming and data set handling classes (tuples) \\
143({\tt Histo Histo2D NTuple DataTable} \ldots)
144\item[] {\bf SkyMap/} Local and full sky maps, and some 3D geometry
145handling utility classes. \\
146({\tt PixelMap<T>, LocalMap<T>, SphericalMap<T>, \ldots})
147\item[] {\bf SUtils/} This module contains few utility classes, such as the
148{\tt DataCard} class, as well as string manipulation functions in C and C++.
149\item[] {\bf SysTools/} This module contains classes implementing
150an interface to various OS specific services, such
151threads and dynamic link/shared library handling.
152
153\end{itemize}
154
155The modules listed below are more tightly related to the
156CMB (Cosmic Microwave Background) data analysis problem:
157\begin{itemize}
158\item[] {\bf SkyT/}
159classes for spectral emission and detector frequency response modelling \\
160({\tt SpectralResponse, RadSpectra, BlackBody} \ldots). This module needs extensive
161checking and corrections.
162\item[] {\bf Samba/} Spherical harmonic analysis, noise generators \ldots
163\end{itemize}
164
165The following modules contain the interface classes with
166external libraries:
167\begin{itemize}
168\item[] {\bf FitsIOServer/} Classes for handling file input-output
169in FITS format using the cfitsio library.
170FITS is maintained by NASA and SAO and is available from: \\
171\href{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html}
172{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html}
173\item[] {\bf LinAlg/} Interface with Lapack linear algebra package.
174Lapack is a linear algebra package and can be downloaded from: \\
175\href{http://www.netlib.org/lapack/}{http://www.netlib.org/lapack/}
176\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a).
177FFTW is a package for performing Fourier transforms, written in C.
178Documentation and source code can be found at: \\
179\href{http://www.fftw.org/}{http://www.fftw.org/}
180\item[] {\bf XAstroPack/} Interface to some common astronomical
181computation libraries. Presently, this module uses an external library
182extracted from the {\bf Xephem } source code. The corresponding source
183code is also available from SOPHYA cvs repository, module {\bf XephemAstroLib}.
184Information on Xephem can be found at : \\
185\href{http://www.clearskyinstitute.com/xephem/}{http://www.clearskyinstitute.com/xephem/}
186
187\end{itemize}
188
189The following modules contain each a set of related programs using the
190SOPHYA library.
191\begin{itemize}
192\item[] {\bf Tests/} Simple test programs. Many of these test programs can be also used
193as examples for using SOPHYA.
194\item[] {\bf PrgUtil/} Various utility programs (runcxx, scanppf, scanfits, \ldots)
195\item[] {\bf PrgMap/} Programs performing operations on skymaps: projections,
196power spectrum in harmonic space, \ldots
197\end{itemize}
198
199As a companion to SOPHYA, the {\bf (s)piapp} interactive data analysis
200program is built on top of SOPHYA and the {\bf PI} GUI class library
201and application framework. The {\bf PI} ({\bf P}eida {\bf Interactive})
202development started in 1995, in the EROS \footnote{EROS: {\bf E}xp\'erience
203de {\bf R}echerche d'{\bf O}bjets {\bf S}ombres - http://eros.in2p3.fr}
204microlensing search collaboration, with PEIDA++ \footnote {PEIDA++:
205The EROS data analysis class library -
206http://www.lal.in2p3.fr/recherche/eros/PeidaDoc/}.
207The {\bf PI} documentation and the {\bf piapp} user's guide are available
208from \href{http://www.sophya.org}{http://www.sophya.org}.
209%\href{http://www.sophya.org}{http://www.sophya.org}.
210The {\bf PI} is organized as the following modules:
211\begin{itemize}
212\item[] {\bf PI/} Portable GUI class library and application development
213framework kernel.
214\item[] {\bf PIGcont/} Contour-plot drawing classes.
215\item[] {\bf PIext/} Specific drawers and adapters for SOPHYA objects,
216and the {\bf piapp} interactive data analysis framework.
217\item[] {\bf ProgPI/} interactive analysis tool main program and pre-loaded
218modules.
219\end{itemize}
220
221Modules containing examples and demo programs and scripts:
222\begin{itemize}
223\item[] {\bf Examples/} Sample SOPHYA codes and example programs and
224makefiles.
225\item[] {\bf DemoPIApp/} Sample scripts and programs for (s)piapp
226interactive analysis tools.
227\end{itemize}
228
229The following modules contains additional programs or libraries, based on SOPHYA :
230\begin{itemize}
231\item[] {\bf AnaLC} contains program files extracted from the EROS/PEIDA software
232and adapted to be compiled with SOPHYA. It can be used in particular to read and analyse
233EROS light curve data files.
234\item[] {\bf LUC} {\bf L}ittle {\bf U}niverse {\bf C}alculator is a module containing classes to
235perform basic computation related to the universe geometry (FRW metric).
236\item[] {\bf PMixer/} skymixer and related programs. This module is {\bf obsolete}.
237\end{itemize}
238
239\newpage
240
241\section{Using Sophya}
242The organisation of SOPHYA directories and some of the associated
243utility programs are described in this section.
244Basic usage of Sophya classes is described in the following sections.
245Complete Sophya documentation can be found at our web site
246{\bf http://www.sophya.org}.
247
248\subsection{Directories, environment variables, configuration files}
249\label{directories}
250The environment variable {\bf SOPHYABASE} is used
251to define the path where the Sophya libraries and binaries are installed.
252\begin{itemize}
253\item \$SOPHYABASE/include : Include (.h) files
254\item \$SOPHYABASE/lib : Path for the archive libraries (.a)
255\item \$SOPHYABASE/slb: Shared library path (.so or .dylib on Darwin/MacOS)
256\item \$SOPHYABASE/exe : Path for binary program files
257\end{itemize}
258
259The directory { \tt \$SOPHYABASE/include/SophyaConfInfo/ } contains files
260describing the installed configuration of SOPHYA software.
261
262The file { \tt \$SOPHYABASE/include/machdefs.h } contains definitions
263(flags, typedef) used in SOPHYA, while some more specific flags,
264are found in { \tt \$SOPHYABASE/include/sspvflags.h }
265
266The file { \tt \$SOPHYABASE/include/sophyamake.inc } contains the
267compilation commands and flags used for building the software.
268Users can use most of compilation and link commands defined in this file:
269 {\tt \$CCOMPILE , \$CXXCOMPILE . \$CXXLINK \ldots}.
270 (See module Example).
271
272The configure script (BuildMgr/configure) creates the directory tree and the
273above files. It also copy (or create symbolic link) for all SOPHYA include
274files as well as symbolic links for external libraries
275include files path in {\tt \$SOPHYABASE/include} (FitsIO, FFTW, XAstro \ldots).
276
277Object files for each module are grouped in a static archive library
278by the build procedure (libXXX.a for module
279XXX, with XXX = BaseTools, TArray, HiStats, FitsIOServer \ldots).
280
281When shared libraries are build, all stand alone SOPHYA modules
282are grouped in {\tt libsophya.so}, {\tt libextsophya.so} contains
283the interface modules with external libraries {\bf (FitsIOServer, LinAlg \ldots)},
284while {\bf PI, PIext, PIGcont} modules are grouped in {\tt libPI.so}.
285Alternatively, it is possible to group all modules in a single shared
286library {\tt libAsophyaextPI.so} (See \ref{build})
287
288In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
289should contain the Sophya shared library path
290({\tt \$SOPHYABASE/slb}).
291On Silicon Graphics machines with IRIX64 operating system,
292the default SOPHYA configuration correspond to the 64 bit architecture.
293The environment variable { \bf LD\_LIBRARY64\_PATH } replace in
294this case the usual {\bf LD\_LIBRARY\_PATH} variable.
295On IBM machines with AIX, the {\bf LIBPATH} environment variables
296contains the shared libraries search path.
297
298When using the dynamic load services in SOPHYA ({\tt PDynLinkMgr}
299class), in runcxx or (s)piapp applications for example, the shared
300library search path must contain the current working directory (
301dot . in unix).
302
303\subsection{Copy constructor and assignment operator}
304\label{memgt}
305In C++, objects can be copied by assignment or by initialisation.
306Copying by initialisation corresponds to creating an object and
307initialising its value through the copy constructor.
308The copy constructor has its first argument as a reference, or
309const reference to the object's class type. It can have
310more arguments, if default values are provided.
311Copying by assignment applies to an existing object and
312is performed through the assignment operator (=).
313The copy constructor implements this for identical type objects:
314\begin{verbatim}
315class MyObject {
316public:
317 MyObject(); // Default constructor
318 MyObject(MyObject const & a); // Copy constructor
319 MyObject & operator = (MyObject const & a) // Assignment operator
320}
321\end{verbatim}
322The copy constructors play an important role, as they are
323called when class objects are passed by value,
324returned by value, or thrown as an exception.
325\begin{verbatim}
326// A function declaration with an argument of type MyObject,
327// passed by value, and returning a MyObject
328MyObject f(MyObject x)
329{
330 MyObject r;
331 ...
332 return(r); // Copy constructor is called here
333}
334// Calling the function :
335MyObject a;
336f(a); // Copy constructor called for a
337\end{verbatim}
338It should be noted that the C++ syntax is ambiguous for the
339assignment operator. {\tt MyObject x; x=y; } and
340{\tt MyObject x=y;} have different meaning.
341\begin{verbatim}
342MyObject a; // default constructor call
343MyObject b(a); // copy constructor call
344MyObject bb = a; // identical to bb(a) : copy constructor call
345MyObject c; // default constructor call
346c = a; // assignment operator call
347\end{verbatim}
348
349As a general rule in SOPHYA, objects which implements
350reference sharing on their data members have a copy constructor
351which shares the data, while the assignment operator copies or
352duplicate the data.
353
354\subsection{Multi-thread programming with SOPHYA}
355Multi-thread programming is usually safe as long as different threads DO NOT access
356the same memory locations. SOPHYA is mainly organized as classes having only
357data members, with very few cases having static data members (memory locations
358common to all instances of a class), or seldom, global variables. \\
359Using different instances of a class in different threads is thus safe for most
360classes / methods in SOPHYA. \\
361A major exception is the reference sharing mechanism, where different instances
362of a class may shared memory locations. This reference sharing mechanism has been
363made thread-safe in SOPHYA, from version V=2.1. \\
364As a consequence, different execution threads can access non overlapping sub-arrays
365and access (read/write) the corresponding elements without the need of
366mutex and thread synchonisation. Moreover, thread safe filling of
367NTuple and DataTable objects can be activated, on an object per object basis. \\
368The {\tt ZThread, ZMutex \ldots} classes in the {\bf SysTools } module offer a relatively
369easy way of writing multi-threaded programs.
370
371\subsection{the runcxx program}
372\index{runcxx} \label{runcxx}
373{\bf runcxx} is a simple program which can be used to compile, link
374and run simple C++ programs. It handles the creation of a
375complete program file, containing the basic set C++ include files,
376the necessary include files for SOPHYA SysTools, TArray, HiStats
377and NTools modules, and the main program with exception handling.
378Other Sophya modules can be included using the {\tt -import} flag.
379Use of additional include files can be specified using the
380{\tt -inc} flag.
381\begin{verbatim}
382csh> runcxx -h
383 PIOPersist::Initialize() Starting Sophya Persistence management service
384SOPHYA Version 1.9 Revision 0 (V_Mai2005) -- May 31 2005 15:11:32 cxx
385 runcxx : compiling and running of a piece of C++ code
386 Usage: runcxx [-compopt CompileOptions] [-linkopt LinkOptions]
387 [-tmpdir TmpDirectory] [-f C++CodeFileName]
388 [-inc includefile] [-inc includefile ...]
389 [-import modulename] [-import modulename ...]
390 [-uarg UserArg1 UserArg2 ...]
391 if no file name is specified, read from standard input
392 modulenames: SkyMap, Samba, SkyT, FitsIOServer,
393 LinAlg, IFFTW, XAstroPack
394\end{verbatim}
395Most examples in this manual can be tested using runcxx.
396The preprocessor macro {\tt KeepObj()} is defined by runcxx and let the user
397to save an objet to an PPF file, with the same name as the corresponding variable.
398The example below shows how to compile, link and run a sample
399code.
400\begin{verbatim}
401## File example.icc :
402
403Matrix mxa(3,3);
404mxa = IdentityMatrix(1.);
405cout << mxa ;
406// Save the object mxa in a PPF file named mxa
407KeepObj(mxa);
408
409## Executing this sample code
410csh> runcxx -f example.icc
411\end{verbatim}
412
413\subsection{the scanppf program}
414\index{scanppf} \label{scanppf}
415{\bf scanppf} is a simple SOPHYA application which can be used to check
416PPF files and list their contents. It can also provide the list of all registered
417PPF handlers.
418\begin{verbatim}
419csh> scanppf -h
420 PIOPersist::Initialize() Starting Sophya Persistence management service
421SOPHYA Version 2.0 Revision 0 (V_Jul2006) -- Jul 17 2006 14:13:27 cxx
422 Usage: scanppf [flags] filename
423 flags = -s -n -a0 -a1 -a2 -a3 -lh -lho -lmod
424 -s[=default} : Sequential reading of objects
425 -n : Object reading at NameTags
426 -a0...a3 : Tag List with PInPersist.AnalyseTags(0...3)
427 -lh : List PPersist handler classes
428 -lho : List PPersist handler and dataobj classes
429 -lmod : List initialized/registered modules
430\end{verbatim}
431
432\subsection{the scanfits program}
433\index{scanfits} \label{scanfits}
434{\bf scanfits} is a SOPHYA program using the FitsIOServer
435\footnote{FitsIOServer module uses the cfitsio library. scanfits has to be linked with
436with FitsIOServer module and cfitsio libraries, or libextsophya.so}
437module which can be used
438to analyse the content of FITS files. It can list the FITS headers, the appropriate
439SOPHYA-FITS handler (implementing {\tt FitsHandlerInterface}) class, and the list of
440all registered FITS handlers.
441\begin{verbatim}
442csh> scanfits -h
443 PIOPersist::Initialize() Starting Sophya Persistence management service
444SOPHYA Version 2.0 Revision 0 (V_Jul2006) -- Jul 17 2006 14:13:27 cxx
445 Usage: scanfits [flags] filename
446 flags = -V1 -lh -rd -header
447 -V1 : Scan using old (V1) code version
448 -lh : Print the list of registered handlers (FitsHandlerInterface)
449 -rd : try to read each HDU data using appropriate handler
450 -header : List header information
451\end{verbatim}
452
453\subsection{the spiapp program}
454\index{spiapp} \label{spiapp}
455{\bf spiapp} is an interactive data analysis program, built on top of the SOPHYA
456library, , the PI GUI library. The interactive data analysis framework is defined
457by the classes in the {\bf PIext} module. spiapp has a c-shell like
458command interpreter and can be used to manipulate SOPHYA and display
459objects. Refer to the piapp user manual for more information.
460\begin{verbatim}
461csh> spiapp -h
462 SophyaInitiator::SophyaInitiator() BaseTools Init
463 PIOPersist::Initialize() Starting Sophya Persistence management service
464SOPHYA Version 2.1 Revision 0 (V_Nov2007) -- Nov 24 2007 13:08:58 gcc 3.3 20030304 (Apple Computer, Inc. build 1495)
465
466 piapp: Interactive data analysis and visualisation program
467 Usage: piapp [-nored] [-doublered] [-termread] [-term]
468 [-hidezswin] [-small] [-nosig] [-nosigfpe] [-nosigsegv]
469 [-tmpdir TmpDirectory] [-help2tex] [-exec file [args]]
470 -nored : Don't redirect stdout/stderr to piapp console
471 -doublered : Redirect stdout/stderr to piapp console AND the terminal
472 -termread : Read commands on terminal (stdin)
473 -term : equivalent to -nored -termread -small
474 -hidezswin : Hide Zoom/Stat/ColMap window
475 -small : Create small size main piapp window
476 -nosig : Don't catch SigFPE, SigSEGV
477 -nosigfpe -nosigsegv: Don t catch SigFPE / SigSEGV
478 -tmpdir TmpDirectory: defines TMDIR for temporary files
479 -help2tex: Create a LaTeX help file (piahelp.tex)
480 -exec file [args] : Execute command file (last option)
481
482\end{verbatim}
483
484
485\newpage
486\section{Module BaseTools}
487
488{\bf BaseTools} contains utility classes such as
489{\tt DVlist}, an hierarchy of exception classes for Sophya, a template
490class {\tcls{NDataBlock}} for handling reference counting on numerical
491arrays, as well as classes providing the services for implementing simple
492serialisation (object persistence services).
493\vspace*{5mm}
494
495\subsection{Initialisation}
496\index{SophyaInitiator}
497A number of actions have to be taken before
498some of the services provided by SOPHYA become operational. This is the case
499of SOPHYA persistence, as well as FITS I/O facilities.
500Initialisation of many SOPHYA modules is performed through an initialiser class,
501which inherits from {\bf SophyaInitiator}.
502\par
503Static instance of each initialiser class exist in the library and the various SOPHYA services
504should be operational when the user code ({\tt main()}) starts, except for
505modules in the second or third shared libraries
506({\tt libextsophya.so libPI.so}). Indeed, a problem related
507to the initialisation of shared libraries arises on some systems
508(Darwin/Mac OS X in particular) causing program crash at start-up,
509if static instance of initialiser class is present in the second shared library.
510The FitsIOServer module should thus be explicitly initialised in the user
511program.
512\par
513In cases where the run time loader does not perform correctly the static
514object initialisation, the initialiser class for the modules used in the
515program must be instanciated in the beginning of your main program: \\
516{\tt TArrayInitiator , HiStatsInitiator , SkyMapInitiator , FitsIOServer \ldots}
517%%%
518\subsection{SOPHYA persistence}
519\label{ppfdesc}
520\index{PPersist} \index{PInPersist} \index{POutPersist}
521\begin{figure}[hbt]
522\dclsa{PPersist}
523\dclsccc{PPFBinarIOStream}{PPFBinaryInputStream}{PInPersist}
524\dclscc{PPFBinaryOutputStream}{POutPersist}
525\caption{partial class diagram for classes handling persistence in Sophya}
526\end{figure}
527A simple persistence mechanism is defined in SOPHYA. Its main
528features are:
529\begin{itemize}
530\item[] Portable file format, containing the description of the data structures
531and object hierarchy. \\
532{\bf PPF} {\bf P}ortable {\bf P}ersistence file {\bf F}ormat.
533\index{PPF}
534\item[] Handling of read/write for multiply referenced objects.
535\item[] All write operations are carried using sequential access only. This
536holds also for read operations, unless positional tags are used.
537SOPHYA persistence services can thus be used to transfer objects
538through network links.
539\item[] The serialisation (reading/writing) for objects for a given class
540is implemented through a handler object. The handler class inherits
541from {\tt PPersist} class.
542\item[] A run time registration mechanism is used in conjunction with
543RTTI (Run Time Type Identification) for identifying handler classes
544when reading {\bf PInPersist} streams, or for associating handlers
545with data objects {\bf AnyDataObject} for write operations.
546\end{itemize}
547The most useful methods for using Sophya
548persistence are listed below. A brief description of the PPF file format
549and some guidelines for writing writing delegate classes for handling
550object persistence can be found in the following paragraphs.
551\begin{itemize}
552\item[] {\tt POutPersist::PutObject(AnyDataObj \& o)} \\
553Writes the data object {\bf o} to the output stream.
554\item[] {\tt POutPersist::PutObject(AnyDataObj \& o, string tagname)} \\
555Writes the data object {\bf o} to the output stream, associated with an
556identification tag {\bf tagname}.
557\item[] {\tt PInPersist::GetObject(AnyDataObj \& o)} \\
558Reads the next object in stream into {\bf o}. An exception is
559generated for incompatible object types.
560\item[] {\tt PInPersist::GetObject(AnyDataObj \& o, string tagname)} \\
561Reads the object associated with the tag {\bf tagname} into {\bf o}.
562An exception is generated for incompatible object types.
563\end{itemize}
564The operators {\tt operator << (POutPersist ...) } and
565{\tt operator >> (PInPersist ...) } are often overloaded
566to perform {\tt PutObject()} and {\tt GetObject()} operations.
567the {\bf PPFNameTag} (ppfnametag.h) class can be used in conjunction with
568{\tt << >> } operators to write objects with a name tag or to retrieve
569an object identified with a name tag. The example below shows the
570usage of these operators:
571\begin{verbatim}
572// Creating and filling a histogram
573Histo hw(0.,10.,100);
574...
575// Writing histogram to a PPF stream
576POutPersist os("hw.ppf");
577os << PPFNameTag("myhisto") << hw;
578
579// Reading a histogram from a PPF stream
580PInPersist is("hr.ppf");
581is >> PPFNameTag("myhisto") >> hr;
582\end{verbatim}
583
584The {\bf scanppf} program can be used to list the content of a PPF file.
585
586\subsubsection{PPF file format}
587The PPF file consist of :
588\begin{itemize}
589\item[\rond] A file header (3x32=96 bytes) , containing an identification string,
590PPF format version (currently V3), the file endiannes {\tt (BIG-ENDIAN , LITTLE-ENDIAN)}
591and the file creation date and time. It should be noted that the present SOPHYA version
592(V=2.2) does not allow updating an existing PPF file.
593\item[\rond] A collection of tagged data items and data structuring tags.
594the PPF tags are one byte (8 bits) long and may be followed by a length information
595for arrays, or strings. The length information might be 4 bytes or 8 bytes long.
596Characters, various length (1,2,4,8 bytes long) signed and unsigned integers,
597simple or double precision (4/8) bytes floating point and complex numbers are
598the basic data types handled by PPF streams.
599\begin{verbatim}
600tag=PPS_SIMPLE_Integer4 + data (4 bytes)
601tag=PPS_SIMPLE_Float4 + data (4 bytes)
602tag=PPS_SIMPLE_Complex8 + data (2x8=16 bytes)
603tag=PPS_STRING + Length (4 bytes) + data (Length bytes)
604tag=PPS_SIMPLE_ARRAY4_UInt2 + Length (4 bytes) + data (2xLength bytes)
605\end{verbatim}
606The two tags {\tt PPS\_OBJECT} and {\tt PPS\_ENDOBJECT} are used for identifying
607objects.
608\item[\rond] The file trailer contains some global statistics such as the total number
609of objects in file, the list of name tags, with their corresponding positions.
610These informations are grouped under the identifier tag {\tt PPS\_NAMETAG\_TABLE }.
611The last byte in file contains the value {\tt PPS\_EOF}
612\end{itemize}
613We show below the result of the analysis of a PPF file using
614 {\tt PPFBinaryInputStream::AnalyseTags()}. This can easily be done using the
615 {\bf runcxx} program.
616The PPF file contains four objects:
617a {\tt TimeStamp} object, two {\tt NDataBlock} objects, the second one sharing its data with
618the first one. The second NDataBlock is only written as a reference to the first object. The fourth
619object is a {\tt RandomGenerator} which has an {\tt NDataBlock} object as data member.
620Notice the corresponding nested structure of object marker tags.
621The first object and the last objects in the file are written associated with a {\bf NameTag}.
622{\small
623\begin{verbatim}
624 ----------------------------------------------------------
625 PPFBinaryInputStream::AnalyseTags(Level= 49)
626 FileName= toto.ppf
627 Version= 3 FileSize= 8884 Creation Date= Fri Dec 7 15:30:58 2007
628
629 NbPosTag=0 NbNameTag=2 NbObjs=4 NbTopLevObjs=3 NbRefs=1 MaxNest=2
630
631<PPS_NAMETAG_MARK> tag at position 60
632<PPS_OBJECT> tag at position 61 ClassId= c09dfe032b341cee ObjectId= 10
633 <PPS_SIMPLE> tag at position 72 DataType=INTEGER x4
634 <PPS_SIMPLE> tag at position 77 DataType=INTEGER x8
635 <PPS_SIMPLE> tag at position 80 DataType=FLOAT x8
636<PPS_ENDOBJECT> tag at position 89 ObjectId= 10
637<PPS_OBJECT> tag at position 92 ClassId= e670300b367585d1 ObjectId= 21
638 <PPS_SIMPLE_ARRAY4> tag at position a3 DataType=UNSIGNED x8 NElts= 3
639 <PPS_SIMPLE_ARRAY4> tag at position c0 DataType=INTEGER x4 NElts= 50
640<PPS_ENDOBJECT> tag at position 18d ObjectId= 21
641<PPS_REFERENCE> tag at position 196 ObjectId= 21 OrigPos=92
642<PPS_NAMETAG_MARK> tag at position 1a7
643<PPS_OBJECT> tag at position 1a8 ClassId= eb6c427a5a30caed ObjectId= 30
644 <PPS_SIMPLE_ARRAY4> tag at position 1b9 DataType=UNSIGNED x4 NElts= 6
645 <PPS_SIMPLE> tag at position 1d6 DataType=UNSIGNED x8
646 <PPS_SIMPLE> tag at position 1df DataType=UNSIGNED x8
647 <PPS_OBJECT> tag at position 1e8 ClassId= 6531a8f47336d4aa ObjectId= 41
648 <PPS_SIMPLE_ARRAY4> tag at position 1f9 DataType=UNSIGNED x8 NElts= 3
649 <PPS_SIMPLE_ARRAY4> tag at position 216 DataType=FLOAT x8 NElts= 1024
650 <PPS_ENDOBJECT> tag at position 221b ObjectId= 41
651<PPS_ENDOBJECT> tag at position 2224 ObjectId= 30
652<PPS_NAMETAG_TABLE> tag at position 222d
653<PPS_NAMETAG_ENTRY> NameTag=rg-randgen NameTagMark Position=1a7
654<PPS_NAMETAG_ENTRY> NameTag=ts-timestamp NameTagMark Position=60
655<PPS_EOF> tag at position 22b3 TagPos=222d
656 PPFBinaryInputStream::AnalyseTags() - End - Total Number of Tags= 23
657 ----------------------------------------------------------
658\end{verbatim}
659} % fin de small
660
661\subsubsection{Writing PPF handlers}
662Here are some guidelines for creating PPF handler classes to make
663a class persistent through the SOPHYA PPersist mechanism. The example
664discussed here can be found in the {\bf Examples} module, in the directory
665{\bf MyPPF}.
666\begin{enumerate}
667\item The class should inherit from the {\bf SOPHYA::AnyDataObj} class.
668 In the example here, we want to make the class {\bf Vfs} persistent
669\begin{verbatim}
670class Vfs : public SOPHYA::AnyDataObj {
671 ...
672}
673\end{verbatim}
674%%%%
675\item The PPF handler class
676should inherit from {\bf SOPHYA::PPersist} . The pure virtual methods
677of PPersist class (DataObj() , SetDataObj() , ReadSelf(), WriteSelf()
678must be implemented.
679\begin{verbatim}
680class PPFHandlerVfs : public SOPHYA::PPersist {
681public:
682 virtual SOPHYA::AnyDataObj* DataObj();
683 virtual void SetDataObj(SOPHYA::AnyDataObj &);
684protected:
685 virtual void ReadSelf(SOPHYA::PInPersist&);
686 virtual void WriteSelf(SOPHYA::POutPersist&) const;
687}
688\end{verbatim}
689%%%
690It is possible to use the template class {\tt SOPHYA::ObjFileIO<T>}, and
691specialize it for the tharget class (here {\tt SOPHYA::ObjFileIO<Vfs>}),
692by defining the two methods : \\
693\hspace*{5mm} {\tt SOPHYA::ObjFileIO<Vfs>::ReadSelf(SOPHYA::PInPersist\&) } \\
694\hspace*{5mm} {\tt SOPHYA::ObjFileIO<Vfs>::WriteSelf(SOPHYA::POutPersist\&) const }
695%%%%
696\item If it is NOT possible to have the target class inherit from {\bf AnyDataObj},
697a wrapper class should be used. The same class can play the role of wrapper
698AND the PPF handler. See the PPF handler / wrapper class for STL vectors : \\
699\hspace*{5mm} {\tt SOPHYA::PPFWrapperSTLVector<T>} in file BaseTools/ppfwrapstlv.h
700%%%
701\item Implement the ReadSelf() and WriteSelf() methods of the PPF handler.
702All the I/O services from the following classes can be used :
703\begin{itemize}
704\item PPFBinaryIOStrem , PPFBinaryInputStream , PInPersist
705\item PPFBinaryIOStrem , PPFBinaryOutputStream , POutPersist
706\end{itemize}
707Writing and reading of the embeded objects for which a handler has been register
708ed can simply be performed by : \\
709\hspace*{5mm} {\tt POutPersist::PutObject() } \\
710\hspace*{5mm} {\tt PInPersist::GetObject() }
711
712{\bf Warning:} The services associated with nametags : \\
713\hspace*{5mm} {\tt PPFBinaryOutputStream::WriteNameTag() } \\
714\hspace*{5mm} {\tt PPFBinaryInputStream::GotoNameTag() } \\
715are {\bf NOT} intented to be used in WriteSelf() , ReadSelf()
716%%%%%
717\item The new PPF handler, as well as the list of classes it can handle has to be
718registered prior to use PPF read/write for the target classes. This must be
719performed during the initialization phase, for example at the beginning of the
720main() program. Another possibility is to use a module initializer
721(See {\bf SophyaInitiator} class in file BaseTools/sophyainit.h )
722and declare a static instance of the class. Notice that this works only if
723the system loader handles correcly the call of constructor
724for the statically declared objects.
725
726The registration can be performed using the CPP macros defined in
727BaseTools/ppersist.h
728\begin{verbatim}
729 // First, register the PPF handler ObjFileIO<Vfs>
730 PPRegister(ObjFileIO<Vfs>);
731 // Register the list of classes which can be handled by ObjFileIO<Vfs>
732 DObjRegister(ObjFileIO<Vfs>, Vfs);
733\end{verbatim}
734%%%%%%%%%%%%%%%%
735\end{enumerate}
736
737%%%%%%%%%%%%
738\subsection{\tcls{NDataBlock}}
739\index{\tcls{NDataBlock}}
740\begin{figure}[hbt]
741\dclsbb{AnyDataObj}{\tcls{NDataBlock}}
742\dclsbb{PPersist}{\tcls{FIO\_NDataBlock}}
743\end{figure}
744The {\bf \tcls{NDataBlock}} is designed to handle reference counting
745and sharing of memory blocs (contiguous arrays) for numerical data
746types. Initialisation, resizing, basic arithmetic operations, as
747well as persistence handling services are provided.
748The persistence handler class ({\tt \tcls{FIO\_NDataBlock}}) insures
749that a single copy of data is written for multiply referenced objects,
750and the data is shared among objects when reading.
751\par
752The example below shows writing of NDataBlock objects through the
753use of overloaded operator $ << $ :
754\begin{verbatim}
755#include "fiondblock.h"
756// ...
757POutPersist pos("aa.ppf");
758NDataBlock<r_4> rdb(40);
759rdb = 567.89;
760pos << rdb;
761// We can also use the PutObject method
762NDataBlock<int_4> idb(20);
763idb = 123;
764pos.PutObject(idb);
765\end{verbatim}
766The following sample programs show the reading of the created PPF file :
767\begin{verbatim}
768PInPersist pis("aa.ppf");
769NDataBlock<r_4> rdb;
770pis >> rdb;
771cout << rdb;
772NDataBlock<int_4> idb;
773cout << idb;
774\end{verbatim}
775
776\subsection{DVList, MuTyV and TimeStamp classes}
777\index{DVList} \index{MuTyV} \index{TimeStamp}
778\begin{figure}[hbt]
779\dclsa{MuTyV}
780\dclsbb{AnyDataObj}{DVList}
781\dclsbb{PPersist}{\tclsc{ObjFileIO}{DVList}}
782\end{figure}
783The {\bf DVList} class objects can be used to create and manage list
784of values, associated with names. A list of pairs of (MuTyV, name(string))
785is maintained by DVList objects. {\bf MuTyV} is a simple class
786capable of holding string, integer, float or complex values,
787providing easy conversion methods between these objects.
788{\bf MuTyV} objects can also hold {\bf TimeStamp } objects.
789\begin{verbatim}
790// Using MuTyV objects
791MuTyV s("hello"); // string type value
792MuTyV x;
793x = "3.14159626"; // string type value, ASCII representation for Pi
794double d = x; // x converted to double = 3.141596
795x = 314; // x contains the integer value = 314
796// Using DVList
797DVList dvl;
798dvl("Pi") = 3.14159626; // float value, named Pi
799dvl("Log2") = 0.30102999; // float value, named Log2
800dvl("FileName") = "myfile.fits"; // string value, named myfile.fits
801// Printing DVList object
802cout << dvl;
803\end{verbatim}
804
805\begin{figure}[hbt]
806\dclsbb{AnyDataObj}{TimeStamp}
807\end{figure}
808%
809The {\bf TimeStamp} class represent date and time and provides
810many standard operations, such as Initialisation from strings,
811conversion to strings and time interval computations. \\
812Usage example:
813\begin{verbatim}
814// Create a object with the current date and time and prints it to cout
815TimeStamp ts;
816cout << ts << endl;
817// Create an object with a specified date and time
818TimeStamp ts2("01/01/1905","00:00:00");
819// Get the number of days since 0 Jan 1901
820cout << ts2.ToDays() << endl;
821
822// Combined use of TimeStamp and MuTyV
823string s;
824TimeStamp ts; // Current date/time
825MuTyV mvt = ts;
826s = mvt; // s contains the current date in string format
827cout << s << endl;
828\end{verbatim}
829
830\subsection{\tcls{SegDataBlock} , \tcls{SwSegDataBlock}}
831\index{\tcls{SegDataBlock}} \index{\tcls{SwSegDataBlock}}
832%%
833\begin{figure}[hbt]
834\dclsccc{AnyDataObj}{\tcls{SegDBInterface}}{ \tcls{SegDataBlock} }
835\dclscc{\tcls{SegDBInterface}}{ \tcls{SwSegDataBlock} }
836\end{figure}
837\begin{itemize}
838\item[\rond] \tcls{SegDataBlock} handles arrays of object of
839type {\bf T} with reference sharing in memory. The array can be extended
840(increase in array size) with fixed segment size. It implements the interface
841defined by \tcls{SegDBInterface}.
842\item[\rond] \tcls{SwSegDataBlock} Implements the same \tcls{SegDBInterface}
843using a data swapper object. Data swappers implement the interface defined in
844(\tcls{DataSwapperInterface} class. \tcls{SwSegDataBlock} can
845thus be used to handle arrays with very large number of objects.
846These classes handles reference sharing.
847\end{itemize}
848
849 \subsection{Pseudo-random number generators}
850 \index{Random numbers}
851\begin{figure}[hbt]
852\dclsbb{ AnyDataObj }{ RandomGeneratorInterface }
853\vspace*{5mm}
854\dclsccc{RandomGeneratorInterface}{DR48RandGen}{ ThSDR48RandGen }
855\dclsbb{ RandomGeneratorInterface }{ FMTRandGen }
856\end{figure}
857The {\bf RandomGeneratorInterface } is a pure virtual class which
858defines the interface for random number generator classes and implements the
859generation algorithm for few probability distribution functions:
860flat, Gaussian, Poisson, Exponential, 2D gaussian \ldots. The pure virtual method
861{\tt Next() } should be implemented by inheriting classes.
862 \index{RandomGeneratorInterface}
863
864Several implementations using different pseud-random generators are provided by SOPHYA:
865\begin{enumerate}
866\item {\bf DR48RandGen} uses the well known {\tt drand48() } pseudo-random
867sequence generator. It implements also the possibility of initializing the generator state, as well saving and
868restoring the generator state through the SOPHYA PPF mechanism. The {\tt drand48() } and thus
869{\bf DR48RandGen} objects are not NOT thread safe. Different instances of {\bf DR48RandGen} share
870the same underlying 48 bit state.
871\item {\bf ThSDR48RandGen} is a thread safe version of DR48RandGen. Using the constructor flag, in can
872be instanciated to reproduce the non thread-safe behaviour of DR48RandGen.
873Several instances of this class can be used in different threads without the risk of
874corrupting the internal state of the drand48() generator. However, in multi-thread applications,
875 there is no guarantee to obtain the same sequence of numbers in each thread.
876\index{ ThSDR48RandGen}
877\item {\bf FMTRandGen } ids an implementation of RandomGeneratorInterface using
878the SIMD-oriented Fast Mersenne Twister (SFMT). It uses the code developed by
879Mutsuo Saito and Makoto Matsumoto from Hiroshima University. For more information on this
880method, refer to \\
881\href{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}
882{ttp://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html} \\
883The instances of FMTRandGen are thread safe, each instance having its own internal state
884and produces an independent sequence of random numbers.
885\index{FMTRandGen }
886\end{enumerate}
887
888A number of c-like functions are defined in the file BaseTools/srandgen.h and can be used
889to generate random sequences: \\
890{\tt drand01() , drandpm1() , GaussianRand() , PoissonRand \ldots} \\
891These functions use a global instance of RandomGeneratorInterface class which can be
892defined and accessed through the static methods:\\
893\hspace*{5mm} { \tt RandomGeneratorInterface* RandomGeneratorInterface::GetGlobalRandGenP() } \\
894\hspace*{5mm} { \tt void SetGlobalRandGenP(RandomGeneratorInterface* rgp) } \\
895
896In multi-threaded programs, an instance of ThSDR48RandGen or FMTRandGen
897should be created in each thread running in parallel and needing to generate random sequences.
898All the three classes discussed here implement the automatic initialization method defined
899{\tt (AutoInit() ) } defined in the interface class and implement a number of other initialization
900methods {\tt (SetSeed() )}. SOPHYA provides also PPF handlers for these classes
901which can be used to save the complete state of the object and the underlying generator to PPF files.
902The generator object state can be restored subsequently from the PPF file. This feature can be used
903to generate very long sequence of random numbers, distributed over several runs.
904The example below illustrates the use of the generator classes and state save/restore
905through PPF files:
906\begin{enumerate}
907\item We create and use a {\bf ThSDR48RandGen} object in a first program. Its sate is saved
908to the PPF file rg.ppf at some point in the program.
909\begin{verbatim}
910// ------ program A
911// A.1- Create a thread safe generator based on drand48()
912ThSDR48RandGen rg;
913// A.2- Auto initilize its state (using the system time)
914rg.AutoInit();
915// A.3- compute and print a smal sequence of random numbers
916int N = 10;
917for(int i=0; i<N; i++)
918 cout << " I=" << i << " rand_flat01= " << rg.Flat01()
919 << " rand.gaussian= " << rg.Gaussian() << endl;
920// A.4- Save the generator state for subsequent use
921POutPersist po("rg.ppf");
922po << rg;
923// A.5- compute and print a second sequence of random numbers
924for(int i=0; i<N; i++)
925 cout << "++ I=" << i+N << " rand_flat01= " << rg.Flat01()
926 << " rand.gaussian= " << rg.Gaussian() << endl;
927\end{verbatim}
928
929\item The pseudo-random generator object state is restored in a second program from the previously created
930PPF file, before being used :
931 \begin{verbatim}
932// ------ program B
933// B.1- Create and initialize the generator from the previously saved state
934ThSDR48RandGen rgr;
935PInPersist pin("rg.ppf");
936pin >> rgr;
937int N = 10;
938// B.2- Compute and print a sequence of random number,
939// should be compared to the sequance A.5
940for(int i=0; i<N; i++)
941 cout << "-- I=" << i << " rand_flat01= " << rgr.Flat01()
942 << " rand.gaussian= " << rgr.Gaussian() << endl;
943\end{verbatim}
944\end{enumerate}
945
946%%%%%%%%%%%%
947\newpage
948\section{Module TArray}
949\index{\tcls{TArray}}
950{\bf TArray} module contains template classes for handling standard
951operations on numerical arrays. Using the class {\tt \tcls{TArray} },
952it is possible to create and manipulate up to 5-dimension numerical
953arrays {\tt (int, float, double, complex, \ldots)}. The include
954file {\tt array.h} declares all the classes and definitions
955in module TArray. {\bf Array} is a typedef for arrays
956with double precision floating value elements. \\
957{\tt typedef TArray$<$r\_8$>$ Array ; }
958
959\begin{figure}[hbt]
960\dclsccc{AnyDataObj}{BaseArray}{\tcls{TArray}}
961\dclsbb{PPersist}{\tcls{FIO\_TArray}}
962\end{figure}
963
964The development of this module started around 1999-2000,
965after evaluation of a number of publicly available
966C++ array hadling packages, including TNT, Lapack++, Blitz++,
967as well as commercial packages from RogueWave (math.h++ \ldots).
968Most of these packages provide interesting functionalities, however,
969not any one package seemed to fulfill most of our requirements.
970\begin{itemize}
971\item Capability to handle {\bf large - multidimensional - dense}
972arrays, for numerical data types. Although we have used templates, for
973data type specialisation, the actual code, apart inline functions is
974not in header files. Instead, we use explicit instanciation, and the
975compiled code for the various numerical types of arrays is the
976library .
977\item The shape and size of the arrays can be defined and changed
978at run time. The classes ensure the memory management of the
979created objets, with reference sharing for the array data.
980The default behaviour of the copy constructor is to share the data,
981avoiding expensive memory copies.
982\item The package provides transparent management of sub-arrays
983and slices, in an intuitive way, somehow similar to what is
984available in Mathlab or Scilab.
985\item The memory organisation for arrays, specially matrices
986(row-major or column major) can be
987controled. This provide compatibility when using existing C or
988Fortran coded numerical libraries.
989\item The classes provide efficient methods to perform basic arithmetic
990and mathematical operations on arrays. In addition, operator overload
991provides intuitive programming for element acces and most basic
992arithmetic operations.
993\item Conversion can be performed between arrays with different
994data types. Copy and arithmetic operations can be done transparently
995between arrays with different memory organisation patterns.
996\item This module does not provide more complex operations
997such as FFT or linear algebra. Additional libraries are used, with interface
998classes for these operations.
999\item ASCII formatted I/O, for printing and read/write operations to/from text files.
1000\item Efficient binary I/O for object persistence (PPF format), or import/export
1001to other data formats, such as FITS are provided by helper or handler classes.
1002\end{itemize}
1003
1004\subsection{Using arrays}
1005\index{Sequence} \index{RandomSequence} \index{RegularSequence}
1006\index{EnumeratedSequence}
1007The example below shows basic usage of arrays, creation, initialisation
1008and arithmetic operations. Different kind of {\bf Sequence} objects
1009can be used for initialising arrays.
1010
1011\begin{figure}[hbt]
1012\dclsbb{Sequence}{RandomSequence}
1013\dclsb{RegularSequence}
1014\dclsb{EnumeratedSequence}
1015\end{figure}
1016
1017The example below shows basic usage of arrays:
1018\index{\tcls{TArray}}
1019\begin{verbatim}
1020// Creating and initialising a 1-D array of integers
1021TArray<int> ia(5);
1022EnumeratedSequence es;
1023es = 24, 35, 46, 57, 68;
1024ia = es;
1025cout << "Array<int> ia = " << ia;
1026// 2-D array of floats
1027TArray<r_4> b(6,4), c(6,4);
1028// Initializing b with a constant
1029b = 2.71828;
1030// Filling c with random numbers
1031c = RandomSequence();
1032// Arithmetic operations
1033TArray<r_4> d = b+0.3f*c;
1034cout << "Array<float> d = " << d;
1035\end{verbatim}
1036
1037The copy constructor shares the array data, while the assignment operator
1038copies the array elements, as illustrated in the following example:
1039\begin{verbatim}
1040TArray<int> a1(4,3);
1041a1 = RegularSequence(0,2);
1042// Array a2 and a1 shares their data
1043TArray<int> a2(a1);
1044// a3 and a1 have the same size and identical elements
1045TArray<int> a3;
1046a3 = a1;
1047// Changing one of the a2 elements
1048a2(1,1,0) = 555;
1049// a1(1,1) is also changed to 555, but not a3(1,1)
1050cout << "Array<int> a1 = " << a1;
1051cout << "Array<int> a3 = " << a3;
1052\end{verbatim}
1053
1054\subsection{Arithmetic operations}
1055The four usual arithmetic operators ({\bf + \, - \, * \, / }) are defined
1056to perform constant addition, subtraction, multiplication and division.
1057The three operators ({\bf + \, - \, / }) between two arrays of the same type
1058are defined to perform element by element addition, subtraction
1059and division. In order to avoid confusion with matrix multiplication,
1060element by element multiplication is defined by overloading the
1061operator {\bf \, \&\& \, }, as shown in the example below:
1062\begin{verbatim}
1063TArray<int_4> a(4,3), b(4,3), c , d, e;
1064a = RegularSequence(1.,1.);
1065b = RegularSequence(10.,10.);
1066cout << a << b ;
1067c = a && b;
1068d = c / a;
1069e = (c / b) - a;
1070cout << c << d << e;
1071\end{verbatim}
1072
1073\subsection{Matrices and vectors}
1074\index{\tcls{TMatrix}} \index{\tcls{TVector}}
1075\begin{figure}[hbt]
1076\dclsccc{\tcls{TArray}}{\tcls{TMatrix}}{\tcls{TVector}}
1077\end{figure}
1078Vectors and matrices are 2 dimensional arrays. The array size
1079along one dimension is equal 1 for vectors. Column vectors
1080have {\tt NCols() = 1} and row vectors have {\tt NRows() = 1}.
1081Mathematical expressions involving matrices and vectors can easily
1082be translated into C++ code using {\tt TMatrix} and
1083{\tt TVector} objects. {\bf Matrix} and {\bf Vector} are
1084typedefs for double precision float matrices and vectors.
1085The operator {\bf *} beteween matrices is redefined to
1086perform matrix multiplication. One can then write: \\
1087\begin{verbatim}
1088 // We create a row vector
1089 Vector v(1000, BaseArray::RowVector);
1090 // Initialize values with a random sequence
1091 v = RandomSequence();
1092 // Compute the vector length (norm)
1093 double norm = (v*v.Transpose()).toScalar();
1094 cout << "Norm(v) = " << norm << endl;
1095\end{verbatim}
1096
1097This module contains basic array and matrix operations
1098such as the Gauss matrix inversion algorithm
1099which can be used to solve linear systems, as illustrated by the
1100example below:
1101\begin{verbatim}
1102#include "sopemtx.h"
1103// ...
1104// Creation of a random 5x5 matrix
1105Matrix A(5,5);
1106A = RandomSequence(RandomSequence::Flat);
1107Vector X0(5);
1108X0 = RandomSequence(RandomSequence::Gaussian);
1109// Computing B = A*X0
1110Vector B = A*X0;
1111// Solving the system A*X = B
1112Vector X;
1113LinSolve(A, B, X);
1114// Checking the result
1115Vector diff = X-X0;
1116cout << "X-X0= " << diff ;
1117double min,max;
1118diff.MinMax(min, max);
1119cout << " Min(X-X0) = " << min << " Max(X-X0) = " << max << endl;
1120\end{verbatim}
1121
1122{\bf Warning: } The operations defined in {\tt sopemtx.h}, such as
1123matrix inversion and linear system solver use a basic Gauss pivot
1124algorithm which are not adapted for large matrices ($>\sim 100x100$).
1125The services provided in other modules, such as {\bf LinAlg} should
1126be preferred in such cases.
1127
1128\subsection{Working with sub-arrays and Ranges}
1129\index{Range}
1130A powerful mechanism is included in array classes for working with
1131sub-arrays. The class {\bf Range} can be used to specify range of array
1132indexes in any of the array dimensions. Any regularly spaced index
1133range can be specified, using the {\tt start} and {\tt end} index
1134and an optional step (or stride). It is also possible to specify
1135the {\tt start} index and the number of elements.
1136\begin{itemize}
1137\item {\bf Range::all()} {\tt = Range(Range::firstIndex(), Range::lastIndex())} \\
1138return a Range objects representing all valid indexes along the
1139corresponding axe.
1140\item {\bf Range::first()} {\tt = Range(Range::firstIndex())} \\
1141return a Range object representing the first valid index
1142\item {\bf Range::last()} {\tt = Range(Range::lastIndex())}
1143return a Range object representing the last valid index
1144\item {\bf Range(idx)} represents a single index ({\bf = idx})
1145\item {\bf Range(first, last)} represents the range of indices
1146{\bf first} $\leq$ index $\leq$ {\bf last}.\\
1147The static method {\tt Range::lastIndex()} can be used
1148to specify the last valid index.
1149\item {\bf Range(first, last, step)} represents the range of index
1150which is equivalent to \\ {\tt for(index=first; index <= last; index += step) }
1151\item { \bf Range (first, last, size, step) } the general form can be used
1152to specify an index range, using the number of elements.
1153It is possible to specify a range of index, ending with the last valid index.
1154For example \\
1155\hspace*{5mm}
1156{\tt Range(Range::lastIndex(), Range::lastIndex(), 3, 2) } \\
1157defines the index range: \hspace*{5mm} last-4, last-2, last.
1158
1159\begin{center}
1160\begin{tabular}{ll}
1161\hline \\
1162\multicolumn{2}{c}{ {\bf Range} {\tt (start, end, size, step) } } \\[2mm]
1163\hline \\
1164{\bf Range} {\tt r(7); } & index range: \hspace{2mm} 7 \\
1165{\bf Range} {\tt r(3,6); } & index range: \hspace{2mm} 3,4,5,6 \\
1166{\bf Range} {\tt r(3,7,2); } & index range: \hspace{2mm} 3,5,7 \\
1167{\bf Range} {\tt r(7,0,3,1); } & index range: \hspace{2mm} 7,8,9 \\
1168{\bf Range} {\tt r(10,0,5,2); } & index range: \hspace{2mm} 10,12,14,16,18 \\[2mm]
1169\hline
1170\end{tabular}
1171\end{center}
1172\end{itemize}
1173
1174The method {\tt TArray<T>SubArray(Range ...)} can be used
1175to extract subarrays and slices. The operator {\tt operator() (Range rx, Range ry ...)}
1176is also overloaded for sub-array extraction.
1177For matrices, {\tt TMatrix<T>::Row()} and {\tt TMatrix<T>::Column()}
1178extract selected matrix rows and columns.
1179
1180The example illustrates the sub-array extraction using Range objects:
1181\begin{verbatim}
1182 // Creating and initialising a 2-D (6 x 4) array of integers
1183 TArray<int> iaa(6, 4);
1184 iaa = RegularSequence(1,2);
1185 cout << "Array<int> iaa = \n" << iaa;
1186 // We extract a sub-array - data is shared with iaa
1187 TArray<int> iae = iaa(Range(1, Range::lastIndex(), 3) ,
1188 Range::all(), Range::first() );
1189 cout << "Array<int> iae=subarray(iaa) = \n" << iae;
1190 // Changing iae elements changes corresponding iaa elements
1191 iae = 0;
1192 cout << "Array<int> iae=0 --> iaa = \n" << iaa;
1193
1194\end{verbatim}
1195
1196In the following example, a simple low-pass filter, on a one
1197dimensional stream (Vector) has been written using sub-arrays:
1198
1199\begin{verbatim}
1200// Input Vector containing a noisy periodic signal
1201 Vector in(1024), out(1024);
1202 in = RandomSequence(RandomSequence::Gaussian, 0., 1.);
1203 for(int kk=0; kk<in.Size(); kk++)
1204 in(kk) += 2*sin(kk*0.05);
1205// Compute the output vector by a simple low pass filter
1206 Vector out(1024);
1207 int w = 2;
1208 for(int k=w; k<in.Size()-w; k++)
1209 out(k) = in(Range(k-w, k+w).Sum()/(2.*w+1.);
1210\end{verbatim}
1211
1212
1213\subsection{Persistence, IO}
1214Arrays can easily be saved to, or restored from files in different formats.
1215SOPHYA library can handle array I/O to ASCII formatted files, to PPF streams,
1216as well as to files in FITS format.
1217FITS format input/output is provided through the classes in
1218{\bf FitsIOServer} module. Only arrays with data types
1219supported by the FITS standard can be handled during
1220I/O operations to and from FITS streams (See the FitsIOServer section
1221for additional details).
1222
1223\subsubsection{PPF streams}
1224
1225SOPHYA persistence (PPF streams) handles reference sharing, and multiply
1226referenced objects are only written once. A hierarchy of arrays and sub-arrays
1227written to a PPF stream is thus completely recovered, when the stream is read.
1228The following example illustrates this point:
1229\begin{verbatim}
1230{
1231// Saving an array with a sub-array into a POutPersist file
1232Matrix A(3,4);
1233A = RegularSequence(10,5);
1234// Create a sub-array of A
1235Matrix AS = A(Range(1,2), Range(2,3));
1236// Save the two arrays to a PPF stream
1237POutPersist pos("aas.ppf");
1238pos << A << AS;
1239}
1240{
1241// Reading arrays from the previously created PPF file aas.ppf
1242PInPersist pis("aas.ppf");
1243Matrix B,BS;
1244pis >> B >> BS;
1245// BS is a sub-array of B, modifying BS changes also B
1246BS(1,1) = 98765.;
1247cout << " B , BS after BS(1,1) = 98765. "
1248 << B << BS << endl;
1249}
1250\end{verbatim}
1251The execution of this sample code creates the file {\tt aas.ppf} and
1252its output is reproduced here. Notice that the array hierarchy is
1253recovered. BS is a sub-array of B, and modifying BS changes also
1254the corresponding element in B.
1255\begin{verbatim}
1256 B , BS after BS(1,1) = 98765.
1257
1258--- TMatrix<double>(NRows=3, NCols=4) ND=2 SizeX*Y*...= 4x3 ---
125910 15 20 25
126030 35 40 45
126150 55 60 98765
1262
1263--- TMatrix<double>(NRows=2, NCols=2) ND=2 SizeX*Y*...= 2x2 ---
126440 45
126560 98765
1266\end{verbatim}
1267
1268{\bf Warning: }
1269There is a drawback in this behaviour: only a single
1270copy of an array is written to a file, even if the array is modified,
1271without being resized and written (dumped) again to the same PPF stream.
1272However, this behavior can be changed using the {\tt RenewObjId()} method,
1273as illustrated below.
1274\begin{verbatim}
1275{
1276POutPersist pos("mca.ppf");
1277TArray<int_4> ia(5,3);
1278ia = 8;
1279pos << ia; // (1)
1280ia = 16;
1281pos << ia; // (2) Only a reference to the previously ia array is written
1282ia = 32;
1283ia.RenewObjId(); // We change the object Id
1284pos << ia; // (3) The complete array is dumped again
1285}
1286\end{verbatim}
1287
1288Only a single copy of the data is effectively written to the output
1289PPF file, corresponding to the value 8 for array elements, for the first two
1290write operations. When we
1291read the three array from the file mca.ppf, the same array elements
1292are obtained two times (all elements equal to 8), and a different array is obtained
1293the third time
1294\begin{verbatim}
1295{
1296PInPersist pis("mca.ppf");
1297TArray<int_4> ib;
1298pis >> ib;
1299cout << " First array read from mca.ppf : " << ib;
1300pis >> ib;
1301cout << " Second array read from mca.ppf : " << ib;
1302pis >> ib;
1303cout << " Third array read from mca.ppf : " << ib;
1304}
1305\end{verbatim}
1306
1307\subsubsection{ASCII streams}
1308
1309The {\bf WriteASCII} method can be used to dump an array to an ASCII
1310formatted file, while the {\bf ReadASCII} method can be used to decode
1311ASCII formatted files. Space or tabs are the possible separators.
1312Complex numbers should be specified as a pair of comma separated
1313real and imaginary parts, enclosed in parenthesis.
1314
1315\begin{verbatim}
1316{
1317// Creating array A and writing it to an ASCII file (aaa.txt)
1318Array A(4,6);
1319A = RegularSequence(0.5, 0.2);
1320ofstream ofs("aaa.txt");
1321A.WriteASCII(ofs);
1322}
1323{
1324// Decoding the ASCII file aaa.txt
1325ifstream ifs("aaa.txt");
1326Array B;
1327sa_size_t nr, nc;
1328B.ReadASCII(ifs,nr,nc);
1329cout << " Array B; B.ReadASCII() from file " << endl;
1330cout << B ;
1331}
1332\end{verbatim}
1333
1334\subsection{Cast without conversion}
1335Data conversion between arrays with different data type is handled transparently,
1336through the copy constructor or the assignment (=) operator . However, in rare cases,
1337one wants to access the same memory locations, without data type conversion.
1338The template functions defined in {\tt arrctcast.h} can be used to access the same
1339memory locations, by arrays with different data types. The SOPHYA/NDataBlock
1340reference sharing mechanism is effective When using these functions.
1341Notice that the array size or stride may change during these cast operations. \\
1342 {\tt arrctcast.h} has been introduced in version V=2.1 (Nov 2007), and has not been
1343 fully tested for non packed arrays.
1344\begin{verbatim}
1345// We define and initialize a real array :
1346TArray<r_4> fa(5);
1347fa = RegularSequence(1.25,0.5);
1348cout << " fa= " << fa;
1349// We construct an integer array from fa, where the floating point values
1350// are converted to integer values
1351TArray<uint_2> ia(fa);
1352cout << " ia= " << ia;
1353cout << " ia (in hex)= " << hex << ia << dec;
1354// We can also access the fa memory locations interpreted as short integers
1355uint_2 ui2;
1356// Note that sfia size is double the fa size
1357TArray<uint_2> sfia = ArrayCast(fa, ui2);
1358cout << " sfia= " << sfia;
1359cout << " sfia (in hex)= " << hex << sfia << dec;
1360\end{verbatim}
1361One of the most useful case of these array type cast without conversion
1362correspond to accessing the real or imaginary part of a complex array.
1363Two specific template functions {\tt SDRealPart()} and {\tt SDImagPart()}
1364are also defined in {\tt arrctcast.h}. Two other functions {\tt ArrCastR2C()}
1365and {\tt ArrCastC2R() } are also defined for real to complex, and
1366complex to real cast.
1367Their usage is shown in the next paragraph on complex arrays.
1368
1369\subsection{Complex arrays}
1370The {\bf TArray} module provides few functions for manipulating
1371arrays of complex numbers (single and double precision).
1372These functions are declared in {\tt matharr.h}.
1373\begin{itemize}
1374\item[\bul] Creating a complex array through the specification of the
1375real and imaginary parts.
1376\item[\bul] Functions returning arrays corresponding to real and imaginary
1377parts of a complex array: {\tt real(za) , imag(za) }
1378({\bf Warning:} Note that the these functions create an array and copies
1379the data from the original complex array. They do not provide
1380shared memory access to real and imaginary parts.
1381For shared memory access, use functions {\tt SDRealPart()} and
1382{\tt SDImagPart() } (see below).
1383\item[\bul] Functions returning arrays corresponding to the module,
1384phase, and module squared of a complex array:
1385 {\tt phase(za) , module(za) , module2(za) }
1386\end{itemize}
1387
1388\begin{verbatim}
1389 TVector<r_4> p_real(10, BaseArray::RowVector);
1390 TVector<r_4> p_imag(10, BaseArray::RowVector);
1391 p_real = RegularSequence(0., 0.5);
1392 p_imag = RegularSequence(0., 0.25);
1393 TVector< complex<r_4> > zvec = ComplexArray(p_real, p_imag);
1394 cout << " :: zvec= " << zvec;
1395 cout << " :: real(zvec) = " << real(zvec) ;
1396 cout << " :::: imag(zvec) = " << imag(zvec) ;
1397 cout << " :::: module2(zvec) = " << module2(zvec) ;
1398 cout << " :::: module(zvec) = " << module(zvec) ;
1399 cout << " :::: phase(zvec) = " << phase(zvec) ;
1400\end{verbatim}
1401
1402The decoding of complex numbers from an ASCII formatted stream
1403is illustrated by the next example. As mentionned already,
1404complex numbers should be specified as a pair of comma separated
1405real and imaginary parts, enclosed in parenthesis.
1406
1407\begin{verbatim}
1408csh> cat zzz.txt
1409(1.,-1) (2., 2.5) -3. 12.
1410-24. (-6.,7.) 14.2 (8.,64.)
1411
1412// Decoding of complex numbers from an ASCII file
1413// Notice that the << operator can be used instead of ReadASCII
1414TArray< complex<r_4> > Z;
1415ifstream ifs("zzz.txt");
1416ifs >> Z;
1417cout << " TArray< complex<r_4> > Z from file zzz.txt " << Z ;
1418\end{verbatim}
1419
1420\noindent {\bf Shared data access :} It is possible to access a complex array
1421elements (real and imaginary parts) through the template functions defined
1422in {\tt arrctcast.h} and discussed above. Four specific template functions are
1423defined for shared data access to complex arrays:
1424\begin{itemize}
1425\item {\bf SDRealPart } return a float/double TArray, when applied to an
1426array with complex elements, corresponding the the {\bf real} part of the complex values.
1427\item {\bf SDImagPart } return a float/double TArray, when applied to an
1428array with complex elements, corresponding the the {\bf imaginary} part of the complex values.
1429\item {\bf ArrCastC2R } return a float/double TArray, when applied to an
1430array with complex elements. The data is shared between the two arrays, the real array
1431containing real and imaginary parts as successive elements.
1432\item {\bf ArrCastR2C } return a complex valued TArray, when applied to a float/double array.
1433\end{itemize}
1434The example below shows how to use these functions:
1435
1436\begin{verbatim}
1437// We define a complex array
1438TArray< complex<r_4> > za(5);
1439cout << " za= " << za;
1440// And two real arrays, corresponding to the real and imaginary parts
1441TArray<r_4> rza = SDRealPart(za);
1442TArray<r_4> iza = SDImagPart(za);
1443// We initialize the real and imaginary parts of the complex array
1444rza = RegularSequence(10.,2.);
1445iza = RegularSequence(5.,0.75);
1446cout << " rza=..., iza=... ----> za = " << za;
1447// The complex array seen as a real array (double size)
1448TArray<r_4> aza = ArrCastC2R(za);
1449cout << " za --> aza= " << aza;
1450TArray< complex<r_4> > zb;
1451zb = ArrCastR2C(aza);
1452KeepObj(zb);
1453\end{verbatim}
1454
1455\subsection{Memory organisation}
1456{\tt \tcls{TArray} } can handle numerical arrays with various memory
1457organisation, as long as the spacing (steps) along each axis is
1458regular. The five axis are labeled X,Y,Z,T,U. The examples below
1459illustrates the memory location for a 2-dimensional, $N_x=4 \times N_y=3$.
1460The first index is along the X axis and the second index along the Y axis.
1461\begin{verbatim}
1462 | (0,0) (1,0) (2,0) (3,0) |
1463 | (0,1) (1,1) (2,1) (3,1) |
1464 | (0,2) (1,2) (2,2) (3,2) |
1465\end{verbatim}
1466In the first case, the array is completely packed
1467($Step_X=1, Step_Y=N_X=4$), with zero offset,
1468while in the second case, $Step_X=2, Step_Y=10, Offset=10$:
1469\begin{verbatim}
1470 | 0 1 2 3 | | 10 12 14 16 |
1471Ex1 | 4 5 6 7 | Ex2 | 20 22 24 26 |
1472 | 8 9 10 11 | | 30 32 34 36 |
1473\end{verbatim}
1474
1475For matrices and vectors, an optional argument ({\tt MemoryMapping})
1476can be used to select the memory mapping, where two basic schemes
1477are available: \\
1478{\tt CMemoryMapping} and {\tt FortranMemoryMapping}. \\
1479In the case where {\tt CMemoryMapping} is used, a given matrix line
1480is packed in memory, while the columns are packed when
1481{\tt FortranMemoryMapping} is used. The first index when addressing
1482the matrix elements (line number index) runs along
1483the Y-axis if {\tt CMemoryMapping} is used, and along the X-axis
1484in the case of {\tt FortranMemoryMapping}.
1485Arithmetic operations between matrices
1486with different memory organisation is allowed as long as
1487the two matrices have the same sizes (Number of rows and columns).
1488The following code example and the corresponding output illustrates
1489these two memory mappings. The {\tt \tcls{TMatrix}::TransposeSelf() }
1490method changes effectively the matrix memory mapping, which is also
1491the case of {\tt \tcls{TMatrix}::Transpose() } method without argument.
1492
1493\begin{verbatim}
1494BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
1495TArray<r_4> X(4,2);
1496X = RegularSequence(1,1);
1497cout << "Array X= " << X << endl;
1498TMatrix<r_4> X_C(X, true);
1499cout << "Matrix X_C (CMemoryMapping) = " << X_C << endl;
1500TMatrix<r_4> X_F(X_C, true);
1501X_F.TransposeSelf(); // we make it a FortranMemoryMapping matrix
1502cout << "Matrix X_F (FortranMemoryMapping) = " << X_F << endl;
1503\end{verbatim}
1504This code would produce the following output (X\_F = Transpose(X\_C)) :
1505\begin{verbatim}
1506Array X=
1507--- TArray<f> ND=2 SizeX*Y*...= 4x2 ---
15081, 2, 3, 4
15095, 6, 7, 8
1510
1511Matrix X_C (CMemoryMapping) =
1512--- TMatrix<f>(NRows=2, NCols=4) ND=2 SizeX*Y*...= 4x2 ---
15131, 2, 3, 4
15145, 6, 7, 8
1515
1516Matrix X_F (FortranMemoryMapping) =
1517--- TMatrix<f>(NRows=4, NCols=2) ND=2 SizeX*Y*...= 4x2 ---
15181, 5
15192, 6
15203, 7
15214, 8
1522\end{verbatim}
1523
1524{\bf Note:} Only the {\tt TMatrix} class is sensitive to the memory organisation. It should also
1525be noted that the first index for a matrix is the row index. For a {\tt TArray} object, the first
1526index correspond always to the X axis. For a matrix in the {\tt CMemoryMapping}, the row
1527index is the Y index of the corresponding array, while for a matrix
1528with {\tt FortranMemoryMapping}, the row index is the X index.
1529\begin{verbatim}
1530# 2D array arr , matrix mx
1531TArray<r_4> arr;
1532...
1533TMatrix<r_4> mx(arr);
1534...
1535### CMemoryMapping :
1536mx( row , col ) -> arr( ix=col, jy= row )
1537### FortranMemoryMapping :
1538mx( row , col ) -> arr( ix=row, jy= col )
1539\end{verbatim}
1540The following code fragment shows the difference between element access
1541index order for the different memory organisation:
1542\begin{verbatim}
1543 TArray<int_4> a(5,3);
1544 a = RegularSequence(10,2);
1545 cout << "Array a = " << a << endl;
1546 TMatrix<r_4> mac(a);
1547 cout << "Matrix mac (CMemoryMapping) = " << mac << endl;
1548 TMatrix<r_4> maf(mac);
1549 maf.TransposeSelf(); // we make it a FortranMemoryMapping matrix
1550 cout << "Matrix maf (FortranMemoryMapping) = " << maf << endl;
1551 // ---- element access
1552 cout << " Array(ix,jy) : a(2,0)= " << a(2,0) << " a(1,2)= " << a(1,2)
1553 << " a(0,2)= " << a(0,2) << endl;
1554 cout << " mac(row,col) : mac(2,0)= " << mac(2,0) << " mac(1,2)= " << mac(1,2)
1555 << " mac(0,2)= " << mac(0,2) << endl;
1556 cout << " maf(row,col) : maf(2,0)= " << maf(2,0) << " maf(1,2)= " << maf(1,2)
1557 << " maf(0,2)= " << maf(0,2) << endl;
1558\end{verbatim}
1559
1560
1561\subsection{Special Square Matrices (SSQM) }
1562SOPHYA V2.2 introduces new classes for handling special square matrices (diagonal, symmetric, triangular \ldots).
1563The figure below shows the corresponding class hierarchy:
1564\begin{figure}[hbt]
1565\dclsccc{AnyDataObj}{\tcls{SpecialSquareMatrix}}{ \tcls{ DiagonalMatrix } }
1566\dclsc{ \tcls{ LowerTriangularMatrix } }
1567\dclsc{ \tcls{ SymmetricMatrix } }
1568\end{figure}
1569
1570\index{ \tcls{DiagonalMatrix} }
1571\index{ \tcls{SymmetricMatrix} }
1572\index{ \tcls{TriangularMatrix} }
1573
1574Theses classes provides methods similar to TArray/TMatrix objects for manipulating these special kind of matrices.
1575However, no sub-matrix extraction method is provided. The following features are currently implemented:
1576\begin{itemize}
1577\item These classes implements reference sharing and behave in a way similar to TArray objects. The copy
1578constructor shares the data, while the equal operator (=) performs an element by element copy.
1579\item The {\bf FIO\_SpecialSquareMatrix$<$T$>$ } manages the persistence (I/O) in PPF format for these classes.
1580\item Constructor from and conversion to general {\bf TMatrix$<$T$>$ } objects.
1581 Non relevant elements are ignored when constructing from a general matrix and the special square matrix can be
1582converted to the general matrix through the {\tt ConvertToStdMatrix () }.
1583\item Definition of element access operator {\tt mx(r,c) } as well as global assignment operator {\tt mx = a; mxb=mxa; }.
1584operator {\tt mx[k] } can be used to access the non zero element k.
1585\item Addition, subtraction of a constant, multiplication or division by a constant, as well as
1586\item Element by element addition, subtraction, multiplication and division methods between same type of SSQM
1587matrices. As in the cases of general matrices, the four binary operators ( + , - , \&\& , /) are overloaded to perform
1588these operations.
1589\item Matrix multiplication operation between same type SSQM objects or between a general matrix and an SSQM
1590object is defined through the methods: \\
1591{\\ Multiply() , MultiplyXG(), MultiplyGX(), X=D,S,T } \\
1592The binary multiplication operator ( * ) is then overloaded to perform matrix multiplication.
1593\end{itemize}
1594The sample code below illustrates the use of these special matrices:
1595\begin{verbatim}
1596// Create and initialize a diagonal matrix
1597DiagonalMatrix<double> diag(3);
1598diag(0,0)=1.; diag(1,1)=2.; diag(2,2)=3.;
1599// use of multiplication by a constant operator
1600diag *= 10.;
1601// check the diagonal matrix
1602cout << diag << diag.ConvertToStdMatrix() << endl;
1603// define and initialize a general matrix
1604TMatrix<double> mx(3,3);
1605mx=RegularSequence();
1606cout << mx;
1607// check the result of a matrix mutiplication
1608cout << mx*diag;
1609\end{verbatim}
1610
1611\newpage
1612
1613\section{Module HiStats}
1614\begin{figure}[hbt]
1615\dclsbb{AnyDataObj}{Histo}
1616\dclscc{Histo}{HProf}
1617\dclsbb{AnyDataObj}{Histo2D}
1618\dclsbb{AnyDataObj}{HistoErr}
1619\dclsbb{AnyDataObj}{Histo2DErr}
1620\caption{partial class diagram for histograms and ntuples}
1621\end{figure}
1622
1623{\bf HiStats} contains classes for creating, filling, printing and
1624doing various operations on one or two dimensional histograms
1625{\tt Histo} and {\tt Histo2D} as well as profile histograms {\tt HProf}. \\
1626This module also contains {\tt NTuple} and {\tt DataTable} which are
1627more or less the same as the binary or ascii FITS tables.
1628
1629\subsection{Histograms}
1630\subsubsection{1D Histograms}
1631\index{Histo}
1632For 1D histograms, various numerical methods are provided such as
1633computing means and sigmas, finding maxima, fitting, rebinning,
1634integrating \dots \\
1635The example below shows creating and filling a one dimensional histogram
1636of 100 bins from $-5.$ to $+5.$ to create a Gaussian normal distribution
1637with errors~:
1638\begin{verbatim}
1639#include "histos.h"
1640// ...
1641Histo H(-0.5,0.5,100);
1642H.Errors();
1643for(int i=0;i<25000;i++) {
1644 double x = NorRand();
1645 H.Add(x);
1646}
1647H.Print(80);
1648\end{verbatim}
1649
1650\subsubsection{2D Histograms}
1651\index{Histo2D}
1652Much of these operations are also valid for 2D histograms. 1D projection
1653or slices can be set~:
1654\begin{verbatim}
1655#include "histos2.h"
1656// ...
1657Histo2D H2(-1.,1.,100,0.,60.,50);
1658H2.SetProjX(); // create the 1D histo for X projection
1659H2.SetBandX(25.,35.); // create 1D histo projection for 25.<y<35.
1660H2.SetBandX(35.,45.); // create 1D histo projection for 35.<y<45.
1661H2.SetBandX(40.,55.); // create 1D histo projection for 40.<y<55.
1662//... fill H2 with what ever you want
1663H2.Print();
1664Histo *hx = H2.HProjX();
1665 hx->Print(80);
1666Histo *hbx2 = HBandX(1); // Get the second X band (35.<y<45.)
1667 hbx2->Print(80);
1668\end{verbatim}
1669
1670\subsubsection{Profile Histograms}
1671\index{HProf}
1672Profiles histograms {\bf HProf} contains the mean and the
1673sigma of the distribution
1674of the values filled in each bin. The sigma can be changed to
1675the error on the mean. When filled, the profile histogram looks
1676like a 1D histogram and much of the operations that can be done on 1D histo
1677may be applied onto profile histograms.
1678
1679\subsubsection{Histograms HistoErr and Histo2DErr}
1680\index{HistoErr}
1681The {\bf HistoErr} are basic histograms where the number of entries for each bin is kept.
1682Methods to compute of the mean and the variance in each bin are provided.
1683The {\bf Histo2DErr} is the same for $2$ dimensions.
1684
1685\subsection{Data tables (tuples)}
1686\begin{figure}[hbt]
1687\dclsbb{AnyDataObj}{NTuple}
1688\dclsccc{AnyDataObj}{BaseDataTable}{DataTable}
1689\dclscc{BaseDataTable}{SwPPFDataTable}
1690\end{figure}
1691
1692\subsubsection{NTuple}
1693\index{NTuple}
1694{\bf NTuple} are memory resident tables of 32 or 64 bits floating values
1695(float/double).They are arranged in columns. Each line is often called an event.
1696These objects are frequently used to analyze data.
1697The piapp graphicals tools can plot a column against an other one
1698with respect to various selection cuts. \\
1699Here is an example of creation and filling~:
1700\begin{verbatim}
1701#include "ntuple.h"
1702#include "srandgen.h"
1703// ...
1704char* nament[4] = {"i","x","y","ey"};
1705r_4 xnt[4];
1706NTuple NT(4,nament);
1707for(i=0;i<5000;i++) {
1708 xnt[0] = i+1;
1709 xnt[1] = 5.*drandpm1(); // a random value between -5 and +5
1710 xnt[2] = 100.*exp(-0.5*xnt[1]*xnt[1]) + 1.;
1711 xnt[3] = sqrt(xnt[2]);
1712 xnt[2] += xnt[3] * NorRand(); // add a random gaussian error
1713 NT.Fill(xnt);
1714}
1715\end{verbatim}
1716
1717{\bf XNTuple} provide additional functionalities, compared to NTuple. However,
1718this class is deprecated and superseded by classes inheriting from BaseDataTable.
1719It is only kept for backward compatibility and should not be used anymore.
1720Use DataTable and SwPPFDataTable instead.
1721Object of type XNTuple handle various types
1722of column values (double,float,int,string,...) and can handle
1723very large data sets, through swap space on disk.
1724
1725\subsubsection{DataTables}
1726\label{datatables}
1727\index{DataTable}
1728The class {\bf DataTable} extends significantly the functionalities provided by
1729NTuple. DataTable is a memory resident implementation of the interface
1730{\bf BaseDataTable } which organizes the data as a 2-D table. User can define
1731the name and data type of each column. Data is added to the table as rows.
1732The table is extended as necessary when adding rows.
1733The sample code below shows an example of DataTable usage :
1734\begin{verbatim}
1735 #include "datatable.h"
1736 // ...
1737 {
1738 DataTable dt(64);
1739 dt.AddFloatColumn("X0_f");
1740 dt.AddFloatColumn("X1_f");
1741 dt.AddDoubleColumn("X0X0pX1X1_d");
1742 double x[5];
1743 for(int i=0; i<63; i++) {
1744 x[0] = (i/9)-4.; x[1] = (i/9)-3.; x[2] = x[0]*x[0]+x[1]*x[1];
1745 dt.AddLine(x);
1746 }
1747 // Printing table info
1748 cout << dt ;
1749 // Saving object into a PPF file
1750 POutPersist po("dtable.ppf");
1751 po << dt ;
1752 }
1753\end{verbatim}
1754
1755
1756\index{SwPPFDataTable}
1757The class {\bf SwPPFDataTable} implements the BaseDataTable interface
1758using segmented data blocks with swap on PPF streams. Very large data sets
1759can be created and manipulated through this class. A similar class
1760SwFitsDataTable (\ref{SwFitsDataTable}), using
1761FITS files as swap space is also provided in the FitsIOServer module.
1762
1763\index{DataTableRow}
1764The class {\bf DataTableRow } is an auxiliary class which simplifies the manipulation
1765of BaseDataTable object rows.
1766The example below show how to create and filling a table, using a PPF stream as
1767swap space. In addition, we have used a {\tt DataTableRow} to prepare data
1768for each table line.
1769\begin{verbatim}
1770 #include "swppfdtable.h"
1771 // ...
1772 {
1773 // ---------- Create an output PPF stream (file)
1774 POutPersist po("swdtable.ppf");
1775 // ------------------
1776 // Create a table with 3 columns, using the above stream as swap space
1777 SwPPFDataTable dtrow(po, 64);
1778 dtrow.AddStringColumn("sline");
1779 dtrow.AddIntegerColumn("line");
1780 dtrow.AddDateTimeColumn("datime");
1781 //
1782 TimeStamp ts, ts2; // Initialize current date and time
1783 string sline;
1784 //---- Create a table row with the required structure
1785 DataTableRow row = dtrow.EmptyRow();
1786 // ----- Fill the table
1787 for(int k = 0; k<2500; k++) {
1788 sline = "L-";
1789 sline += (string)MuTyV(k);
1790 row["sline"] = sline;
1791 row[1] = k;
1792 ts2.Set(ts.ToDays()+(double)k);
1793 row["datime"] = ts2;
1794 dtrow.AddRow(row);
1795 }
1796 //------ Write the table itself to the stream, before closing the file
1797 po << PPFNameTag("SwTable") << dtrow;
1798 }
1799\end{verbatim}
1800%%
1801The previously created table can easily be read in, as shown below:
1802%%
1803\begin{verbatim}
1804 #include "swppfdtable.h"
1805 // ...
1806 {
1807 // ------ Create the input PPF stream (file)
1808 PInPersist pin("swdtable.ppf");
1809 // ------ Read in the SwPPFDataTable object
1810 SwPPFDataTable dtr;
1811 pin >> PPFNameTag("SwTable") >> dtr;
1812 // ---- Create a table row with the required structure
1813 DataTableRow row = dtr.EmptyRow();
1814 // ---- Acces and print two of the table rows :
1815 cout << dtr.GetRow(6, row) << endl;
1816 cout << dtr.GetRow(17, row) << endl;
1817 }
1818\end{verbatim}
1819
1820\subsection{Writing, viewing \dots }
1821%%
1822Histogram and NTuple/DataTable objects have PPF handlers and
1823can be written to or read from PPF files. Sophya graphical tool (spiapp) can
1824automatically display and operate on all these objects.
1825The following example shows how to write the previously created objects
1826into such a file~:
1827\begin{verbatim}
1828{
1829char *fileout = "myfile.ppf";
1830string tag;
1831POutPersist outppf(fileout);
1832tag = "H"; outppf.PutObject(H,tag);
1833tag = "H2"; outppf.PutObject(H2,tag);
1834tag = "NT"; outppf.PutObject(NT,tag);
1835} // closing ``}'' destroy ``outppf'' and automatically close the file !
1836\end{verbatim}
1837
1838\newpage
1839\section{Module NTools}
1840
1841This module provides elementary numerical tools for numerical integration,
1842fitting, sorting and ODE solving. FFTs are also provided (Mayer,FFTPack).
1843
1844\subsection{Fitting}
1845\index{Fitting} \index{Minimisation}
1846Fitting is done with two classes {\tt GeneralFit} and {\tt GeneralFitData}
1847and is based on the Levenberg-Marquardt method.
1848\index{GeneralFit} \index{GeneralFitData}
1849GeneralFitData is a class which provide a description of the data
1850to be fitted. GeneralFit is the fitter class. Parametrized functions
1851can be given as classes which inherit {\tt GeneralFunction}
1852or as simple C functions. Classes of pre-defined functions are provided
1853(see files fct1dfit.h and fct2dfit.h). The user interface is very close
1854from that of the CERN {\tt Minuit} fitter.
1855Number of objects (Histo, HProf \dots ) are interfaced with GeneralFit
1856and can be easily fitted. \\
1857Here is a very simple example for fitting the previously created NTuple
1858with a Gaussian~:
1859\begin{verbatim}
1860#include "fct1dfit.h"
1861// ...
1862
1863// Read from ppf file
1864NTuple nt;
1865{
1866PInPersist pis("myfile.ppf");
1867string tag = "NT"; pis.GetObject(nt,tag);
1868}
1869
1870// Fill GeneralData
1871GeneralData mGdata(nt.NEntry());
1872for(int i=0; i<nt.NEntry(); i++)
1873 mGdata.AddData1(xnt[1],xnt[2],xnt[3]); // Fill x, y and error on y
1874mGData.PrintStatus();
1875
1876// Function for fitting : y = f(x) + noise
1877Gauss1DPol mFunction; // gaussian + constant
1878
1879// Prepare for fit
1880GeneralFit mFit(&mFunction); // create a fitter for the choosen function
1881mFit.SetData(&mGData); // connect data to the fitter
1882
1883// Set and initialise the parameters (that's non-linear fitting!)
1884// (num par, name, guess start, step, [limits min and max])
1885mFit.SetParam(0,"high",90.,1..);
1886mFit.SetParam(1,"xcenter",0.05,0.01);
1887mFit.SetParam(2,"sigma",sig,0.05,0.01,10.);
1888 // Give limits to avoid division by zero
1889mFit.SetParam(3,"constant",0.,1.);
1890
1891// Fit and print result
1892int rcfit = mFit.Fit();
1893mFit.PrintFit();
1894if(rcfit>0) {)
1895 cout<<"Reduce_Chisquare = "<<mFit.GetChi2Red()
1896 <<" nstep="<<mFit.GetNStep()<<" rc="<<rcfit<<endl;
1897} else {
1898 cout<<"Fit_Error, rc = "<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
1899 mFit.PrintFitErr(rcfit);
1900}
1901
1902// Get the result for further use
1903TVector<r_8> ParResult = mFit.GetParm();
1904cout<<ParResult;
1905\end{verbatim}
1906
1907Much more usefull possibilities and detailed informations might be found
1908in the HTML pages of the Sophya manual.
1909
1910\subsection{Simplex method}
1911The class {\bf MinZSimplex} implements the simplex method for non linear
1912optimization / minimization. See the SOPHYA manual for more information.
1913
1914\subsection{Polynomial}
1915\index{Polynomial} \index{Poly} \index{Poly2}
1916Polynomials of 1 or 2 variables are supported ({\tt Poly} and {\tt Poly2}).
1917Various operations are supported~:
1918\begin{itemize}
1919\item elementary operations between polynomials $(+,-,*,/) $
1920\item setting or getting coefficients
1921\item computing the value of the polynomial for a given value
1922 of the variable(s),
1923\item derivating
1924\item computing roots (degre 1 or 2)
1925\item fitting the polynomial to vectors of data.
1926\end{itemize}
1927Here is an example of polynomial fitting~:
1928\begin{verbatim}
1929#include "poly.h"
1930// ...
1931Poly pol(2);
1932pol[0] = 100.; pol[1] = 0.; pol[2] = 0.01; // Setting coefficients
1933TVector<r_8> x(100);
1934TVector<r_8> y(100);
1935TVector<r_8> ey(100);
1936for(int i=0;i<100;i++) {
1937 x(i) = i;
1938 ey(i) = 10.;
1939 y(i) = pol((double) i) + ey(i)*NorRand();
1940 ey(i) *= ey(i)
1941}
1942
1943TVector<r_8> errcoef;
1944Poly polfit;
1945polfit.Fit(x,y,ey,2,errcoef);
1946
1947cout<<"Fit Result"<<polfit<<endl;
1948cout<<"Errors :"<<errcoef;
1949\end{verbatim}
1950
1951Similar operations can be done on polynomials with 2 variables.
1952
1953\subsection{Integration, Differential equations}
1954\index{Integration}
1955The NTools module provide also simple classes for numerical integration
1956of functions and differential equations.
1957\begin{figure}[hbt]
1958\dclsbb{Integrator}{GLInteg}
1959\dclsb{TrpzInteg}
1960\end{figure}
1961
1962\index{GLInteg} \index{TrpzInteg}
1963{\bf GLInteg} implements the integration through Gauss-Legendre method
1964and {\bf TrpzInteg} implements trapeze integration. For {\bf TrpzInteg},
1965number of steps specify the number of trapeze, and integration step,
1966their width.
1967The sample code below illustrates the use of TrpzInteg class:
1968\begin{verbatim}
1969#include "integ.h"
1970// ......................................................
1971// Function to be integrated
1972double myf(double x)
1973{
1974// Simple a x + b x^2 (a=2 b=3)
1975return (x*(2.+3.*x));
1976}
1977// ......................................................
1978
1979// Compute Integral(myf, 2., 5.) between xmin=2., xmax=5.
1980TrpzInteg trpz(myf, 2., 5.);
1981// We specify an integration step
1982trpz.DX(0.01);
1983// The integral can be computed as trpz.Value()
1984double myf_integral = trpz.Value();
1985// We could have used the cast operator :
1986cout << "Integral[myf, 2., 5.]= " << (double)trpz << endl;
1987// Limits can be specified through ValueBetween() method
1988cout << "Integral[myf, 0., 4.]= " << trpz.ValueBetween(0.,4.) << endl;
1989\end{verbatim}
1990
1991\subsection{Fourier transform (FFT)}
1992\index{FFT} \index{FFTPackServer}
1993An abstract interface for performing FFT operations is defined by the
1994{\bf FFTServerInterface} class. The {\bf FFTPackSever} class implements
1995one dimensional FFT, on real and complex data. FFTPackServer uses an
1996adapted and extended version of FFTPack (available from netlib),
1997translated in C, and can operate on single and double precision
1998({\tt float, double}) data.
1999
2000The sample code below illustrates the use of FFTServers:
2001\begin{verbatim}
2002#include "fftpserver.h"
2003 // ...
2004TVector<r_8> in(32);
2005TVector< complex<r_8> > out;
2006in = RandomSequence();
2007FFTPackServer ffts;
2008ffts.setNormalize(true); // To have normalized transforms
2009cout << " FFTServer info string= " << ffts.getInfo() << endl;
2010cout << "in= " << in << endl;
2011cout << " Calling ffts.FFTForward(in, out) : " << endl;
2012ffts.FFTForward(in, out);
2013cout << "out= " << out << endl;
2014\end{verbatim}
2015
2016% \newpage
2017\section{Module SUtils}
2018Some utility classes and C/C++ string manipulation functions are gathered
2019in {\bf SUtils} module.
2020\subsection{Using DataCards}
2021\index{DataCards}
2022The {\bf DataCards} class can be used to read parameters from a file.
2023Each line in the file starting with \@ defines a set of values
2024associated with a keyword. In the example below, we read the
2025parameters corresponding with the keyword {\tt SIZE} from the
2026file {\tt ex.d}. We suppose that {\tt ex.d} contains the line: \\
2027{\tt @SIZE 400 250} \\
2028\begin{verbatim}
2029#include "datacards.h"
2030// ...
2031// Initialising DataCards object dc from file ex.d
2032DataCards dc( "ex.d" );
2033// Getting the first and second parameters for keyword size
2034// We define a default value 100
2035int size_x = dc.IParam("SIZE", 0, 100);
2036int size_y = dc.IParam("SIZE", 1, 100);
2037cout << " size_x= " << size_x << " size_y= " << size_y << endl;
2038\end{verbatim}
2039
2040\section{Module SysTools}
2041The {\bf SysTools} module contains classes implementing interface to some
2042OS specific services, such as thread creation and management, dynamic loading and
2043resource usage information. For example, yhe class {\bf Periodic} provides the
2044necessary services needed to implement the execution of a periodic action.
2045
2046\subsection{Resource usage (CPU, memory \ldots) }
2047 The class {\bf ResourceUsage} \index{ResourceUsage}
2048and {\bf Timer} \index{Timer} provides access to information
2049about various resource usage (memory, CPU, ...).
2050The class {\bf Timer} \index{time (CPU, elapsed)} and c-functions
2051{\tt InitTim() , PrtTim(const char * Comm) } can be used to print
2052the amount of CPU and elapsed time in programs.
2053
2054The following sample code illustrates the use of {\bf ResourceUsage} :
2055\begin{verbatim}
2056 // How to check resource usage for a given part of the program
2057 ResourceUsage res;
2058 // --- Part of the program to be checked : Start
2059 // ...
2060 res.Update();
2061 cout << " Memory size increase (KB):" << res.getDeltaMemorySize() << endl;
2062 cout << " Resource usage info : \n" << res << endl;
2063\end{verbatim}
2064
2065\subsection{Thread management classes}
2066\index{ZThread} \index{ZMutex}
2067A basic interface to POSIX threads is also provided
2068through the \index{threads} {\bf ZThread}, {\bf ZMutex} and {\bf ZSync}
2069classes. The best way to use thread management classes is by inheriting
2070from {\bf ZThread} and redefining the {\tt run() } method.
2071It is also possible to use the default {\tt run() } implementation and associate
2072a function to perform the action, as in the example below :
2073\begin{verbatim}
2074 // The functions to perform computing
2075 void fun1(void *arg) { }
2076 void fun2(void *arg) { }
2077 // ...
2078 ZThread zt1;
2079 zt1.setAction(fun1, arg[1]);
2080 ZThread zt2;
2081 zt2.setAction(fun2, arg[1]);
2082 cout << " Starting threads ... " << endl;
2083 zt1.start();
2084 zt2.start();
2085 cout << " Waiting for threads to end ... " << endl;
2086 zt1.join();
2087 zt2.join();
2088\end{verbatim}
2089The classes {\bf ZMutex} \index{mutex} and {\bf ZSync} can be used
2090to perform synchronisation and signaling between threads.
2091Example multithread programs using these classes can be found in
2092the {\bf Tests} module : \\
2093\hspace{10mm} {\tt zthr.cc , tmtdt.cc , tmtrnd.cc }
2094
2095\begin{figure}[hbt]
2096\dclsbb{ZThread}{ParalExThread}
2097\dclsa{ParallelExecutor}
2098\end{figure}
2099%% CHECK
2100{\bf \large DECRIRE l'utilisation des nouvelles classes ParallelExecutor ... }
2101
2102\subsection{Dynamic linker and C++ compiler classes}
2103\index{PDynLinkMgr}
2104The class {\bf PDynLinkMgr} can be used for managing shared libraries
2105at run time. The example below shows the run time linking of a function:\\
2106{\tt extern "C" { void myfunc(); } } \\
2107\begin{verbatim}
2108#include "pdlmgr.h"
2109// ...
2110string soname = "mylib.so";
2111string funcname = "myfunc";
2112PDynLinkMgr dyl(soname);
2113DlFunction f = dyl.GetFunction(funcname);
2114if (f != NULL) {
2115// Calling the function
2116 f();
2117}
2118\end{verbatim}
2119
2120\index{CxxCompilerLinker}
2121The {\bf CxxCompilerLinker} class provides the services to compile C++ code and building
2122shared libraries, using the same compiler and options which have
2123been used to create the SOPHYA shared library.
2124The sample program below illustrates using this class to build
2125the shared library (myfunc.so) from the source file myfunc.cc :
2126\begin{verbatim}
2127#include "cxxcmplnk.h"
2128// ...
2129string flnm = "myfunc.cc";
2130string oname, soname;
2131int rc;
2132CxxCompilerLinker cxx;
2133// The Compile method provides a default object file name
2134rc = cxx.Compile(flnm, oname);
2135if (rc != 0 ) { // Error when compiling ... }
2136// The BuildSO method provides a default shared object file name
2137rc = cxx.BuildSO(oname, soname);
2138if (rc != 0 ) { // Error when creating shared object ... }
2139\end{verbatim}
2140
2141\subsection{Command interpreter}
2142\index{Commander}
2143\index{Expression evaluation}
2144The class {\bf Commander} can be used in interactive programs to provide
2145c-shell like command interpreter and scripting capabilties.
2146Arithmetic expression evaluation is implemented through the {\bf CExpressionEvaluator}
2147and {\bf RPNExpressionEvaluator} classes. The latter can be used to perform evaluation
2148in reverse polish notation (old HP calculators),
2149while the former uses the usual algebraic (c-language like) expressions.
2150The command language provides variable manipulation through the usual
2151{\tt \$varname} vector variable and arithmetic expression extensions, as well
2152as the control and test blocs.
2153\begin{verbatim}
2154#include "commander.h"
2155...
2156Commander cmd;
2157char* ss[3] = {"foreach f ( AA bbb CCCC ddddd )", "echo $f" , "end"};
2158for(int k=0; k<3; k++) {
2159 string line = ss[k];
2160 cmd.Interpret(line);
2161}
2162\end{verbatim}
2163
2164\newpage
2165\section{Module SkyMap}
2166\begin{figure}[hbt]
2167\dclsbb{AnyDataObj}{PixelMap}
2168\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
2169\dclsc{SphereThetaPhi}
2170\dclsc{SphereECP}
2171\dclsb{LocalMap}
2172\caption{partial class diagram for spherical map classes in Sophya}
2173\end{figure}
2174The {\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.
2175\subsection{3D geometry}
2176Some of the classes in this module simplify the angle manipulation, or the 3D geometry
2177and rotation calculations, such as the {\bf Vector3d} and {\bf UnitVector} classes.
2178\index{Vector3d}
2179The classes {\bf Vector3d} and {\bf UnitVector} are useful for angle conversions.
2180\index{Angle}
2181{\bf LongLat} class and {\bf Angle} class (defined in file vector3d.h) can be used for
2182angle conversions :
2183\begin{verbatim}
2184 // Example to convert 0.035 radians to arcsec
2185 double vr = 0.035;
2186 cout << "Angle rad= " << vr << " ->arcsec= " << Angle(vr).ToArcSec() << endl;
2187 // Example to convert 2.3 arcmin to radian - we use the conversion operator
2188 double vam = 2.3;
2189 cout << "Angle arcmin= " << vam << " ->rad= "
2190 << (double)Angle(vam, Angle::ArcMin) << endl;
2191\end{verbatim}
2192%%%
2193\subsection {Spherical maps}
2194SkyMap module provides three classes for representing data with pixels distributed over complete ($4 \pi$ steradians) spheres. These classes implements the common interface defined
2195in the base class \tcls{SphericalMap} with three different algorithms or
2196pixelization scheme.
2197The spherical maps can be instanciated for the followind data types: \\
2198\hspace*{5mm} int\_4 , r\_4 (float) , r\_8 (double) , complex$<$r\_4$>$ , complex$<$r\_8$>$. \\
2199The SphereHEALPix can in addition be instanciated for T=uint\_2.
2200
2201\begin{enumerate}
2202\item \index{\tcls{SphereHEALPix}}
2203{\bf \tcls{SphereHEALPix}} implements the HEALPix
2204({\bf H}ierarchical {\bf E}qual {\bf A}rea iso{\bf L}atitude {\bf Pix}elization) scheme,
2205developped originaly by K. Gorski \& E. Hivon. Refer to
2206\href{http://healpix.jpl.nasa.gov/}{HEALPix home page} for
2207detailed information about this pixelisation scheme and related software
2208\footnote{HEALPix home page: http://healpix.jpl.nasa.gov/ }.
2209FITS read/write for SphereHEALPix objects is handled by the \tcls{FITS\_SphereHEALPix}
2210class in module FitsIOServer.
2211\item \index{\tcls{SphereThetaPhi}}
2212{\bf \tcls{SphereThetaPhi}} represents spheres pixelized following an algorithm
2213developed at LAL-ORSAY, for SOPHYA.
2214The sphere is divided into a number of rings or slices
2215along the parallels, corresponding to different values of the angle $\theta$.
2216Each slice is then divided into a number of pixels, with an aspect ratio close
2217to one (square pixels). Pixels are exactly iso-latitude and very uniform in surface,
2218over the sphere.
2219
2220\item \index{\tcls{SphereECP}}
2221
2222The {\bf \tcls{SphereECP}} class correspond to the cylindrical projection.
2223Like SphereThetaPhi class, the sphere is divided into a number of rings
2224or slices, and each ring is divided into a constant number of pixels
2225along the $\varphi$ direction. Although the \tcls{SphereECP} does not have
2226equal area pixels when used for complete spheres, it can be used for
2227representing partial or full spherical maps.
2228\end{enumerate}
2229
2230The example below shows creating and filling of a
2231SphereHEALPix with nside = 8
2232(The sphere will have $12 \times 8 \times 8= 768$ pixels) :
2233
2234\begin{verbatim}
2235#include "spherehealpix.h"
2236// ...
2237SphereHEALPix<double> sph(8);
2238for (int k=0; k< sph.NbPixels(); k++) sph(k) = (double)(10*k);
2239\end{verbatim}
2240
2241Pixels at an angular posistion can be directly accessed through the operator \\
2242\hspace*{15mm} T \tcls{SphericalMap}::operator()($\theta, \varphi$) :
2243\begin{verbatim}
2244#include "vector3d.h"
2245#include "spherethetaphi.h"
2246...
2247// Create a sphere with 40 arcmin resolution
2248int M = SphereThetaPhi<r_4>::ResolToSizeIndex(
2249 Angle(40., Angle::ArcMin).ToRadian() );
2250SphereThetaPhi<r_4> sphtp(M);
2251double tet, phi;
2252for (int k=0; k< sphtp.NbPixels(); k++) {
2253 sphtp.PixThetaPhi(k, tet, phi);
2254 sphtp(k) = cos(5.*tet)*sin(7.*phi);
2255}
2256cout << sphtp;
2257// To save sphtp to file sphtp (if executed through runcxx)
2258KeepObj(sphtp);
2259\end{verbatim}
2260
2261\index{\tcls{SphereThetaPhi}}
2262
2263\subsection {Local maps}
2264\index{\tcls{LocalMap}}
2265A 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).
2266
2267Internally, 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(...))
2268
2269The 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).
2270\begin{verbatim}
2271#include "localmap.h"
2272//..............
2273 LocalMap<r_4> locmap(4,5);
2274 for (int k=0; k<locmap.NbPixels();k++) locmap(k)=10.*k;
2275 locmap.SetOrigin();
2276 locmap.SetSize(30.,30.);
2277\end{verbatim}
2278
2279\subsection{Writing, viewing \dots }
2280
2281All these objects have been design to be written to or read from a persistant file.
2282The following example shows how to write the previously created objects
2283into such a file~:
2284\begin{verbatim}
2285//-- Writing
2286
2287#include "fiospherehealpix.h"
2288//................
2289
2290char *fileout = "myfile.ppf";
2291POutPersist outppf(fileout);
2292FIO_SphereHEALPix<r_8> outsph(sph);
2293outsph.Write(outppf);
2294FIO_LocalMap<r_8> outloc(locmap);
2295outloc.Write(outppf);
2296// It is also possible to use the << operator
2297POutPersist os("sph.ppf");
2298os << outsph;
2299os << outloc;
2300\end{verbatim}
2301
2302Sophya graphical tools (spiapp) can automatically display and operate
2303all these objects.
2304
2305\newpage
2306\section{Samba and SkyT}
2307\subsection{Samba}
2308\index{Spherical Harmonics}
2309\index{SphericalTransformServer}
2310The 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...
2311\begin{verbatim}
2312#include "skymap.h"
2313#include "samba.h"
2314....................
2315
2316// Generate input spectra a + b* l + c * gaussienne(l, 50, 20)
2317int lmax = 92;
2318Vector clin(lmax);
2319for(int l=0; l<lmax; l++) {
2320 double xx = (l-50.)/10.;
2321 clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
2322}
2323
2324// Compute map from spectra
2325SphericalTransformServer<r_8> ylmserver;
2326int m = 128; // HealPix pixelisation parameter
2327SphereHEALPix<r_8> map(m);
2328ylmserver.GenerateFromCl(map, m, clin, 0.);
2329// Compute power spectrum from map
2330Vector clout = ylmserver.DecomposeToCl(map, lmax, 0.);
2331\end{verbatim}
2332
2333\subsection{Module SkyT}
2334\index{RadSpectra} \index{SpectralResponse}
2335The SkyT module is composed of two types of classes:
2336\begin{itemize}
2337\item{} one which corresponds to an emission spectrum of
2338radiation, which is called RadSpectra
2339\item{} one which corresponds to the spectral response
2340of a given detector (i.e. corresponding to a detector
2341filter in a given frequency domain), which is called
2342SpectralResponse.
2343\end{itemize}
2344\begin{figure}[hbt]
2345\dclsbb{RadSpectra}{RadSpectraVec}
2346\dclsb{BlackBody}
2347\dclsccc{AnyDataObj}{SpectralResponse}{SpecRespVec}
2348\dclsc{GaussianFilter}
2349\caption{partial class for SkyT module}
2350\end{figure}
2351
2352\begin{verbatim}
2353#include "skyt.h"
2354// ....
2355// Compute the flux from a blackbody at 2.73 K through a square filter
2356BlackBody myBB(2.73);
2357// We define a square filter from 100 - 200 GHz
2358SquareFilter mySF(100,200);
2359// Compute the filtered integrated flux :
2360double flux = myBB.filteredIntegratedFlux(mySF);
2361\end{verbatim}
2362
2363A more detailed description of SkyT module can be found in:
2364{\it The SkyMixer (SkyT and PMixer modules) - Sophya Note No 2. }
2365available also from Sophya Web site.
2366
2367\newpage
2368\section{Module FitsIOServer}
2369This module provides classes for handling file input-output in FITS
2370\footnote{http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html}
2371format using the cfitsio library. Its
2372design is similar to the SOPHYA persistence (see Module BaseTools).
2373Delegate classes or handlers perform the actual read/write from/to fits files.
2374\par
2375Compared to the SOPHYA native persistence (PPF format),
2376FITS format has the advantage of being used extensively, and handled
2377by a many different software tools. It is a de facto standard in
2378astronomy and astrophysics.
2379However, FITS lacks some of the features present in SOPHYA PPF, and although
2380many SOPHYA objects can be saved in FITS files, FITS persistence has
2381some limitations. For example, FITS does not currently handle complex arrays.
2382\subsection{FITS streams}
2383\index{FITS} \index{FitsInOutFile}
2384%%
2385The class {\bf FitsInOutFile} can be seen as a wrapper class for the cfitsio library functions.
2386This class has been introduced in 2005 (V=1.9), when the module has been
2387extensively changed. In order to keep backward compatibility, the old fits wrapper
2388classes ({\bf FitsFile, FitsInFile, FitsOutFile}) has been changed to inherit from
2389 {\bf FitsInOutFile}. The use of class FitsFile and specific services of these old classes
2390 should be avoided, but FitsInFile, FitsOutFile can be safely considered a specialisation
2391 of FitsInOutFile for read/input and write/output operations respectively.
2392 Most c-fitsio errors are converted to an exception: {\bf FitsIOException}.
2393 \par
2394 File names are passed to cfitsio library. It is thus possible to use cfitsio file name conventions,
2395 such as {\bf ! } as the first character, for overwriting existing files (when creating files).
2396 The diagram below shows the class hierarchy for cfitsio wrapper classes.
2397\begin{figure}[hbt]
2398\dclsa{FitsInOutFile}
2399\dclscc{FitsFile}{FitsInFile}
2400\dclscc{FitsFile}{FitsOutFile}
2401\end{figure}
2402%%%%
2403\subsection{FITS handlers and I/O operators}
2404\index{FitsManager}
2405Handlers classes inheriting from {\bf FitsHandlerInterface} perform write/read operations
2406for AnyDataObj objects to/from FitsInOutFile streams. The {\bf FitsManager} class provides
2407top level services for object read/write in FITS files.
2408\par In most cases,
2409\hspace{5mm} {\tt FitsInOutFile\& $<<$ } \, and \, {\tt FitsInOutFile\& $>>$ } \hspace{5mm}
2410operators can be used to write and read objects.
2411When reading objects from a fits file using the {\tt FitsInOutFile\& $>>$ } operator,
2412the fits file is positioned on the next HDU, after reading. Also, if the {\bf FitsInOutFile}
2413object is positioned on a first empty HDU (without data, naxis=0), reading in objects
2414corresponding to a binary or ascii table using the operator $>>$ will skip automatically
2415the empty HDU and position the fits file on the second HDU, before trying to read in
2416the object.
2417\par
2418The two main types of fits data structures, images and tables
2419{\tt (IMAGE\_HDU , BINARY\_TBL , ASCII\_TBL)} are handled by the generic handlers: \\
2420{\bf \tcls{FitsArrayHandler}} and {\bf FitsHandler$<$BaseDataTable$>$}.
2421\par
2422A number of more specific handlers are also available, in particular for NTuple,
2423\tcls{SphereHealPix} and \tcls{SphereThetaPhi}. \\[2mm]
2424{\bf Warning:} Some handlers were written with the old FitsIOServer classes.
2425They inherit from the intermediate class {\bf FitsIOHandler} and
2426have been adapted to the new scheme. \\[2mm]
2427%%%
2428The examples below illustrates the usage of FitsIOServer classes. They can be compiled
2429and executed using runcxx, without the {\tt include} lines: \\[1mm]
2430\hspace*{5mm} {\tt csh> runcxx -import SkyMap -import FitsIOServer -inc fiosinit.h }
2431%%%
2432\begin{enumerate}
2433\item Saving an array and a HealPix map to a Fits file
2434\begin{verbatim}
2435#include "fitsioserver.h"
2436#include "fiosinit.h"
2437// ....
2438{
2439// Make sure FitsIOServer module is initialised :
2440FitsIOServerInit();
2441// Create and open a fits file named myfile.fits
2442FitsInOutFile fos("myfile.fits", FitsInOutFile ::Fits_Create);
2443// Create and save a 15x11 matrix of integers
2444TMatrix<int_4> mxi(15, 11);
2445mxi = RegularSequence(10.,10.);
2446fos << mxi;
2447// Save a HEALPix spherical map using FitsManager services
2448SphereHEALPix<r_8> sph(16);
2449sph = 48.3;
2450FitsManager::Write(fos, sph);
2451// --- The << operator could have been used instead : fos << sph;
2452}
2453\end{verbatim}
2454%%%%
2455%%%%
2456\item Reading objects and the header from the previously created fits file:
2457\begin{verbatim}
2458{
2459FitsIOServerInit(); // Initialisation
2460// ---- Open the fits file named myfile.fits
2461FitsInFile fis("myfile.fits");
2462//---- print file information on cout
2463cout << fis << endl;
2464//--- Read in the array
2465TArray<int_4> arr;
2466fis >> arr;
2467arr.Show();
2468//--- Position on second HDU
2469fis.MoveAbsToHDU(2);
2470//--- read and display header information
2471DVList hdu2;
2472fis.GetHeaderRecords(hdu2, true, true);
2473cout << hdu2;
2474//--- read in the HEALPix map
2475SphereHEALPix<r_8> sph;
2476FitsManager::Read(fis, sph);
2477// --- The >> operator could have been used instead : fis >> sph;
2478sph.Show();
2479}
2480\end{verbatim}
2481%%%%%%%
2482%%%
2483\item DataTable objects can be read from and written to FITS files as ASCII or
2484binary tables. The example belo show reading the DataTable created in the example
2485in section \ref{datatables} from a PPF file and saving it to a fits file.
2486\begin{verbatim}
2487#include "swfitsdtable.h"
2488// ....
2489{
2490FitsIOServerInit(); // FitsIOServer Initialisation
2491//--- Reading in DataTable object from PPF file
2492PInPersist pin("dtable.ppf");
2493DataTable dt;
2494pin >> dt;
2495dt.Show();
2496//--- Saving table to FITS
2497FitsInOutFile fos("!dtable.fits", FitsInOutFile ::Fits_Create);
2498fos << dt;
2499}
2500\end{verbatim}
2501%%%%
2502\end{enumerate}
2503%%%
2504A partial class diagram of FITS persistence handling classes is shown below. The
2505class {\tt FitsIOhandler} conforms to the old FitsIOServer module design and
2506should not be used anymore.
2507\begin{figure}[hbt]
2508\dclsbb{FitsHandlerInterface}{FitsArrayHandler$<$T$>$}
2509\dclsb{\tcls{FitsHandler}}
2510\dclscc{FitsIOhandler}{FITS\_NTuple}
2511\dclsc{FITS\_SphereHEALPix}
2512% \dclsb{FITS\_LocalMap}
2513\end{figure}
2514
2515\subsection{SwFitsDataTable and other classes}
2516\label{SwFitsDataTable}
2517\index{SwFitsDataTable}
2518The {\bf SwFitsDataTable} class implements the BaseDataTable interface
2519using a FITS file as swap space. Compared to SwPPFDataTable, they can be
2520used in R/W mode (reading from the table, when it is being created / filled).
2521They can be used in a way similar to DataTable and SwPPFDataTable.
2522When creating the table, a {\tt FitsInOutFile } stream, opened for writing has
2523to be passed to the creator. No further operation is needed.
2524\begin{verbatim}
2525// ....
2526FitsInOutFile so("!myswtable.fits", FitsInOutFile::Fits_Create);
2527SwFitsDataTable dt(so, 16);
2528// define table columns
2529dt.AddIntegerColumn("X0_i");
2530dt.AddFloatColumn("X1_f");
2531// ... Fill the table
2532r_8 x[5];
2533for(int i=0; i<63; i++) {
2534 x[0] = (i%9)-4.; x[1] = (i/9)-3.;
2535 dt.AddLine(x);
2536}
2537\end{verbatim}
2538The class {\bf FitsBTNtuIntf } provide an alternative tool to read FITS tables.
2539{\bf FitsABTColRd} , {\bf FitsABTWriter } and {\bf FitsImg2DWriter } can also
2540be used to manipulate FITS files.
2541\par
2542The {\bf scanfits} program can be used to check FITS files and analyse their
2543content (See \ref{scanfits}).
2544
2545%%%%
2546\newpage
2547\section{LinAlg and IFFTW modules}
2548An interface to use LAPACK library (available from {\tt http://www.netlib.org})
2549is implemented by the {\bf LapackServer} class, in module LinAlg.
2550\index{LapackServer}.
2551The sample code below shows how to use SVD (Singular Value Decomposition)
2552through LapackServer:
2553\begin{verbatim}
2554#include "intflapack.h"
2555// ...
2556// Use FortranMemoryMapping as default
2557BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
2558// Create an fill the arrays A and its copy AA
2559int n = 20;
2560Matrix A(n , n), AA;
2561A = RandomSequence(RandomSequence::Gaussian, 0., 4.);
2562AA = A; // AA is a copy of A
2563// Compute the SVD decomposition
2564Vector S; // Vector of singular values
2565Matrix U, VT;
2566LapackServer<r_8> lpks;
2567lpks.SVD(AA, S, U, VT);
2568// We create a diagonal matrix using S
2569Matrix SM(n, n);
2570for(int k=0; k<n; k++) SM(k,k) = S(k);
2571// Check the result : A = U*SM*VT
2572Matrix diff = U*(SM*VT) - A;
2573double min, max;
2574diff.MinMax(min, max);
2575cout << " Min/Max difference Matrix (?=0) , Min= " << min
2576 << " Max= " << max << endl;
2577\end{verbatim}
2578
2579\index{FFTWServer}
2580The {\bf FFTWServer} class (in module FFTW) implements FFTServerInterface class
2581methods, for one dimensional and multi-dimensional Fourier
2582transforms using the FFTW package (available from {\tt http://www.fftw.org}).
2583Depending on the configuration options during library build, only
2584transforms on double precision data might be available.
2585
2586\section{Multi-thread and parallel programming with Sophya}
2587
2588\newpage
2589\section{Building and installing Sophya}
2590\subsection{Supported platforms}
2591Presently, the SOPHYA and PI libraries and tools has been compiled and tested with
2592the following compiler/platform pairs:
2593
2594\begin{center}
2595\begin{tabular}{|l|l|}
2596\hline
2597OS & compiler \\
2598\hline
2599Linux (SL5, 64) & g++ 4.1 \\
2600Linux (SL5, 64) & icc - Intel compiler (10.1) \\
2601Linux Ubuntu & g++ 4.4 \\
2602MacOSX/Darwin 10.3 (PowerPC) \, 10.4 (Intel) & Apple g++ 3.3 \, 4.0 \\
2603MacOSX/Darwin 10.5 (Intel) & Apple g++ 4.0 \\
2604MacOSX/Darwin 10.6 (Intel) Snow Leopard & Apple g++ 4.2 \\
2605IBM AIX 5.3 & xlC 8.0 \\
2606HP/Compaq/DEC Tru64 ( OSF1) & cxx 6.1 \, 6.3 \\
2607\hline
2608\end{tabular}
2609\end{center}
2610
2611The SOPHYA and PI library and associated tools (spiapp, runcxx \ldots) were ported and tested
2612on the Silicon Graphics (SGI) IRIX system and the corresponding native c++ compiler ( IRIX64, CC 7.3 ),
2613up to SOPHYA V2.1 (V\_Nov2007).
2614
2615\subsection{Library and makefile structure}
2616%
2617The object files from a given Sophya module are grouped in an archive library
2618with the module's name ({\tt libmodulename.a}). All Sophya modules
2619 are grouped in a single shared library ({\tt libsophya.so}), while the
2620modules with reference to external libraries are grouped in
2621({\tt libextsophya.so}). The {\bf PI} and {\bf PIext} modules are
2622grouped in ({\tt libPI.so}).
2623Alternatively, it is possible to group all modules in a single shared
2624library {\tt libAsophyaextPI.so}.
2625\par
2626Each library module has a {\tt Makefile} which compiles the source files
2627and build the correspond static (archive) library using the compilation
2628rules and flags defined in \\
2629\hspace*{5mm} {\tt \$SOPHYABASE/include/sophyamake.inc}. \\
2630Each program module has a {\tt Makefile} which compiles and link the
2631corresponding programs using the compilation rules and libraries
2632defined in {\$SOPHYABASE/include/sophyamake.inc}.
2633The top level Makefile in BuildMgr/ compiles each library modules
2634and builds shared libraries.
2635\par
2636Some of the modules in the Sophya package uses external libraries. The
2637{\bf FitsIOServer} is an example of such a module, where the {\tt libcfitsio.a}
2638is used. The list of all Sophya modules using external libraries is
2639presented in section \ref{sopmodules}.
2640The external libraries should be installed before the configure step
2641(see below) and the compilation of the corresponding Sophya modules.
2642\par
2643The series of Makefiles use the link to {\tt sophyamake.inc} in BuildMgr.
2644There are also the {\tt smakefile} series which uses the explicit path, using
2645{\tt \$SOPHYABASE} environment variable.
2646
2647\subsection{Build instructions}
2648\label{build}
2649The build procedure has two main steps:
2650\begin{enumerate}
2651\item The configure step (BuildMgr/configure) setup the directory structure and
2652the necessary configuration file. Refer to section \ref{directories} for
2653the description of SOPHYA directory tree and files.
2654\item The make step compiles the different sources files, create the library and optionaly
2655builds all or some of the associated executables.
2656\end{enumerate}
2657
2658{\tt BuildMgr/configure } is a c-shell script with a number of arguments:
2659\begin{verbatim}
2660csh> ./configure -h
2661configure [-sbase SOPHYABASE] [-scxx SOPHYACXX] [-incln]
2662 [-minc mymake.inc] [-compopt 'cc/cxxOptions']
2663 [-arch64] [-arch32] [-sasz64] [-ldble128] [-nofpic] [-nothsafe] [-boundcheck] [-sodebug]
2664 [-extp dir1 -extp dir2 ...] [-extip dir1 -extip dir2 ... ] [-extlp dir1 -extlp dir2 ... ]
2665 [-noextlib] [-noext fits] [-noext fftw] [-noext lapack] [-noext astro]
2666 [-noPI] [-slballinone]
2667 [-alsofftwfloat] [-usefftw2] [-uselapack2]
2668 (See SOPHYA manual/web pages for a detailed description of configure options)
2669\end{verbatim}
2670\begin{itemize}
2671\item[] -sbase : define SOPHYA installation base directory. \$SOPHYABASE is used
2672if not specified.
2673\item[] -scxx : selects the C++ compiler. \$SOPHYACXX s used
2674if not specified.
2675\item[] -incln : creates symbolic link for include files, instead of copying them.
2676\item[] -minc : give an explicit name for the file used to generate
2677\$SOPHYABASE/include/sophyamake.inc. If {\tt -minc} is not specified, one of
2678the files in BuildMgr/ directory is selected, based on the system name and the
2679compiler {\tt Linux\_g++\_make.inc , OSF1\_cxx\_make.inc , AIX\_xlC\_make.inc \ldots}
2680\item[] -compopt : additional compiler options \\[2mm]
2681\item[] -arch64 : select/force 64 bits compilation mode. This is the default on 64 bits Linux systems and AIX.
2682\item[] -arch32 : select/force 32 bits compilation mode. useful on 64 bits Linux systems or AIX
2683\item[] -sasz64 : select 8 byte (64 bits) size for array indices, useful if you need to allocate an manipulate
2684large arrays, with more than $2^32$ elements along a single array dimension.
2685\item[] -ldble128 : set the flags activating {\tt long double} (128 bits) arrays. \\[2mm]
2686\item[] -nofpic : disable -fPIC (Position Independent Code) generation flag
2687\item[] -nothsafe : disable thread-safety procedures in SOPHYA : reference sharing \ldots.
2688 ( activate the conditional compilation flag {\tt SO\_NOTHSAFE } )
2689\item[] -boundcheck : compile SOPHYA and bound checking activated when accessing array elements
2690( activate conditional compilation flag {\tt SO\_BOUNDCHECKING } )
2691\item[] -sodebug : activate conditional compilation flag {\tt SOPHYA\_DEBUG } \\[2mm]
2692\item[] -extp : Adds the specied path to the search path of the external libraries
2693include files and archive library.
2694\item[] -extip : Adds the specied path to the search path of the external libraries
2695include files.
2696\item[] -extp : Adds the specied path to the search path of the external libraries
2697archive (libxxx.a).
2698\item[] -noextlib : Disable compiling of modules referencing external libraries.
2699\item[] -noext : Disable compiling of the specified module (with reference to external
2700library : {\tt -noext fits , -noext fftw -noext lapack -noext astro } \\[2mm]
2701\item[] -usefftw2: Use FFTW V2 instead of the default FFTW V3 - A preprocessor
2702flag will be defined in sspvflags.h
2703\item[] -uselapack2: Lapack V2 is being used (defaulr V3) - A preprocessor
2704flag will be defined in sspvflags.h
2705\item[] -alsofftwfloat : compile single precision (float) version of the Fourier
2706transform methods (module IFFTW, class FFTWServer). A preprocessor
2707flag will be defined in sspvflags.h and float version of the FFTW library
2708(libfftw3f.a) will be linked with SOPHYA, in addition to the default double
2709precision library (libfftw3.a).
2710\item[] -noPI : has currently no effect
2711\item[] -slballinone: A single shared library for all SOPHYA, PI and external library interface
2712modules will be build. A compilation flag
2713will be defined in sspvflags.h . \\ See also target {\tt slballinone} below.
2714\end{itemize}
2715
2716{\large \bf configure steps } \\[1mm]
2717The configure script performs the following actions :
2718\begin{enumerate}
2719\item Creates directory tree under {\tt \$SOPHYABASE }
2720\item Copy include files to {\tt \$SOPHYABASE/include/ } (or creates symbolic link)
2721\item Search for external libraries include files and create the necessary links
2722in {\tt \$SOPHYABASE/include}
2723\item Search for external libraries (-lfits \ldots) and add the corresponding
2724directories to the library search path, in {\tt sophyamake.inc}
2725\item Creates the file { \tt \$SOPHYABASE/include/sophyamake.inc }
2726\item Creates the file {\tt machdefs.h} from { BaseTools/machdefs\_mkmf.h } and
2727{\tt sspvflags.h }
2728\item Creates the object list files for shared library creation
2729\end{enumerate}
2730
2731{\large \bf Example } \\[1mm]
2732
2733In the example below, we assume that we want to install Sophya from a
2734released (tagged) version in the source directory {\tt \$SRC} in the
2735{\tt /usr/local/Sophya} directory, using {\tt g++}. We assume that
2736the external libraries can be found in {\tt /usr/local/ExtLibs/}.
2737We disable the compilation of the XAstroPack package.
2738
2739\vspace*{3mm}
2740\begin{verbatim}
2741# Create the top level directory
2742csh> mkdir /usr/local/Sophya/
2743csh> cd $SRC/BuildMgr
2744# Step 1.a : Run the configuration script
2745csh> ./configure -sbase /usr/local/Sophya -scxx g++ -extp /usr/local/ExtLibs/ \
2746-noext astro
2747# Step 1.b : Check the generated file $SOPHYABASE/include/sophyamake.inc
2748csh> ls -lt *.inc
2749csh> more sophyamake.inc
2750\end{verbatim}
2751If necessary, edit the generated file {\tt sophyamake.inc } in order to modify
2752compilation flags, library list. The file is rather short and self documented.
2753\begin{verbatim}
2754# Step 2.a: Compile the modules without external library reference
2755csh> make libs
2756# Step 2.b: Compile the modules WITH external library reference (optional)
2757csh> make extlibs
2758# Step 2.c: Build libsophya.so
2759csh> make slb
2760# Step 2.d: Build libextsophya.so (optional)
2761csh> make slbext
2762# Step 2.e: Compile the PI and PIext modules (optional)
2763csh> make PI
2764# Step 2.f: Build the corresponding shared library libPI.so (optional)
2765csh> make slbpi
2766\end{verbatim}
2767
2768To compile all modules and build the shared libraries, it is possible
2769to perform the steps 2.a to 2.f using the targets {\tt all} and {\tt slball}
2770defined in the Makefile
2771\begin{verbatim}
2772# Step 2.a ... 2.f
2773csh> make all slball
2774\end{verbatim}
2775
2776It is also possible to group all modules in a single shared library using
2777the target {\tt slballinone}.
2778\begin{verbatim}
2779# Step 2.a ... 2.f
2780csh> make all slballinone
2781\end{verbatim}
2782
2783At this step, all libraries should have been made. Programs using
2784Sophya libraries can now be built:
2785\begin{verbatim}
2786# To compile some of the test programs
2787csh> make basetests
2788# To compile runcxx , scanppf , scanfits
2789csh> make prgutil
2790# To build (s)piapp (libPI.so is needed)
2791csh> make piapp
2792\end{verbatim}
2793
2794If no further modification or update of source files is foreseen, it is possible
2795to remove all .o files:
2796\begin{verbatim}
2797# To clean $SOPHYABASE/obj directory :
2798csh> make cleanobj
2799\end{verbatim}
2800
2801
2802\subsection{Notes}
2803\begin{itemize}
2804\item[{\bf Makefile}] List of top level Makefile build targets
2805\begin{verbatim}
2806> libs extlibs PI = all
2807> slb slbext slbpi = slball (OR = slballinone)
2808> clean cleanobj
2809> tests prgutil prgmap progpi = prgall
2810> basetests piapp (ou progpi) pmixer
2811\end{verbatim}
2812\item[{\bf MacOS X}] A high performance mathematic and signal processing
2813library, including LAPACK and BLAS is packaged in Darwin/MacOS X (10.3, 10.4) : \\
2814\hspace*{5mm} {\bf -framework Accelerate}
2815\item[{\bf Tru64/OSF}] An optimised math library with LAPACK and BLAS might
2816optionaly be installed {\bf (-lcxlm) }. On our system, this libray contained Lapack V2.
2817So we used the LAPACK, as compiled from the public sources, and linked with
2818the Tru64 native BLAS.
2819\item[{\bf IRIX64}] We used the math library with LAPACK V2 and BLAS
2820from SGI : {\bf -lcomplib.sgimath}
2821\item[{\bf AIX}] There seem to be a problem on AIX when several shared
2822libraries are used. We have been able to run SOPHYA programs either
2823using static libraries, or a single shared library (libAsophyaextPI.so)
2824if extlibs and PI are needed, in addition to stand alone SOPHYA modules.
2825It has not been possible to link SOPHYA with fortran libraries
2826\item[{\bf Mgr}] This module contains makefiles and build scripts
2827that were used in SOPHYA up to version 1.7 (2004) : OBSOLETE.
2828\end{itemize}
2829
2830\subsection{Files and scripts in BuildMgr/ }
2831\begin{itemize}
2832\item[] {\bf Makefile:} Top level Makefile for building SOPHYA.
2833{\tt smakefile} is similar to Makefile, except that it uses
2834the smakefiles in each module.
2835\item[] {\bf mkmflib:} c-shell script for creation of library module
2836Makefile / smakefile. \\
2837\hspace*{5mm} {\tt ./mkmflib -sbase /tmp/sbase SUtils }
2838\item[] {\b mkmfprog:}
2839c-shell script for creation of programs module Makefile / smakefile \\
2840\hspace*{5mm} {\tt ./mkmfprog -sbase /tmp/sbase ProgPI }
2841\item[] {\bf domkmf:} c-shell script - calls mkmflib for all modules \\
2842\hspace*{5mm} {\tt ./domkmf -sbase /tmp/sbase}
2843\item[] {\bf xxx\_make.inc:} Configuration files for different compilers and OS
2844{\tt ( Linux\_g++\_make.inc , OSF1\_cxx\_make.inc \ldots )}.
2845These files are used to generate {\tt sophyamake.inc}
2846\end{itemize}
2847
2848
2849
2850\newpage
2851\appendix
2852\section{SOPHYA Exceptions}
2853\index{Exception classes} \index{PThrowable} \index{PError} \index{PException}
2854SOPHYA library defines a set of exceptions which are used
2855for signalling error conditions. From version V2.2 , all SOPHYA exceptions inherit from
2856the standard C++ exception class ( {\bf std::std::exception}). The method {\tt const char* what()} is
2857thus implemented and returns the corresponding error message.
2858The figure below shows a partial class diagram for exception classes in SOPHYA.
2859\begin{figure}[hbt]
2860\dclsbb{PThrowable}{PError}
2861\dclscc{PError}{AllocationError}
2862\dclscc{PError}{NullPtrError}
2863\dclscc{PError}{ForbiddenError}
2864\dclscc{PError}{AssertionFailedError}
2865\dclsbb{PThrowable}{PException}
2866\dclscc{PException}{IOExc}
2867\dclscc{PException}{SzMismatchError}
2868\dclscc{PException}{RangeCheckError}
2869\dclscc{PException}{ParmError}
2870\dclscc{PException}{TypeMismatchExc}
2871\dclscc{PException}{MathExc}
2872\dclscc{PException}{CaughtSignalExc}
2873\caption{partial class diagram for exception handling in Sophya}
2874\end{figure}
2875
2876For simple programs, it is a good practice to handle
2877the exceptions at least at high level, in the {\tt main()} function.
2878The example below shows the exception handling and the usage
2879of Sophya persistence.
2880
2881\input{ex1.inc}
2882
2883\newpage
2884\addcontentsline{toc}{section}{Index}
2885\printindex
2886\end{document}
2887
Note: See TracBrowser for help on using the repository browser.