Changeset 817

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

Files:

• SVD/matrix_tools.cpp (modified) (3 diffs)
• SVD/matrix_tools.h (modified) (1 diff)
• simil_model.cpp (modified) (7 diffs)
• simil_model.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);
 };
 

cpp/frams/model/similarity/simil_model.cpp

@@ -35 +35 @@
 
 static ParamEntry MSparam_tab[] = {
-                { "Creature: Similarity", 1, 6, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016", },
+                { "Creature: Similarity", 1, 7, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016", },
                 { "simil_parts", 0, 0, "Weight of parts count", "f 0 100 0", FIELD(m_adFactors[0]), "Differing number of parts is also handled by the 'part degree' similarity component.", },
                 { "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
@@ -41 +41 @@
                 { "simil_partgeom", 0, 0, "Weight of parts' geometric distances", "f 0 100 0", FIELD(m_adFactors[3]), "", },
                 { "simil_fixedZaxis", 0, 0, "Fix 'z' (vertical) axis?", "d 0 1 0", FIELD(fixedZaxis), "", },
+                { "simil_weightedMDS", 0, 0, "Should weighted MDS be used?", "d 0 1 0", FIELD(wMDS), "", },
                 { "evaluateDistance", 0, PARAM_DONTSAVE | PARAM_USERHIDDEN, "evaluate model dissimilarity", "p f(oGeno,oGeno)", PROCEDURE(p_evaldistance), "Calculates dissimilarity between two models created from Geno objects.", },
                 { 0, },
@@ -75 +76 @@
         isFuzzy = 0;
         fuzzyDepth = 10;
+        
+        //Determines whether weighted MDS should be used.
+        wMDS = 0;
 }
 
@@ -1981 +1985 @@
                 int nSize = m_Mod[iMod]->getPartCount();
 
-                double *pDistances = (double *)malloc(nSize * nSize * sizeof(double));
+                double *pDistances =  new double[nSize * nSize];
 
                 for (int i = 0; i < nSize; i++)
@@ -1993 +1997 @@
                 Pt3D P1Pos, P2Pos; // positions of parts 
                 double dDistance; // the distance between Parts
+                
+                double *weights = new double[nSize];
+                for (int i = 0; i < nSize; i++)
+                {
+                        if (wMDS==1)
+                                weights[i] = 0;
+                        else
+                                weights[i] = 1;
+                }
+                                
+                if (wMDS==1)
+                        for (int i = 0; i < pModel->getJointCount(); i++)
+                        {
+                                weights[pModel->getJoint(i)->p1_refno]++;
+                                weights[pModel->getJoint(i)->p2_refno]++; 
+                        }
+                
                 for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
                 {
@@ -2020 +2041 @@
                 } // for (iP1)
 
-                MatrixTools::SVD(vEigenvalues, nSize, pDistances, m_aPositions[iMod]);
+                MatrixTools::weightedMDS(vEigenvalues, nSize, pDistances, m_aPositions[iMod], weights);
                 if (fixedZaxis == 1) //restore the original vertical coordinate of each Part
                 {
@@ -2029 +2050 @@
                 }
 
-                free(pDistances);
+                delete[] pDistances;
+                delete[] weights;
         }
 

cpp/frams/model/similarity/simil_model.h

@@ -88 +88 @@
         int fuzzyDepth;
         int isFuzzy;
+        
+        //For wMDS = 1 weighted MDS with vertex degrees as weights is used for the alignment.
+        int wMDS;
 
         /// Interface to local parameters