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

Last change on this file since 2278 was 2278, checked in by ansari, 23 years ago

MAJ sophya.tex pour V1.4 (liste des modules, etc ...) Reza 26/11/2002

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