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

Last change on this file since 2452 was 2437, checked in by cmv, 22 years ago

chgt DPC...REP en SOPHYA...REP et EROSCXX en SOPHYACXX cmv 17/09/2003

File size: 59.0 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.5 (V\_Dec2002) }
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.
61We have decided to use as much as possible available
62numerical analysis libraries, encapsulating them whenever
63possible. Although some the modules in SOPHYA have been
64designed with the specific goal of providing some of the tools for
65the Planck-HFI data processing software, most of the
66packages presented here have a more general scope than the
67CMB analysis and Planck mission problem.
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/CVSPlanck}
80is organised into a number of modules.
81
82\begin{itemize}
83\item[] {\bf Mgr/} 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 XNTuple} \ldots)
95\item[] {\bf SkyMap/} Local and full sky maps, and few geometry
96handling utility classes. \\
97({\tt PixelMap<T>, LocalMap<T>, SphericalMap<T>, \ldots})
98\item[] {\bf SUtils/} This module contains
99few utility classes, such as the
100{\tt DataCard} class, as well as string manipulation functions in C and C++.
101\item[] {\bf SysTools/} This module contains classes implementing
102an interface to various OS specific services, such
103as shared object and dynamic link handling.
104
105\end{itemize}
106
107The modules listed below are more tightly related to the
108CMB (Cosmic Microwave Background) data analysis problem:
109\begin{itemize}
110\item[] {\bf SkyT/}
111classes for spectral emission and detector frequency response modelling \\
112({\tt SpectralResponse, RadSpectra, BlackBody} \ldots)
113\item[] {\bf Samba/} Spherical harmonic analysis, noise generators \ldots
114\end{itemize}
115
116The following modules contain the interface classes with
117external libraries:
118\begin{itemize}
119\item[] {\bf FitsIOServer/} Classes for handling file input-output
120in FITS format using the cfitsio library.
121\item[] {\bf LinAlg/} Interface with Lapack linear algebra package
122\item[] {\bf IFFTW/} Interface with FFTW package (libfftw.a)
123\item[] {\bf XAstroPack/} Interface to some common astronomical
124computation libraries. Presently, this module uses an external library
125extracted from the {\bf Xephem } source code. The corresponding source
126code is also available from SOPHYA cvs repository, module {\bf XephemAstroLib}.
127\end{itemize}
128
129The following modules contain each a set of related programs using the
130SOPHYA library.
131\begin{itemize}
132\item[] {\bf Tests/} Simple test programs
133\item[] {\bf PrgUtil/} Various utility programs (runcxx, scanppf, scanfits, \ldots)
134\item[] {\bf PrgMap/} Programs performing operations on skymaps: projections,
135power spectrum in harmonic space, \ldots
136\item[] {\bf PMixer/} skymixer and related programs
137\end{itemize}
138
139As a companion to SOPHYA, the {\bf (s)piapp} interactive data analysis
140program is built on top of SOPHYA and the {\bf PI} GUI class library
141and application framework. The {\bf PI} ({\bf P}eida {\bf Interactive})
142development started in 1995, in the EROS \footnote{EROS: {\bf E}xp\'erience
143de {\bf R}echerche d'{\bf O}bjets {\bf S}ombres - http://eros.in2p3.fr}
144microlensing search collaboration, with PEIDA++ \footnote {PEIDA++:
145The EROS data analysis class library -
146http://www.lal.in2p3.fr/recherche/eros/PeidaDoc/}.
147The {\bf PI} documentation and the {\bf piapp} user's guide are available
148from \href{http://www.sophya.org}{http://www.sophya.org}.
149%\href{http://www.sophya.org}{http://www.sophya.org}.
150The {\bf PI} is organized as the following modules:
151\begin{itemize}
152\item[] {\bf PI/} Portable GUI class library and application development
153framework kernel.
154\item[] {\bf PIGcont/} Contour-plot drawing classes.
155\item[] {\bf PIext/} Specific drawers and adapters for SOPHYA objects,
156and the {\bf piapp} interactive data analysis framework.
157\item[] {\bf ProgPI/} interactive analysis tool main program and pre-loaded
158modules.
159\end{itemize}
160
161Modules containing examples and demo programs:
162\begin{itemize}
163\item[] {\bf Examples/} Sample SOPHYA codes and example programs and
164makefiles (auto\_makefile and ex\_makefile).
165\item[] {\bf DemoPIApp/} Sample exripts and programs for (s)piapp
166interactive analysis tools.
167\end{itemize}
168\newpage
169
170\section{Using Sophya}
171Basic usage of Sophya classes are described in in the following sections.
172Complete Sophya documentation can be found at our web site
173{\bf http://www.sophya.org}.
174
175\subsection{Environment variables}
176Two environment variables {\bf SOPHYABASEREP} and {\bf SOPHYACXX} are used
177to define the path where the Sophya libraries and executable are installed.
178{\bf SOPHYABASEREP} defines the base directory path and {\bf SOPHYACXX} the
179name of the C++ compiler. The complete path is built using {\bf SOPHYABASEREP},
180the operating system name (as obtained by the {\tt uname} command), and
181the compiler name. In the example below, we show the complete path
182for a {\tt Linux} system, using the GNU g++ compiler:
183
184\begin{itemize}
185\item \$SOPHYABASEREP/Include : Include (.h) files
186\item \$SOPHYABASEREP/Linux-g++/Libs : Path for the archive libraries (.a)
187\item \$SOPHYABASEREP/Linux-g++/ShLibs : Shared library path (.so)
188\item \$SOPHYABASEREP/Linux-g++/Exec : Executable file path
189\end{itemize}
190
191In order to use the shared libraries, the {\bf LD\_LIBRARY\_PATH} variable
192should contain the Sophya shared library path
193({\tt \$SOPHYABASEREP/Linux-g++/ShLibs } when using g++ compiler on Linux)
194
195For modules using external libraries, the {\bf EXTLIBDIR}
196environment variable should contain the path to these libraries
197and corresponding include files.
198C-FitsIO anf FFTW include files should be accessible through: \\
199{\tt \$EXTLIBDIR/Include/FitsIO } \\
200{\tt \$EXTLIBDIR/Include/FFTW } \\
201The corresponding libraries are expected to be found in: \\
202{\tt \$EXTLIBDIR/Linux-g++/Libs} \\
203
204\subsection{User makefiles}
205The file {\tt \$SOPHYABASEREP/Include/MakefileUser.h} defines the compilation
206flags and the list of Sophya libraries. It should be included in the
207user's makefile. The default compilation rules assumes that the object (.o)
208and executable files would be put in the following directories: \\
209{\tt \$HOME/`uname`-\$SOPHYACXX/Objs} \\
210{\tt \$HOME/`uname`-\$SOPHYACXX/Exec}.
211In the case of a {\tt Linux} system and using {\tt g++} as the C++ compiler,
212these two directories would be translated to \\
213{\tt \$HOME/Linux-g++/Objs} and {\tt \$HOME/Linux-g++/Exec}.
214The GNU make program should be used.
215\par
216The file {\tt Examples/auto\_makefile} defines the rules to compile
217a given source program, and link it against the Sophya libraries to produce
218an executable. The example below shows the steps to compile a program named
219{\tt trivial.cc }.
220\begin{verbatim}
221csh> cp Examples/auto_makefile makefile
222csh> cp Examples/ex1.cc trivial.cc
223csh> make trivial
224\end{verbatim}
225This command should compile the {\tt trivial.cc} file,
226and link it against the sophya libraries.
227The file {\tt Examples/ex\_makefile} provides another
228example makefile.
229
230\subsection{the runcxx program}
231\index{runcxx}
232{\bf runcxx} is a simple program which can be used to compile, link
233and run simple C++ programs. It handles the creation of a
234complete program file, containing the basic set C++ include files,
235the necessary include files for SOPHYA SysTools, TArray, HiStats
236and NTools modules, and the main program with exception handling.
237Other Sophya modules can be included using the {\tt -import} flag.
238Use of additional include files can be specified using the
239{\tt -inc} flag.
240\begin{verbatim}
241csh> runcxx -h
242SOPHYA Version 1.1 Revision 0 (V_Fev2001) -- Feb 28 2001 11:19:17 cxx
243 runcxx : compiling and running of a piece of C++ code
244 Usage: runcxx [-compopt CompileOptions] [-linkopt LinkOptions]
245 [-tmpdir TmpDirectory] [-f C++CodeFileName]
246 [-inc includefile] [-inc includefile ...]
247 [-import modulename] [-import modulename ...]
248 [-uarg UserArg1 UserArg2 ...]
249 if no file name is specified, read from standard input
250 modulenames: SkyMap, Samba, SkyT, FitsIOServer, LinAlg, IFFTW
251\end{verbatim}
252Most examples in this manual can be tested using runcxx. The
253example below shows how to compile, link and run a sample
254code.
255\begin{verbatim}
256// File example.icc
257Matrix a(3,3);
258a = IdentityMatrix(1.);
259cout << a ;
260// Executing this sample code
261csh> runcxx -f example.icc
262\end{verbatim}
263
264\subsection{the scanppf program}
265{\bf scanppf} is a simple SOPHYA application which can be used to check
266PPF files and list their contents.
267\begin{verbatim}
268csh> scanppf -h
269 PIOPersist::Initialize() Starting Sophya Persistence management service
270SOPHYA Version 1.4 Revision 0 (V_Nov2002) -- Nov 15 2002 10:32:12 cxx
271 Usage: scanppf filename [s/n/a0/a1/a2/a3]
272 s[=default} : Sequential reading of objects
273 n : Object reading at NameTags
274 a0...a3 : Tag List with PInPersist.AnalyseTags(0...3)
275\end{verbatim}
276
277
278\newpage
279
280\section{Copy constructor and assignment operator}
281In C++, objects can be copied by assignment or by initialisation.
282Copying by initialisation corresponds to creating an object and
283initialising its value through the copy constructor.
284The copy constructor has its first argument as a reference, or
285const reference to the object's class type. It can have
286more arguments, if default values are provided.
287Copying by assignment applies to an existing object and
288is performed through the assignment operator (=).
289The copy constructor implements this for identical type objects:
290\begin{verbatim}
291class MyObject {
292public:
293 MyObject(); // Default constructor
294 MyObject(MyObject const & a); // Copy constructor
295 MyObject & operator = (MyObject const & a) // Assignment operator
296}
297\end{verbatim}
298The copy constructors play an important role, as they are
299called when class objects are passed by value,
300returned by value, or thrown as an exception.
301\begin{verbatim}
302// A function declaration with an argument of type MyObject,
303// passed by value, and returning a MyObject
304MyObject f(MyObject x)
305{
306 MyObject r;
307 ...
308 return(r); // Copy constructor is called here
309}
310// Calling the function :
311MyObject a;
312f(a); // Copy constructor called for a
313\end{verbatim}
314It should be noted that the C++ syntax is ambiguous for the
315assignment operator. {\tt MyObject x; x=y; } and
316{\tt MyObject x=y;} have different meaning.
317\begin{verbatim}
318MyObject a; // default constructor call
319MyObject b(a); // copy constructor call
320MyObject bb = a; // identical to bb(a) : copy constructor call
321MyObject c; // default constructor call
322c = a; // assignment operator call
323\end{verbatim}
324
325As a general rule in SOPHYA, objects which implements
326reference sharing on their data members have a copy constructor
327which shares the data, while the assignment operator copies or
328duplicate the data.
329
330\newpage
331\section{Module BaseTools}
332
333{\bf BaseTools} contains utility classes such as
334{\tt DVlist}, an hierarchy of exception classes for Sophya, a template
335class {\tcls{NDataBlock}} for handling reference counting on numerical
336arrays, as well as classes providing the services for implementing simple
337serialisation.
338\vspace*{5mm}
339
340\subsection{SOPHYA persistence}
341\index{PPersist} \index{PInPersist} \index{POutPersist}
342\begin{figure}[hbt]
343\dclsa{PPersist}
344\dclsbb{PIOPersist}{PInPersist}
345\dclsb{POutPersist}
346\caption{partial class diagram for classes handling persistence in Sophya}
347\end{figure}
348A simple persistence mechanism is defined in SOPHYA. Its main
349features are:
350\begin{itemize}
351\item[] Portable file format, containing the description of the data structures
352and object hierarchy. \\
353{\bf PPF} {\bf P}ortable {\bf P}ersistence file {\bf F}ormat.
354\index{PPF}
355\item[] Handling of read/write for multiply referenced objects.
356\item[] All write operations are carried using sequential access only. This
357holds also for read operations, unless positional tags are used.
358SOPHYA persistence services can thus be used to transfer objects
359through network links.
360\item[] The serialisation (reading/writing) for objects for a given class
361is implemented through a handler object. The handler class inherits
362from {\tt PPersist} class.
363\item[] A run time registration mechanism is used in conjunction with
364RTTI (Run Time Type Identification) for identifying handler classes
365when reading {\bf PInPersist} streams, or for associating handlers
366with data objects {\bf AnyDataObject} for write operations.
367\end{itemize}
368A complete description of SOPHYA persistence mechanism and guidelines
369for writing delegate classes for handling object persistence is beyond
370the scope of this document. The most useful methods for using Sophya
371persistence are listed below:
372\begin{itemize}
373\item[] {\tt POutPersist::PutObject(AnyDataObj \& o)} \\
374Writes the data object {\bf o} to the output stream.
375\item[] {\tt POutPersist::PutObject(AnyDataObj \& o, string tagname)} \\
376Writes the data object {\bf o} to the output stream, associated with an
377identification tag {\bf tagname}.
378\item[] {\tt PInPersist::GetObject(AnyDataObj \& o)} \\
379Reads the next object in stream into {\bf o}. An exception is
380generated for incompatible object types.
381\item[] {\tt PInPersist::GetObject(AnyDataObj \& o, string tagname)} \\
382Reads the object associated with the tag {\bf tagname} into {\bf o}.
383An exception is generated for incompatible object types.
384\end{itemize}
385The operators {\tt operator << (POutPersist ...) } and
386{\tt operator >> (PInPersist ...) } are often overloaded
387to perform {\tt PutObject()} and {\tt GetObject()} operations,
388as illustrated in the example below:
389\begin{verbatim}
390// Creating and filling a histogram
391Histo hw(0.,10.,100);
392...
393// Writing histogram to a PPF stream
394POutPersist os("hw.ppf");
395os << hw;
396// Reading a histogram from a PPF stream
397PInPersist is("hr.ppf");
398is >> hr;
399\end{verbatim}
400
401The {\bf scanppf} program can be used to list the content of a PPF file.
402\index{scanppf}
403\begin{verbatim}
404csh> scanppf -h
405SOPHYA Version 1.1 Revision 0 (V_Fev2001) -- Feb 28 2001 11:19:17 cxx
406 Usage: scanppf filename [s/n/a0/a1/a2/a3]
407 s[=default} : Sequential reading of objects
408 n : Object reading at NameTags
409 a0...a3 : Tag List with PInPersist.AnalyseTags(0...3)
410\end{verbatim}
411
412
413\subsection{\tcls{NDataBlock}}
414\index{\tcls{NDataBlock}}
415\begin{figure}[hbt]
416\dclsbb{AnyDataObj}{\tcls{NDataBlock}}
417\dclsbb{PPersist}{\tcls{FIO\_NDataBlock}}
418\end{figure}
419The {\bf \tcls{NDataBlock}} is designed to handle reference counting
420and sharing of memory blocs (contiguous arrays) for numerical data
421types. Initialisation, resizing, basic arithmetic operations, as
422well as persistence handling services are provided.
423The persistence handler class ({\tt \tcls{FIO\_NDataBlock}}) insures
424that a single copy of data is written for multiply referenced objects,
425and the data is shared among objects when reading.
426\par
427The example below shows writing of NDataBlock objects through the
428use of overloaded operator $ << $ :
429\begin{verbatim}
430#include "fiondblock.h"
431// ...
432POutPersist pos("aa.ppf");
433NDataBlock<r_4> rdb(40);
434rdb = 567.89;
435pos << rdb;
436// We can also use the PutObject method
437NDataBlock<int_4> idb(20);
438idb = 123;
439pos.PutObject(idb);
440\end{verbatim}
441The following sample programs show the reading of the created PPF file :
442\begin{verbatim}
443PInPersist pis("aa.ppf");
444NDataBlock<r_4> rdb;
445pis >> rdb;
446cout << rdb;
447NDataBlock<int_4> idb;
448cout << idb;
449\end{verbatim}
450
451\subsection{Using DVList}
452\index{DVList} \index{MuTyV}
453\begin{figure}[hbt]
454\dclsbb{AnyDataObj}{DVList}
455\dclsbb{PPersist}{\tclsc{ObjFileIO}{DVList}}
456\end{figure}
457The {\bf DVList} class objects can be used to create and manage list
458of values, associated with names. A list of pairs of (MuTyV, name(string))
459is maintained by DVList objects. {\bf MuTyV} is a simple class
460capable of holding string, integer, float or complex values,
461providing easy conversion methods between these objects.
462\begin{verbatim}
463// Using MuTyV objects
464MuTyV s("hello"); // string type value
465MuTyV x;
466x = "3.14159626"; // string type value, ASCII representation for Pi
467double d = x; // x converted to double = 3.141596
468x = 314; // x contains the integer value = 314
469// Using DVList
470DVList dvl;
471dvl("Pi") = 3.14159626; // float value, named Pi
472dvl("Log2") = 0.30102999; // float value, named Log2
473dvl("FileName") = "myfile.fits"; // string value, named myfile.fits
474// Printing DVList object
475cout << dvl;
476\end{verbatim}
477
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\dclsb{XNtuple}
873\caption{partial class diagram for histograms and ntuples}
874\end{figure}
875
876{\bf HiStats} contains classes for creating, filling, printing and
877doing various operations on one or two dimensional histograms
878{\tt Histo} and {\tt Histo2D} as well as profile histograms {\tt HProf}. \\
879This module also contains {\tt NTuple} and {\tt XNTuple} which are
880more or less the same that the binary FITS tables.
881
882\subsection{1D Histograms}
883\index{Histo}
884For 1D histograms, various numerical methods are provided such as
885computing means and sigmas, finding maxima, fitting, rebinning,
886integrating \dots \\
887The example below shows creating and filling a one dimensional histogram
888of 100 bins from $-5.$ to $+5.$ to create a Gaussian normal distribution
889with errors~:
890\begin{verbatim}
891#include "histos.h"
892// ...
893Histo H(-0.5,0.5,100);
894H.Errors();
895for(int i=0;i<25000;i++) {
896 double x = NorRand();
897 H.Add(x);
898}
899H.Print(80);
900\end{verbatim}
901
902\subsection{2D Histograms}
903\index{Histo2D}
904Much of these operations are also valid for 2D histograms. 1D projection
905or slices can be set~:
906\begin{verbatim}
907#include "histos2.h"
908// ...
909Histo2D H2(-1.,1.,100,0.,60.,50);
910H2.SetProjX(); // create the 1D histo for X projection
911H2.SetBandX(25.,35.); // create 1D histo projection for 25.<y<35.
912H2.SetBandX(35.,45.); // create 1D histo projection for 35.<y<45.
913H2.SetBandX(40.,55.); // create 1D histo projection for 40.<y<55.
914//... fill H2 with what ever you want
915H2.Print();
916Histo *hx = H2.HProjX();
917 hx->Print(80);
918Histo *hbx2 = HBandX(1); // Get the second X band (35.<y<45.)
919 hbx2->Print(80);
920\end{verbatim}
921
922\subsection{Profile Histograms}
923\index{HProf}
924Profiles histograms {\bf HProf} contains the mean and the
925sigma of the distribution
926of the values filled in each bin. The sigma can be changed to
927the error on the mean. When filled, the profile histogram looks
928like a 1D histogram and much of the operations that can be done on 1D histo
929may be applied onto profile histograms.
930
931\subsection{Data tables (tuples)}
932\index{NTuple}
933NTuple are memory resident tables of 32 bits floating values (float).
934They are arranged in columns. Each line is often called an event.
935These objects are frequently used to analyze data.
936Graphicals tools (spiapp) can plot a column against an other one
937with respect to various selection cuts. \\
938Here is an example of creation and filling~:
939\begin{verbatim}
940#include "ntuple.h"
941#include "srandgen.h"
942// ...
943char* nament[4] = {"i","x","y","ey"};
944r_4 xnt[4];
945NTuple NT(4,nament);
946for(i=0;i<5000;i++) {
947 xnt[0] = i+1;
948 xnt[1] = 5.*drandpm1(); // a random value between -5 and +5
949 xnt[2] = 100.*exp(-0.5*xnt[1]*xnt[1]) + 1.;
950 xnt[3] = sqrt(xnt[2]);
951 xnt[2] += xnt[3] * NorRand(); // add a random gaussian error
952 NT.Fill(xnt);
953}
954\end{verbatim}
955
956XNTuple are sophisticated NTuple : they accept various types
957of column values (double,float,int,string,...) and can handle
958very large data sets, through swap space on disk.
959\index{XNTuple}
960In the sample code below we show how to create a XNTuple
961object with four columns (double, double, int, string).
962Several entries (lines) are then appended to the table,
963which is saved to a PPF file.
964\begin{verbatim}
965#include "xntuple.h"
966// ...
967char * names[4] = {"X", "X2", "XInt","XStr"};
968// XNTuple (Table) creation with 4 columns, of integer,
969// double(2) and string type
970XNTuple xnt(2,0,1,1, names);
971// Filling the NTuple
972r_8 xd[2];
973int_4 xi[2];
974char xss[2][32];
975char * xs[2] = {xss[0], xss[1]} ;
976for(int i=0; i<50; i++) {
977 xi[0] = i; xd[0] = i+0.5; xd[1] = xd[0]*xd[0];
978 sprintf(xs[0],"X=%g", xd[0]);
979 xnt.Fill(xd, NULL, xi, xs);
980}
981// Printing table info
982cout << xnt ;
983// Saving object into a PPF file
984POutPersist po("xnt.ppf");
985po << xnt ;
986\end{verbatim}
987
988\subsection{Writing, viewing \dots }
989
990All these objects have been design to be written to or read from a persistent file.
991The following example shows how to write the previously created objects
992into such a file~:
993\begin{verbatim}
994//-- Writing
995{
996char *fileout = "myfile.ppf";
997string tag;
998POutPersist outppf(fileout);
999tag = "H"; outppf.PutObject(H,tag);
1000tag = "H2"; outppf.PutObject(H2,tag);
1001tag = "NT"; outppf.PutObject(NT,tag);
1002} // closing ``}'' destroy ``outppf'' and automatically close the file !
1003\end{verbatim}
1004
1005Sophya graphical tools (spiapp) can automatically display and operate
1006all these objects.
1007
1008\newpage
1009\section{Module NTools}
1010
1011This module provides elementary numerical tools for numerical integration,
1012fitting, sorting and ODE solving. FFTs are also provided (Mayer,FFTPack).
1013
1014\subsection{Fitting}
1015\index{Fitting} \index{Minimisation}
1016Fitting is done with two classes {\tt GeneralFit} and {\tt GeneralFitData}
1017and is based on the Levenberg-Marquardt method.
1018\index{GeneralFit} \index{GeneralFitData}
1019GeneralFitData is a class which provide a description of the data
1020to be fitted. GeneralFit is the fitter class. Parametrized functions
1021can be given as classes which inherit {\tt GeneralFunction}
1022or as simple C functions. Classes of pre-defined functions are provided
1023(see files fct1dfit.h and fct2dfit.h). The user interface is very close
1024from that of the CERN {\tt Minuit} fitter.
1025Number of objects (Histo, HProf \dots ) are interfaced with GeneralFit
1026and can be easily fitted. \\
1027Here is a very simple example for fitting the previously created NTuple
1028with a Gaussian~:
1029\begin{verbatim}
1030#include "fct1dfit.h"
1031// ...
1032
1033// Read from ppf file
1034NTuple nt;
1035{
1036PInPersist pis("myfile.ppf");
1037string tag = "NT"; pis.GetObject(nt,tag);
1038}
1039
1040// Fill GeneralData
1041GeneralData mGdata(nt.NEntry());
1042for(int i=0; i<nt.NEntry(); i++)
1043 mGdata.AddData1(xnt[1],xnt[2],xnt[3]); // Fill x, y and error on y
1044mGData.PrintStatus();
1045
1046// Function for fitting : y = f(x) + noise
1047Gauss1DPol mFunction; // gaussian + constant
1048
1049// Prepare for fit
1050GeneralFit mFit(&mFunction); // create a fitter for the choosen function
1051mFit.SetData(&mGData); // connect data to the fitter
1052
1053// Set and initialise the parameters (that's non-linear fitting!)
1054// (num par, name, guess start, step, [limits min and max])
1055mFit.SetParam(0,"high",90.,1..);
1056mFit.SetParam(1,"xcenter",0.05,0.01);
1057mFit.SetParam(2,"sigma",sig,0.05,0.01,10.);
1058 // Give limits to avoid division by zero
1059mFit.SetParam(3,"constant",0.,1.);
1060
1061// Fit and print result
1062int rcfit = mFit.Fit();
1063mFit.PrintFit();
1064if(rcfit>0) {)
1065 cout<<"Reduce_Chisquare = "<<mFit.GetChi2Red()
1066 <<" nstep="<<mFit.GetNStep()<<" rc="<<rcfit<<endl;
1067} else {
1068 cout<<"Fit_Error, rc = "<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
1069 mFit.PrintFitErr(rcfit);
1070}
1071
1072// Get the result for further use
1073TVector<r_8> ParResult = mFit.GetParm();
1074cout<<ParResult;
1075\end{verbatim}
1076
1077Much more usefull possibilities and detailed informations might be found
1078in the HTML pages of the Sophya manual.
1079
1080\subsection{Polynomial}
1081\index{Polynomial} \index{Poly} \index{Poly2}
1082Polynomials of 1 or 2 variables are supported ({\tt Poly} and {\tt Poly2}).
1083Various operations are supported~:
1084\begin{itemize}
1085\item elementary operations between polynomials $(+,-,*,/) $
1086\item setting or getting coefficients
1087\item computing the value of the polynomial for a given value
1088 of the variable(s),
1089\item derivating
1090\item computing roots (degre 1 or 2)
1091\item fitting the polynomial to vectors of data.
1092\end{itemize}
1093Here is an example of polynomial fitting~:
1094\begin{verbatim}
1095#include "poly.h"
1096// ...
1097Poly pol(2);
1098pol[0] = 100.; pol[1] = 0.; pol[2] = 0.01; // Setting coefficients
1099TVector<r_8> x(100);
1100TVector<r_8> y(100);
1101TVector<r_8> ey(100);
1102for(int i=0;i<100;i++) {
1103 x(i) = i;
1104 ey(i) = 10.;
1105 y(i) = pol((double) i) + ey(i)*NorRand();
1106 ey(i) *= ey(i)
1107}
1108
1109TVector<r_8> errcoef;
1110Poly polfit;
1111polfit.Fit(x,y,ey,2,errcoef);
1112
1113cout<<"Fit Result"<<polfit<<endl;
1114cout<<"Errors :"<<errcoef;
1115\end{verbatim}
1116
1117Similar operations can be done on polynomials with 2 variables.
1118
1119\subsection{Integration, Differential equations}
1120\index{Integration}
1121The NTools module provide also simple classes for numerical integration
1122of functions and differential equations.
1123\begin{figure}[hbt]
1124\dclsbb{Integrator}{GLInteg}
1125\dclsb{TrpzInteg}
1126\end{figure}
1127
1128\index{GLInteg} \index{TrpzInteg}
1129{\bf GLInteg} implements the integration through Gauss-Legendre method
1130and {\bf TrpzInteg} implements trapeze integration. For {\bf TrpzInteg},
1131number of steps specify the number of trapeze, and integration step,
1132their width.
1133The sample code below illustrates the use of TrpzInteg class:
1134\begin{verbatim}
1135#include "integ.h"
1136// ......................................................
1137// Function to be integrated
1138double myf(double x)
1139{
1140// Simple a x + b x^2 (a=2 b=3)
1141return (x*(2.+3.*x));
1142}
1143// ......................................................
1144
1145// Compute Integral(myf, 2., 5.) between xmin=2., xmax=5.
1146TrpzInteg trpz(myf, 2., 5.);
1147// We specify an integration step
1148trpz.DX(0.01);
1149// The integral can be computed as trpz.Value()
1150double myf_integral = trpz.Value();
1151// We could have used the cast operator :
1152cout << "Integral[myf, 2., 5.]= " << (double)trpz << endl;
1153// Limits can be specified through ValueBetween() method
1154cout << "Integral[myf, 0., 4.]= " << trpz.ValueBetween(0.,4.) << endl;
1155\end{verbatim}
1156
1157\subsection{Fourier transform (FFT)}
1158\index{FFT} \index{FFTPackServer}
1159An abstract interface for performing FFT operations is defined by the
1160{\bf FFTServerInterface} class. The {\bf FFTPackSever} class implements
1161one dimensional FFT, on real and complex data. FFTPackServer uses an
1162adapted and extended version of FFTPack (available from netlib),
1163translated in C, and can operate on single and double precision
1164({\tt float, double}) data.
1165
1166The sample code below illustrates the use of FFTServers:
1167\begin{verbatim}
1168#include "fftpserver.h"
1169 // ...
1170TVector<r_8> in(32);
1171TVector< complex<r_8> > out;
1172in = RandomSequence();
1173FFTPackServer ffts;
1174ffts.setNormalize(true); // To have normalized transforms
1175cout << " FFTServer info string= " << ffts.getInfo() << endl;
1176cout << "in= " << in << endl;
1177cout << " Calling ffts.FFTForward(in, out) : " << endl;
1178ffts.FFTForward(in, out);
1179cout << "out= " << out << endl;
1180\end{verbatim}
1181
1182% \newpage
1183\section{Module SUtils}
1184Some utility classes and C/C++ string manipulation functions are gathered
1185in {\bf SUtils} module.
1186\subsection{Using DataCards}
1187\index{DataCards}
1188The {\bf DataCards} class can be used to read parameters from a file.
1189Each line in the file starting with \@ defines a set of values
1190associated with a keyword. In the example below, we read the
1191parameters corresponding with the keyword {\tt SIZE} from the
1192file {\tt ex.d}. We suppose that {\tt ex.d} contains the line: \\
1193{\tt @SIZE 400 250} \\
1194\begin{verbatim}
1195#include "datacards.h"
1196// ...
1197// Initialising DataCards object dc from file ex.d
1198DataCards dc( "ex.d" );
1199// Getting the first and second parameters for keyword size
1200// We define a default value 100
1201int size_x = dc.IParam("SIZE", 0, 100);
1202int size_y = dc.IParam("SIZE", 1, 100);
1203cout << " size_x= " << size_x << " size_y= " << size_y << endl;
1204\end{verbatim}
1205
1206\section{Module SysTools}
1207The {\bf SysTools} module contains classes implementing interface to some
1208OS specific services. The class {\bf ResourceUsage} \index{ResourceUsage}
1209and {\bf Timer} {\index{Timer} provides access to information
1210about various resource usage (memory, CPU, ...).
1211The class {\bf Periodic} provides the necessary services needed to
1212implement the execution of a periodic action.
1213A basic interface to POSIX threads \index{thread} is also provided
1214through the \index{ZThread} {\bf ZThread}, {\bf ZMutex} and {\bf ZSync}
1215classes.
1216
1217\subsection{Dynamic linker}
1218\index{PDynLinkMgr}
1219The class {\bf PDynLinkMgr} can be used for managing shared libraries
1220at run time. The example below shows the run time linking of a function:\\
1221{\tt extern "C" { void myfunc(); } } \\
1222\begin{verbatim}
1223#include "pdlmgr.h"
1224// ...
1225string soname = "mylib.so";
1226string funcname = "myfunc";
1227PDynLinkMgr dyl(soname);
1228DlFunction f = dyl.GetFunction(funcname);
1229if (f != NULL) {
1230// Calling the function
1231 f();
1232}
1233\end{verbatim}
1234
1235\subsection{CxxCompilerLinker class}
1236\index{CxxCompilerLinker}
1237This class provides the services to compile C++ code and building
1238shared libraries, using the same compiler and options which have
1239been used to create the SOPHYA shared library.
1240The sample program below illustrates using this class to build
1241the shared library (myfunc.so) from the source file myfunc.cc :
1242\begin{verbatim}
1243#include "cxxcmplnk.h"
1244// ...
1245string flnm = "myfunc.cc";
1246string oname, soname;
1247int rc;
1248CxxCompilerLinker cxx;
1249// The Compile method provides a default object file name
1250rc = cxx.Compile(flnm, oname);
1251if (rc != 0 ) { // Error when compiling ... }
1252// The BuildSO method provides a default shared object file name
1253rc = cxx.BuildSO(oname, soname);
1254if (rc != 0 ) { // Error when creating shared object ... }
1255\end{verbatim}
1256
1257\newpage
1258\section{Module SkyMap}
1259\begin{figure}[hbt]
1260\dclsbb{AnyDataObj}{PixelMap}
1261\dclsccc{PixelMap}{Sphericalmap}{SphereHEALPix}
1262\dclsc{SphereThetaPhi}
1263\dclsb{LocalMap}
1264\caption{partial class diagram for pixelization classes in Sophya}
1265\end{figure}
1266The {\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.
1267\subsection {Spherical maps}
1268There are two kinds of spherical maps according pixelization algorithms. SphereHEALPix represents spheres pixelized following the HEALPIix algorithm (E. Hivon, K. Gorski)
1269\footnote{see the HEALPix Homepage: http://www.eso.org/kgorski/healpix/ }
1270, 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) :
1271\index{\tcls{SphereHEALPix}}
1272\index{\tcls{SphereThetaPhi}}
1273
1274\begin{verbatim}
1275#include "spherehealpix.h"
1276// ...
1277SphereHEALPix<double> sph(8);
1278for (int k=0; k< sph.NbPixels(); k++) sph(k) = (double)(10*k);
1279\end{verbatim}
1280
1281SphereThetaPhi is used in a similar way with an argument representing number of slices in theta (Euler angle) for an hemisphere.
1282\index{\tcls{SphereThetaPhi}}
1283
1284\subsection {Local maps}
1285\index{\tcls{LocalMap}}
1286A 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).
1287
1288Internally, 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(...))
1289
1290The 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).
1291\begin{verbatim}
1292#include "localmap.h"
1293//..............
1294 LocalMap<r_4> locmap(4,5);
1295 for (int k=0; k<locmap.NbPixels();k++) locmap(k)=10.*k;
1296 locmap.SetOrigin();
1297 locmap.SetSize(30.,30.);
1298\end{verbatim}
1299
1300\subsection{Writing, viewing \dots }
1301
1302All these objects have been design to be written to or read from a persistant file.
1303The following example shows how to write the previously created objects
1304into such a file~:
1305\begin{verbatim}
1306//-- Writing
1307
1308#include "fiospherehealpix.h"
1309//................
1310
1311char *fileout = "myfile.ppf";
1312POutPersist outppf(fileout);
1313FIO_SphereHEALPix<r_8> outsph(sph);
1314outsph.Write(outppf);
1315FIO_LocalMap<r_8> outloc(locmap);
1316outloc.Write(outppf);
1317// It is also possible to use the << operator
1318POutPersist os("sph.ppf");
1319os << outsph;
1320os << outloc;
1321\end{verbatim}
1322
1323Sophya graphical tools (spiapp) can automatically display and operate
1324all these objects.
1325
1326\newpage
1327\section{Module Samba}
1328\index{Spherical Harmonics}
1329\index{SphericalTransformServer}
1330The 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...
1331\begin{verbatim}
1332#include "skymap.h"
1333#include "samba.h"
1334....................
1335
1336// Generate input spectra a + b* l + c * gaussienne(l, 50, 20)
1337int lmax = 92;
1338Vector clin(lmax);
1339for(int l=0; l<lmax; l++) {
1340 double xx = (l-50.)/10.;
1341 clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
1342}
1343
1344// Compute map from spectra
1345SphericalTransformServer<r_8> ylmserver;
1346int m = 128; // HealPix pixelisation parameter
1347SphereHEALPix<r_8> map(m);
1348ylmserver.GenerateFromCl(map, m, clin, 0.);
1349// Compute power spectrum from map
1350Vector clout = ylmserver.DecomposeToCl(map, lmax, 0.);
1351\end{verbatim}
1352
1353\newpage
1354\section{Module SkyT}
1355\index{RadSpectra} \index{SpectralResponse}
1356The SkyT module is composed of two types of classes:
1357\begin{itemize}
1358\item{} one which corresponds to an emission spectrum of
1359radiation, which is called RadSpectra
1360\item{} one which corresponds to the spectral response
1361of a given detector (i.e. corresponding to a detector
1362filter in a given frequency domain), which is called
1363SpectralResponse.
1364\end{itemize}
1365\begin{figure}[hbt]
1366\dclsbb{RadSpectra}{RadSpectraVec}
1367\dclsb{BlackBody}
1368\dclsccc{AnyDataObj}{SpectralResponse}{SpecRespVec}
1369\dclsc{GaussianFilter}
1370\caption{partial class for SkyT module}
1371\end{figure}
1372
1373\begin{verbatim}
1374#include "skyt.h"
1375// ....
1376// Compute the flux from a blackbody at 2.73 K through a square filter
1377BlackBody myBB(2.73);
1378// We define a square filter from 100 - 200 GHz
1379SquareFilter mySF(100,200);
1380// Compute the filtered integrated flux :
1381double flux = myBB.filteredIntegratedFlux(mySF);
1382\end{verbatim}
1383
1384A more detailed description of SkyT module can be found in:
1385{\it The SkyMixer (SkyT and PMixer modules) - Sophya Note No 2. }
1386available also from Sophya Web site.
1387
1388\newpage
1389\section{Module FitsIOServer}
1390\begin{figure}[hbt]
1391\dclsbb{FitsFile}{FitsInFile}
1392\dclsb{FitsOutFile}
1393\end{figure}
1394\index{FITS} \index{FitsInFile} \index{FitsOutFile}
1395This 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 :
1396\begin{verbatim}
1397#include "spherehealpix.h"
1398#include "fitsspherehealpix.h"
1399#include "fitstarray.h"
1400#include "tmatrix.h"
1401//...........................
1402
1403int m=...;
1404SphereHEALPix<r_8> sph(m);
1405................
1406int dim1=...;
1407int dim2=...;
1408TMatrix<r_8> mat(dim1,dim2);
1409............
1410
1411FITS_SphereHEALPix<r_8> sph_temp(sph);
1412FITS_TArray<r_8> mat_temp(mat);
1413// writing
1414
1415FitsOutFile os("myfile.fits");
1416sph_temp.Write(os);
1417mat_temp.Write(os);
1418
1419// reading
1420FitsInFile is("myfile.fits");
1421sph_temp.Read(is);
1422mat_temp.Read(is);
1423SphereHEALPix<r_8> new_sph=(SphereHEALPix<r_8>)sph_temp;
1424TMatrix<r_8> new_mat=(TMatrix<r_8>)mat_temp;
1425................
1426
1427\end{verbatim}
1428
1429The operators {\tt operator << (FitsOutFile ...)} and
1430{\tt operator >> (FitsInFile ...)} are defined in order
1431to facilitate the FITS file operations:
1432\begin{verbatim}
1433// Writing an array object to a FITS file
1434#include "fitstarray.h"
1435FitsOutFile fio("arr.fits");
1436Matrix m(20,30);
1437m = 12345.;
1438fio << m;
1439// .....
1440// Reading a binary table to a XNTuple
1441#include "fitsxntuple.h"
1442XNTuple xn;
1443FitsInFile fii("table.fits");
1444fii >> xn;
1445\end{verbatim}
1446
1447The class {\bf FITS\_AutoReader} provides a limited FITS files reading
1448and decoding capabilities. A partial class diagram of FITS persistence
1449handling classes is shown below:
1450\begin{figure}[hbt]
1451\dclsbb{FitsIOhandler}{FITS\_TArray}
1452\dclsb{FITS\_NTuple}
1453% \dclsb{FITS\_XNTuple}
1454\dclsb{FITS\_SphereHEALPix}
1455% \dclsb{FITS\_LocalMap}
1456\end{figure}
1457
1458\newpage
1459\section{LinAlg and IFFTW modules}
1460An interface to use LAPACK library (available from {\tt http://www.netlib.org})
1461is implemented by the {\bf LapackServer} class, in module LinAlg.
1462\index{LapackServer}.
1463The sample code below shows how to use SVD (Singular Value Decomposition)
1464through LapackServer:
1465\begin{verbatim}
1466#include "intflapack.h"
1467// ...
1468// Use FortranMemoryMapping as default
1469BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
1470// Create an fill the arrays A and its copy AA
1471int n = 20;
1472Matrix A(n , n), AA;
1473A = RandomSequence(RandomSequence::Gaussian, 0., 4.);
1474AA = A; // AA is a copy of A
1475// Compute the SVD decomposition
1476Vector S; // Vector of singular values
1477Matrix U, VT;
1478LapackServer<r_8> lpks;
1479lpks.SVD(AA, S, U, VT);
1480// We create a diagonal matrix using S
1481Matrix SM(n, n);
1482for(int k=0; k<n; k++) SM(k,k) = S(k);
1483// Check the result : A = U*SM*VT
1484Matrix diff = U*(SM*VT) - A;
1485double min, max;
1486diff.MinMax(min, max);
1487cout << " Min/Max difference Matrix (?=0) , Min= " << min
1488 << " Max= " << max << endl;
1489\end{verbatim}
1490
1491\index{FFTWServer}
1492The {\bf FFTWServer} class (in module FFTW) implements FFTServerInterface class
1493methods, for one dimensional and multi-dimensional Fourier
1494transforms on double precision data using the FFTW package
1495(available from {\tt http://www.fftw.org}).
1496
1497\newpage
1498\section{Building and installing Sophya}
1499Presently, the Sophya library has been tested with the following
1500compiler/platform pairs:
1501
1502\begin{center}
1503\begin{tabular}{ll}
1504Compaq/DEC OSF1 & cxx (6.0 , 6.2) \\
1505Compaq/DEC OSF1 & g++ (2.95) \\
1506Linux & g++ (2.91 , 2.95) \\
1507Linux & KCC (3.4) \\
1508Solaris & g++ (2.95) \\
1509SGI IRIX64 & CC (7.3) \\
1510MacOSX/Darwin 10.1 & c++/gcc 2.95 \\
1511\end{tabular}
1512\end{center}
1513
1514Some of the modules in the Sophya package uses external libraries. The
1515{\bf FitsIOServer} is the example of such a module, where the {\tt libcfitsio.a}
1516is used.
1517The build procedure expects to find the include files and the libraries in: \\
1518{\tt \$EXTLIBDIR/Include/FitsIO } \\
1519{\tt \$EXTLIBDIR/`uname`-\$SOPHYACXX/Libs} \\
1520The complete directory tree content for the various external libraries used
1521by this version of SOPHYA is listed below:
1522\begin{verbatim}
1523csh> ls $EXTLIBDIR/Include/
1524FFTW FitsIO XAstro
1525csh> ls $EXTLIBDIR/Include/FFTW/
1526fftw.h rfftw.h
1527csh> ls $EXTLIBDIR/Include/XAstro/
1528P_.h circum.h satlib.h vector.h
1529astro.h deepconst.h satspec.h vsop87.h
1530chap95.h preferences.h sattypes.h
1531csh> ls $EXTLIBDIR/Include/FitsIO/
1532cfortran.h eval_defs.h fitsio.h grparser.h ricecomp.h
1533compress.h eval_tab.h fitsio2.h longnam.h
1534drvrsmem.h f77_wrap.h group.h region.h
1535csh> ls $EXTLIBDIR/Linux-g++/Libs/
1536libcfitsio.a libfftw.la librfftw.a libxastro.a
1537libblas.a libfftw.a liblapack.a librfftw.la
1538\end{verbatim}
1539The object files from a given Sophya module are grouped in an archive library
1540with the module's name ({\tt libmodulename.a}). All Sophya modules
1541 are grouped in a single shared library ({\tt libsophya.so}), while the
1542modules with reference to external libraries are grouped in
1543({\tt libextsophya.so}). The {\bf PI} and {\bf PIext} modules are
1544grouped in ({\tt libPI.so}).
1545
1546The environment variables {\bf SOPHYADEVREP}, {\bf EXTLIBDIR} and {\bf SOPHYACXX}
1547must be defined in order to install the Sophya package.
1548In the example below, we assume that we want to install Sophya from a
1549released (tagged) version in the source directory {\tt \$SRC} in the
1550{\tt /usr/local/Sophya} diretory, using {\tt g++}. We assume that
1551the external libraries directory tree has been set up in
1552{\tt /usr/local/ExtLibs/}. \\[3mm]
1553\centerline{
1554 \rule{20mm}{0.5mm}
1555 {\bf \large the use of GNU make is mandatory}
1556 \rule{20mm}{0.5mm} }
1557
1558\vspace*{3mm}
1559\begin{verbatim}
1560# We select our C++ compiler
1561csh> setenv SOPHYACXX g++
1562# Setup the build directory
1563csh> mkdir /usr/local/Sophya/
1564csh> setenv SOPHYADEVREP /usr/local/Sophya/
1565csh> setenv EXTLIBDIR /usr/local/ExtLibs/
1566# Use the top level makefile in Mgr/
1567csh> cd \$SRC
1568csh> cp Mgr/Makefile Makefile
1569# Step 1: Create the directory tree and copy the include files (.h)
1570csh> make depend
1571# Step 2: Compile the modules without external library reference
1572csh> make libs
1573# Step 3: Compile the modules WITH external library reference (optional)
1574csh> make extlibs
1575# Step 4: Build libsophya.so
1576csh> make slb
1577# Step 5: Build libextsophya.so (optional)
1578csh> make slbext
1579# Step 6: Compile the PI and PIext modules (optional)
1580csh> make PI
1581# Step 7: Build the corresponding shared library libPI.so (optional)
1582csh> make slbpi
1583\end{verbatim}
1584
1585To compile all modules and build the shared libraries, it is possible
1586to use:
1587\begin{verbatim}
1588# Step 2,3,6
1589csh> make all
1590# Step 4,5,7
1591csh> make slball
1592\end{verbatim}
1593
1594At this step, all libraries should have been made. Programs using
1595Sophya libraries can now be built:
1596\begin{verbatim}
1597# To compile test programs
1598csh> cd Tests
1599csh> make arrt ...
1600csh> cd ..
1601# To compile other programs, for example from the PMixer module
1602csh> cd PMixer
1603csh> make
1604csh> cd ..
1605# To build (s)piapp (libPI.so is needed)
1606csh> cd ProgPI
1607csh> make
1608csh> cd ..
1609\end{verbatim}
1610
1611\subsection{Mgr module}
1612This module contains scripts which can be used for generating the
1613makefiles for each module.
1614\begin{itemize}
1615\item {\bf Makefile} Top level Makefile for building the libraries.
1616\item {\bf Makefile.h} contains the definition of compilation flags for the
1617different compilers and systems. This file is used for building the
1618library and generating {\bf MakefileUser.h} (to be included in makefiles).
1619\item {\bf Makefile.slb} contains the rules for building shared libraries
1620for the different compilers and systems. (to be included in makefiles)
1621\item {\bf crerep\_sophya} c-shell script for creating the directory tree
1622under {\tt \$SOPHYABASEREP} and {\tt \$SOPHYADEVREP}
1623\item {\bf install\_sophya} c-shell script for installing the Sophya package.
1624Usually from {\tt \$SOPHYADEVREP} to {\tt \$SOPHYABASEREP}
1625\item {\bf mkmflien} c-shell script for making symbolic links or copying
1626include files to {\tt \$SOPHYADEVREP/Include} or {\tt \$SOPHYABASEREP/Include}
1627\item {\bf mkmf} c-shell script for generating module makefiles and the
1628top level makefile (named GNUmakefile)
1629\item {\bf mkmflib} c-shell script for generating each library module
1630makefile (named GNUmakefile)
1631\item {\bf mkmfprog} c-shell script for generating makefile for a module
1632containing the source for executable programs (named GNUmakefile).
1633\item {\bf mkmfPI} c-shell script for generating makefile for PI and PIext
1634modules (named GNUmakefile)
1635\item {\bf libdirs} List of Sophya modules without reference to external
1636libraries.
1637\item {\bf extlibdirs} List of Sophya modules with reference to external
1638libraries.
1639
1640\end{itemize}
1641
1642\newpage
1643\appendix
1644\section{SOPHYA Exceptions}
1645\index{Exception classes} \index{PThrowable} \index{PError} \index{PException}
1646SOPHYA library defines a set of exceptions which are used
1647for signalling error conditions. The figure below shows a partial
1648class diagram for exception classes in SOPHYA.
1649\begin{figure}[hbt]
1650\dclsbb{PThrowable}{PError}
1651\dclscc{PError}{AllocationError}
1652\dclscc{PError}{NullPtrError}
1653\dclscc{PError}{ForbiddenError}
1654\dclscc{PError}{AssertionFailedError}
1655\dclsbb{PThrowable}{PException}
1656\dclscc{PException}{IOExc}
1657\dclscc{PException}{SzMismatchError}
1658\dclscc{PException}{RangeCheckError}
1659\dclscc{PException}{ParmError}
1660\dclscc{PException}{TypeMismatchExc}
1661\dclscc{PException}{MathExc}
1662\dclscc{PException}{CaughtSignalExc}
1663\caption{partial class diagram for exception handling in Sophya}
1664\end{figure}
1665
1666For simple programs, it is a good practice to handle
1667the exceptions at least at high level, in the {\tt main()} function.
1668The example below shows the exception handling and the usage
1669of Sophya persistence.
1670
1671\input{ex1.inc}
1672
1673
1674\newpage
1675\addcontentsline{toc}{section}{Index}
1676\printindex
1677\end{document}
1678
Note: See TracBrowser for help on using the repository browser.