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

Last change on this file since 1414 was 1414, checked in by ansari, 25 years ago

MAJ documentation - Reza 22/2/2001

File size: 41.8 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\makeindex
20
21\begin{document}
22
23\begin{titlepage}
24% The title page - top of the page with the title of the paper
25\titrehp{Sophya \\ An overview }
26% Authors list
27\auteurs{
28R. Ansari & ansari@lal.in2p3.fr \\
29E. Aubourg & aubourg@hep.saclay.cea.fr \\
30G. Le Meur & lemeur@lal.in2p3.fr \\
31C. Magneville & cmv@hep.saclay.cea.fr \\
32S. Henrot-Versille & versille@in2p3.fr
33}
34% \auteursall
35% The title page - bottom of the page with the paper number
36\vspace{1cm}
37\begin{center}
38{\bf \Large Document being updated !}
39\end{center}
40\titrebp{1}
41\end{titlepage}
42
43\tableofcontents
44
45\newpage
46
47\section{Introduction}
48
49{\bf SOPHYA} ({\bf SO}ftware for {\bf PHY}sics {\bf A}nalysis)
50is a collection of C++ classes designed for numerical and
51physics analysis software development. Our goal is to provide
52easy to use, yet powerful classes which can be used by scientists.
53We have decided to use as much as possible available
54numerical analysis libraries, encapsulating them whenever
55possible.
56
57 The SOPHYA design and implementation has been carried out
58with the specific goal of providing the general framework for
59the Planck-HFI data processing software. However, most of the
60packages presented here have a more general scope than the CMB analysis
61and Planck mission problem.
62 The source directory tree
63\footnote{ CVS: cvsserver.lal.in2p3.fr:/exp/eros/CVSPlanck}
64is organised into a number of modules.
65
66\begin{itemize}
67\item[] {\bf Mgr/} Scripts for code management,
68makefile generation and software installation
69\item[] {\bf SysTools/} General architecture support classes such
70as {\tt PPersist, NDataBlock<T>}, and few utility classes
71({\tt DataCard, DVList} \ldots).
72\item[] {\bf TArray/} template numerical arrays, vectors and matrices \\
73({\tt PixelMap<T> SphericalMap<T>} \ldots)
74\item[] {\bf NTools/} Some standard numerical analysis tools
75(linear, and non linear parameter fitting, FFT, \ldots)
76\item[] {\bf HiStats/} Histogram-ming and data set handling classes \\
77({\tt Histo Histo2D NTuple XNTuple} \ldots)
78\end{itemize}
79
80The modules listed below are more tightly related to the
81CMB (Cosmic Microwave Background) data analysis problem:
82\begin{itemize}
83\item[] {\bf SkyMap/} Local and full sky maps, and few geometry
84handling utility classes. \\
85({\tt PixelMap<T>, LocalMap<T>, SphericalMap<T>, \ldots})
86\item[] {\bf SkyT/}
87classes for spectral emission and detector frequency response modelling \\
88({\tt SpectralResponse, RadSpectra, BlackBody} \ldots)
89\item[] {\bf Samba/} Spherical harmonic analysis.
90\end{itemize}
91
92The following modules contain the interface classes with
93external libraries:
94\begin{itemize}
95\item[] {\bf FitsIOServer/} Classes for handling file input-output
96in FITS format using the cfitsio library.
97\item[] {\bf LinAlg/} Interface with Lapack linear algebra package
98\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a)
99\end{itemize}
100
101Other modules:
102\begin{itemize}
103\item[] {\bf Tests/} Simple test programs
104\item[] {\bf PrgUtil/} Various utility programs (runcxx, scanppf, scanfits, \ldots)
105\item[] {\bf PMixer/} skymixer and related programs
106\item[] {\bf ProgPI/} interactive analysis tool - It should be noted that
107this module uses the SOPHYA class library and is based on {\bf PI}
108which is a C++ library defining a complete GUI program
109architecture. An additional module (PIext) define the interactive
110analysis program framework and the interfaces with the objects
111in SOPHYA. The {\bf PI/} \footnote{the PI package documentation
112is available from {\bf http://www.lal.in2p3.fr/recherche/eros/PeidaDoc/} }
113and {\bf PIext/} modules are not currently part
114of the SOPHYA CVS structure.
115\end{itemize}
116
117\newpage
118
119\section{Using Sophya}
120Basic usage of Sophya classes are described in in the following sections.
121Complete Sophya documentation can be found at our web site: \\
122{\bf http://hfi-l2.in2p3.fr}.
123
124\subsection{Environment variables}
125Two environment variables {\bf DPCBASEREP} and {\bf EROSCXX} are used
126to define the path where the Sophya libraries and executable are installed.
127{\bf DPCBASEREP} defines the base directory path and {\bf EROSCXX} the
128name of the C++ compiler. The complete path is built using {\bf DPCBASEREP},
129the operating system name (as obtained by the {\tt uname} command), and
130the compiler name. In the example below, we show the complete path
131for a {\tt Linux} system, using the GNU g++ compiler:
132
133\begin{itemize}
134\item \$DPCBASEREP/Include : Include (.h) files
135\item \$DPCBASEREP/Linux-g++/Libs : Path for the archive libraries (.a)
136\item \$DPCBASEREP/Linux-g++/ShLibs : Shared library path (.so)
137\item \$DPCBASEREP/Linux-g++/Exec : Executable file path
138\end{itemize}
139
140In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
141should contain the Sophya shared library path
142({\tt \$DPCBASEREP/Linux-g++/ShLibs } when using g++ compiler on Linux)
143
144For modules using external libraries, the {\bf EXTLIBDIR}
145environment variable should contain the path to these libraries
146and corresponding include files.
147C-FitsIO anf FFTW include files should be accessible through: \\
148{\tt \$EXTLIBDIR/Include/FitsIO } \\
149{\tt \$EXTLIBDIR/Include/FFTW } \\
150The corresponding libraries are expected to be found in: \\
151{\tt \$EXTLIBDIR/Linux-g++/Libs} \\
152
153\subsection{User makefiles}
154The file {\tt \$DPCBASEREP/Include/MakefileUser.h} defines the compilation
155flags and the list of Sophya libraries. It should be included in the
156user's makefile. The default compilation rules assumes that the object (.o)
157and executable files would be put in the following diretories: \\
158{\tt \$HOME/`uname`-\$EROSCXX/Objs} \\
159{\tt \$HOME/`uname`-\$EROSCXX/Exec}.
160In the case of a {\tt Linux} system and using {\tt g++} as the C++ compiler,
161these two directories would be translated to \\
162{\tt \$HOME/Linux-g++/Objs} and {\tt \$HOME/Linux-g++/Exec}.
163The GNU make program should be used.
164\par
165The file {\tt \$DPCBASEREP/Include/makefile\_auto} defines the rules to compile
166a given source program, and link it against the Sophya libraries to produce
167an executable. The example below shows the steps to compile a program named
168{\tt trivial.cc }.
169\begin{verbatim}
170csh> cp \$DPCBASEREP/Include/makefile_auto makefile
171csh> make trivial
172\end{verbatim}
173This command should compile the {\tt trivial.cc} file,
174and link it against the sophya libraries. The object and executable
175file names are: \\
176{\tt \$HOME/`uname`-\$EROSCXX/Objs/trivial.o} \\
177{\tt \$HOME/`uname`-\$EROSCXX/Exec/trivial}.
178\par
179The file {\tt \$DPCBASEREP/Include/makefile\_example} provides another
180example makefile.
181
182\subsection{the runcxx program}
183\index{runcxx}
184{\bf runcxx} is a simple program which can be used to compile, link
185and run simple C++ programs. It handles the creation of a
186complete program file, containing the basic set C++ include files,
187the necessary include files for SOPHYA SysTools, TArray, HiStats
188and NTools modules, and the main program with exception handling.
189Other Sophya modules can be included using the {\tt -import} flag.
190\begin{verbatim}
191csh> runcxx -h
192SOPHYA Version 0.9 Revision 97 (V_Oct2000) -- Nov 9 2000 16:20:52 cxx
193 runcxx : compiling and running of a piece of C++ code
194 Usage: runcxx [-compopt CompileOptions] [-linkopt LinkOptions]
195 [-tmpdir TmpDirectory] [-f C++CodeFileName]
196 [-inc includefile] [-inc includefile ...]
197 [-import modulename] [-import modulename ...]
198 [-uarg UserArg1 UserArg2 ...]
199 if no file name is specified, read from standard input
200 modulenames: SkyMap, Samba, SkyT, FitsIOServer, LinAlg, IFFTW
201\end{verbatim}
202Most examples in this manual can be tested using runcxx. The
203example below shows how to compile, link and run a sample
204code.
205\begin{verbatim}
206// File example.icc
207Matrix a(3,3);
208a = IdentityMatrix(1.);
209cout << a ;
210// Executing this sample code
211csh> runcxx -f example.icc
212\end{verbatim}
213
214
215\section{Copy constructor and assignment operator}
216In C++, objects can be copied by assignment or by initialization.
217Copying by initialization corresponds to creating an object and
218initializing its value through the copy constructor.
219The copy constructor has its first argument as a reference, or
220const reference to the object's class type. It can have
221more arguments, if default values are provided.
222Copying by assignment applies to an existing object and
223is performed through the assignment operator (=).
224The copy constructor implements this for identical type objects:
225\begin{verbatim}
226class MyObject {
227public:
228 MyObject(); // Default constructor
229 MyObject(MyObject const & a); // Copy constructor
230 MyObject & operator = (MyObject const & a) // Assignment operator
231}
232\end{verbatim}
233The copy constructors play an important role, as they are
234called when class objects are passed by value,
235returned by value, or thrown as an exception.
236\begin{verbatim}
237// A function declaration with an argument of type MyObject,
238// passed by value, and returning a MyObject
239MyObject f(MyObject x)
240{
241 MyObject r;
242 ...
243 return(r); // Copy constructor is called here
244}
245// Calling the function :
246MyObject a;
247f(a); // Copy constructor called for a
248\end{verbatim}
249It should be noted that the C++ syntax is ambiguous for the
250assignment operator. {\tt MyObject x; x=y; } and
251{\tt MyObject x=y;} have different meaning.
252\begin{verbatim}
253MyObject a; // default constructor call
254MyObject b(a); // copy constructor call
255MyObject bb = a; // identical to bb(a) : copy constructor call
256MyObject c; // default constructor call
257c = a; // assignment operator call
258\end{verbatim}
259
260As a general rule in SOPHYA, objects which implements
261reference sharing on their data members have a copy constructor
262which shares the data, while the assignment operator copies or
263duplicate the data.
264
265\section{Module SysTools}
266
267{\bf SysTools} contains utility classes such as {\tt DataCards} or
268{\tt DVlist}, an hierarchy of exception classes for Sophya, a template
269class {\tcls{NDataBlock}} for handling reference counting on numerical
270arrays, as well as classes providing the services for implementing simple
271serialization.
272\vspace*{5mm}
273
274\subsection{SOPHYA persistence}
275\index{PPersist} \index{PInPersist} \index{POutPersist}
276\begin{figure}[hbt]
277\dclsa{PPersist}
278\dclsbb{PIOPersist}{PInPersist}
279\dclsb{POutPersist}
280\caption{partial class diagram for classes handling persistence in Sophya}
281\end{figure}
282A simple persistence mechanism is defined in SOPHYA. Its main
283features are:
284\begin{itemize}
285\item[] Portable file format, containing the description of the data structures
286and object hierarchy. \\
287{\bf PPF} {\bf P}ortable {\bf P}ersistence file {\bf F}ormat.
288\item[] Handling of read/write for mutiply referenced objects.
289\item[] All write operations are carried using sequential access only. This
290holds also for read operations, unless positional tags are used.
291SOPHYA persistence services can thus be used to transfer objects
292through network links.
293\item[] The serialisation (reading/writing) for objects for a given class
294is implemented through a delegate object. The delegate class inherits
295from {\tt PPersist} class.
296\end{itemize}
297A complete description of SOPHYA persistence mechanism and guidelines
298for writing delegate classes for handling object persistence is beyond
299the scope of this document. The example in the next paragraph shows
300simple use of SOPHYA persistence.
301
302\subsection{\tcls{NDataBlock}}
303\index{\tcls{NDataBlock}}
304\begin{figure}[hbt]
305\dclsbb{AnyDataObj}{\tcls{NDataBlock}}
306\dclsbb{PPersist}{\tcls{FIO\_NDataBlock}}
307\end{figure}
308The {\bf \tcls{NDataBlock}} is designed to handle reference counting
309and sharing of memory blocs (contiguous arrays) for numerical data
310types. Initialisation, resizing, basic arithmetic operations, as
311well as persistence handling services are provided.
312The persistence handler class ({\tt \tcls{FIO\_NDataBlock}}) insures
313that a single copy of data is written for multiply referenced objects,
314and the data is shared among objects when reading.
315\par
316The example below shows writing of NDataBlock objects through the
317use of overloaded operator $ << $ :
318\begin{verbatim}
319#include "fiondblock.h"
320// ...
321POutPersist pos("aa.ppf");
322NDataBlock<r_4> rdb(40);
323rdb = 567.89;
324pos << rdb;
325// We can also use the PutObject method
326NDataBlock<int_4> idb(20);
327idb = 123;
328pos.PutObject(idb);
329\end{verbatim}
330The following sample programs show the reading of the created PPF file :
331\begin{verbatim}
332PInPersist pis("aa.ppf");
333NDataBlock<r_4> rdb;
334pis >> rdb;
335cout << rdb;
336NDataBlock<int_4> idb;
337cout << idb;
338\end{verbatim}
339
340\subsection{Using DVList}
341\index{DVList} \index{MuTyV}
342\begin{figure}[hbt]
343\dclsbb{AnyDataObj}{DVList}
344\dclsbb{PPersist}{\tclsc{ObjFileIO}{DVList}}
345\end{figure}
346The {\bf DVList} class objects can be used to create and manage list
347of values, associated with names. A list of pairs of (MuTyV, name(string))
348is maintained by DVList objects. {\bf MuTyV} is a simple class
349capable of holding string, integer, float or complex values,
350providing easy conversion methods between these objects.
351\begin{verbatim}
352// Using MuTyV objects
353MuTyV s("hello"); // string type value
354MuTyV x;
355x = "3.14159626"; // string type value, ascii representation for Pi
356double d = x; // x converted to double = 3.141596
357x = 314; // x contains the integer value = 314
358// Using DVList
359DVList dvl;
360dvl("Pi") = 3.14159626; // float value, named Pi
361dvl("Log2") = 0.30102999; // float value, named Log2
362dvl("FileName") = "myfile.fits"; // string value, named myfile.fits
363// Printing DVList object
364cout << dvl;
365\end{verbatim}
366
367\subsection{Using DataCards}
368\index{DataCards}
369The {\bf DataCards} class can be used to read parameters from a file.
370Each line in the file starting with \@ defines a set of values
371associated with a keyword. In the example below, we read the
372parameters corresponding with the keyword {\tt SIZE} from the
373file {\tt ex.d}. We suppose that {\tt ex.d} contains the line: \\
374{\tt @SIZE 400 250} \\
375\begin{verbatim}
376#include "datacards.h"
377// ...
378// Initializing DataCards object dc from file ex.d
379DataCards dc( "ex.d" );
380// Getting the first and second parameters for keyword size
381// We define a default value 100
382int size_x = dc.IParam("SIZE", 0, 100);
383int size_y = dc.IParam("SIZE", 1, 100);
384cout << " size_x= " << size_x << " size_y= " << size_y << endl;
385\end{verbatim}
386
387\subsection{Dynamic linker}
388\index{PDynLinkMgr}
389The class {\bf PDynLinkMgr} can be used for managing shared libraries
390at run time. The example below shows the run time linking of a function:\\
391{\tt extern "C" { void myfunc(); } } \\
392\begin{verbatim}
393#include "pdlmgr.h"
394// ...
395string soname = "mylib.so";
396string funcname = "myfunc";
397PDynLinkMgr dyl(soname);
398DlFunction f = dyl.GetFunction(funcname);
399if (f != NULL) {
400// Calling the function
401 f();
402}
403\end{verbatim}
404
405\subsection{CxxCompilerLinker class}
406\index{CxxCompilerLinker}
407This class provides the services to compile C++ code and building
408shared libraries, using the same compiler and options which have
409been used to create the SOPHYA shared library.
410The sample program below illustrates using this class to build
411the shared library (myfunc.so) from the source file myfunc.cc :
412\begin{verbatim}
413#include "cxxcmplnk.h"
414// ...
415string flnm = "myfunc.cc";
416string oname, soname;
417int rc;
418CxxCompilerLinker cxx;
419// The Compile method provides a default object file name
420rc = cxx.Compile(flnm, oname);
421if (rc != 0 ) { // Error when compiling ... }
422// The BuildSO method provides a default shared object file name
423rc = cxx.BuildSO(oname, soname);
424if (rc != 0 ) { // Error when creating shared object ... }
425\end{verbatim}
426
427\section{Module TArray}
428\index{\tcls{TArray}}
429{\bf TArray} module contains template classes for handling standard
430operations on numerical arrays. Using the class {\tt \tcls{TArray} },
431it is possible to create and manipulate up to 5-dimension numerical
432arrays {\tt (int, float, double, complex, \ldots)}. The include
433file {\tt array.h} declares all the classes and definitions
434in module TArray. {\bf Array} is a typdef for arrays
435with double precision floating value elements. \\
436{\tt typedef TArray$<$r\_8$>$ Array ; }
437
438\begin{figure}[hbt]
439\dclsccc{AnyDataObj}{BaseArray}{\tcls{TArray}}
440\dclsbb{PPersist}{\tcls{FIO\_TArray}}
441\end{figure}
442
443
444\subsection{Using arrays}
445\index{Sequence} \index{RandomSequence} \index{RegularSequence}
446\index{EnumeratedSequence}
447The example below shows basic usage of arrays, creation, initialization
448and arithmetic operations. Different kind of {\bf Sequence} objects
449can be used for initializing arrays.
450
451\begin{figure}[hbt]
452\dclsbb{Sequence}{RandomSequence}
453\dclsb{RegularSequence}
454\dclsb{EnumeratedSequence}
455\end{figure}
456
457The example below shows basic usage of arrays:
458\index{\tcls{TArray}}
459\begin{verbatim}
460// Creating and initializing a 1-D array of integers
461TArray<int> ia(5);
462EnumeratedSequence es;
463es = 24, 35, 46, 57, 68;
464ia = es;
465cout << "Array<int> ia = " << ia;
466// 2-D array of floats
467TArray<r_4> b(6,4), c(6,4);
468// Initializing b with a constant
469b = 2.71828;
470// Filling c with random numbers
471c = RandomSequence();
472// Arithmetic operations
473TArray<r_4> d = b+0.3f*c;
474cout << "Array<float> d = " << d;
475\end{verbatim}
476
477The copy constructor shares the array data, while the assignment operator
478copies the array elements, as illustrated in the following example:
479\begin{verbatim}
480TArray<int> a1(4,3);
481a1 = RegularSequence(0,2);
482// Array a2 and a1 shares their data
483TArray<int> a2(a1);
484// a3 and a1 have the same size and identical elements
485TArray<int> a3;
486a3 = a1;
487// Changing one of the a2 elements
488a2(1,1,0) = 555;
489// a1(1,1) is also changed to 555, but not a3(1,1)
490cout << "Array<int> a1 = " << a1;
491cout << "Array<int> a3 = " << a3;
492\end{verbatim}
493
494\subsection{Matrices and vectors}
495\index{\tcls{TMatrix}} \index{\tcls{TVector}}
496\begin{figure}[hbt]
497\dclsccc{\tcls{TArray}}{\tcls{TMatrix}}{\tcls{TVector}}
498\end{figure}
499Vectors and matrices are 2 dimensional arrays. The array size
500along one dimension is equal 1 for vectors. Column vectors
501have {\tt NCols() = 1} and row vectors have {\tt NRows() = 1}.
502Mathematical expressions involving matrices and vectors can easily
503be translated into C++ code using {\tt TMatrix} and
504{\tt TVector} objects. {\bf Matrix} and {\bf Vector} are
505typedefs for double precision float matrices and vectors.
506The operator {\bf *} beteween matrices is redefined to
507perform matrix multiplication. One can then write: \\
508\begin{verbatim}
509 // We create a row vector
510 Vector v(1000, BaseArray::RowVector);
511 // Initialize values with a random sequence
512 v = RandomSequence();
513 // Compute the vector length (norm)
514 double norm = (v*v.Transpose()).toScalar();
515 cout << "Norm(v) = " << norm << endl;
516\end{verbatim}
517
518This module contains basic array and matrix operations
519such as the Gauss matrix inversion algorithm
520which can be used to solve linear systems, as illustrated by the
521example below:
522\begin{verbatim}
523#include "sopemtx.h"
524// ...
525// Creation of a random 5x5 matrix
526Matrix A(5,5);
527A = RandomSequence(RandomSequence::Flat);
528Vector X0(5);
529X0 = RandomSequence(RandomSequence::Gaussian);
530// Computing B = A*X0
531Vector B = A*X0;
532// Solving the system A*X = B
533Vector X;
534LinSolve(A, B, X);
535// Checking the result
536Vector diff = X-X0;
537cout << "X-X0= " << diff ;
538double min,max;
539diff.MinMax(min, max);
540cout << " Min(X-X0) = " << min << " Max(X-X0) = " << max << endl;
541\end{verbatim}
542
543\subsection{Working with sub-arrays and Ranges}
544\index{Range}
545A powerful mechanism is included in array classes for working with
546sub-arrays. The class {\bf Range} can be used to specify range of array
547indexes in any of the array dimensions. Any regularly spaced index
548range can be specified, using the {\tt start} and {\tt end} index
549and an optional step (or stride). It is also possible to specify
550the {\tt start} index and the number of elements.
551In the following example, a simple low-pas filter, on a one
552dimensional stream (Vector) has been written using sub-arrays:
553
554\begin{verbatim}
555// Input Vector containing a noisy periodic signal
556 Vector in(1024), out(1024);
557 in = RandomSequence(RandomSequence::Gaussian, 0., 1.);
558 for(int kk=0; kk<in.Size(); kk++)
559 in(kk) += 2*sin(kk*0.05);
560// Compute the output vector by a simple low pass filter
561 Vector out(1024);
562 int w = 2;
563 for(int k=w; k<in.Size()-w; k++)
564 out(k) = in(Range(k-w, k+w).Sum()/(2.*w+1.);
565\end{verbatim}
566
567\newpage
568\subsection{Memory organisation}
569{\tt \tcls{TArray} } can handle numerical arrays with various memory
570organisation, as long as the spacing (steps) along each axis is
571regular. The five axis are labeled X,Y,Z,T,U. The examples below
572illustrates the memory location for a 2-dimensional, $N_x=4 \times N_y=3$.
573The first index is along the X axis and the second index along the Y axis.
574\begin{verbatim}
575 | (0,0) (0,1) (0,2) (0,3) |
576 | (1,0) (1,1) (1,2) (1,3) |
577 | (2,0) (2,1) (2,2) (2,3) |
578\end{verbatim}
579In the first case, the array is completely packed
580($Step_X=1, Step_Y=N_X=4$), with zero offset,
581while in the second case, $Step_X=2, Step_Y=10, Offset=10$:
582\begin{verbatim}
583 | 0 1 2 3 | | 10 12 14 16 |
584Ex1 | 4 5 6 7 | Ex2 | 20 22 24 26 |
585 | 8 9 10 11 | | 30 32 34 36 |
586\end{verbatim}
587
588For matrices and vectors, an optional argument ({\tt MemoryMapping})
589can be used to select the memory mapping, where two basic schemes
590are available: \\
591{\tt CMemoryMapping} and {\tt FortranMemoryMapping}. \\
592In the case where {\tt CMemoryMapping} is used, a given matrix line
593is packed in memory, while the columns are packed when
594{\tt FortranMemoryMapping} is used. The first index when addressing
595the matrix elements (line number index) runs along
596the Y-axis if {\tt CMemoryMapping} is used, and along the X-axis
597in the case of {\tt FortranMemoryMapping}.
598Arithmetic operations between matrices
599with different memory organisation is allowed as long as
600the two matrices have the same sizes (Number of rows and columns).
601The following code example and the corresponding output illustrates
602these two memory mappings. The {\tt \tcls{TMatrix}::TransposeSelf() }
603method changes effectively the matrix memory mapping, which is also
604the case of {\tt \tcls{TMatrix}::Transpose() } method without argument.
605
606\begin{verbatim}
607TArray<r_4> X(4,2);
608X = RegularSequence(1,1);
609cout << "Array X= " << X << endl;
610TMatrix<r_4> X_C(X, true, BaseArray::CMemoryMapping);
611cout << "Matrix X_C (CMemoryMapping) = " << X_C << endl;
612TMatrix<r_4> X_F(X, true, BaseArray::FortranMemoryMapping);
613cout << "Matrix X_F (FortranMemoryMapping) = " << X_F << endl;
614\end{verbatim}
615\newpage
616This code would produce the following output (X\_F = Transpose(X\_C)) :
617\begin{verbatim}
618Array X=
619--- TArray<f> ND=2 SizeX*Y*...= 4x2 ---
6201, 2, 3, 4
6215, 6, 7, 8
622
623Matrix X_C (CMemoryMapping) =
624--- TMatrix<f>(NRows=2, NCols=4) ND=2 SizeX*Y*...= 4x2 ---
6251, 2, 3, 4
6265, 6, 7, 8
627
628Matrix X_F (FortranMemoryMapping) =
629--- TMatrix<f>(NRows=4, NCols=2) ND=2 SizeX*Y*...= 4x2 ---
6301, 5
6312, 6
6323, 7
6334, 8
634\end{verbatim}
635
636\newpage
637
638\section{Module HiStats}
639\begin{figure}[hbt]
640\dclsccc{AnyDataObj}{Histo}{HProf}
641\dclsbb{AnyDataObj}{Histo2D}
642\dclsbb{AnyDataObj}{Ntuple}
643\dclsb{XNtuple}
644\caption{partial class diagram for histograms and ntuples}
645\end{figure}
646
647{\bf HiStats} contains classes for creating, filling, printing and
648doing various operations on one or two dimensional histograms
649{\tt Histo} and {\tt Histo2D} as well as profile histograms {\tt HProf}. \\
650This module also contains {\tt NTuple} and {\tt XNTuple} which are
651more or less the same that the binary FITS tables.
652
653\subsection{1D Histograms}
654\index{Histo}
655For 1D histograms, various numerical methods are provided such as
656computing means and sigmas, finding maxima, fitting, rebinning,
657integrating \dots \\
658The example below shows creating and filling a one dimensionnal histogram
659of 100 bins from $-5.$ to $+5.$ to create a gaussian normal distribution
660with errors~:
661\begin{verbatim}
662#include "histos.h"
663// ...
664Histo H(-0.5,0.5,100);
665H.Errors();
666for(int i=0;i<25000;i++) {
667 double x = NorRand();
668 H.Add(x);
669}
670H.Print(80);
671\end{verbatim}
672
673\subsection{2D Histograms}
674\index{Histo2D}
675Much of these operations are also valid for 2D histograms. 1D projection
676or slices can be set~:
677\begin{verbatim}
678#include "histos2.h"
679// ...
680Histo2D H2(-1.,1.,100,0.,60.,50);
681H2.SetProjX(); // create the 1D histo for X projection
682H2.SetBandX(25.,35.); // create 1D histo projection for 25.<y<35.
683H2.SetBandX(35.,45.); // create 1D histo projection for 35.<y<45.
684H2.SetBandX(40.,55.); // create 1D histo projection for 40.<y<55.
685//... fill H2 with what ever you want
686H2.Print();
687Histo *hx = H2.HProjX();
688 hx->Print(80);
689Histo *hbx2 = HBandX(1); // Get the second X band (35.<y<45.)
690 hbx2->Print(80);
691\end{verbatim}
692
693\subsection{Profile Histograms}
694\index{HProf}
695Profiles histograms {\bf HProf} contains the mean and the
696sigma of the distribution
697of the values filled in each bin. The sigma can be changed to
698the error on the mean. When filled, the profile histogram looks
699like a 1D histogram and much of the operations that can be done on 1D histo
700may be applied onto profile histograms.
701
702\subsection{NTuples}
703\index{NTuple}
704NTuple are memory resident tables of 32 bits floating values (float).
705They are arranged in columns. Each line is often called an event.
706These objects are frequently used to analyze data.
707Graphicals tools (spiapp) can plot a column against an other one
708with respect to various selection cuts. \\
709Here is an example of creation and filling~:
710\begin{verbatim}
711#include "ntuple.h"
712#include "srandgen.h"
713// ...
714char* nament[4] = {"i","x","y","ey"};
715r_4 xnt[4];
716NTuple NT(4,nament);
717for(i=0;i<5000;i++) {
718 xnt[0] = i+1;
719 xnt[1] = 5.*drandpm1(); // a random value between -5 and +5
720 xnt[2] = 100.*exp(-0.5*xnt[1]*xnt[1]) + 1.;
721 xnt[3] = sqrt(xnt[2]);
722 xnt[2] += xnt[3] * NorRand(); // add a random gaussian error
723 NT.Fill(xnt);
724}
725\end{verbatim}
726
727XNTuple are sophisticated NTuple : they accept various types
728of column values (double,float,int,string,...) and can handle
729very large data sets, through swap space on disk.
730\index{XNTuple}
731In the sample code below we show how to create a XNTuple
732object with four columns (double, double, int, string).
733Several entries (lines) are then appended to the table,
734which is saved to a PPF file.
735\begin{verbatim}
736#include "xntuple.h"
737// ...
738char * names[4] = {"X", "X2", "XInt","XStr"};
739// XNTuple (Table) creation with 4 columns, of integer,
740// double(2) and string type
741XNTuple xnt(2,0,1,1, names);
742// Filling the NTuple
743r_8 xd[2];
744int_4 xi[2];
745char xss[2][32];
746char * xs[2] = {xss[0], xss[1]} ;
747for(int i=0; i<50; i++) {
748 xi[0] = i; xd[0] = i+0.5; xd[1] = xd[0]*xd[0];
749 sprintf(xs[0],"X=%g", xd[0]);
750 xnt.Fill(xd, NULL, xi, xs);
751}
752// Printing table info
753cout << xnt ;
754// Saving object into a PPF file
755POutPersist po("xnt.ppf");
756po << xnt ;
757\end{verbatim}
758
759\subsection{Writing, viewing \dots }
760
761All these objects have been design to be written to or read from a persistant file.
762The following example shows how to write the previously created objects
763into such a file~:
764\begin{verbatim}
765//-- Writing
766{
767char *fileout = "myfile.ppf";
768string tag;
769POutPersist outppf(fileout);
770tag = "H"; outppf.PutObject(H,tag);
771tag = "H2"; outppf.PutObject(H2,tag);
772tag = "NT"; outppf.PutObject(NT,tag);
773} // closing ``}'' destroy ``outppf'' and automatically close the file !
774\end{verbatim}
775
776Sophya graphical tools (spiapp) can automatically display and operate
777all these objects.
778
779\subsection{Fourier transform (FFT)}
780\index{FFT} \index{FFTServer}
781
782\newpage
783\section{Module SkyMap}
784\begin{figure}[hbt]
785\dclsbb{AnyDataObj}{PixelMap}
786\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
787\dclsc{SphereThetaPhi}
788\dclsb{LocalMap}
789\caption{partial class diagram for pixelization classes in Sophya}
790\end{figure}
791The {\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.
792\subsection {Spherical maps}
793There are two kinds of spherical maps according pixelization algorithms. SphereHEALPix represents spheres pixelized following the HEALPIix algorithm (E. Yvon, K. Gorski), 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) :
794\index{\tcls{SphereHEALPix}}
795
796\begin{verbatim}
797#include ``spherehealpix.h''
798// ...
799SphereHEALPix<double> sph(8);
800for (int k=0; k< sph.NbPixels(); k++) sph(k) = (double)(10*k);
801\end{verbatim}
802
803SphereThetaPhi is used in a similar way with an argument representing number of slices in theta (Euler angle) for an hemisphere.
804\index{\tcls{SphereThetaPhi}}
805
806\subsection {Local maps}
807\index{\tcls{LocalMap}}
808A 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).
809
810Internally, 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(...))
811
812The 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).
813\begin{verbatim}
814#include "localmap.h"
815//..............
816 LocalMap<r_4> locmap(4,5);
817 for (int k=0; k<locmap.NbPixels();k++) locmap(k)=10.*k;
818 locmap.SetOrigin();
819 locmap.SetSize(30.,30.);
820\end{verbatim}
821
822\subsection{Writing, viewing \dots }
823
824All these objects have been design to be written to or read from a persistant file.
825The following example shows how to write the previously created objects
826into such a file~:
827\begin{verbatim}
828//-- Writing
829
830#include "fiospherehealpix.h"
831//................
832
833char *fileout = "myfile.ppf";
834POutPersist outppf(fileout);
835FIO_SphereHEALPix<r_8> outsph(sph);
836outsph.Write(outppf);
837FIO_LocalMap<r_8> outloc(locmap);
838outloc.Write(outppf);
839
840\end{verbatim}
841
842Sophya graphical tools (spiapp) can automatically display and operate
843all these objects.
844
845
846\section{Module NTools}
847
848This module provides elementary numerical tools for numerical integration,
849fitting, sorting and ODE solving. FFTs are also provided (Mayer,FFTPack).
850
851\subsection{Fitting}
852\index{Fitting} \index{Minimization}
853Fitting is done with two classes {\tt GeneralFit} and {\tt GeneralFitData}
854and is based on the Levenberg-Marquardt method.
855\index{GeneralFit} \index{GeneralFitData}
856GeneralFitData is a class which provide a description of the data
857to be fitted. GeneralFit is the fitter class. Parametrized functions
858can be given as classes which inherit {\tt GeneralFunction}
859or as simple C functions. Classes of pre-defined functions are provided
860(see files fct1dfit.h and fct2dfit.h). The user interface is very close
861from that of the CERN {\tt Minuit} fitter.
862Number of objects (Histo, HProf \dots ) are interfaced with GeneralFit
863and can be easily fitted. \\
864Here is a very simple example for fitting the previously created NTuple
865with a gaussian~:
866\begin{verbatim}
867#include "fct1dfit.h"
868// ...
869
870// Read from ppf file
871NTuple nt;
872{
873PInPersist pis("myfile.ppf");
874string tag = "NT"; pis.GetObject(nt,tag);
875}
876
877// Fill GeneralData
878GeneralData mGdata(nt.NEntry());
879for(int i=0; i<nt.NEntry(); i++)
880 mGdata.AddData1(xnt[1],xnt[2],xnt[3]); // Fill x, y and error on y
881mGData.PrintStatus();
882
883// Function for fitting : y = f(x) + noise
884Gauss1DPol mFunction; // gaussian + constant
885
886// Prepare for fit
887GeneralFit mFit(&mFunction); // create a fitter for the choosen function
888mFit.SetData(&mGData); // connect data to the fitter
889
890// Set and initialize the parameters (that's non-linear fitting!)
891// (num par, name, guess start, step, [limits min and max])
892mFit.SetParam(0,"high",90.,1..);
893mFit.SetParam(1,"xcenter",0.05,0.01);
894mFit.SetParam(2,"sigma",sig,0.05,0.01,10.);
895 // Give limits to avoid division by zero
896mFit.SetParam(3,"constant",0.,1.);
897
898// Fit and print result
899int rcfit = mFit.Fit();
900mFit.PrintFit();
901if(rcfit>0) {)
902 cout<<"Reduce_Chisquare = "<<mFit.GetChi2Red()
903 <<" nstep="<<mFit.GetNStep()<<" rc="<<rcfit<<endl;
904} else {
905 cout<<"Fit_Error, rc = "<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
906 mFit.PrintFitErr(rcfit);
907}
908
909// Get the result for further use
910TVector<r_8> ParResult = mFit.GetParm();
911cout<<ParResult;
912\end{verbatim}
913
914Much more usefull possibilities and detailed informations might be found
915in the HTML pages of the Sophya manual.
916
917\subsection{Polynomial}
918\index{Polynomial} \index{Poly} \index{Poly2}
919Polynomials of 1 or 2 variables are supported ({\tt Poly} and {\tt Poly2}).
920Various operations are supported~:
921\begin{itemize}
922\item elementary operations between polynomials $(+,-,*,/) $
923\item setting or getting coefficients
924\item computing the value of the polynomial for a given value
925 of the variable(s),
926\item derivating
927\item computing roots (degre 1 or 2)
928\item fitting the polynomial to vectors of data.
929\end{itemize}
930Here is an example of polynomial fitting~:
931\begin{verbatim}
932#include "poly.h"
933// ...
934Poly pol(2);
935pol[0] = 100.; pol[1] = 0.; pol[2] = 0.01; // Setting coefficients
936TVector<r_8> x(100);
937TVector<r_8> y(100);
938TVector<r_8> ey(100);
939for(int i=0;i<100;i++) {
940 x(i) = i;
941 ey(i) = 10.;
942 y(i) = pol((double) i) + ey(i)*NorRand();
943 ey(i) *= ey(i)
944}
945
946TVector<r_8> errcoef;
947Poly polfit;
948polfit.Fit(x,y,ey,2,errcoef);
949
950cout<<"Fit Result"<<polfit<<endl;
951cout<<"Errors :"<<errcoef;
952\end{verbatim}
953
954Similar operations can be done on polynomials with 2 variables.
955
956\section{Module Samba}
957\index{Spherical Harmonics}
958\index{SphericalTransformServer}
959The 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 analyzed back to Cl's...
960\begin{verbatim}
961#include "skymap.h"
962#include "samba.h"
963....................
964
965// Generate input spectra a + b* l + c * gaussienne(l, 50, 20)
966int lmax = 92;
967Vector clin(lmax);
968for(int l=0; l<lmax; l++) {
969 double xx = (l-50.)/10.;
970 clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
971}
972
973// Compute map from spectra
974SphericalTransformServer<r_8> ylmserver;
975int m = 128; // HealPix pixelisation parameter
976SphereHEALPix<r_8> map(m);
977ylmserver.GenerateFromCl(map, m, clin, 0.);
978// Compute power spectrum from map
979Vector clout = ylmserver.DecomposeToCl(map, lmax, 0.);
980\end{verbatim}
981
982\section{Module SkyT}
983
984\section{Module FitsIOServer}
985\begin{figure}[hbt]
986\dclsbb{FitsFile}{FitsInFile}
987\dclsb{FitsOutFile}
988\end{figure}
989\index{FITS}
990This 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 :
991\begin{verbatim}
992#include "spherehealpix.h"
993#include "fitsspherehealpix.h"
994#include "fitstarray.h"
995#include "tmatrix.h"
996//...........................
997
998int m=...;
999SphereHEALPix<r_8> sph(m);
1000................
1001int dim1=...;
1002int dim2=...;
1003TMatrix<r_8> mat(dim1,dim2);
1004............
1005
1006FITS_SphereHEALPix<r_8> sph_temp(sph);
1007FITS_TArray<r_8> mat_temp(mat);
1008// writing
1009
1010FitsOutFile os("myfile.fits");
1011sph_temp.Write(os);
1012mat_temp.Write(os);
1013
1014// reading
1015FitsInFile is("myfile.fits");
1016sph_temp.Read(is);
1017mat_temp.Read(is);
1018SphereHEALPix<r_8> new_sph=(SphereHEALPix<r_8>)sph_temp;
1019TMatrix<r_8> new_mat=(TMatrix<r_8>)mat_temp;
1020................
1021
1022\end{verbatim}
1023
1024
1025
1026\newpage
1027\section{Building and installing Sophya}
1028Presently, the Sophya library compilation has been tested with the following
1029compiler/platform pairs:
1030
1031\begin{center}
1032\begin{tabular}{ll}
1033Compaq/DEC OSF1 & cxx (6.0 , 6.2) \\
1034Linux & g++ (2.91 , 2.95) \\
1035Linux & KCC (3.4) \\
1036Solaris & g++ (2.95) \\
1037SGI IRIX64 & CC \\
1038\end{tabular}
1039\end{center}
1040
1041Some of the modules in the Sophya package uses external libraries. The
1042{\bf FitsIOServer} is the example of such a module, where the {\tt libcfitsio.a}
1043is used.
1044The build procedure expects to find the include files and the libraries in: \\
1045{\tt \$EXTLIBDIR/Include/FitsIO } \\
1046{\tt \$EXTLIBDIR/`uname`-\$EROSCXX/Libs} \\
1047
1048The object files from a given Sophya module are grouped in an archive library
1049with the module's name ({\tt libmodulename.a}). All Sophya modules
1050 are grouped in a single shared library ({\tt libsophya.so}), while the
1051modules with reference to external libraries are grouped in
1052({\tt libextsophya.so}). The {\bf PI} and {\bf PIext} modules are
1053grouped in ({\tt libPI.so}).
1054
1055The environment variables {\bf DPCDEVREP}, {\bf EXTLIBDIR} and {\bf EROSCXX}
1056must be defined in order to install the Sophya package.
1057In the example below, we assume that we want to install Sophya from a
1058released (tagged) version in the source directory {\tt \$SRC} in the
1059{\tt /usr/local/Sophya} diretory, using {\tt g++}. We assume that
1060the external libraries directory tree has been set up in
1061{\tt /usr/local/ExtLibs/}. \\[3mm]
1062\centerline{
1063 \rule{20mm}{0.5mm}
1064 {\bf \large the use of GNU make is mandatory}
1065 \rule{20mm}{0.5mm} }
1066
1067\vspace*{3mm}
1068\begin{verbatim}
1069# We select our C++ compiler
1070csh> setenv EROSCXX g++
1071# Setup the build directory
1072csh> mkdir /usr/local/Sophya/
1073csh> setenv DPCDEVREP /usr/local/Sophya/
1074csh> setenv EXTLIBDIR /usr/local/ExtLibs/
1075# Use the top level makefile in Mgr/
1076csh> cd \$SRC
1077csh> cp Mgr/Makefile Makefile
1078# Step 1: Create the directory tree and copy the include files (.h)
1079csh> make depend
1080# Step 2: Compile the modules without external library reference
1081csh> make libs
1082# Step 3: Compile the modules WITH external library reference (optional)
1083csh> make extlibs
1084# Step 4: Build libsophya.so
1085csh> make slb
1086# Step 5: Build libextsophya.so (optional)
1087csh> make slbext
1088# Step 6: Compile the PI and PIext modules (optional)
1089csh> make PI
1090# Step 7: Build the corresponding shared library libPI.so (optional)
1091csh> make slbpi
1092\end{verbatim}
1093
1094To compile all modules and build the shared libraries, it is possible
1095to use:
1096\begin{verbatim}
1097# Step 2,3,6
1098csh> make all
1099# Step 4,5,7
1100csh> make slball
1101\end{verbatim}
1102
1103At this step, all libraries sould have been made. Programs using
1104Sophya libraries can now be built:
1105\begin{verbatim}
1106# To compile test programs
1107csh> cd Tests
1108csh> make arrt ...
1109csh> cd ..
1110# To compile other programs, for example from the PMixer module
1111csh> cd PMixer
1112csh> make
1113csh> cd ..
1114# To build (s)piapp (libPI.so is needed)
1115csh> cd ProgPI
1116csh> make
1117csh> cd ..
1118\end{verbatim}
1119
1120\subsection{Mgr module}
1121This module contains scripts which can be used for generating the
1122makefiles for each module.
1123\begin{itemize}
1124\item {\bf Makefile} Top level Makefile for builiding the libraries.
1125\item {\bf Makefile.h} contains the defintion of compilation flags for the
1126different compilers and systems. This file is used for building the
1127library and generating {\bf MakefileUser.h} (to be included in makefiles).
1128\item {\bf Makefile.slb} contains the rules for building shared libraries
1129for the different compilers and systems. (to be included in makefiles)
1130\item {\bf crerep\_sophya} c-shell script for creating the directory tree
1131under {\tt \$DPCBASEREP} and {\tt \$DPCDEVREP}
1132\item {\bf install\_sophya} c-shell script for installing the Sophya package.
1133Usually from {\tt \$DPCDEVREP} to {\tt \$DPCBASEREP}
1134\item {\bf mkmflien} c-shell script for making symbolic links or copying
1135include files to {\tt \$DPCDEVREP/Include} or {\tt \$DPCBASEREP/Include}
1136\item {\bf mkmf} c-shell script for generating module makefiles and the
1137top level makefile (named GNUmakefile)
1138\item {\bf mkmflib} c-shell script for generating each library module
1139makefile (named GNUmakefile)
1140\item {\bf mkmfprog} c-shell script for generating makefile for a module
1141containing the source for executable programs (named GNUmakefile).
1142\item {\bf mkmfPI} c-shell script for generating makefile for PI and PIext
1143modules (named GNUmakefile)
1144\item {\bf libdirs} List of Sophya modules without reference to external
1145libraries.
1146\item {\bf extlibdirs} List of Sophya modules with reference to external
1147libraries.
1148
1149\end{itemize}
1150
1151\newpage
1152\appendix
1153\section{SOPHYA Exceptions}
1154\index{Exception classes} \index{PThrowable} \index{PError} \index{PException}
1155SOPHYA library defines a set of exceptions which are used
1156for signaling error conditions. The figure below shows a partial
1157class diagram for exception classes in SOPHYA.
1158\begin{figure}[hbt]
1159\dclsbb{PThrowable}{PError}
1160\dclscc{PError}{AllocationError}
1161\dclscc{PError}{NullPtrError}
1162\dclscc{PError}{ForbiddenError}
1163\dclscc{PError}{AssertionFailedError}
1164\dclsbb{PThrowable}{PException}
1165\dclscc{PException}{IOExc}
1166\dclscc{PException}{SzMismatchError}
1167\dclscc{PException}{RangeCheckError}
1168\dclscc{PException}{ParmError}
1169\dclscc{PException}{TypeMismatchExc}
1170\dclscc{PException}{MathExc}
1171\dclscc{PException}{CaughtSignalExc}
1172\caption{partial class diagram for exception handling in Sophya}
1173\end{figure}
1174
1175For simple programs, it is a good practice to handle
1176the exceptions at least at high level, in the {\tt main()} function.
1177The example below shows the exception handling and the usage
1178of Sophya persistence.
1179
1180\input{ex1.inc}
1181
1182
1183\newpage
1184\addcontentsline{toc}{section}{Index}
1185\printindex
1186\end{document}
1187
Note: See TracBrowser for help on using the repository browser.