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

Last change on this file since 2911 was 2835, checked in by ansari, 20 years ago

MAJ fichier liste de modifs - suite modifs TimeStamp et DataTable - Reza 8/11/2005

File size: 60.4 KB
Line 
1\documentclass[twoside,11pt]{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\begin{document}
29
30\begin{titlepage}
31% The title page - top of the page with the title of the paper
32\titrehp{Sophya \\ An overview }
33% Authors list
34\auteurs{
35R. Ansari & ansari@lal.in2p3.fr \\
36E. Aubourg & aubourg@hep.saclay.cea.fr \\
37G. Le Meur & lemeur@lal.in2p3.fr \\
38C. Magneville & cmv@hep.saclay.cea.fr \\
39S. Henrot-Versille & versille@in2p3.fr
40}
41% \auteursall
42% The title page - bottom of the page with the paper number
43\vspace{1cm}
44\begin{center}
45{\bf \Large Sophya Version: 1.9 (V\_Mai2005) }
46% Document revision 1.0 }
47\end{center}
48\titrebp{1}
49\end{titlepage}
50
51\tableofcontents
52
53\newpage
54
55\section{Introduction}
56
57{\bf SOPHYA} ({\bf SO}ftware for {\bf PHY}sics {\bf A}nalysis)
58is a collection of C++ classes designed for numerical and
59physics analysis software development. Our goal is to provide
60easy to use, yet powerful classes which can be used by scientists.
61Although some of the SOPHYA modules (SkyMap, Samba, SkyT)
62have been designed with the specific goal CMB data analysis, most
63modules presented here have a much broader scope and can be
64used in scientific data analysis and modeling/simulation.
65Whenever possible, we use existing numerical package and libraries,
66encapsulating them whenever
67possible.
68\par
69\vspace*{2mm}
70This documents
71presents only a brief overview of the class library,
72mainly from the user's point of view. A more complete description
73can be found in the reference manual, available from the SOPHYA
74web site: % {\bf http://www.sophya.org}.
75\href{http://www.sophya.org}{http://www.sophya.org}.
76\par
77\vspace*{2mm}
78The source directory tree
79\footnote{ CVS: cvsserver.lal.in2p3.fr:/exp/eros/CVSSophya}
80is organised into a number of modules.
81
82\begin{itemize}
83\item[] {\bf BuildMgr/} Scripts for code management,
84makefile generation and software installation
85\item[] {\bf BaseTools/} General architecture support classes such
86as {\tt PPersist, NDataBlock<T>}, and few utility classes
87such as the dynamic variable list manager ({\tt DVList}) as well
88as the basic set of exception classes used in SOPHYA.
89\item[] {\bf TArray/} template numerical arrays, vectors and matrices \\
90({\tt TArray<T> TMatrix<T> TVector<T> } \ldots)
91\item[] {\bf NTools/} Some standard numerical analysis tools
92(linear, and non linear parameter fitting, FFT, \ldots)
93\item[] {\bf HiStats/} Histogram-ming and data set handling classes (tuples) \\
94({\tt Histo Histo2D NTuple DataTable} \ldots)
95\item[] {\bf SkyMap/} Local and full sky maps, and some 3D geometry
96handling utility classes. \\
97({\tt PixelMap<T>, LocalMap<T>, SphericalMap<T>, \ldots})
98\item[] {\bf SUtils/} This module contains few utility classes, such as the
99{\tt DataCard} class, as well as string manipulation functions in C and C++.
100\item[] {\bf SysTools/} This module contains classes implementing
101an interface to various OS specific services, such
102threads and dynamic link/shared library handling.
103
104\end{itemize}
105
106The modules listed below are more tightly related to the
107CMB (Cosmic Microwave Background) data analysis problem:
108\begin{itemize}
109\item[] {\bf SkyT/}
110classes for spectral emission and detector frequency response modelling \\
111({\tt SpectralResponse, RadSpectra, BlackBody} \ldots)
112\item[] {\bf Samba/} Spherical harmonic analysis, noise generators \ldots
113\end{itemize}
114
115The following modules contain the interface classes with
116external libraries:
117\begin{itemize}
118\item[] {\bf FitsIOServer/} Classes for handling file input-output
119in FITS format using the cfitsio library.
120\item[] {\bf LinAlg/} Interface with Lapack linear algebra package
121\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a)
122\item[] {\bf XAstroPack/} Interface to some common astronomical
123computation libraries. Presently, this module uses an external library
124extracted from the {\bf Xephem } source code. The corresponding source
125code is also available from SOPHYA cvs repository, module {\bf XephemAstroLib}.
126\item[] {\bf MinuitAdapt/} Wrapper classes to CERN minimization routines (Minuit).
127
128\end{itemize}
129
130The following modules contain each a set of related programs using the
131SOPHYA library.
132\begin{itemize}
133\item[] {\bf Tests/} Simple test programs
134\item[] {\bf PrgUtil/} Various utility programs (runcxx, scanppf, scanfits, \ldots)
135\item[] {\bf PrgMap/} Programs performing operations on skymaps: projections,
136power spectrum in harmonic space, \ldots
137\item[] {\bf PMixer/} skymixer and related programs
138\end{itemize}
139
140As a companion to SOPHYA, the {\bf (s)piapp} interactive data analysis
141program is built on top of SOPHYA and the {\bf PI} GUI class library
142and application framework. The {\bf PI} ({\bf P}eida {\bf Interactive})
143development started in 1995, in the EROS \footnote{EROS: {\bf E}xp\'erience
144de {\bf R}echerche d'{\bf O}bjets {\bf S}ombres - http://eros.in2p3.fr}
145microlensing search collaboration, with PEIDA++ \footnote {PEIDA++:
146The EROS data analysis class library -
147http://www.lal.in2p3.fr/recherche/eros/PeidaDoc/}.
148The {\bf PI} documentation and the {\bf piapp} user's guide are available
149from \href{http://www.sophya.org}{http://www.sophya.org}.
150%\href{http://www.sophya.org}{http://www.sophya.org}.
151The {\bf PI} is organized as the following modules:
152\begin{itemize}
153\item[] {\bf PI/} Portable GUI class library and application development
154framework kernel.
155\item[] {\bf PIGcont/} Contour-plot drawing classes.
156\item[] {\bf PIext/} Specific drawers and adapters for SOPHYA objects,
157and the {\bf piapp} interactive data analysis framework.
158\item[] {\bf ProgPI/} interactive analysis tool main program and pre-loaded
159modules.
160\end{itemize}
161
162Modules containing examples and demo programs:
163\begin{itemize}
164\item[] {\bf Examples/} Sample SOPHYA codes and example programs and
165makefiles.
166\item[] {\bf DemoPIApp/} Sample scripts and programs for (s)piapp
167interactive analysis tools.
168\end{itemize}
169\newpage
170
171\section{Using Sophya}
172Basic usage of Sophya classes are described in in the following sections.
173Complete Sophya documentation can be found at our web site
174{\bf http://www.sophya.org}.
175
176\subsection{Directories, environment variables, configuration files}
177The environment variable {\bf SOPHYABASE} is used
178to define the path where the Sophya libraries and executable are installed.
179\begin{itemize}
180\item \$SOPHYABASE/include : Include (.h) files
181\item \$SOPHYABASE/lib : Path for the archive libraries (.a)
182\item \$SOPHYABASE/slb: Shared library path (.so)
183\item \$SOPHYABASE/exe : Executable file path
184\end{itemize}
185
186In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
187should contain the Sophya shared library path
188({\tt \$SOPHYABASE/slb }). \\
189On Silocon Graphics machines (IRIX64) , the default SOPHYA configuration
190correspond to the 64 bit architecture.
191The environment variable { \bf LD\_LIBRARY64\_PATH } defines
192the shared library path in this case and should contain ({\tt \$SOPHYABASE/slb }.
193
194The configure script creates links for external libraries include files in :
195\begin{itemize}
196\item \$SOPHYABASE/include/FitsIO : c-fitsio library include files
197\item \$SOPHYABASE/include/FFTW : FFTW library include files
198\item \$SOPHYABASE/include/XAstro : XEphem library include files
199\end{itemize}
200
201The directory { \tt \$SOPHYABASE/include/SophyaConfInfo/ } contains files
202describing the installed configuration of SOPHYA software.
203
204The file { \tt \$SOPHYABASE/include/machdefs.h } contains definitions
205(flags, typedef) used in SOPHYA.
206
207The file { \tt \$SOPHYABASE/include/sophyamake.inc } contains the
208compilation commands and flags used for building the software.
209Users can use most of compilation and link commands defined in this file:
210 {\tt \$CCOMPILE , \$CXXCOMPILE . \$CXXLINK \ldots}.
211 (See module Example).
212
213The configure script (BuildMgr/configure) creates the directory tree and the
214above files.
215
216\subsection{the runcxx program}
217\index{runcxx}
218{\bf runcxx} is a simple program which can be used to compile, link
219and run simple C++ programs. It handles the creation of a
220complete program file, containing the basic set C++ include files,
221the necessary include files for SOPHYA SysTools, TArray, HiStats
222and NTools modules, and the main program with exception handling.
223Other Sophya modules can be included using the {\tt -import} flag.
224Use of additional include files can be specified using the
225{\tt -inc} flag.
226\begin{verbatim}
227csh> runcxx -h
228 PIOPersist::Initialize() Starting Sophya Persistence management service
229SOPHYA Version 1.9 Revision 0 (V_Mai2005) -- May 31 2005 15:11:32 cxx
230 runcxx : compiling and running of a piece of C++ code
231 Usage: runcxx [-compopt CompileOptions] [-linkopt LinkOptions]
232 [-tmpdir TmpDirectory] [-f C++CodeFileName]
233 [-inc includefile] [-inc includefile ...]
234 [-import modulename] [-import modulename ...]
235 [-uarg UserArg1 UserArg2 ...]
236 if no file name is specified, read from standard input
237 modulenames: SkyMap, Samba, SkyT, FitsIOServer,
238 LinAlg, IFFTW, XAstroPack
239\end{verbatim}
240Most examples in this manual can be tested using runcxx. The
241example below shows how to compile, link and run a sample
242code.
243\begin{verbatim}
244// File example.icc
245Matrix a(3,3);
246a = IdentityMatrix(1.);
247cout << a ;
248// Executing this sample code
249csh> runcxx -f example.icc
250\end{verbatim}
251
252\subsection{the scanppf program}
253{\bf scanppf} is a simple SOPHYA application which can be used to check
254PPF files and list their contents. It can also provide the list of all registered
255PPF handlers.
256\begin{verbatim}
257csh> scanppf -h
258 PIOPersist::Initialize() Starting Sophya Persistence management service
259SOPHYA Version 1.9 Revision 0 (V_Mai2005) -- May 31 2005 15:11:32 cxx
260 Usage: scanppf [flags] filename
261 flags = -s -n -a0 -a1 -a2 -a3 -lh -lho
262 -s[=default} : Sequential reading of objects
263 -n : Object reading at NameTags
264 -a0...a3 : Tag List with PInPersist.AnalyseTags(0...3)
265 -lh : List PPersist handler classes
266 -lho : List PPersist handler and dataobj classes
267\end{verbatim}
268
269
270\newpage
271
272\section{Copy constructor and assignment operator}
273In C++, objects can be copied by assignment or by initialisation.
274Copying by initialisation corresponds to creating an object and
275initialising its value through the copy constructor.
276The copy constructor has its first argument as a reference, or
277const reference to the object's class type. It can have
278more arguments, if default values are provided.
279Copying by assignment applies to an existing object and
280is performed through the assignment operator (=).
281The copy constructor implements this for identical type objects:
282\begin{verbatim}
283class MyObject {
284public:
285 MyObject(); // Default constructor
286 MyObject(MyObject const & a); // Copy constructor
287 MyObject & operator = (MyObject const & a) // Assignment operator
288}
289\end{verbatim}
290The copy constructors play an important role, as they are
291called when class objects are passed by value,
292returned by value, or thrown as an exception.
293\begin{verbatim}
294// A function declaration with an argument of type MyObject,
295// passed by value, and returning a MyObject
296MyObject f(MyObject x)
297{
298 MyObject r;
299 ...
300 return(r); // Copy constructor is called here
301}
302// Calling the function :
303MyObject a;
304f(a); // Copy constructor called for a
305\end{verbatim}
306It should be noted that the C++ syntax is ambiguous for the
307assignment operator. {\tt MyObject x; x=y; } and
308{\tt MyObject x=y;} have different meaning.
309\begin{verbatim}
310MyObject a; // default constructor call
311MyObject b(a); // copy constructor call
312MyObject bb = a; // identical to bb(a) : copy constructor call
313MyObject c; // default constructor call
314c = a; // assignment operator call
315\end{verbatim}
316
317As a general rule in SOPHYA, objects which implements
318reference sharing on their data members have a copy constructor
319which shares the data, while the assignment operator copies or
320duplicate the data.
321
322\newpage
323\section{Module BaseTools}
324
325{\bf BaseTools} contains utility classes such as
326{\tt DVlist}, an hierarchy of exception classes for Sophya, a template
327class {\tcls{NDataBlock}} for handling reference counting on numerical
328arrays, as well as classes providing the services for implementing simple
329serialisation.
330\vspace*{5mm}
331
332\subsection{SOPHYA persistence}
333\index{PPersist} \index{PInPersist} \index{POutPersist}
334\begin{figure}[hbt]
335\dclsa{PPersist}
336\dclsccc{PPFBinarIOStream}{PPFBinaryInputStream}{PInPersist}
337\dclscc{PPFBinaryOutputStream}{POutPersist}
338\caption{partial class diagram for classes handling persistence in Sophya}
339\end{figure}
340A simple persistence mechanism is defined in SOPHYA. Its main
341features are:
342\begin{itemize}
343\item[] Portable file format, containing the description of the data structures
344and object hierarchy. \\
345{\bf PPF} {\bf P}ortable {\bf P}ersistence file {\bf F}ormat.
346\index{PPF}
347\item[] Handling of read/write for multiply referenced objects.
348\item[] All write operations are carried using sequential access only. This
349holds also for read operations, unless positional tags are used.
350SOPHYA persistence services can thus be used to transfer objects
351through network links.
352\item[] The serialisation (reading/writing) for objects for a given class
353is implemented through a handler object. The handler class inherits
354from {\tt PPersist} class.
355\item[] A run time registration mechanism is used in conjunction with
356RTTI (Run Time Type Identification) for identifying handler classes
357when reading {\bf PInPersist} streams, or for associating handlers
358with data objects {\bf AnyDataObject} for write operations.
359\end{itemize}
360A complete description of SOPHYA persistence mechanism and guidelines
361for writing delegate classes for handling object persistence is beyond
362the scope of this document. The most useful methods for using Sophya
363persistence are listed below:
364\begin{itemize}
365\item[] {\tt POutPersist::PutObject(AnyDataObj \& o)} \\
366Writes the data object {\bf o} to the output stream.
367\item[] {\tt POutPersist::PutObject(AnyDataObj \& o, string tagname)} \\
368Writes the data object {\bf o} to the output stream, associated with an
369identification tag {\bf tagname}.
370\item[] {\tt PInPersist::GetObject(AnyDataObj \& o)} \\
371Reads the next object in stream into {\bf o}. An exception is
372generated for incompatible object types.
373\item[] {\tt PInPersist::GetObject(AnyDataObj \& o, string tagname)} \\
374Reads the object associated with the tag {\bf tagname} into {\bf o}.
375An exception is generated for incompatible object types.
376\end{itemize}
377The operators {\tt operator << (POutPersist ...) } and
378{\tt operator >> (PInPersist ...) } are often overloaded
379to perform {\tt PutObject()} and {\tt GetObject()} operations.
380the {\bf PPFNameTag} (ppfnametag.h) class can be used in conjunction with
381{\tt << >> } operators to write objects with a name tag or to retrieve
382an object identified with a name tag. The example below shows the
383usage of these operators:
384\begin{verbatim}
385// Creating and filling a histogram
386Histo hw(0.,10.,100);
387...
388// Writing histogram to a PPF stream
389POutPersist os("hw.ppf");
390os << PPFNameTag("myhisto") << hw;
391
392// Reading a histogram from a PPF stream
393PInPersist is("hr.ppf");
394is >> PPFNameTag("myhisto") >> hr;
395\end{verbatim}
396
397The {\bf scanppf} program can be used to list the content of a PPF file.
398
399\subsection{\tcls{NDataBlock}}
400\index{\tcls{NDataBlock}}
401\begin{figure}[hbt]
402\dclsbb{AnyDataObj}{\tcls{NDataBlock}}
403\dclsbb{PPersist}{\tcls{FIO\_NDataBlock}}
404\end{figure}
405The {\bf \tcls{NDataBlock}} is designed to handle reference counting
406and sharing of memory blocs (contiguous arrays) for numerical data
407types. Initialisation, resizing, basic arithmetic operations, as
408well as persistence handling services are provided.
409The persistence handler class ({\tt \tcls{FIO\_NDataBlock}}) insures
410that a single copy of data is written for multiply referenced objects,
411and the data is shared among objects when reading.
412\par
413The example below shows writing of NDataBlock objects through the
414use of overloaded operator $ << $ :
415\begin{verbatim}
416#include "fiondblock.h"
417// ...
418POutPersist pos("aa.ppf");
419NDataBlock<r_4> rdb(40);
420rdb = 567.89;
421pos << rdb;
422// We can also use the PutObject method
423NDataBlock<int_4> idb(20);
424idb = 123;
425pos.PutObject(idb);
426\end{verbatim}
427The following sample programs show the reading of the created PPF file :
428\begin{verbatim}
429PInPersist pis("aa.ppf");
430NDataBlock<r_4> rdb;
431pis >> rdb;
432cout << rdb;
433NDataBlock<int_4> idb;
434cout << idb;
435\end{verbatim}
436
437\subsection{Using DVList}
438\index{DVList} \index{MuTyV}
439\begin{figure}[hbt]
440\dclsbb{AnyDataObj}{DVList}
441\dclsbb{PPersist}{\tclsc{ObjFileIO}{DVList}}
442\end{figure}
443The {\bf DVList} class objects can be used to create and manage list
444of values, associated with names. A list of pairs of (MuTyV, name(string))
445is maintained by DVList objects. {\bf MuTyV} is a simple class
446capable of holding string, integer, float or complex values,
447providing easy conversion methods between these objects.
448\begin{verbatim}
449// Using MuTyV objects
450MuTyV s("hello"); // string type value
451MuTyV x;
452x = "3.14159626"; // string type value, ASCII representation for Pi
453double d = x; // x converted to double = 3.141596
454x = 314; // x contains the integer value = 314
455// Using DVList
456DVList dvl;
457dvl("Pi") = 3.14159626; // float value, named Pi
458dvl("Log2") = 0.30102999; // float value, named Log2
459dvl("FileName") = "myfile.fits"; // string value, named myfile.fits
460// Printing DVList object
461cout << dvl;
462\end{verbatim}
463
464\subsection{\tcls{SegDataBlock} , \tcls{SwSegDataBlock}}
465\begin{figure}[hbt]
466\dclsccc{AnyDataObj}{\tcls{SegDBInterface}}{ \tcls{SegDataBlock} }
467\dclscc{\tcls{SegDBInterface}}{ \tcls{SwSegDataBlock} }
468\end{figure}
469\begin{itemize}
470\item[] \tcls{SegDataBlock} handles arrays of object of
471type {\bf T} with reference sharing in memory. The array can be extended
472(increase in array size) with fixed segment size. It implements the interface
473defined by tcls{SegDBInterface}.
474\item[] \tcls{SwSegDataBlock} Implements the same tcls{SegDBInterface}
475using a data swapper (\tcls{DataSwapperInterface} .
476Can thus be used for very large objects.
477\end{itemize}
478
479\newpage
480\section{Module TArray}
481\index{\tcls{TArray}}
482{\bf TArray} module contains template classes for handling standard
483operations on numerical arrays. Using the class {\tt \tcls{TArray} },
484it is possible to create and manipulate up to 5-dimension numerical
485arrays {\tt (int, float, double, complex, \ldots)}. The include
486file {\tt array.h} declares all the classes and definitions
487in module TArray. {\bf Array} is a typedef for arrays
488with double precision floating value elements. \\
489{\tt typedef TArray$<$r\_8$>$ Array ; }
490
491\begin{figure}[hbt]
492\dclsccc{AnyDataObj}{BaseArray}{\tcls{TArray}}
493\dclsbb{PPersist}{\tcls{FIO\_TArray}}
494\end{figure}
495
496
497\subsection{Using arrays}
498\index{Sequence} \index{RandomSequence} \index{RegularSequence}
499\index{EnumeratedSequence}
500The example below shows basic usage of arrays, creation, initialisation
501and arithmetic operations. Different kind of {\bf Sequence} objects
502can be used for initialising arrays.
503
504\begin{figure}[hbt]
505\dclsbb{Sequence}{RandomSequence}
506\dclsb{RegularSequence}
507\dclsb{EnumeratedSequence}
508\end{figure}
509
510The example below shows basic usage of arrays:
511\index{\tcls{TArray}}
512\begin{verbatim}
513// Creating and initialising a 1-D array of integers
514TArray<int> ia(5);
515EnumeratedSequence es;
516es = 24, 35, 46, 57, 68;
517ia = es;
518cout << "Array<int> ia = " << ia;
519// 2-D array of floats
520TArray<r_4> b(6,4), c(6,4);
521// Initializing b with a constant
522b = 2.71828;
523// Filling c with random numbers
524c = RandomSequence();
525// Arithmetic operations
526TArray<r_4> d = b+0.3f*c;
527cout << "Array<float> d = " << d;
528\end{verbatim}
529
530The copy constructor shares the array data, while the assignment operator
531copies the array elements, as illustrated in the following example:
532\begin{verbatim}
533TArray<int> a1(4,3);
534a1 = RegularSequence(0,2);
535// Array a2 and a1 shares their data
536TArray<int> a2(a1);
537// a3 and a1 have the same size and identical elements
538TArray<int> a3;
539a3 = a1;
540// Changing one of the a2 elements
541a2(1,1,0) = 555;
542// a1(1,1) is also changed to 555, but not a3(1,1)
543cout << "Array<int> a1 = " << a1;
544cout << "Array<int> a3 = " << a3;
545\end{verbatim}
546
547\subsection{Matrices and vectors}
548\index{\tcls{TMatrix}} \index{\tcls{TVector}}
549\begin{figure}[hbt]
550\dclsccc{\tcls{TArray}}{\tcls{TMatrix}}{\tcls{TVector}}
551\end{figure}
552Vectors and matrices are 2 dimensional arrays. The array size
553along one dimension is equal 1 for vectors. Column vectors
554have {\tt NCols() = 1} and row vectors have {\tt NRows() = 1}.
555Mathematical expressions involving matrices and vectors can easily
556be translated into C++ code using {\tt TMatrix} and
557{\tt TVector} objects. {\bf Matrix} and {\bf Vector} are
558typedefs for double precision float matrices and vectors.
559The operator {\bf *} beteween matrices is redefined to
560perform matrix multiplication. One can then write: \\
561\begin{verbatim}
562 // We create a row vector
563 Vector v(1000, BaseArray::RowVector);
564 // Initialize values with a random sequence
565 v = RandomSequence();
566 // Compute the vector length (norm)
567 double norm = (v*v.Transpose()).toScalar();
568 cout << "Norm(v) = " << norm << endl;
569\end{verbatim}
570
571This module contains basic array and matrix operations
572such as the Gauss matrix inversion algorithm
573which can be used to solve linear systems, as illustrated by the
574example below:
575\begin{verbatim}
576#include "sopemtx.h"
577// ...
578// Creation of a random 5x5 matrix
579Matrix A(5,5);
580A = RandomSequence(RandomSequence::Flat);
581Vector X0(5);
582X0 = RandomSequence(RandomSequence::Gaussian);
583// Computing B = A*X0
584Vector B = A*X0;
585// Solving the system A*X = B
586Vector X;
587LinSolve(A, B, X);
588// Checking the result
589Vector diff = X-X0;
590cout << "X-X0= " << diff ;
591double min,max;
592diff.MinMax(min, max);
593cout << " Min(X-X0) = " << min << " Max(X-X0) = " << max << endl;
594\end{verbatim}
595
596\subsection{Working with sub-arrays and Ranges}
597\index{Range}
598A powerful mechanism is included in array classes for working with
599sub-arrays. The class {\bf Range} can be used to specify range of array
600indexes in any of the array dimensions. Any regularly spaced index
601range can be specified, using the {\tt start} and {\tt end} index
602and an optional step (or stride). It is also possible to specify
603the {\tt start} index and the number of elements:
604\begin{center}
605
606\begin{tabular}{ll}
607\multicolumn{2}{c}{ {\bf Range} {\tt (start=0, end=0, size=1, step=1) } } \\[2mm]
608\hline \\
609{\bf Range} {\tt r(3,6); } & index range 3,4,5,6 \\
610{\bf Range} {\tt r(7,0,3); } & index range 7,8,9 \\
611{\bf Range} {\tt r(10,0,3,5); } & index range 10,12,14,16,18 \\
612\end{tabular}
613\end{center}
614
615In the following example, a simple low-pass filter, on a one
616dimensional stream (Vector) has been written using sub-arrays:
617
618\begin{verbatim}
619// Input Vector containing a noisy periodic signal
620 Vector in(1024), out(1024);
621 in = RandomSequence(RandomSequence::Gaussian, 0., 1.);
622 for(int kk=0; kk<in.Size(); kk++)
623 in(kk) += 2*sin(kk*0.05);
624// Compute the output vector by a simple low pass filter
625 Vector out(1024);
626 int w = 2;
627 for(int k=w; k<in.Size()-w; k++)
628 out(k) = in(Range(k-w, k+w).Sum()/(2.*w+1.);
629\end{verbatim}
630
631\subsection{Input, Output}
632Arrays can easily be saved to, or restored from files in different formats.
633SOPHYA library can handle array I/O to ASCII formatted files, to PPF streams,
634as well as to files in FITS format.
635FITS format input/output is provided through the classes in
636{\bf FitsIOServer} module. Onnly arrays with data types
637supported by the FITS standard can be handled during
638I/O operations to and from FITS streams (See the FitsIOServer section
639for additional details).
640
641\subsubsection{PPF streams}
642
643SOPHYA persistence (PPF streams) handles reference sharing, and multiply
644referenced objects are only written once. A hierarchy of arrays and sub-arrays
645written to a PPF stream is thus completely recovered, when the stream is read.
646The following example illustrates this point:
647\begin{verbatim}
648{
649// Saving an array with a sub-array into a POutPersist file
650Matrix A(3,4);
651A = RegularSequence(10,5);
652// Create a sub-array of A
653Matrix AS = A(Range(1,2), Range(2,3));
654// Save the two arrays to a PPF stream
655POutPersist pos("aas.ppf");
656pos << A << AS;
657}
658{
659// Reading arrays from the previously created PPF file aas.ppf
660PInPersist pis("aas.ppf");
661Matrix B,BS;
662pis >> B >> BS;
663// BS is a sub-array of B, modifying BS changes also B
664BS(1,1) = 98765.;
665cout << " B , BS after BS(1,1) = 98765. "
666 << B << BS << endl;
667}
668\end{verbatim}
669The execution of this sample code creates the file {\tt aas.ppf} and
670its output is reproduced here. Notice that the array hierarch is
671recovered. BS is a sub-array of B, and modifying BS changes also
672the corresponding element in B.
673\begin{verbatim}
674 B , BS after BS(1,1) = 98765.
675
676--- TMatrix<double>(NRows=3, NCols=4) ND=2 SizeX*Y*...= 4x3 ---
67710 15 20 25
67830 35 40 45
67950 55 60 98765
680
681--- TMatrix<double>(NRows=2, NCols=2) ND=2 SizeX*Y*...= 2x2 ---
68240 45
68360 98765
684\end{verbatim}
685
686\centerline{\bf Warning: }
687
688There is a drawback in this behaviour: only a single
689copy of an array is written to a file, even if the array is modified,
690without being resized and written to a PPF stream.
691\begin{verbatim}
692{
693POutPersist pos("mca.ppf");
694TArray<int_4> ia(5,3);
695ia = 8;
696pos << ia;
697ia = 16;
698pos << ia;
699ia = 32;
700pos << ia;
701}
702\end{verbatim}
703
704Only a single copy of the data is effectively written to the output
705PPF file, corresponding to the value 8 for array elements. When we
706read the three array from the file mca.ppf, the same array elements
707are obtained three times (all elements equal to 8):
708\begin{verbatim}
709{
710PInPersist pis("mca.ppf");
711TArray<int_4> ib;
712pis >> ib;
713cout << " First array read from mca.ppf : " << ib;
714pis >> ib;
715cout << " Second array read from mca.ppf : " << ib;
716pis >> ib;
717cout << " Third array read from mca.ppf : " << ib;
718}
719\end{verbatim}
720
721\subsubsection{ASCII streams}
722
723The {\bf WriteASCII} method can be used to dump an array to an ASCII
724formatted file, while the {\bf ReadASCII} method can be used to decode
725ASCII formatted files. Space or tabs are the possible separators.
726Complex numbers should be specified as a pair of comma separated
727real and imaginary parts, enclosed in parenthesis.
728
729\begin{verbatim}
730{
731// Creating array A and writing it to an ASCII file (aaa.txt)
732Array A(4,6);
733A = RegularSequence(0.5, 0.2);
734ofstream ofs("aaa.txt");
735A.WriteASCII(ofs);
736}
737{
738// Decoding the ASCII file aaa.txt
739ifstream ifs("aaa.txt");
740Array B;
741sa_size_t nr, nc;
742B.ReadASCII(ifs,nr,nc);
743cout << " Array B; B.ReadASCII() from file " << endl;
744cout << B ;
745}
746\end{verbatim}
747
748
749\subsection{Complex arrays}
750The {\bf TArray} module provides few functions for manipulating
751arrays of complex numbers (single and double precision).
752These functions are declared in {\tt matharr.h}.
753\begin{itemize}
754\item[\bul] Creating a complex array through the specification of the
755real and imaginary parts.
756\item[\bul] Functions returning arrays corresponding to real and imaginary
757parts of a complex array: {\tt real(za) , imag(za) }
758({\bf Warning:} Note that the present implementation does not provide
759shared memory access to real and imaginary parts.)
760\item[\bul] Functions returning arrays corresponding to the module,
761phase, and module squared of a complex array:
762 {\tt phase(za) , module(za) , module2(za) }
763\end{itemize}
764
765\begin{verbatim}
766 TVector<r_4> p_real(10, BaseArray::RowVector);
767 TVector<r_4> p_imag(10, BaseArray::RowVector);
768 p_real = RegularSequence(0., 0.5);
769 p_imag = RegularSequence(0., 0.25);
770 TVector< complex<r_4> > zvec = ComplexArray(p_real, p_imag);
771 cout << " :: zvec= " << zvec;
772 cout << " :: real(zvec) = " << real(zvec) ;
773 cout << " :::: imag(zvec) = " << imag(zvec) ;
774 cout << " :::: module2(zvec) = " << module2(zvec) ;
775 cout << " :::: module(zvec) = " << module(zvec) ;
776 cout << " :::: phase(zvec) = " << phase(zvec) ;
777\end{verbatim}
778
779The decoding of complex numbers from an ASCII formatted stream
780is illustrated by the next example. As mentionned already,
781complex numbers should be specified as a pair of comma separated
782real and imaginary parts, enclosed in parenthesis.
783
784\begin{verbatim}
785csh> cat zzz.txt
786(1.,-1) (2., 2.5) -3. 12.
787-24. (-6.,7.) 14.2 (8.,64.)
788
789// Decoding of complex numbers from an ASCII file
790// Notice that the << operator can be used instead of ReadASCII
791TArray< complex<r_4> > Z;
792ifstream ifs("zzz.txt");
793ifs >> Z;
794cout << " TArray< complex<r_4> > Z from file zzz.txt " << Z ;
795\end{verbatim}
796
797
798\subsection{Memory organisation}
799{\tt \tcls{TArray} } can handle numerical arrays with various memory
800organisation, as long as the spacing (steps) along each axis is
801regular. The five axis are labeled X,Y,Z,T,U. The examples below
802illustrates the memory location for a 2-dimensional, $N_x=4 \times N_y=3$.
803The first index is along the X axis and the second index along the Y axis.
804\begin{verbatim}
805 | (0,0) (0,1) (0,2) (0,3) |
806 | (1,0) (1,1) (1,2) (1,3) |
807 | (2,0) (2,1) (2,2) (2,3) |
808\end{verbatim}
809In the first case, the array is completely packed
810($Step_X=1, Step_Y=N_X=4$), with zero offset,
811while in the second case, $Step_X=2, Step_Y=10, Offset=10$:
812\begin{verbatim}
813 | 0 1 2 3 | | 10 12 14 16 |
814Ex1 | 4 5 6 7 | Ex2 | 20 22 24 26 |
815 | 8 9 10 11 | | 30 32 34 36 |
816\end{verbatim}
817
818For matrices and vectors, an optional argument ({\tt MemoryMapping})
819can be used to select the memory mapping, where two basic schemes
820are available: \\
821{\tt CMemoryMapping} and {\tt FortranMemoryMapping}. \\
822In the case where {\tt CMemoryMapping} is used, a given matrix line
823is packed in memory, while the columns are packed when
824{\tt FortranMemoryMapping} is used. The first index when addressing
825the matrix elements (line number index) runs along
826the Y-axis if {\tt CMemoryMapping} is used, and along the X-axis
827in the case of {\tt FortranMemoryMapping}.
828Arithmetic operations between matrices
829with different memory organisation is allowed as long as
830the two matrices have the same sizes (Number of rows and columns).
831The following code example and the corresponding output illustrates
832these two memory mappings. The {\tt \tcls{TMatrix}::TransposeSelf() }
833method changes effectively the matrix memory mapping, which is also
834the case of {\tt \tcls{TMatrix}::Transpose() } method without argument.
835
836\begin{verbatim}
837TArray<r_4> X(4,2);
838X = RegularSequence(1,1);
839cout << "Array X= " << X << endl;
840TMatrix<r_4> X_C(X, true, BaseArray::CMemoryMapping);
841cout << "Matrix X_C (CMemoryMapping) = " << X_C << endl;
842TMatrix<r_4> X_F(X, true, BaseArray::FortranMemoryMapping);
843cout << "Matrix X_F (FortranMemoryMapping) = " << X_F << endl;
844\end{verbatim}
845This code would produce the following output (X\_F = Transpose(X\_C)) :
846\begin{verbatim}
847Array X=
848--- TArray<f> ND=2 SizeX*Y*...= 4x2 ---
8491, 2, 3, 4
8505, 6, 7, 8
851
852Matrix X_C (CMemoryMapping) =
853--- TMatrix<f>(NRows=2, NCols=4) ND=2 SizeX*Y*...= 4x2 ---
8541, 2, 3, 4
8555, 6, 7, 8
856
857Matrix X_F (FortranMemoryMapping) =
858--- TMatrix<f>(NRows=4, NCols=2) ND=2 SizeX*Y*...= 4x2 ---
8591, 5
8602, 6
8613, 7
8624, 8
863\end{verbatim}
864
865\newpage
866
867\section{Module HiStats}
868\begin{figure}[hbt]
869\dclsccc{AnyDataObj}{Histo}{HProf}
870\dclsbb{AnyDataObj}{Histo2D}
871\dclsbb{AnyDataObj}{Ntuple}
872\caption{partial class diagram for histograms and ntuples}
873\end{figure}
874
875{\bf HiStats} contains classes for creating, filling, printing and
876doing various operations on one or two dimensional histograms
877{\tt Histo} and {\tt Histo2D} as well as profile histograms {\tt HProf}. \\
878This module also contains {\tt NTuple} and {\tt XNTuple} which are
879more or less the same that the binary FITS tables.
880
881\subsection{1D Histograms}
882\index{Histo}
883For 1D histograms, various numerical methods are provided such as
884computing means and sigmas, finding maxima, fitting, rebinning,
885integrating \dots \\
886The example below shows creating and filling a one dimensional histogram
887of 100 bins from $-5.$ to $+5.$ to create a Gaussian normal distribution
888with errors~:
889\begin{verbatim}
890#include "histos.h"
891// ...
892Histo H(-0.5,0.5,100);
893H.Errors();
894for(int i=0;i<25000;i++) {
895 double x = NorRand();
896 H.Add(x);
897}
898H.Print(80);
899\end{verbatim}
900
901\subsection{2D Histograms}
902\index{Histo2D}
903Much of these operations are also valid for 2D histograms. 1D projection
904or slices can be set~:
905\begin{verbatim}
906#include "histos2.h"
907// ...
908Histo2D H2(-1.,1.,100,0.,60.,50);
909H2.SetProjX(); // create the 1D histo for X projection
910H2.SetBandX(25.,35.); // create 1D histo projection for 25.<y<35.
911H2.SetBandX(35.,45.); // create 1D histo projection for 35.<y<45.
912H2.SetBandX(40.,55.); // create 1D histo projection for 40.<y<55.
913//... fill H2 with what ever you want
914H2.Print();
915Histo *hx = H2.HProjX();
916 hx->Print(80);
917Histo *hbx2 = HBandX(1); // Get the second X band (35.<y<45.)
918 hbx2->Print(80);
919\end{verbatim}
920
921\subsection{Profile Histograms}
922\index{HProf}
923Profiles histograms {\bf HProf} contains the mean and the
924sigma of the distribution
925of the values filled in each bin. The sigma can be changed to
926the error on the mean. When filled, the profile histogram looks
927like a 1D histogram and much of the operations that can be done on 1D histo
928may be applied onto profile histograms.
929
930\subsection{Data tables (tuples)}
931\index{NTuple}
932NTuple are memory resident tables of 32 or 64 bits floating values
933(float/double).They are arranged in columns. Each line is often called an event.
934These objects are frequently used to analyze data.
935The piapp graphicals tools can plot a column against an other one
936with respect to various selection cuts. \\
937Here is an example of creation and filling~:
938\begin{verbatim}
939#include "ntuple.h"
940#include "srandgen.h"
941// ...
942char* nament[4] = {"i","x","y","ey"};
943r_4 xnt[4];
944NTuple NT(4,nament);
945for(i=0;i<5000;i++) {
946 xnt[0] = i+1;
947 xnt[1] = 5.*drandpm1(); // a random value between -5 and +5
948 xnt[2] = 100.*exp(-0.5*xnt[1]*xnt[1]) + 1.;
949 xnt[3] = sqrt(xnt[2]);
950 xnt[2] += xnt[3] * NorRand(); // add a random gaussian error
951 NT.Fill(xnt);
952}
953\end{verbatim}
954
955XNTuple provide additional functionalities, compared to NTuple.
956They are deprecated and are only kept for backward compatibility
957and should not be used anymore. Use DataTable and
958SwPPFDataTable instead.
959Object of type XNTuple handle various types
960of column values (double,float,int,string,...) and can handle
961very large data sets, through swap space on disk.
962
963\index{DataTable}
964
965The class {\bf DataTable} extends significantly the functionalities provided by
966NTuple. DataTable is a memory resident implementation of the interface
967{\bf BaseDataTable } which organizes the data as a 2-D table. User can define
968the name and data type of each column. Data is added to the table as rows.
969The table is extended as necessary when adding rows.
970The sample code below shows an example of DataTable usage :
971\begin{verbatim}
972 #include "datatable.h"
973 // ...
974 DataTable dt(64);
975 dt.AddFloatColumn("X0_f");
976 dt.AddFloatColumn("X1_f");
977 dt.AddDoubleColumn("X0X0pX1X1_d");
978 double x[5];
979 for(int i=0; i<63; i++) {
980 x[0] = (i/9)-4.; x[1] = (i/9)-3.; x[2] = x[0]*x[0]+x[1]*x[1];
981 dt.AddLine(x);
982 }
983 // Printing table info
984 cout << dt ;
985 // Saving object into a PPF file
986 POutPersist po("dtable.ppf");
987 po << dt ;
988
989\end{verbatim}
990
991\begin{figure}[hbt]
992\dclsccc{AnyDataObj}{BaseDataTable}{DataTable}
993\dclscc{BaseDataTable}{SwPPFDataTable}
994\end{figure}
995
996\index{SwPPFDataTable}
997The class {\bf SwPPFDataTable} implements the BaseDataTable interface
998using segmented data blocks with swap on PPF streams. Very large data sets
999can be created and manipulated through tis class
1000
1001\subsection{Writing, viewing \dots }
1002
1003All these objects have been design to be written to or read from a persistent file.
1004The following example shows how to write the previously created objects
1005into such a file~:
1006\begin{verbatim}
1007//-- Writing
1008{
1009char *fileout = "myfile.ppf";
1010string tag;
1011POutPersist outppf(fileout);
1012tag = "H"; outppf.PutObject(H,tag);
1013tag = "H2"; outppf.PutObject(H2,tag);
1014tag = "NT"; outppf.PutObject(NT,tag);
1015} // closing ``}'' destroy ``outppf'' and automatically close the file !
1016\end{verbatim}
1017
1018Sophya graphical tools (spiapp) can automatically display and operate
1019all these objects.
1020
1021\newpage
1022\section{Module NTools}
1023
1024This module provides elementary numerical tools for numerical integration,
1025fitting, sorting and ODE solving. FFTs are also provided (Mayer,FFTPack).
1026
1027\subsection{Fitting}
1028\index{Fitting} \index{Minimisation}
1029Fitting is done with two classes {\tt GeneralFit} and {\tt GeneralFitData}
1030and is based on the Levenberg-Marquardt method.
1031\index{GeneralFit} \index{GeneralFitData}
1032GeneralFitData is a class which provide a description of the data
1033to be fitted. GeneralFit is the fitter class. Parametrized functions
1034can be given as classes which inherit {\tt GeneralFunction}
1035or as simple C functions. Classes of pre-defined functions are provided
1036(see files fct1dfit.h and fct2dfit.h). The user interface is very close
1037from that of the CERN {\tt Minuit} fitter.
1038Number of objects (Histo, HProf \dots ) are interfaced with GeneralFit
1039and can be easily fitted. \\
1040Here is a very simple example for fitting the previously created NTuple
1041with a Gaussian~:
1042\begin{verbatim}
1043#include "fct1dfit.h"
1044// ...
1045
1046// Read from ppf file
1047NTuple nt;
1048{
1049PInPersist pis("myfile.ppf");
1050string tag = "NT"; pis.GetObject(nt,tag);
1051}
1052
1053// Fill GeneralData
1054GeneralData mGdata(nt.NEntry());
1055for(int i=0; i<nt.NEntry(); i++)
1056 mGdata.AddData1(xnt[1],xnt[2],xnt[3]); // Fill x, y and error on y
1057mGData.PrintStatus();
1058
1059// Function for fitting : y = f(x) + noise
1060Gauss1DPol mFunction; // gaussian + constant
1061
1062// Prepare for fit
1063GeneralFit mFit(&mFunction); // create a fitter for the choosen function
1064mFit.SetData(&mGData); // connect data to the fitter
1065
1066// Set and initialise the parameters (that's non-linear fitting!)
1067// (num par, name, guess start, step, [limits min and max])
1068mFit.SetParam(0,"high",90.,1..);
1069mFit.SetParam(1,"xcenter",0.05,0.01);
1070mFit.SetParam(2,"sigma",sig,0.05,0.01,10.);
1071 // Give limits to avoid division by zero
1072mFit.SetParam(3,"constant",0.,1.);
1073
1074// Fit and print result
1075int rcfit = mFit.Fit();
1076mFit.PrintFit();
1077if(rcfit>0) {)
1078 cout<<"Reduce_Chisquare = "<<mFit.GetChi2Red()
1079 <<" nstep="<<mFit.GetNStep()<<" rc="<<rcfit<<endl;
1080} else {
1081 cout<<"Fit_Error, rc = "<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
1082 mFit.PrintFitErr(rcfit);
1083}
1084
1085// Get the result for further use
1086TVector<r_8> ParResult = mFit.GetParm();
1087cout<<ParResult;
1088\end{verbatim}
1089
1090Much more usefull possibilities and detailed informations might be found
1091in the HTML pages of the Sophya manual.
1092
1093\subsection{Polynomial}
1094\index{Polynomial} \index{Poly} \index{Poly2}
1095Polynomials of 1 or 2 variables are supported ({\tt Poly} and {\tt Poly2}).
1096Various operations are supported~:
1097\begin{itemize}
1098\item elementary operations between polynomials $(+,-,*,/) $
1099\item setting or getting coefficients
1100\item computing the value of the polynomial for a given value
1101 of the variable(s),
1102\item derivating
1103\item computing roots (degre 1 or 2)
1104\item fitting the polynomial to vectors of data.
1105\end{itemize}
1106Here is an example of polynomial fitting~:
1107\begin{verbatim}
1108#include "poly.h"
1109// ...
1110Poly pol(2);
1111pol[0] = 100.; pol[1] = 0.; pol[2] = 0.01; // Setting coefficients
1112TVector<r_8> x(100);
1113TVector<r_8> y(100);
1114TVector<r_8> ey(100);
1115for(int i=0;i<100;i++) {
1116 x(i) = i;
1117 ey(i) = 10.;
1118 y(i) = pol((double) i) + ey(i)*NorRand();
1119 ey(i) *= ey(i)
1120}
1121
1122TVector<r_8> errcoef;
1123Poly polfit;
1124polfit.Fit(x,y,ey,2,errcoef);
1125
1126cout<<"Fit Result"<<polfit<<endl;
1127cout<<"Errors :"<<errcoef;
1128\end{verbatim}
1129
1130Similar operations can be done on polynomials with 2 variables.
1131
1132\subsection{Integration, Differential equations}
1133\index{Integration}
1134The NTools module provide also simple classes for numerical integration
1135of functions and differential equations.
1136\begin{figure}[hbt]
1137\dclsbb{Integrator}{GLInteg}
1138\dclsb{TrpzInteg}
1139\end{figure}
1140
1141\index{GLInteg} \index{TrpzInteg}
1142{\bf GLInteg} implements the integration through Gauss-Legendre method
1143and {\bf TrpzInteg} implements trapeze integration. For {\bf TrpzInteg},
1144number of steps specify the number of trapeze, and integration step,
1145their width.
1146The sample code below illustrates the use of TrpzInteg class:
1147\begin{verbatim}
1148#include "integ.h"
1149// ......................................................
1150// Function to be integrated
1151double myf(double x)
1152{
1153// Simple a x + b x^2 (a=2 b=3)
1154return (x*(2.+3.*x));
1155}
1156// ......................................................
1157
1158// Compute Integral(myf, 2., 5.) between xmin=2., xmax=5.
1159TrpzInteg trpz(myf, 2., 5.);
1160// We specify an integration step
1161trpz.DX(0.01);
1162// The integral can be computed as trpz.Value()
1163double myf_integral = trpz.Value();
1164// We could have used the cast operator :
1165cout << "Integral[myf, 2., 5.]= " << (double)trpz << endl;
1166// Limits can be specified through ValueBetween() method
1167cout << "Integral[myf, 0., 4.]= " << trpz.ValueBetween(0.,4.) << endl;
1168\end{verbatim}
1169
1170\subsection{Fourier transform (FFT)}
1171\index{FFT} \index{FFTPackServer}
1172An abstract interface for performing FFT operations is defined by the
1173{\bf FFTServerInterface} class. The {\bf FFTPackSever} class implements
1174one dimensional FFT, on real and complex data. FFTPackServer uses an
1175adapted and extended version of FFTPack (available from netlib),
1176translated in C, and can operate on single and double precision
1177({\tt float, double}) data.
1178
1179The sample code below illustrates the use of FFTServers:
1180\begin{verbatim}
1181#include "fftpserver.h"
1182 // ...
1183TVector<r_8> in(32);
1184TVector< complex<r_8> > out;
1185in = RandomSequence();
1186FFTPackServer ffts;
1187ffts.setNormalize(true); // To have normalized transforms
1188cout << " FFTServer info string= " << ffts.getInfo() << endl;
1189cout << "in= " << in << endl;
1190cout << " Calling ffts.FFTForward(in, out) : " << endl;
1191ffts.FFTForward(in, out);
1192cout << "out= " << out << endl;
1193\end{verbatim}
1194
1195% \newpage
1196\section{Module SUtils}
1197Some utility classes and C/C++ string manipulation functions are gathered
1198in {\bf SUtils} module.
1199\subsection{Using DataCards}
1200\index{DataCards}
1201The {\bf DataCards} class can be used to read parameters from a file.
1202Each line in the file starting with \@ defines a set of values
1203associated with a keyword. In the example below, we read the
1204parameters corresponding with the keyword {\tt SIZE} from the
1205file {\tt ex.d}. We suppose that {\tt ex.d} contains the line: \\
1206{\tt @SIZE 400 250} \\
1207\begin{verbatim}
1208#include "datacards.h"
1209// ...
1210// Initialising DataCards object dc from file ex.d
1211DataCards dc( "ex.d" );
1212// Getting the first and second parameters for keyword size
1213// We define a default value 100
1214int size_x = dc.IParam("SIZE", 0, 100);
1215int size_y = dc.IParam("SIZE", 1, 100);
1216cout << " size_x= " << size_x << " size_y= " << size_y << endl;
1217\end{verbatim}
1218
1219\section{Module SysTools}
1220The {\bf SysTools} module contains classes implementing interface to some
1221OS specific services. The class {\bf ResourceUsage} \index{ResourceUsage}
1222and {\bf Timer} {\index{Timer} provides access to information
1223about various resource usage (memory, CPU, ...).
1224The class {\bf Periodic} provides the necessary services needed to
1225implement the execution of a periodic action.
1226
1227\subsection{Thread management classes}
1228A basic interface to POSIX threads \index{thread} is also provided
1229through the \index{ZThread} {\bf ZThread}, {\bf ZMutex} and {\bf ZSync}
1230classes.
1231
1232\subsection{Dynamic linker and C++ compiler classes}
1233\index{PDynLinkMgr}
1234The class {\bf PDynLinkMgr} can be used for managing shared libraries
1235at run time. The example below shows the run time linking of a function:\\
1236{\tt extern "C" { void myfunc(); } } \\
1237\begin{verbatim}
1238#include "pdlmgr.h"
1239// ...
1240string soname = "mylib.so";
1241string funcname = "myfunc";
1242PDynLinkMgr dyl(soname);
1243DlFunction f = dyl.GetFunction(funcname);
1244if (f != NULL) {
1245// Calling the function
1246 f();
1247}
1248\end{verbatim}
1249
1250\index{CxxCompilerLinker}
1251The {\bf CxxCompilerLinker} class provides the services to compile C++ code and building
1252shared libraries, using the same compiler and options which have
1253been used to create the SOPHYA shared library.
1254The sample program below illustrates using this class to build
1255the shared library (myfunc.so) from the source file myfunc.cc :
1256\begin{verbatim}
1257#include "cxxcmplnk.h"
1258// ...
1259string flnm = "myfunc.cc";
1260string oname, soname;
1261int rc;
1262CxxCompilerLinker cxx;
1263// The Compile method provides a default object file name
1264rc = cxx.Compile(flnm, oname);
1265if (rc != 0 ) { // Error when compiling ... }
1266// The BuildSO method provides a default shared object file name
1267rc = cxx.BuildSO(oname, soname);
1268if (rc != 0 ) { // Error when creating shared object ... }
1269\end{verbatim}
1270
1271\subsection{Command interpreter}
1272The class {\bf Commander} can be used in interactive programs to provide
1273c-shell like command interpreter and scripting capabilties.
1274Arithmetic expression evaluation is implemented through the {\bf CExpressionEvaluator}
1275and {\bf RPNExpressionEvaluator} classes.
1276
1277\newpage
1278\section{Module SkyMap}
1279\begin{figure}[hbt]
1280\dclsbb{AnyDataObj}{PixelMap}
1281\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
1282\dclsc{SphereThetaPhi}
1283\dclsb{LocalMap}
1284\caption{partial class diagram for pixelization classes in Sophya}
1285\end{figure}
1286The {\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.
1287\subsection {Spherical maps}
1288There are two kinds of spherical maps according pixelization algorithms. SphereHEALPix represents spheres pixelized following the HEALPIix algorithm (E. Hivon, K. Gorski)
1289\footnote{see the HEALPix Homepage: http://www.eso.org/kgorski/healpix/ }
1290, SphereThetaPhi represents spheres pixelized following an algorithm developed at LAL-ORSAY. The example below shows creating and filling of a SphereHEALPix with nside = 8 (it will be 12*8*8= 768 pixels) :
1291\index{\tcls{SphereHEALPix}}
1292\index{\tcls{SphereThetaPhi}}
1293
1294\begin{verbatim}
1295#include "spherehealpix.h"
1296// ...
1297SphereHEALPix<double> sph(8);
1298for (int k=0; k< sph.NbPixels(); k++) sph(k) = (double)(10*k);
1299\end{verbatim}
1300
1301SphereThetaPhi is used in a similar way with an argument representing number of slices in theta (Euler angle) for an hemisphere.
1302\index{\tcls{SphereThetaPhi}}
1303
1304\subsection {Local maps}
1305\index{\tcls{LocalMap}}
1306A 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).
1307
1308Internally, 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(...))
1309
1310The 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).
1311\begin{verbatim}
1312#include "localmap.h"
1313//..............
1314 LocalMap<r_4> locmap(4,5);
1315 for (int k=0; k<locmap.NbPixels();k++) locmap(k)=10.*k;
1316 locmap.SetOrigin();
1317 locmap.SetSize(30.,30.);
1318\end{verbatim}
1319
1320\subsection{Writing, viewing \dots }
1321
1322All these objects have been design to be written to or read from a persistant file.
1323The following example shows how to write the previously created objects
1324into such a file~:
1325\begin{verbatim}
1326//-- Writing
1327
1328#include "fiospherehealpix.h"
1329//................
1330
1331char *fileout = "myfile.ppf";
1332POutPersist outppf(fileout);
1333FIO_SphereHEALPix<r_8> outsph(sph);
1334outsph.Write(outppf);
1335FIO_LocalMap<r_8> outloc(locmap);
1336outloc.Write(outppf);
1337// It is also possible to use the << operator
1338POutPersist os("sph.ppf");
1339os << outsph;
1340os << outloc;
1341\end{verbatim}
1342
1343Sophya graphical tools (spiapp) can automatically display and operate
1344all these objects.
1345
1346\newpage
1347\section{Module Samba}
1348\index{Spherical Harmonics}
1349\index{SphericalTransformServer}
1350The 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...
1351\begin{verbatim}
1352#include "skymap.h"
1353#include "samba.h"
1354....................
1355
1356// Generate input spectra a + b* l + c * gaussienne(l, 50, 20)
1357int lmax = 92;
1358Vector clin(lmax);
1359for(int l=0; l<lmax; l++) {
1360 double xx = (l-50.)/10.;
1361 clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
1362}
1363
1364// Compute map from spectra
1365SphericalTransformServer<r_8> ylmserver;
1366int m = 128; // HealPix pixelisation parameter
1367SphereHEALPix<r_8> map(m);
1368ylmserver.GenerateFromCl(map, m, clin, 0.);
1369// Compute power spectrum from map
1370Vector clout = ylmserver.DecomposeToCl(map, lmax, 0.);
1371\end{verbatim}
1372
1373\newpage
1374\section{Module SkyT}
1375\index{RadSpectra} \index{SpectralResponse}
1376The SkyT module is composed of two types of classes:
1377\begin{itemize}
1378\item{} one which corresponds to an emission spectrum of
1379radiation, which is called RadSpectra
1380\item{} one which corresponds to the spectral response
1381of a given detector (i.e. corresponding to a detector
1382filter in a given frequency domain), which is called
1383SpectralResponse.
1384\end{itemize}
1385\begin{figure}[hbt]
1386\dclsbb{RadSpectra}{RadSpectraVec}
1387\dclsb{BlackBody}
1388\dclsccc{AnyDataObj}{SpectralResponse}{SpecRespVec}
1389\dclsc{GaussianFilter}
1390\caption{partial class for SkyT module}
1391\end{figure}
1392
1393\begin{verbatim}
1394#include "skyt.h"
1395// ....
1396// Compute the flux from a blackbody at 2.73 K through a square filter
1397BlackBody myBB(2.73);
1398// We define a square filter from 100 - 200 GHz
1399SquareFilter mySF(100,200);
1400// Compute the filtered integrated flux :
1401double flux = myBB.filteredIntegratedFlux(mySF);
1402\end{verbatim}
1403
1404A more detailed description of SkyT module can be found in:
1405{\it The SkyMixer (SkyT and PMixer modules) - Sophya Note No 2. }
1406available also from Sophya Web site.
1407
1408\newpage
1409\section{Module FitsIOServer}
1410\begin{figure}[hbt]
1411\dclsbb{FitsFile}{FitsInFile}
1412\dclsb{FitsOutFile}
1413\end{figure}
1414\index{FITS} \index{FitsInFile} \index{FitsOutFile}
1415This module provides classes for handling file input-output in FITS format using the cfitsio library. It works like the SOPHYA persistence (see Module SysTools), using delegate objects, but its design is simpler. The following example writes a matrix (see module TArray) and a spherical map (see module SkyMap) on a FITS file and reads back from FITS file and creates new objects :
1416\begin{verbatim}
1417#include "spherehealpix.h"
1418#include "fitsspherehealpix.h"
1419#include "fitstarray.h"
1420#include "tmatrix.h"
1421//...........................
1422
1423int m=...;
1424SphereHEALPix<r_8> sph(m);
1425................
1426int dim1=...;
1427int dim2=...;
1428TMatrix<r_8> mat(dim1,dim2);
1429............
1430
1431FITS_SphereHEALPix<r_8> sph_temp(sph);
1432FITS_TArray<r_8> mat_temp(mat);
1433// writing
1434
1435FitsOutFile os("myfile.fits");
1436sph_temp.Write(os);
1437mat_temp.Write(os);
1438
1439// reading
1440FitsInFile is("myfile.fits");
1441sph_temp.Read(is);
1442mat_temp.Read(is);
1443SphereHEALPix<r_8> new_sph=(SphereHEALPix<r_8>)sph_temp;
1444TMatrix<r_8> new_mat=(TMatrix<r_8>)mat_temp;
1445................
1446
1447\end{verbatim}
1448
1449The operators {\tt operator << (FitsOutFile ...)} and
1450{\tt operator >> (FitsInFile ...)} are defined in order
1451to facilitate the FITS file operations:
1452\begin{verbatim}
1453// Writing an array object to a FITS file
1454#include "fitstarray.h"
1455FitsOutFile fio("arr.fits");
1456Matrix m(20,30);
1457m = 12345.;
1458fio << m;
1459// .....
1460// Reading a binary table to a XNTuple
1461#include "fitsxntuple.h"
1462XNTuple xn;
1463FitsInFile fii("table.fits");
1464fii >> xn;
1465\end{verbatim}
1466
1467The class {\bf FITS\_AutoReader} provides a limited FITS files reading
1468and decoding capabilities. A partial class diagram of FITS persistence
1469handling classes is shown below:
1470\begin{figure}[hbt]
1471\dclsbb{FitsIOhandler}{FITS\_TArray}
1472\dclsb{FITS\_NTuple}
1473% \dclsb{FITS\_XNTuple}
1474\dclsb{FITS\_SphereHEALPix}
1475% \dclsb{FITS\_LocalMap}
1476\end{figure}
1477
1478\newpage
1479\section{LinAlg and IFFTW modules}
1480An interface to use LAPACK library (available from {\tt http://www.netlib.org})
1481is implemented by the {\bf LapackServer} class, in module LinAlg.
1482\index{LapackServer}.
1483The sample code below shows how to use SVD (Singular Value Decomposition)
1484through LapackServer:
1485\begin{verbatim}
1486#include "intflapack.h"
1487// ...
1488// Use FortranMemoryMapping as default
1489BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
1490// Create an fill the arrays A and its copy AA
1491int n = 20;
1492Matrix A(n , n), AA;
1493A = RandomSequence(RandomSequence::Gaussian, 0., 4.);
1494AA = A; // AA is a copy of A
1495// Compute the SVD decomposition
1496Vector S; // Vector of singular values
1497Matrix U, VT;
1498LapackServer<r_8> lpks;
1499lpks.SVD(AA, S, U, VT);
1500// We create a diagonal matrix using S
1501Matrix SM(n, n);
1502for(int k=0; k<n; k++) SM(k,k) = S(k);
1503// Check the result : A = U*SM*VT
1504Matrix diff = U*(SM*VT) - A;
1505double min, max;
1506diff.MinMax(min, max);
1507cout << " Min/Max difference Matrix (?=0) , Min= " << min
1508 << " Max= " << max << endl;
1509\end{verbatim}
1510
1511\index{FFTWServer}
1512The {\bf FFTWServer} class (in module FFTW) implements FFTServerInterface class
1513methods, for one dimensional and multi-dimensional Fourier
1514transforms on double precision data using the FFTW package
1515(available from {\tt http://www.fftw.org}).
1516
1517\newpage
1518\section{Building and installing Sophya}
1519\subsection{supported platforms}
1520Presently, the Sophya library has been tested with the following
1521compiler/platform pairs:
1522
1523\begin{center}
1524\begin{tabular}{|l|l|}
1525\hline
1526OS & compiler \\
1527\hline
1528HP/Compaq/DEC Tru64 ( OSF1) & cxx (6.1 , 6.3) \\
1529Linux (RH) & g++ (3.2) \\
1530Linux (SCL) & icc (8.1) (Intel compiler) \\
1531SGI IRIX64 & CC (7.3) \\
1532MacOSX/Darwin 10.3 & g++ 3.3 \\
1533\hline
1534\end{tabular}
1535\end{center}
1536
1537Some of the modules in the Sophya package uses external libraries. The
1538{\bf FitsIOServer} is the example of such a module, where the {\tt libcfitsio.a}
1539is used.
1540par
1541The object files from a given Sophya module are grouped in an archive library
1542with the module's name ({\tt libmodulename.a}). All Sophya modules
1543 are grouped in a single shared library ({\tt libsophya.so}), while the
1544modules with reference to external libraries are grouped in
1545({\tt libextsophya.so}). The {\bf PI} and {\bf PIext} modules are
1546grouped in ({\tt libPI.so}).
1547
1548\subsection{Installation}
1549
1550The build procedure has two main steps: \\
1551- The configure step (BuildMgr/configure) setup the directory structure and
1552the necessary configuration file. \\
1553- The make step compiles the different sources files, create the library and optionaly
1554builds all or some of the associated executables.
1555
1556\par
1557{\tt BuildMgr/configure } is a c-shell script with a number of arguments:
1558\begin{verbatim}
1559csh> ./configure -h
1560configure [-sbase SOPHYABASE] [-scxx SOPHYACXX] [-incln]
1561 [-minc mymake.inc]
1562 [-extp dir1 -extp dir2 ...] [-extip dir1 -extip dir2 ... ]
1563 [-extlp dir1 -extlp dir2 ... ]
1564 [-noextlib -noext fits -noext fftw -noext lapack ]
1565 [-noext astro -noext minuit]
1566\end{verbatim}
1567\begin{itemize}
1568\item[] -sbase : define SOPHYA installation base directory. \$SOPHYABASE is used
1569if not specified.
1570\item[] -scxx : selects the C++ compiler. \$SOPHYACXX s used
1571if not specified.
1572\item[] -incln : creates symbolic link for include files, instead of copying them.
1573\item[] -minc : give an explicit name for the file used to generate
1574\$SOPHYABASE/include/sophyamake.inc.
1575\item[] -extp : Adds the specied path to the search path of the external libraries
1576include files and archive library.
1577\item[] -extip : Adds the specied path to the search path of the external libraries
1578include files.
1579\item[] -extp : Adds the specied path to the search path of the external libraries
1580archive (libxxx.a).
1581\item[] -noextlib : Disable compiling of modules referencing external libraries.
1582\item[] -noext : Disable compiling of the specified module (with reference to external
1583library.
1584\end{itemize}
1585
1586In the example below, we assume that we want to install Sophya from a
1587released (tagged) version in the source directory {\tt \$SRC} in the
1588{\tt /usr/local/Sophya} directory, using {\tt g++}. We assume that
1589the external libraries can be found in {\tt /usr/local/ExtLibs/}.
1590We disable the compilation of the MinuitAdapt and XAstrPack packages.
1591
1592\vspace*{3mm}
1593\begin{verbatim}
1594# Create the top level directory
1595csh> mkdir /usr/local/Sophya/
1596csh> cd $SRC/BuildMgr
1597# Step 1.a : Run the configuration script
1598csh> ./configure -sbase /usr/local/Sophya -scxx g++ -extp /usr/local/ExtLibs/ \
1599-noext astro -noext minuit
1600# Step 1.b : Check the generated file $SOPHYABASE/include/
1601# Step 2.a: Compile the modules without external library reference
1602csh> make libs
1603# Step 2.b: Compile the modules WITH external library reference (optional)
1604csh> make extlibs
1605# Step 2.c: Build libsophya.so
1606csh> make slb
1607# Step 2.d: Build libextsophya.so (optional)
1608csh> make slbext
1609# Step 2.e: Compile the PI and PIext modules (optional)
1610csh> make PI
1611# Step 2.f: Build the corresponding shared library libPI.so (optional)
1612csh> make slbpi
1613\end{verbatim}
1614
1615To compile all modules and build the shared libraries, it is possible
1616to use:
1617\begin{verbatim}
1618# Step 2.a ... 2.f
1619csh> make all slball
1620\end{verbatim}
1621
1622At this step, all libraries should have been made. Programs using
1623Sophya libraries can now be built:
1624\begin{verbatim}
1625# To compile test programs
1626csh> cd ../PrgUtil
1627csh> make
1628# To build (s)piapp (libPI.so is needed)
1629csh> cd ../ProgPI
1630csh> make
1631csh> cd ..
1632\end{verbatim}
1633
1634\subsection{Mgr module}
1635This module contains scripts which can be used for generating the
1636makefiles for each module.
1637\begin{itemize}
1638\item {\bf Makefile} Top level Makefile for building the libraries.
1639\item {\bf Makefile.h} contains the definition of compilation flags for the
1640different compilers and systems. This file is used for building the
1641library and generating {\bf MakefileUser.h} (to be included in makefiles).
1642\item {\bf Makefile.slb} contains the rules for building shared libraries
1643for the different compilers and systems. (to be included in makefiles)
1644\item {\bf crerep\_sophya} c-shell script for creating the directory tree
1645under {\tt \$SOPHYABASEREP} and {\tt \$SOPHYADEVREP}
1646\item {\bf install\_sophya} c-shell script for installing the Sophya package.
1647Usually from {\tt \$SOPHYADEVREP} to {\tt \$SOPHYABASEREP}
1648\item {\bf mkmflien} c-shell script for making symbolic links or copying
1649include files to {\tt \$SOPHYADEVREP/Include} or {\tt \$SOPHYABASEREP/Include}
1650\item {\bf mkmf} c-shell script for generating module makefiles and the
1651top level makefile (named GNUmakefile)
1652\item {\bf mkmflib} c-shell script for generating each library module
1653makefile (named GNUmakefile)
1654\item {\bf mkmfprog} c-shell script for generating makefile for a module
1655containing the source for executable programs (named GNUmakefile).
1656\item {\bf mkmfPI} c-shell script for generating makefile for PI and PIext
1657modules (named GNUmakefile)
1658\item {\bf libdirs} List of Sophya modules without reference to external
1659libraries.
1660\item {\bf extlibdirs} List of Sophya modules with reference to external
1661libraries.
1662
1663\end{itemize}
1664
1665\newpage
1666\appendix
1667\section{SOPHYA Exceptions}
1668\index{Exception classes} \index{PThrowable} \index{PError} \index{PException}
1669SOPHYA library defines a set of exceptions which are used
1670for signalling error conditions. The figure below shows a partial
1671class diagram for exception classes in SOPHYA.
1672\begin{figure}[hbt]
1673\dclsbb{PThrowable}{PError}
1674\dclscc{PError}{AllocationError}
1675\dclscc{PError}{NullPtrError}
1676\dclscc{PError}{ForbiddenError}
1677\dclscc{PError}{AssertionFailedError}
1678\dclsbb{PThrowable}{PException}
1679\dclscc{PException}{IOExc}
1680\dclscc{PException}{SzMismatchError}
1681\dclscc{PException}{RangeCheckError}
1682\dclscc{PException}{ParmError}
1683\dclscc{PException}{TypeMismatchExc}
1684\dclscc{PException}{MathExc}
1685\dclscc{PException}{CaughtSignalExc}
1686\caption{partial class diagram for exception handling in Sophya}
1687\end{figure}
1688
1689For simple programs, it is a good practice to handle
1690the exceptions at least at high level, in the {\tt main()} function.
1691The example below shows the exception handling and the usage
1692of Sophya persistence.
1693
1694\input{ex1.inc}
1695
1696
1697\newpage
1698\addcontentsline{toc}{section}{Index}
1699\printindex
1700\end{document}
1701
Note: See TracBrowser for help on using the repository browser.