Changeset 817 for cpp/frams/model/similarity/SVD

Timestamp:

10/07/18 13:23:11 (6 years ago)

Author:

oriona

Message:

MDS procedure replaced with weighted MDS procedure.

Location:

cpp/frams/model/similarity/SVD

Files:

• matrix_tools.cpp (modified) (3 diffs)
• matrix_tools.h (modified) (1 diff)

cpp/frams/model/similarity/SVD/matrix_tools.cpp

@@ -84 +84 @@
 }
 
-double *Transpose(double *&A, int arow, int acol)
+double *Transpose(double *&A, int arow, int acol, double *&toDel, int delSize)
 {
         double *result = Create(acol * arow);
@@ -94 +94 @@
                 }
 
+        if (delSize != 0)
+                delete[] toDel;
+        
         return result;
-
-}
-
-/** Computes the SVD of the nSize x nSize distance matrix
+}
+
+//Weighted centring of a matrix.
+//https://github.com/vegandevs/vegan/blob/master/src/goffactor.c
+void wcentre(double *x, double *w, int *nr, int *nc)
+{
+        int i, j, ij;
+        double sw, swx;
+
+        for (i = 0, sw = 0.0; i < (*nr); i++)
+                sw += w[i];
+
+        for (j = 0; j < (*nc) ; j++)
+        {
+                for (i = 0, swx = 0.0, ij = (*nr)*j; i < (*nr); i++, ij++)
+                {
+                        swx += w[i] * x[ij];
+                }
+                swx /= sw;
+                for (i = 0,  ij = (*nr)*j; i < (*nr); i++, ij++)
+                {
+                        x[ij] -= swx;
+                        x[ij] *= sqrt(w[i]);
+                }
+        }
+}
+
+/** Computes the weighted MDS of the nSize x nSize distance matrix
                 @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the
                 decomposed distance matrix (or rather, to be strict, of the matrix of scalar products
@@ -106 +133 @@
                 @param pDistances [IN] matrix of distances between parts.
                 @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix.
+                @param weights [IN] vector of row weights.
                 */
-void MatrixTools::SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates)
-{
-        // compute squares of elements of this array
-        // compute the matrix B that is the parameter of SVD
-        double *B = Create(nSize * nSize);
-        {
-                // use additional scope to delete temporary matrices
-                double *Ones, *Eye, *Z, *D;
-
-                D = Create(nSize * nSize);
-                D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
-
-                Ones = Create(nSize * nSize);
-                for (int i = 0; i < nSize; i++)
-                        for (int j = 0; j < nSize; j++)
-                        {
-                        Ones[i * nSize + j] = 1;
-                        }
-
-                Eye = Create(nSize * nSize);
-                for (int i = 0; i < nSize; i++)
-                {
-                        for (int j = 0; j < nSize; j++)
-                        {
-                                if (i == j)
-                                {
-                                        Eye[i * nSize + j] = 1;
-                                }
-                                else
-                                {
-                                        Eye[i * nSize + j] = 0;
-                                }
-                        }
-                }
-
-                Z = Create(nSize * nSize);
-                for (int i = 0; i < nSize; i++)
-                {
-                        for (int j = 0; j < nSize; j++)
-                        {
-                                Z[i * nSize + j] = 1.0 / ((double)nSize) * Ones[i * nSize + j];
-                        }
-                }
-
-                for (int i = 0; i < nSize; i++)
-                {
-                        for (int j = 0; j < nSize; j++)
-                        {
-                                Z[i * nSize + j] = Eye[i * nSize + j] - Z[i * nSize + j];
-                        }
-                }
-
-                for (int i = 0; i < nSize; i++)
-                {
-                        for (int j = 0; j < nSize; j++)
-                        {
-                                B[i * nSize + j] = Z[i * nSize + j] * -0.5;
-                        }
-                }
-
-                B = Multiply(B, D, nSize, nSize, nSize, B, nSize);
-                B = Multiply(B, Z, nSize, nSize, nSize, B, nSize);
-
-                delete[] Ones;
-                delete[] Eye;
-                delete[] Z;
-                delete[] D;
-        }
+void MatrixTools::weightedMDS(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates, double *weights)
+{
+        // compute the matrix D that is the parameter of SVD
+        double *D = Create(nSize * nSize);
+        D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
+
+        for (int i = 0; i < 2; i++)
+        {
+                wcentre(D, weights, &nSize, &nSize);
+                D = Transpose(D, nSize, nSize, D, nSize);
+        }
+                
+        for (int i = 0; i < nSize; i++)
+                for (int j = 0; j < nSize; j++)
+                {
+                        D[i * nSize + j] *= -0.5;
+                }
 
         double *Eigenvalues = Create(nSize);
         double *S = Create(nSize * nSize);
 
-        // call SVD function
+        // call the SVD function
         double *Vt = Create(nSize * nSize);
         size_t astep = nSize * sizeof(double);
-        Lapack::JacobiSVD(B, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
-
-        double *W = Transpose(Vt, nSize, nSize);
-
-        delete[] B;
+        Lapack::JacobiSVD(D, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
+
+        double *W = Transpose(Vt, nSize, nSize, W, 0);
+
+        delete[] D;
         delete[] Vt;
+        
+        // deweight
+        double row_weight = 1;
+        for (int i = 0; i < nSize; i++)
+        {
+                row_weight = weights[i];
+                for (int j = 0; j < nSize; j++)
+                {
+                        W[i * nSize + j] /= sqrt(row_weight);
+                }
+        }
 
         for (int i = 0; i < nSize; i++)

cpp/frams/model/similarity/SVD/matrix_tools.h

@@ -12 +12 @@
 {
 public:
-        static void SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates);
+        static void weightedMDS(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates, double *weights);
 };