Ignore:
Timestamp:
10/07/18 13:23:11 (4 years ago)
Author:
oriona
Message:

MDS procedure replaced with weighted MDS procedure.

File:
1 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++)
Note: See TracChangeset for help on using the changeset viewer.