Changeset 817 for cpp/frams/model/similarity/SVD/matrix_tools.cpp
 Timestamp:
 10/07/18 13:23:11 (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

cpp/frams/model/similarity/SVD/matrix_tools.cpp
r455 r817 84 84 } 85 85 86 double *Transpose(double *&A, int arow, int acol )86 double *Transpose(double *&A, int arow, int acol, double *&toDel, int delSize) 87 87 { 88 88 double *result = Create(acol * arow); … … 94 94 } 95 95 96 if (delSize != 0) 97 delete[] toDel; 98 96 99 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 104 void 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 101 128 @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the 102 129 decomposed distance matrix (or rather, to be strict, of the matrix of scalar products … … 106 133 @param pDistances [IN] matrix of distances between parts. 107 134 @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix. 135 @param weights [IN] vector of row weights. 108 136 */ 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 } 137 void 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 } 177 154 178 155 double *Eigenvalues = Create(nSize); 179 156 double *S = Create(nSize * nSize); 180 157 181 // call SVD function158 // call the SVD function 182 159 double *Vt = Create(nSize * nSize); 183 160 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; 189 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 } 190 178 191 179 for (int i = 0; i < nSize; i++)
Note: See TracChangeset
for help on using the changeset viewer.