source: cpp/frams/model/similarity/SVD/lapack.cpp @ 468

Last change on this file since 468 was 389, checked in by Maciej Komosinski, 10 years ago

Code formatting

  • Property svn:eol-style set to native
File size: 10.6 KB
Line 
1// This file is a part of Framsticks SDK.  http://www.framsticks.com/
2// Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
3// See LICENSE.txt for details.
4
5////////////////////////////////////////////////////////////////////////////////////////
6//
7//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
8//
9//  By downloading, copying, installing or using the software you agree to this license.
10//  If you do not agree to this license, do not download, install,
11//  copy or use the software.
12//
13//
14//                           License Agreement
15//                For Open Source Computer Vision Library
16//
17// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
18// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
19// Third party copyrights are property of their respective owners.
20//
21// Redistribution and use in source and binary forms, with or without modification,
22// are permitted provided that the following conditions are met:
23//
24//   * Redistribution's of source code must retain the above copyright notice,
25//     this list of conditions and the following disclaimer.
26//
27//   * Redistribution's in binary form must reproduce the above copyright notice,
28//     this list of conditions and the following disclaimer in the documentation
29//     and/or other materials provided with the distribution.
30//
31//   * The name of the copyright holders may not be used to endorse or promote products
32//     derived from this software without specific prior written permission.
33//
34// This software is provided by the copyright holders and contributors "as is" and
35// any express or implied warranties, including, but not limited to, the implied
36// warranties of merchantability and fitness for a particular purpose are disclaimed.
37// In no event shall the Intel Corporation or contributors be liable for any direct,
38// indirect, incidental, special, exemplary, or consequential damages
39// (including, but not limited to, procurement of substitute goods or services;
40// loss of use, data, or profits; or business interruption) however caused
41// and on any theory of liability, whether in contract, strict liability,
42// or tort (including negligence or otherwise) arising in any way out of
43// the use of this software, even if advised of the possibility of such damage.
44//
45//
46/*
47 * AutoBuffer from https://github.com/Itseez/opencv/blob/master/modules/core/include/opencv2/core/utility.hpp
48 * VBLAS, JacobiSVDImpl_ from https://github.com/Itseez/opencv/blob/master/modules/core/src/lapack.cpp
49 * changes MK, May 2015:
50 * - "RNG rng(0x12345678)" and "rng.next()" replaced with a local PRBS-7 implementation so that this file does not depend on external random number generators.
51 */
52#include <cstdlib>
53#include <algorithm>
54#include <cmath>
55#include <limits>
56#include <cfloat>
57//#include <assert.h>
58#include <math.h> //hypot(), embarcadero
59#include <stdint.h> //uint8_t
60#include "lapack.h"
61
62
63#if defined _M_IX86 && defined _MSC_VER && _MSC_VER < 1700
64#pragma float_control(precise, on)
65#endif
66
67template<typename _Tp, size_t fixed_size = 1024 / sizeof(_Tp) + 8 > class AutoBuffer
68{
69public:
70        typedef _Tp value_type;
71
72        //! the default constructor
73        AutoBuffer();
74        //! constructor taking the real buffer size
75        AutoBuffer(size_t _size);
76
77        //! the copy constructor
78        AutoBuffer(const AutoBuffer<_Tp, fixed_size>& buf);
79        //! the assignment operator
80        AutoBuffer<_Tp, fixed_size>& operator=(const AutoBuffer<_Tp, fixed_size>& buf);
81
82        //! destructor. calls deallocate()
83        ~AutoBuffer();
84
85        //! allocates the new buffer of size _size. if the _size is small enough, stack-allocated buffer is used
86        void allocate(size_t _size);
87        //! deallocates the buffer if it was dynamically allocated
88        void deallocate();
89        //! resizes the buffer and preserves the content
90        void resize(size_t _size);
91        //! returns the current buffer size
92        size_t size() const;
93        //! returns pointer to the real buffer, stack-allocated or head-allocated
94        operator _Tp* ();
95        //! returns read-only pointer to the real buffer, stack-allocated or head-allocated
96        operator const _Tp* () const;
97
98protected:
99        //! pointer to the real buffer, can point to buf if the buffer is small enough
100        _Tp* ptr;
101        //! size of the real buffer
102        size_t sz;
103        //! pre-allocated buffer. At least 1 element to confirm C++ standard reqirements
104        _Tp buf[(fixed_size > 0) ? fixed_size : 1];
105};
106
107/////////////////////////////// AutoBuffer implementation ////////////////////////////////////////
108
109template<typename _Tp, size_t fixed_size> inline
110AutoBuffer<_Tp, fixed_size>::AutoBuffer()
111{
112        ptr = buf;
113        sz = fixed_size;
114}
115
116template<typename _Tp, size_t fixed_size> inline
117AutoBuffer<_Tp, fixed_size>::AutoBuffer(size_t _size)
118{
119        ptr = buf;
120        sz = fixed_size;
121        allocate(_size);
122}
123
124template<typename _Tp, size_t fixed_size> inline
125AutoBuffer<_Tp, fixed_size>::AutoBuffer(const AutoBuffer<_Tp, fixed_size>& abuf)
126{
127        ptr = buf;
128        sz = fixed_size;
129        allocate(abuf.size());
130        for (size_t i = 0; i < sz; i++)
131                ptr[i] = abuf.ptr[i];
132}
133
134template<typename _Tp, size_t fixed_size> inline AutoBuffer<_Tp, fixed_size>&
135AutoBuffer<_Tp, fixed_size>::operator=(const AutoBuffer<_Tp, fixed_size>& abuf)
136{
137        if (this != &abuf)
138        {
139                deallocate();
140                allocate(abuf.size());
141                for (size_t i = 0; i < sz; i++)
142                        ptr[i] = abuf.ptr[i];
143        }
144        return *this;
145}
146
147template<typename _Tp, size_t fixed_size> inline
148AutoBuffer<_Tp, fixed_size>::~AutoBuffer()
149{
150        deallocate();
151}
152
153template<typename _Tp, size_t fixed_size> inline void
154AutoBuffer<_Tp, fixed_size>::allocate(size_t _size)
155{
156        if (_size <= sz)
157        {
158                sz = _size;
159                return;
160        }
161        deallocate();
162        if (_size > fixed_size)
163        {
164                ptr = new _Tp[_size];
165                sz = _size;
166        }
167}
168
169template<typename _Tp, size_t fixed_size> inline void
170AutoBuffer<_Tp, fixed_size>::deallocate()
171{
172        if (ptr != buf)
173        {
174                delete[] ptr;
175                ptr = buf;
176                sz = fixed_size;
177        }
178}
179
180template<typename _Tp, size_t fixed_size> inline void
181AutoBuffer<_Tp, fixed_size>::resize(size_t _size)
182{
183        if (_size <= sz)
184        {
185                sz = _size;
186                return;
187        }
188        //size_t i, prevsize = sz, minsize = MIN(prevsize, _size);
189        size_t i, prevsize = sz, minsize = prevsize < _size ? prevsize : _size;
190        _Tp* prevptr = ptr;
191
192        ptr = _size > fixed_size ? new _Tp[_size] : buf;
193        sz = _size;
194
195        if (ptr != prevptr)
196                for (i = 0; i < minsize; i++)
197                        ptr[i] = prevptr[i];
198        for (i = prevsize; i < _size; i++)
199                ptr[i] = _Tp();
200
201        if (prevptr != buf)
202                delete[] prevptr;
203}
204
205template<typename _Tp, size_t fixed_size> inline size_t
206AutoBuffer<_Tp, fixed_size>::size() const
207{
208        return sz;
209}
210
211template<typename _Tp, size_t fixed_size> inline
212AutoBuffer<_Tp, fixed_size>::operator _Tp* ()
213{
214        return ptr;
215}
216
217template<typename _Tp, size_t fixed_size> inline
218AutoBuffer<_Tp, fixed_size>::operator const _Tp* () const
219{
220        return ptr;
221}
222
223template<typename T> struct VBLAS
224{
225
226        int dot(const T*, const T*, int, T*) const
227        {
228                return 0;
229        }
230
231        int givens(T*, T*, int, T, T) const
232        {
233                return 0;
234        }
235
236        int givensx(T*, T*, int, T, T, T*, T*) const
237        {
238                return 0;
239        }
240};
241
242template<typename _Tp> void
243JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep,
244int m, int n, int n1, double minval, _Tp eps)
245{
246        VBLAS<_Tp> vblas;
247        AutoBuffer<double> Wbuf(n);
248        double* W = Wbuf;
249        int i, j, k, iter, max_iter = std::max(m, 30);
250        _Tp c, s;
251        double sd;
252        astep /= sizeof(At[0]);
253        vstep /= sizeof(Vt[0]);
254
255        for (i = 0; i < n; i++)
256        {
257                for (k = 0, sd = 0; k < m; k++)
258                {
259                        _Tp t = At[i * astep + k];
260                        sd += (double)t*t;
261                }
262                W[i] = sd;
263
264                if (Vt)
265                {
266                        for (k = 0; k < n; k++)
267                                Vt[i * vstep + k] = 0;
268                        Vt[i * vstep + i] = 1;
269                }
270        }
271
272        for (iter = 0; iter < max_iter; iter++)
273        {
274                bool changed = false;
275
276                for (i = 0; i < n - 1; i++)
277                        for (j = i + 1; j < n; j++)
278                        {
279                        _Tp *Ai = At + i*astep, *Aj = At + j*astep;
280                        double a = W[i], p = 0, b = W[j];
281
282                        for (k = 0; k < m; k++)
283                                p += (double)Ai[k] * Aj[k];
284
285                        if (std::abs(p) <= eps * std::sqrt((double)a * b))
286                                continue;
287
288                        p *= 2;
289                        double beta = a - b, gamma = hypot((double)p, beta);
290                        if (beta < 0)
291                        {
292                                double delta = (gamma - beta)*0.5;
293                                s = (_Tp)std::sqrt(delta / gamma);
294                                c = (_Tp)(p / (gamma * s * 2));
295                        }
296                        else
297                        {
298                                c = (_Tp)std::sqrt((gamma + beta) / (gamma * 2));
299                                s = (_Tp)(p / (gamma * c * 2));
300                        }
301
302                        a = b = 0;
303                        for (k = 0; k < m; k++)
304                        {
305                                _Tp t0 = c * Ai[k] + s * Aj[k];
306                                _Tp t1 = -s * Ai[k] + c * Aj[k];
307                                Ai[k] = t0;
308                                Aj[k] = t1;
309
310                                a += (double)t0*t0;
311                                b += (double)t1*t1;
312                        }
313                        W[i] = a;
314                        W[j] = b;
315
316                        changed = true;
317
318                        if (Vt)
319                        {
320                                _Tp *Vi = Vt + i*vstep, *Vj = Vt + j*vstep;
321                                k = vblas.givens(Vi, Vj, n, c, s);
322
323                                for (; k < n; k++)
324                                {
325                                        _Tp t0 = c * Vi[k] + s * Vj[k];
326                                        _Tp t1 = -s * Vi[k] + c * Vj[k];
327                                        Vi[k] = t0;
328                                        Vj[k] = t1;
329                                }
330                        }
331                        }
332                if (!changed)
333                        break;
334        }
335
336        for (i = 0; i < n; i++)
337        {
338                for (k = 0, sd = 0; k < m; k++)
339                {
340                        _Tp t = At[i * astep + k];
341                        sd += (double)t*t;
342                }
343                W[i] = std::sqrt(sd);
344        }
345
346        for (i = 0; i < n - 1; i++)
347        {
348                j = i;
349                for (k = i + 1; k < n; k++)
350                {
351                        if (W[j] < W[k])
352                                j = k;
353                }
354                if (i != j)
355                {
356                        std::swap(W[i], W[j]);
357                        if (Vt)
358                        {
359                                for (k = 0; k < m; k++)
360                                        std::swap(At[i * astep + k], At[j * astep + k]);
361
362                                for (k = 0; k < n; k++)
363                                        std::swap(Vt[i * vstep + k], Vt[j * vstep + k]);
364                        }
365                }
366        }
367
368        for (i = 0; i < n; i++)
369                _W[i] = (_Tp)W[i];
370
371        if (!Vt)
372                return;
373
374        uint8_t rndstate = 0x02; //PRBS-7 from http://en.wikipedia.org/wiki/Pseudorandom_binary_sequence
375        for (i = 0; i < n1; i++)
376        {
377                sd = i < n ? W[i] : 0;
378
379                while (sd <= minval)
380                {
381                        // if we got a zero singular value, then in order to get the corresponding left singular vector
382                        // we generate a random vector, project it to the previously computed left singular vectors,
383                        // subtract the projection and normalize the difference.
384                        const _Tp val0 = (_Tp)(1. / m);
385                        for (k = 0; k < m; k++)
386                        {
387                                int rndbit = (((rndstate >> 6) ^ (rndstate >> 5)) & 1);
388                                rndstate = ((rndstate << 1) | rndbit) & 0x7f;
389                                _Tp val = rndbit == 0 ? val0 : -val0;
390                                At[i * astep + k] = val;
391                        }
392                        for (iter = 0; iter < 2; iter++)
393                        {
394                                for (j = 0; j < i; j++)
395                                {
396                                        sd = 0;
397                                        for (k = 0; k < m; k++)
398                                                sd += At[i * astep + k] * At[j * astep + k];
399                                        _Tp asum = 0;
400                                        for (k = 0; k < m; k++)
401                                        {
402                                                _Tp t = (_Tp)(At[i * astep + k] - sd * At[j * astep + k]);
403                                                At[i * astep + k] = t;
404                                                asum += std::abs(t);
405                                        }
406                                        asum = asum ? 1 / asum : 0;
407                                        for (k = 0; k < m; k++)
408                                                At[i * astep + k] *= asum;
409                                }
410                        }
411                        sd = 0;
412                        for (k = 0; k < m; k++)
413                        {
414                                _Tp t = At[i * astep + k];
415                                sd += (double)t*t;
416                        }
417                        sd = std::sqrt(sd);
418                }
419
420                s = (_Tp)(1 / sd);
421                for (k = 0; k < m; k++)
422                        At[i * astep + k] *= s;
423        }
424}
425
426void Lapack::JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1 = -1)
427{
428        JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, DBL_MIN, DBL_EPSILON * 10);
429}
Note: See TracBrowser for help on using the repository browser.