root/library/bdm/itpp_ext.cpp @ 662

Revision 662, 10.6 kB (checked in by suzdalev, 15 years ago)

mac compilation

  • Property svn:eol-style set to native
Line 
1//
2// C++ Implementation: itpp_ext
3//
4// Description:
5//
6//
7// Author: smidl <smidl@utia.cas.cz>, (C) 2008
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include "itpp_ext.h"
14
15#ifndef M_PI
16#define M_PI        3.14159265358979323846
17#endif
18// from  algebra/lapack.h
19
20extern "C"    /* QR factorization of a general matrix A  */
21{
22        void dgeqrf_ (int *m, int *n, double *a, int *lda, double *tau, double *work,
23                      int *lwork, int *info);
24};
25
26namespace itpp
27{
28Array<int> to_Arr (const ivec &indices)
29{
30        Array<int> a (indices.size());
31        for (int i = 0; i < a.size(); i++) {
32                a (i) = indices (i);
33        }
34        return a;
35}
36
37ivec linspace (int from, int to)
38{
39        int n = to - from + 1;
40        int i;
41        it_assert_debug (n > 0, "wrong linspace");
42        ivec iv (n);
43        for (i = 0; i < n; i++) iv (i) = from + i;
44        return iv;
45};
46
47void set_subvector (vec &ov, const ivec &iv, const vec &v)
48{
49        it_assert_debug ( (iv.length() <= v.length()),
50                          "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
51                          "of range of v");
52        for (int i = 0; i < iv.length(); i++) {
53                it_assert_debug (iv (i) < ov.length(),
54                                 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
55                                 "of range of v");
56                ov (iv (i)) = v (i);
57        }
58}
59
60vec get_vec (const vec &v, const ivec &indexlist)
61{
62        int size = indexlist.size();
63        vec temp (size);
64        for (int i = 0; i < size; ++i) {
65                temp (i) = v._data() [indexlist (i) ];
66        }
67        return temp;
68}
69
70// Gamma
71#define log std::log
72#define exp std::exp
73#define sqrt std::sqrt
74#define R_FINITE std::isfinite
75
76bvec operator> (const vec &t1, const vec &t2)
77{
78        it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
79        bvec temp (t1.length());
80        for (int i = 0; i < t1.length(); i++)
81                temp (i) = (t1[i] > t2[i]);
82        return temp;
83}
84
85bvec operator< (const vec &t1, const vec &t2)
86{
87        it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
88        bvec temp (t1.length());
89        for (int i = 0; i < t1.length(); i++)
90                temp (i) = (t1[i] < t2[i]);
91        return temp;
92}
93
94
95bvec operator& (const bvec &a, const bvec &b)
96{
97        it_assert_debug (b.size() == a.size(), "operator&(): Vectors of different lengths");
98
99        bvec temp (a.size());
100        for (int i = 0; i < a.size(); i++) {
101                temp (i) = a (i) & b (i);
102        }
103        return temp;
104}
105
106bvec operator| (const bvec &a, const bvec &b)
107{
108        it_assert_debug (b.size() != a.size(), "operator&(): Vectors of different lengths");
109
110        bvec temp (a.size());
111        for (int i = 0; i < a.size(); i++) {
112                temp (i) = a (i) | b (i);
113        }
114        return temp;
115}
116
117//#if 0
118Gamma_RNG::Gamma_RNG (double a, double b)
119{
120        setup (a, b);
121}
122double  Gamma_RNG::sample()
123{
124        //A copy of rgamma code from the R package!!
125        //
126
127        /* Constants : */
128        const static double sqrt32 = 5.656854;
129        const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
130
131        /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
132         * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k)
133         * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
134         */
135        const static double q1 = 0.04166669;
136        const static double q2 = 0.02083148;
137        const static double q3 = 0.00801191;
138        const static double q4 = 0.00144121;
139        const static double q5 = -7.388e-5;
140        const static double q6 = 2.4511e-4;
141        const static double q7 = 2.424e-4;
142
143        const static double a1 = 0.3333333;
144        const static double a2 = -0.250003;
145        const static double a3 = 0.2000062;
146        const static double a4 = -0.1662921;
147        const static double a5 = 0.1423657;
148        const static double a6 = -0.1367177;
149        const static double a7 = 0.1233795;
150
151        /* State variables [FIXME for threading!] :*/
152        static double aa = 0.;
153        static double aaa = 0.;
154        static double s, s2, d;    /* no. 1 (step 1) */
155        static double q0, b, si, c;/* no. 2 (step 4) */
156
157        double e, p, q, r, t, u, v, w, x, ret_val;
158        double a = alpha;
159        double scale = 1.0 / beta;
160
161        if (!R_FINITE (a) || !R_FINITE (scale) || a < 0.0 || scale <= 0.0) {
162                it_error ("Gamma_RNG wrong parameters");
163        }
164
165        if (a < 1.) {   /* GS algorithm for parameters a < 1 */
166                if (a == 0)
167                        return 0.;
168                e = 1.0 + exp_m1 * a;
169                for (;;) {    //VS repeat
170                        p = e * unif_rand();
171                        if (p >= 1.0) {
172                                x = -log ( (e - p) / a);
173                                if (exp_rand() >= (1.0 - a) * log (x))
174                                        break;
175                        } else {
176                                x = exp (log (p) / a);
177                                if (exp_rand() >= x)
178                                        break;
179                        }
180                }
181                return scale * x;
182        }
183
184        /* --- a >= 1 : GD algorithm --- */
185
186        /* Step 1: Recalculations of s2, s, d if a has changed */
187        if (a != aa) {
188                aa = a;
189                s2 = a - 0.5;
190                s = sqrt (s2);
191                d = sqrt32 - s * 12.0;
192        }
193        /* Step 2: t = standard normal deviate,
194                   x = (s,1/2) -normal deviate. */
195
196        /* immediate acceptance (i) */
197        t = norm_rand();
198        x = s + 0.5 * t;
199        ret_val = x * x;
200        if (t >= 0.0)
201                return scale * ret_val;
202
203        /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
204        u = unif_rand();
205        if ( (d * u) <= (t * t * t))
206                return scale * ret_val;
207
208        /* Step 4: recalculations of q0, b, si, c if necessary */
209
210        if (a != aaa) {
211                aaa = a;
212                r = 1.0 / a;
213                q0 = ( ( ( ( ( (q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
214                         + q2) * r + q1) * r;
215
216                /* Approximation depending on size of parameter a */
217                /* The constants in the expressions for b, si and c */
218                /* were established by numerical experiments */
219
220                if (a <= 3.686) {
221                        b = 0.463 + s + 0.178 * s2;
222                        si = 1.235;
223                        c = 0.195 / s - 0.079 + 0.16 * s;
224                } else if (a <= 13.022) {
225                        b = 1.654 + 0.0076 * s2;
226                        si = 1.68 / s + 0.275;
227                        c = 0.062 / s + 0.024;
228                } else {
229                        b = 1.77;
230                        si = 0.75;
231                        c = 0.1515 / s;
232                }
233        }
234        /* Step 5: no quotient test if x not positive */
235
236        if (x > 0.0) {
237                /* Step 6: calculation of v and quotient q */
238                v = t / (s + s);
239                if (fabs (v) <= 0.25)
240                        q = q0 + 0.5 * t * t * ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v
241                                                     + a3) * v + a2) * v + a1) * v;
242                else
243                        q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v);
244
245
246                /* Step 7: quotient acceptance (q) */
247                if (log (1.0 - u) <= q)
248                        return scale * ret_val;
249        }
250
251        for (;;) {   //VS repeat
252                /* Step 8: e = standard exponential deviate
253                 *      u =  0,1 -uniform deviate
254                 *      t = (b,si)-double exponential (laplace) sample */
255                e = exp_rand();
256                u = unif_rand();
257                u = u + u - 1.0;
258                if (u < 0.0)
259                        t = b - si * e;
260                else
261                        t = b + si * e;
262                /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */
263                if (t >= -0.71874483771719) {
264                        /* Step 10:      calculation of v and quotient q */
265                        v = t / (s + s);
266                        if (fabs (v) <= 0.25)
267                                q = q0 + 0.5 * t * t *
268                                    ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
269                                        + a2) * v + a1) * v;
270                        else
271                                q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v);
272                        /* Step 11:      hat acceptance (h) */
273                        /* (if q not positive go to step 8) */
274                        if (q > 0.0) {
275                                // TODO: w = expm1(q);
276                                w = exp (q) - 1;
277                                /*  ^^^^^ original code had approximation with rel.err < 2e-7 */
278                                /* if t is rejected sample again at step 8 */
279                                if ( (c * fabs (u)) <= (w * exp (e - 0.5 * t * t)))
280                                        break;
281                        }
282                }
283        } /* repeat .. until  `t' is accepted */
284        x = s + 0.5 * t;
285        return scale * x * x;
286}
287
288
289bool qr (const mat &A, mat &R)
290{
291        int info;
292        int m = A.rows();
293        int n = A.cols();
294        int lwork = n;
295        int k = std::min (m, n);
296        vec tau (k);
297        vec work (lwork);
298
299        R = A;
300
301        // perform workspace query for optimum lwork value
302        int lwork_tmp = -1;
303        dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
304                 &info);
305        if (info == 0) {
306                lwork = static_cast<int> (work (0));
307                work.set_size (lwork, false);
308        }
309        dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
310
311        // construct R
312        for (int i = 0; i < m; i++)
313                for (int j = 0; j < std::min (i, n); j++)
314                        R (i, j) = 0;
315
316        return (info == 0);
317}
318
319//#endif
320std::string num2str (double d)
321{
322        char tmp[20];//that should do
323        sprintf (tmp, "%f", d);
324        return std::string (tmp);
325};
326std::string num2str (int i)
327{
328        char tmp[10];//that should do
329        sprintf (tmp, "%d", i);
330        return std::string (tmp);
331};
332
333// digamma
334// copied from C. Bonds' source
335#include <math.h>
336#define el 0.5772156649015329
337
338double psi (double x)
339{
340        double s, ps, xa, x2;
341        int n, k;
342        static double a[] = {
343                -0.8333333333333e-01,
344                0.83333333333333333e-02,
345                -0.39682539682539683e-02,
346                0.41666666666666667e-02,
347                -0.75757575757575758e-02,
348                0.21092796092796093e-01,
349                -0.83333333333333333e-01,
350                0.4432598039215686
351        };
352
353        xa = fabs (x);
354        s = 0.0;
355        if ( (x == (int) x) && (x <= 0.0)) {
356                ps = 1e308;
357                return ps;
358        }
359        if (xa == (int) xa) {
360                n = xa;
361                for (k = 1; k < n; k++) {
362                        s += 1.0 / k;
363                }
364                ps =  s - el;
365        } else if ( (xa + 0.5) == ( (int) (xa + 0.5))) {
366                n = xa - 0.5;
367                for (k = 1; k <= n; k++) {
368                        s += 1.0 / (2.0 * k - 1.0);
369                }
370                ps = 2.0 * s - el - 1.386294361119891;
371        } else {
372                if (xa < 10.0) {
373                        n = 10 - (int) xa;
374                        for (k = 0; k < n; k++) {
375                                s += 1.0 / (xa + k);
376                        }
377                        xa += n;
378                }
379                x2 = 1.0 / (xa * xa);
380                ps = log (xa) - 0.5 / xa + x2 * ( ( ( ( ( ( (a[7] * x2 + a[6]) * x2 + a[5]) * x2 +
381                                                        a[4]) * x2 + a[3]) * x2 + a[2]) * x2 + a[1]) * x2 + a[0]);
382                ps -= s;
383        }
384        if (x < 0.0)
385                ps = ps - M_PI * std::cos (M_PI * x) / std::sin (M_PI * x) - 1.0 / x;
386        return ps;
387}
388
389void triu (mat &A)
390{
391        for (int i = 1;i < A.rows();i++) { // row cycle
392                for (int j = 0; j < i; j++) {A (i, j) = 0;}
393        }
394}
395
396//! Storage of randun() internals
397class RandunStorage
398{
399                const int A;
400                const int M;
401                static double seed;
402                static int counter;
403        public:
404                RandunStorage() : A (16807), M (2147483647) {};
405                //!set seed of the randun() generator
406                void set_seed (double seed0) {seed = seed0;}
407                //! generate randun() sample
408                double get() {
409                        long long tmp = A * seed;
410                        tmp = tmp % M;
411                        seed = tmp;
412                        counter++;
413                        return seed / M;
414                }
415};
416static RandunStorage randun_global_storage;
417double RandunStorage::seed = 1111111;
418int RandunStorage::counter = 0;
419double randun() {return randun_global_storage.get();};
420vec randun (int n) {vec res (n); for (int i = 0;i < n;i++) {res (i) = randun();}; return res;};
421mat randun (int n, int m) {mat res (n, m); for (int i = 0;i < n*m;i++) {res (i) = randun();}; return res;};
422
423ivec unique (const ivec &in)
424{
425        ivec uniq (0);
426        int j = 0;
427        bool found = false;
428        for (int i = 0;i < in.length(); i++) {
429                found = false;
430                j = 0;
431                while ( (!found) && (j < uniq.length())) {
432                        if (in (i) == uniq (j)) found = true;
433                        j++;
434                }
435                if (!found) uniq = concat (uniq, in (i));
436        }
437        return uniq;
438}
439
440ivec unique_complement (const ivec &in, const ivec &base)
441{
442        // almost a copy of unique
443        ivec uniq (0);
444        int j = 0;
445        bool found = false;
446        for (int i = 0;i < in.length(); i++) {
447                found = false;
448                j = 0;
449                while ( (!found) && (j < uniq.length())) {
450                        if (in (i) == uniq (j)) found = true;
451                        j++;
452                }
453                j=0;
454                while ( (!found) && (j < base.length())) {
455                        if (in (i) == base (j)) found = true;
456                        j++;
457                }
458                if (!found) uniq = concat (uniq, in (i));
459        }
460        return uniq;
461}
462
463}
Note: See TracBrowser for help on using the browser.