Changeset 817 for cpp


Ignore:
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:
4 edited

Legend:

Unmodified
Added
Removed

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

    r455 r817  
    8484}
    8585
    86 double *Transpose(double *&A, int arow, int acol)
     86double *Transpose(double *&A, int arow, int acol, double *&toDel, int delSize)
    8787{
    8888        double *result = Create(acol * arow);
    … …  
    9494                }
    9595
     96        if (delSize != 0)
     97                delete[] toDel;
     98       
    9699        return result;
    97 
    98 }
    99 
    100 /** Computes the SVD of the nSize x nSize distance matrix
     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
    101128                @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the
    102129                decomposed distance matrix (or rather, to be strict, of the matrix of scalar products
    … …  
    106133                @param pDistances [IN] matrix of distances between parts.
    107134                @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix.
     135                @param weights [IN] vector of row weights.
    108136                */
    109 void MatrixTools::SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates)
    110 {
    111         // compute squares of elements of this array
    112         // compute the matrix B that is the parameter of SVD
    113         double *B = Create(nSize * nSize);
    114         {
    115                 // use additional scope to delete temporary matrices
    116                 double *Ones, *Eye, *Z, *D;
    117 
    118                 D = Create(nSize * nSize);
    119                 D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
    120 
    121                 Ones = Create(nSize * nSize);
    122                 for (int i = 0; i < nSize; i++)
    123                         for (int j = 0; j < nSize; j++)
    124                         {
    125                         Ones[i * nSize + j] = 1;
    126                         }
    127 
    128                 Eye = Create(nSize * nSize);
    129                 for (int i = 0; i < nSize; i++)
    130                 {
    131                         for (int j = 0; j < nSize; j++)
    132                         {
    133                                 if (i == j)
    134                                 {
    135                                         Eye[i * nSize + j] = 1;
    136                                 }
    137                                 else
    138                                 {
    139                                         Eye[i * nSize + j] = 0;
    140                                 }
    141                         }
    142                 }
    143 
    144                 Z = Create(nSize * nSize);
    145                 for (int i = 0; i < nSize; i++)
    146                 {
    147                         for (int j = 0; j < nSize; j++)
    148                         {
    149                                 Z[i * nSize + j] = 1.0 / ((double)nSize) * Ones[i * nSize + j];
    150                         }
    151                 }
    152 
    153                 for (int i = 0; i < nSize; i++)
    154                 {
    155                         for (int j = 0; j < nSize; j++)
    156                         {
    157                                 Z[i * nSize + j] = Eye[i * nSize + j] - Z[i * nSize + j];
    158                         }
    159                 }
    160 
    161                 for (int i = 0; i < nSize; i++)
    162                 {
    163                         for (int j = 0; j < nSize; j++)
    164                         {
    165                                 B[i * nSize + j] = Z[i * nSize + j] * -0.5;
    166                         }
    167                 }
    168 
    169                 B = Multiply(B, D, nSize, nSize, nSize, B, nSize);
    170                 B = Multiply(B, Z, nSize, nSize, nSize, B, nSize);
    171 
    172                 delete[] Ones;
    173                 delete[] Eye;
    174                 delete[] Z;
    175                 delete[] D;
    176         }
     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                }
    177154
    178155        double *Eigenvalues = Create(nSize);
    179156        double *S = Create(nSize * nSize);
    180157
    181         // call SVD function
     158        // call the SVD function
    182159        double *Vt = Create(nSize * nSize);
    183160        size_t astep = nSize * sizeof(double);
    184         Lapack::JacobiSVD(B, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
    185 
    186         double *W = Transpose(Vt, nSize, nSize);
    187 
    188         delete[] B;
     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;
    189166        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        }
    190178
    191179        for (int i = 0; i < nSize; i++)

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

    r389 r817  
    1212{
    1313public:
    14         static void SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates);
     14        static void weightedMDS(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates, double *weights);
    1515};
    1616

  • cpp/frams/model/similarity/simil_model.cpp

    r782 r817  
    3535
    3636static ParamEntry MSparam_tab[] = {
    37                 { "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", },
     37                { "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", },
    3838                { "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.", },
    3939                { "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
    … …  
    4141                { "simil_partgeom", 0, 0, "Weight of parts' geometric distances", "f 0 100 0", FIELD(m_adFactors[3]), "", },
    4242                { "simil_fixedZaxis", 0, 0, "Fix 'z' (vertical) axis?", "d 0 1 0", FIELD(fixedZaxis), "", },
     43                { "simil_weightedMDS", 0, 0, "Should weighted MDS be used?", "d 0 1 0", FIELD(wMDS), "", },
    4344                { "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.", },
    4445                { 0, },
    … …  
    7576        isFuzzy = 0;
    7677        fuzzyDepth = 10;
     78       
     79        //Determines whether weighted MDS should be used.
     80        wMDS = 0;
    7781}
    7882
    … …  
    19811985                int nSize = m_Mod[iMod]->getPartCount();
    19821986
    1983                 double *pDistances = (double *)malloc(nSize * nSize * sizeof(double));
     1987                double *pDistances =  new double[nSize * nSize];
    19841988
    19851989                for (int i = 0; i < nSize; i++)
    … …  
    19931997                Pt3D P1Pos, P2Pos; // positions of parts
    19941998                double dDistance; // the distance between Parts
     1999               
     2000                double *weights = new double[nSize];
     2001                for (int i = 0; i < nSize; i++)
     2002                {
     2003                        if (wMDS==1)
     2004                                weights[i] = 0;
     2005                        else
     2006                                weights[i] = 1;
     2007                }
     2008                               
     2009                if (wMDS==1)
     2010                        for (int i = 0; i < pModel->getJointCount(); i++)
     2011                        {
     2012                                weights[pModel->getJoint(i)->p1_refno]++;
     2013                                weights[pModel->getJoint(i)->p2_refno]++;
     2014                        }
     2015               
    19952016                for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
    19962017                {
    … …  
    20202041                } // for (iP1)
    20212042
    2022                 MatrixTools::SVD(vEigenvalues, nSize, pDistances, m_aPositions[iMod]);
     2043                MatrixTools::weightedMDS(vEigenvalues, nSize, pDistances, m_aPositions[iMod], weights);
    20232044                if (fixedZaxis == 1) //restore the original vertical coordinate of each Part
    20242045                {
    … …  
    20292050                }
    20302051
    2031                 free(pDistances);
     2052                delete[] pDistances;
     2053                delete[] weights;
    20322054        }
    20332055

  • cpp/frams/model/similarity/simil_model.h

    r606 r817  
    8888        int fuzzyDepth;
    8989        int isFuzzy;
     90       
     91        //For wMDS = 1 weighted MDS with vertex degrees as weights is used for the alignment.
     92        int wMDS;
    9093
    9194        /// Interface to local parameters
Note: See TracChangeset for help on using the changeset viewer.