1 | // Utility classes for template numerical arrays
|
---|
2 | // R. Ansari, C.Magneville 03/2000
|
---|
3 |
|
---|
4 | #include "machdefs.h"
|
---|
5 | #include "utilarr.h"
|
---|
6 | #include "srandgen.h"
|
---|
7 |
|
---|
8 | // Classe utilitaires
|
---|
9 |
|
---|
10 | Sequence::~Sequence()
|
---|
11 | {
|
---|
12 | }
|
---|
13 |
|
---|
14 | //////////////////////////////////////////////////////////
|
---|
15 | /*!
|
---|
16 | \class SOPHYA::RandomSequence
|
---|
17 | \ingroup TArray
|
---|
18 | Base class to generate a sequence of random values
|
---|
19 | */
|
---|
20 |
|
---|
21 | //! Constructor
|
---|
22 | /*!
|
---|
23 | \param typ : generator type
|
---|
24 | \param m : mean parameter of the generator (if needed)
|
---|
25 | \param s : sigma parameter of the generator (if needed)
|
---|
26 | */
|
---|
27 | RandomSequence::RandomSequence(int typ, double m, double s)
|
---|
28 | {
|
---|
29 | typ_ = (typ == Flat) ? Flat : Gaussian;
|
---|
30 | mean_ = m;
|
---|
31 | sig_ = s;
|
---|
32 | }
|
---|
33 | RandomSequence::~RandomSequence()
|
---|
34 | {
|
---|
35 | }
|
---|
36 |
|
---|
37 | //! Return random sequence values.
|
---|
38 | /*!
|
---|
39 | \return If typ = Flat : return [-1,+1]*sig + mean
|
---|
40 | \return If typ = Gaussian : return gaussian distributed
|
---|
41 | with \b mean mean and sigma \b sig
|
---|
42 | */
|
---|
43 | double RandomSequence::Rand()
|
---|
44 | {
|
---|
45 | if (typ_ == Flat)
|
---|
46 | return(drandpm1()*sig_ + mean_);
|
---|
47 | else return(GauRnd(mean_, sig_));
|
---|
48 | }
|
---|
49 |
|
---|
50 | MuTyV & RandomSequence::Value(sa_size_t k) const
|
---|
51 | {
|
---|
52 | if (typ_ == Flat) retv_ = drandpm1()*sig_ + mean_;
|
---|
53 | else retv_ = GauRnd(mean_, sig_);
|
---|
54 | return retv_;
|
---|
55 | }
|
---|
56 |
|
---|
57 |
|
---|
58 | //////////////////////////////////////////////////////////
|
---|
59 | /*!
|
---|
60 | \class SOPHYA::RegularSequence
|
---|
61 | \ingroup TArray
|
---|
62 | Class to generate a sequence of values
|
---|
63 | */
|
---|
64 |
|
---|
65 | //! Constructor
|
---|
66 | /*!
|
---|
67 | \param start : start value
|
---|
68 | \param step : step value
|
---|
69 | \param f : pointer to the sequence function (default = NULL, f(x)=x )
|
---|
70 |
|
---|
71 | See \ref RegularSequenceOperat "operator()"
|
---|
72 | */
|
---|
73 | RegularSequence::RegularSequence(double start, double step, Arr_DoubleFunctionOfX f)
|
---|
74 | {
|
---|
75 | start_ = start;
|
---|
76 | step_ = step;
|
---|
77 | myf_ = f;
|
---|
78 | }
|
---|
79 |
|
---|
80 | RegularSequence::~RegularSequence()
|
---|
81 | {
|
---|
82 | }
|
---|
83 |
|
---|
84 | //! Get the \b k th value
|
---|
85 | /*!
|
---|
86 | \param k : index of the value
|
---|
87 | \anchor RegularSequenceOperat
|
---|
88 |
|
---|
89 | \return f(start+k*step)
|
---|
90 |
|
---|
91 | */
|
---|
92 |
|
---|
93 | MuTyV & RegularSequence::Value (sa_size_t k) const
|
---|
94 | {
|
---|
95 | double x = start_+(double)k*step_;
|
---|
96 | if (myf_) x = myf_(x);
|
---|
97 | retv_ = x;
|
---|
98 | return(retv_);
|
---|
99 | }
|
---|
100 |
|
---|
101 | /*!
|
---|
102 | \class SOPHYA::EnumeratedSequence
|
---|
103 | \ingroup TArray
|
---|
104 | Explicitly defined sequence of values. The comma operator has
|
---|
105 | been redefined to let an easy definition of sequences.
|
---|
106 |
|
---|
107 | \code
|
---|
108 | // Initializing a sequence
|
---|
109 | EnumeratedSequence es;
|
---|
110 | es = 11, 22, 33, 44, 55, 66;
|
---|
111 |
|
---|
112 | for(int k=0; k<8; k++)
|
---|
113 | cout << " k= " << k << " es(k)= " << es(k) << endl;
|
---|
114 |
|
---|
115 | // Decoding a sequence from a string
|
---|
116 | EnumeratedSequence ess;
|
---|
117 | int nbad;
|
---|
118 | ess.Append("56.5 (1.,-1.) 4 8 16", nbad);
|
---|
119 | cout << ess;
|
---|
120 | \endcode
|
---|
121 | */
|
---|
122 |
|
---|
123 | //! Default constructor
|
---|
124 | EnumeratedSequence::EnumeratedSequence()
|
---|
125 | {
|
---|
126 | }
|
---|
127 |
|
---|
128 | //! Copy constructor
|
---|
129 | EnumeratedSequence::EnumeratedSequence(EnumeratedSequence const & es)
|
---|
130 | {
|
---|
131 | Append(es);
|
---|
132 | }
|
---|
133 |
|
---|
134 | EnumeratedSequence::~EnumeratedSequence()
|
---|
135 | {
|
---|
136 | }
|
---|
137 |
|
---|
138 | //! Return the k th value in the sequence (default = 0)
|
---|
139 | MuTyV & EnumeratedSequence::Value (sa_size_t k) const
|
---|
140 | {
|
---|
141 | if (k >= vecv_.size()) retv_ = 0;
|
---|
142 | else retv_ = vecv_[k];
|
---|
143 | return(retv_);
|
---|
144 | }
|
---|
145 |
|
---|
146 | //! Appends a new value to the sequence
|
---|
147 | EnumeratedSequence & EnumeratedSequence::operator , (MuTyV const & v)
|
---|
148 | {
|
---|
149 | vecv_.push_back(v);
|
---|
150 | return(*this);
|
---|
151 | }
|
---|
152 |
|
---|
153 | //! Initialize the sequence with a single value \b v
|
---|
154 | EnumeratedSequence & EnumeratedSequence::operator = (MuTyV const & v)
|
---|
155 | {
|
---|
156 | vecv_.clear();
|
---|
157 | vecv_.push_back(v);
|
---|
158 | return(*this);
|
---|
159 | }
|
---|
160 |
|
---|
161 | //! Copy operator
|
---|
162 | EnumeratedSequence & EnumeratedSequence::operator = (EnumeratedSequence const & seq)
|
---|
163 | {
|
---|
164 | Clear();
|
---|
165 | Append(seq);
|
---|
166 | return(*this);
|
---|
167 | }
|
---|
168 |
|
---|
169 |
|
---|
170 | //! Prints the list to the output stream \b os
|
---|
171 | void EnumeratedSequence::Print(ostream& os) const
|
---|
172 | {
|
---|
173 | os << " EnumeratedSequence::Print() - Size()= " << Size() << endl;
|
---|
174 | for(int k=0; k<vecv_.size(); k++) {
|
---|
175 | os << vecv_[k];
|
---|
176 | if ((k > 0) && (k%10 == 0)) os << endl;
|
---|
177 | else os << " " ;
|
---|
178 | }
|
---|
179 | os << endl;
|
---|
180 | return;
|
---|
181 | }
|
---|
182 |
|
---|
183 | //! Append the \b seq to the end of the sequence.
|
---|
184 | /*!
|
---|
185 | \return the number of added elements.
|
---|
186 | */
|
---|
187 | sa_size_t EnumeratedSequence::Append(EnumeratedSequence const & seq)
|
---|
188 | {
|
---|
189 | for(int k=0; k<seq.vecv_.size(); k++)
|
---|
190 | vecv_.push_back(seq.vecv_[k]);
|
---|
191 | return(seq.vecv_.size());
|
---|
192 | }
|
---|
193 |
|
---|
194 | //! Decodes the string, appending values to the end of the sequence.
|
---|
195 | /*!
|
---|
196 | \param str : string to be decoded
|
---|
197 | \param nbad : number of unmatched quotes or parenthesis (returned value)
|
---|
198 | \param sep : word separator in string. Each word is decoded as a MuTyV value.
|
---|
199 | \return the number of added elements.
|
---|
200 | */
|
---|
201 | sa_size_t EnumeratedSequence::Append(string const & str, int& nbad, const char* sep)
|
---|
202 | {
|
---|
203 | nbad = 0;
|
---|
204 | sa_size_t n = 0;
|
---|
205 | size_t l = str.length();
|
---|
206 | if (l < 1) return(0);
|
---|
207 | // if ((str[0] == '#') || (str[0] == '*')) return(0);
|
---|
208 | size_t q = 0;
|
---|
209 | size_t p = 0;
|
---|
210 | /*
|
---|
211 | size_t p = str.find_first_not_of(sep);
|
---|
212 | if ((str[p] == '+') || (str[p] == '-')) {
|
---|
213 | if (p == l-1) return(0);
|
---|
214 | if (!isdigit(str[p+1])) return(0);
|
---|
215 | }
|
---|
216 | else if (!isdigit(str[p]) && (str[p] != '\'') && (str[p] != '(') ) return(0);
|
---|
217 | */
|
---|
218 | while(q < l) {
|
---|
219 | p = str.find_first_not_of(sep,q);
|
---|
220 | if (p >= l) break;
|
---|
221 | if (str[p] == '\'') { // Decodage d'un string
|
---|
222 | q = str.find('\'',p+1);
|
---|
223 | if (q < l) {
|
---|
224 | vecv_.push_back(MuTyV(str.substr(p+1,q-p-1)));
|
---|
225 | n++; q++;
|
---|
226 | }
|
---|
227 | else nbad++;
|
---|
228 | }
|
---|
229 | else if (str[p] == '(') { // Decodage d'un complex
|
---|
230 | q = str.find(')',p);
|
---|
231 | if (q < l) {
|
---|
232 | q++;
|
---|
233 | MuTyV mtv(str.substr(p,q-p));
|
---|
234 | complex<r_8> z = mtv.operator complex<r_8>();
|
---|
235 | vecv_.push_back(MuTyV(z));
|
---|
236 | n++;
|
---|
237 | }
|
---|
238 | else nbad++;
|
---|
239 | }
|
---|
240 | else {
|
---|
241 | q = str.find_first_of(sep,p);
|
---|
242 | if ( !isdigit(str[p]) && !(str[p] == '+')
|
---|
243 | && !(str[p] == '-') && !(str[p] == '.') ) { // une chaine
|
---|
244 | vecv_.push_back(MuTyV(str.substr(p,q-p)));
|
---|
245 | n++;
|
---|
246 | }
|
---|
247 | else { // C'est un nombre
|
---|
248 | if (str.find_first_of(".eE",p) < q) { // c'est un flottant
|
---|
249 | r_8 x = atof(str.substr(p,q-p).c_str());
|
---|
250 | vecv_.push_back(MuTyV(x));
|
---|
251 | }
|
---|
252 | else { // un entier
|
---|
253 | int_8 l = atol(str.substr(p,q-p).c_str());
|
---|
254 | vecv_.push_back(MuTyV(l));
|
---|
255 | }
|
---|
256 | n++;
|
---|
257 | }
|
---|
258 | }
|
---|
259 | }
|
---|
260 | return (n);
|
---|
261 | }
|
---|
262 |
|
---|
263 | //! Decodes the input ASCII stream, creating a sequence of values
|
---|
264 | /*! \param is : Input ASCII stream
|
---|
265 | \param nr : Number of non empty (or comment) lines in stream (return value)
|
---|
266 | \param nc : Number of columns (= ntot/nlines) (return value)
|
---|
267 | \param clm : Lines starting with clm character are treated as comment lines
|
---|
268 | \param sep : word separator in lines
|
---|
269 | \return Number of decoded elements
|
---|
270 | */
|
---|
271 | sa_size_t EnumeratedSequence::FillFromFile(istream& is, sa_size_t& nr, sa_size_t& nc,
|
---|
272 | char clm, const char* sep)
|
---|
273 | {
|
---|
274 | nr = 0;
|
---|
275 | nc = 0;
|
---|
276 | sa_size_t n = 0;
|
---|
277 | char buff[256];
|
---|
278 | string line;
|
---|
279 | int nbad, nbadtot, nel;
|
---|
280 | nbadtot = nbad = 0;
|
---|
281 | while (!is.eof()) {
|
---|
282 | is.clear();
|
---|
283 | is.getline(buff, 256);
|
---|
284 | // cout << " DBG : buff=" << buff << " :state=" << is.rdstate() << endl;
|
---|
285 | line += buff; nel = 0;
|
---|
286 | if (is.good()) {
|
---|
287 | if ((line.length() > 0) && (line[0]!=clm)) {
|
---|
288 | nel = Append(line, nbad);
|
---|
289 | if (nel > 0) {
|
---|
290 | nr++; n += nel;
|
---|
291 | }
|
---|
292 | nbadtot += nbad;
|
---|
293 | }
|
---|
294 | // cout << " Decoding line = " << line << " Nel= " << nel << endl;
|
---|
295 | line = "";
|
---|
296 | }
|
---|
297 | }
|
---|
298 | if ((line.length() > 0) && (line[0]!=clm)) {
|
---|
299 | nel = Append(line, nbad);
|
---|
300 | // cout << " Decoding Eline = " << line << " Nel= " << nel << endl;
|
---|
301 | if (nel > 0) {
|
---|
302 | nr++; n += nel;
|
---|
303 | }
|
---|
304 | nbadtot += nbad;
|
---|
305 | line = "";
|
---|
306 | }
|
---|
307 | if (nbadtot > 0)
|
---|
308 | cout << "EnumeratedSequence::FillFromFile()/Warning " << nbadtot
|
---|
309 | << " bad match (quotes or parenthesis) in stream " << endl;
|
---|
310 | nc = n/nr;
|
---|
311 | return (n);
|
---|
312 | }
|
---|
313 |
|
---|
314 | //////////////////////////////////////////////////////////
|
---|
315 | /*!
|
---|
316 | \class SOPHYA::Range
|
---|
317 | \ingroup TArray
|
---|
318 | Class to define a range of indexes
|
---|
319 | */
|
---|
320 |
|
---|
321 | //! Constructor
|
---|
322 | /*!
|
---|
323 | Define a range of indexes
|
---|
324 | \param start : start index (inclusive)
|
---|
325 | \param end : end index (inclusive)
|
---|
326 | \param size : size (number of elements, used if end \<= start)
|
---|
327 | \param step : step (or stride)
|
---|
328 |
|
---|
329 | \warning If \b end \> \b start, \b size is computed automatically
|
---|
330 | \warning If not \b size is fixed and \b end recomputed
|
---|
331 | */
|
---|
332 | Range::Range(sa_size_t start, sa_size_t end, sa_size_t size, sa_size_t step)
|
---|
333 | {
|
---|
334 | start_ = start;
|
---|
335 | step_ = (step > 0) ? step : 1;
|
---|
336 | if (end > start) { // Taille calcule automatiquement
|
---|
337 | end_ = end;
|
---|
338 | if (step_ > ((end_-start_)+1)) size_ = 1;
|
---|
339 | else size_ = ((end-start)+1)/step_;
|
---|
340 | }
|
---|
341 | else { // Taille fixee
|
---|
342 | size_ = size;
|
---|
343 | end_ = start_+size_*step_;
|
---|
344 | }
|
---|
345 | }
|
---|
346 |
|
---|
347 | /*
|
---|
348 | Range & Range::operator = (sa_size_t start)
|
---|
349 | {
|
---|
350 | start_ = start;
|
---|
351 | size_ = 1;
|
---|
352 | step_ = 1;
|
---|
353 | return (*this);
|
---|
354 | }
|
---|
355 | */
|
---|
356 |
|
---|
357 |
|
---|
358 | //////////////////////////////////////////////////////////
|
---|
359 | /*!
|
---|
360 | \class SOPHYA::IdentityMatrix
|
---|
361 | \ingroup TArray
|
---|
362 | Class to define an identity matrix
|
---|
363 | */
|
---|
364 |
|
---|
365 | //! Constructor of a (n,n) diagonal matrix with value diag on the diagonal
|
---|
366 | IdentityMatrix::IdentityMatrix(double diag, sa_size_t n)
|
---|
367 | {
|
---|
368 | size_ = n;
|
---|
369 | diag_ = diag;
|
---|
370 | }
|
---|