5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
12 /*********************************************************************************************************************************/
14 inline double SIGN(const double a, const double b)
16 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
18 /*********************************************************************************************************************************/
20 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
22 vector<vector<double> > product;
24 int first_rows = first.size();
25 int first_cols = first[0].size();
26 int second_cols = second[0].size();
28 product.resize(first_rows);
29 for(int i=0;i<first_rows;i++){
30 product[i].resize(second_cols);
33 for(int i=0;i<first_rows;i++){
34 for(int j=0;j<second_cols;j++){
36 if (m->control_pressed) { return product; }
39 for(int k=0;k<first_cols;k++){
40 product[i][j] += first[i][k] * second[k][j];
48 m->errorOut(e, "LinearAlgebra", "matrix_mult");
54 /*********************************************************************************************************************************/
56 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
58 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
60 double scale, hh, h, g, f;
67 for(int i=n-1;i>0;i--){
71 for(int k=0;k<l+1;k++){
72 scale += fabs(a[i][k]);
78 for(int k=0;k<l+1;k++){
80 h += a[i][k] * a[i][k];
83 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
88 for(int j=0;j<l+1;j++){
89 a[j][i] = a[i][j] / h;
91 for(int k=0;k<j+1;k++){
92 g += a[j][k] * a[i][k];
94 for(int k=j+1;k<l+1;k++){
95 g += a[k][j] * a[i][k];
101 for(int j=0;j<l+1;j++){
103 e[j] = g = e[j] - hh * f;
104 for(int k=0;k<j+1;k++){
105 a[j][k] -= (f * e[k] + g * a[i][k]);
120 for(int i=0;i<n;i++){
123 for(int j=0;j<l;j++){
125 for(int k=0;k<l;k++){
126 g += a[i][k] * a[k][j];
128 for(int k=0;k<l;k++){
129 a[k][j] -= g * a[k][i];
135 for(int j=0;j<l;j++){
136 a[j][i] = a[i][j] = 0.0;
142 catch(exception& e) {
143 m->errorOut(e, "LinearAlgebra", "tred2");
148 /*********************************************************************************************************************************/
150 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
152 /*********************************************************************************************************************************/
154 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
156 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
159 double s, r, p, g, f, dd, c, b;
162 for(int i=1;i<=n;i++){
167 for(int l=0;l<n;l++){
170 for(myM=l;myM<n-1;myM++){
171 dd = fabs(d[myM]) + fabs(d[myM+1]);
172 if(fabs(e[myM])+dd == dd) break;
175 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
176 g = (d[l+1]-d[l]) / (2.0 * e[l]);
178 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
181 for(i=myM-1;i>=l;i--){
184 e[i+1] = (r=pythag(f,g));
193 r = (d[i] - g) * s + 2.0 * c * b;
194 d[i+1] = g + ( p = s * r);
196 for(int k=0;k<n;k++){
198 z[k][i+1] = s * z[k][i] + c * f;
199 z[k][i] = c * z[k][i] - s * f;
202 if(r == 0.00 && i >= l) continue;
211 for(int i=0;i<n;i++){
213 for(int j=i;j<n;j++){
221 for(int j=0;j<n;j++){
231 catch(exception& e) {
232 m->errorOut(e, "LinearAlgebra", "qtli");
236 /*********************************************************************************************************************************/
237 //groups by dimension
238 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
241 vector< vector<double> > dists; dists.resize(axes.size());
242 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes.size(), 0.0); }
244 if (dimensions == 1) { //one dimension calc = abs(x-y)
246 for (int i = 0; i < dists.size(); i++) {
248 if (m->control_pressed) { return dists; }
250 for (int j = 0; j < i; j++) {
251 dists[i][j] = abs(axes[i][0] - axes[j][0]);
252 dists[j][i] = dists[i][j];
256 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
258 for (int i = 0; i < dists.size(); i++) {
260 if (m->control_pressed) { return dists; }
262 for (int j = 0; j < i; j++) {
264 for (int k = 0; k < dimensions; k++) {
265 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
268 dists[i][j] = sqrt(sum);
269 dists[j][i] = dists[i][j];
277 catch(exception& e) {
278 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
282 /*********************************************************************************************************************************/
283 //returns groups by dimensions from dimensions by groups
284 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
287 vector< vector<double> > dists; dists.resize(axes[0].size());
288 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); }
290 if (axes.size() == 1) { //one dimension calc = abs(x-y)
292 for (int i = 0; i < dists.size(); i++) {
294 if (m->control_pressed) { return dists; }
296 for (int j = 0; j < i; j++) {
297 dists[i][j] = abs(axes[0][i] - axes[0][j]);
298 dists[j][i] = dists[i][j];
302 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
304 for (int i = 0; i < dists[0].size(); i++) {
306 if (m->control_pressed) { return dists; }
308 for (int j = 0; j < i; j++) {
310 for (int k = 0; k < axes.size(); k++) {
311 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
314 dists[i][j] = sqrt(sum);
315 dists[j][i] = dists[i][j];
323 catch(exception& e) {
324 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
328 /*********************************************************************************************************************************/
329 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
332 //find average for - X
333 vector<float> averageEuclid; averageEuclid.resize(euclidDists.size(), 0.0);
334 for (int i = 0; i < euclidDists.size(); i++) {
335 for (int j = 0; j < euclidDists[i].size(); j++) {
336 averageEuclid[i] += euclidDists[i][j];
339 for (int i = 0; i < averageEuclid.size(); i++) { averageEuclid[i] = averageEuclid[i] / (float) euclidDists.size(); }
341 //find average for - Y
342 vector<float> averageUser; averageUser.resize(userDists.size(), 0.0);
343 for (int i = 0; i < userDists.size(); i++) {
344 for (int j = 0; j < userDists[i].size(); j++) {
345 averageUser[i] += userDists[i][j];
348 for (int i = 0; i < averageUser.size(); i++) { averageUser[i] = averageUser[i] / (float) userDists.size(); }
350 double numerator = 0.0;
351 double denomTerm1 = 0.0;
352 double denomTerm2 = 0.0;
354 for (int i = 0; i < euclidDists.size(); i++) {
356 for (int k = 0; k < i; k++) {
358 float Yi = userDists[i][k];
359 float Xi = euclidDists[i][k];
361 numerator += ((Xi - averageEuclid[k]) * (Yi - averageUser[k]));
362 denomTerm1 += ((Xi - averageEuclid[k]) * (Xi - averageEuclid[k]));
363 denomTerm2 += ((Yi - averageUser[k]) * (Yi - averageUser[k]));
367 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
368 double r = numerator / denom;
372 catch(exception& e) {
373 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
377 /*********************************************************************************************************************************/