source: cpp/frams/model/similarity/SVD/matrix_tools.cpp @ 872

Last change on this file since 872 was 817, checked in by oriona, 5 years ago

MDS procedure replaced with weighted MDS procedure.

  • Property svn:eol-style set to native
File size: 4.7 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#include "matrix_tools.h"
7#include "lapack.h"
8#include <cstdlib>
9#include <cmath>
10#include <cstdio>
11#include <stdlib.h> //malloc(), embarcadero
12#include <math.h> //sqrt(), embarcadero
13
14
15double *Create(int nSize)
16{
17        double *matrix = new double[nSize];
18
19        for (int i = 0; i < nSize; i++)
20        {
21                matrix[i] = 0;
22        }
23
24        return matrix;
25}
26
27double *Multiply(double *&a, double *&b, int nrow, int ncol, double ncol2, double *&toDel, int delSize)
28{
29        double *c = Create(nrow * ncol2);
30        int i = 0, j = 0, k = 0;
31
32        for (i = 0; i < nrow; i++)
33        {
34                for (j = 0; j < ncol2; j++)
35                {
36                        for (k = 0; k < ncol; k++)
37                                c[i * nrow + j] += a[i * nrow + k] * b[k * ncol + j];
38                }
39        }
40
41        if (delSize != 0)
42                delete[] toDel;
43        return c;
44}
45
46double *Power(double *&array, int nrow, int ncol, double pow, double *&toDel, int delSize)
47{
48        double *m_Power = Create(nrow * ncol);
49        if (pow == 2)
50        {
51                for (int i = 0; i < nrow; i++)
52                {
53                        for (int j = 0; j < ncol; j++)
54                        {
55                                m_Power[i * nrow + j] = array[i * nrow + j] * array[i * nrow + j];
56                        }
57
58                }
59        }
60        else
61        {
62                for (int i = 0; i < nrow; i++)
63                {
64                        for (int j = 0; j < ncol; j++)
65                        {
66                                m_Power[i * nrow + j] = sqrt(array[i * nrow + j]);
67                        }
68
69                }
70        }
71
72        if (delSize != 0)
73                delete[] toDel;
74
75        return m_Power;
76}
77
78void Print(double *&mat, int nelems)
79{
80        for (int i = 0; i < nelems; i++)
81                printf("%6.2f ", mat[i]);
82        printf("\n");
83
84}
85
86double *Transpose(double *&A, int arow, int acol, double *&toDel, int delSize)
87{
88        double *result = Create(acol * arow);
89
90        for (int i = 0; i < acol; i++)
91                for (int j = 0; j < arow; j++)
92                {
93                result[i * arow + j] = A[j * acol + i];
94                }
95
96        if (delSize != 0)
97                delete[] toDel;
98       
99        return result;
100}
101
102//Weighted centring of a matrix.
103//https://github.com/vegandevs/vegan/blob/master/src/goffactor.c
104void wcentre(double *x, double *w, int *nr, int *nc)
105{
106        int i, j, ij;
107        double sw, swx;
108
109        for (i = 0, sw = 0.0; i < (*nr); i++)
110                sw += w[i];
111
112        for (j = 0; j < (*nc) ; j++)
113        {
114                for (i = 0, swx = 0.0, ij = (*nr)*j; i < (*nr); i++, ij++)
115                {
116                        swx += w[i] * x[ij];
117                }
118                swx /= sw;
119                for (i = 0,  ij = (*nr)*j; i < (*nr); i++, ij++)
120                {
121                        x[ij] -= swx;
122                        x[ij] *= sqrt(w[i]);
123                }
124        }
125}
126
127/** Computes the weighted MDS of the nSize x nSize distance matrix
128                @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the
129                decomposed distance matrix (or rather, to be strict, of the matrix of scalar products
130                created from the matrix of distances). The vector is assumed to be empty before the function call and
131                all variance percentages are pushed at the end of it.
132                @param nSize size of the matrix of distances.
133                @param pDistances [IN] matrix of distances between parts.
134                @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix.
135                @param weights [IN] vector of row weights.
136                */
137void MatrixTools::weightedMDS(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates, double *weights)
138{
139        // compute the matrix D that is the parameter of SVD
140        double *D = Create(nSize * nSize);
141        D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
142
143        for (int i = 0; i < 2; i++)
144        {
145                wcentre(D, weights, &nSize, &nSize);
146                D = Transpose(D, nSize, nSize, D, nSize);
147        }
148               
149        for (int i = 0; i < nSize; i++)
150                for (int j = 0; j < nSize; j++)
151                {
152                        D[i * nSize + j] *= -0.5;
153                }
154
155        double *Eigenvalues = Create(nSize);
156        double *S = Create(nSize * nSize);
157
158        // call the SVD function
159        double *Vt = Create(nSize * nSize);
160        size_t astep = nSize * sizeof(double);
161        Lapack::JacobiSVD(D, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
162
163        double *W = Transpose(Vt, nSize, nSize, W, 0);
164
165        delete[] D;
166        delete[] Vt;
167       
168        // deweight
169        double row_weight = 1;
170        for (int i = 0; i < nSize; i++)
171        {
172                row_weight = weights[i];
173                for (int j = 0; j < nSize; j++)
174                {
175                        W[i * nSize + j] /= sqrt(row_weight);
176                }
177        }
178
179        for (int i = 0; i < nSize; i++)
180                for (int j = 0; j < nSize; j++)
181                {
182                if (i == j)
183                        S[i * nSize + j] = Eigenvalues[i];
184                else
185                        S[i * nSize + j] = 0;
186                }
187
188        // compute coordinates of points
189        double *sqS, *dCoordinates;
190        sqS = Power(S, nSize, nSize, 0.5, S, nSize);
191        dCoordinates = Multiply(W, sqS, nSize, nSize, nSize, W, nSize);
192        delete[] sqS;
193
194        for (int i = 0; i < nSize; i++)
195        {
196                // set coordinate from the SVD solution
197                Coordinates[i].x = dCoordinates[i * nSize];
198                Coordinates[i].y = dCoordinates[i * nSize + 1];
199                if (nSize > 2)
200                        Coordinates[i].z = dCoordinates[i * nSize + 2];
201                else
202                        Coordinates[i].z = 0;
203        }
204
205        // store the eigenvalues in the output vector
206        for (int i = 0; i < nSize; i++)
207        {
208                double dElement = Eigenvalues[i];
209                vdEigenvalues.push_back(dElement);
210        }
211
212        delete[] Eigenvalues;
213        delete[] dCoordinates;
214}
Note: See TracBrowser for help on using the repository browser.