SNAP Library 6.0, Developer Reference  2020-12-09 16:24:20
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
linalg.cpp
Go to the documentation of this file.
1 // Sparse-Column-Matrix
3 void TSparseColMatrix::PMultiply(const TFltVV& B, int ColId, TFltV& Result) const {
4  Assert(B.GetRows() >= ColN && Result.Len() >= RowN);
5  int i, j; TFlt *ResV = Result.BegI();
6  for (i = 0; i < RowN; i++) ResV[i] = 0.0;
7  for (j = 0; j < ColN; j++) {
8  const TIntFltKdV& ColV = ColSpVV[j]; int len = ColV.Len();
9  for (i = 0; i < len; i++) {
10  ResV[ColV[i].Key] += ColV[i].Dat * B(j,ColId);
11  }
12  }
13 }
14 
15 void TSparseColMatrix::PMultiply(const TFltV& Vec, TFltV& Result) const {
16  Assert(Vec.Len() >= ColN && Result.Len() >= RowN);
17  int i, j; TFlt *ResV = Result.BegI();
18  for (i = 0; i < RowN; i++) ResV[i] = 0.0;
19  for (j = 0; j < ColN; j++) {
20  const TIntFltKdV& ColV = ColSpVV[j]; int len = ColV.Len();
21  for (i = 0; i < len; i++) {
22  ResV[ColV[i].Key] += ColV[i].Dat * Vec[j];
23  }
24  }
25 }
26 
27 void TSparseColMatrix::PMultiplyT(const TFltVV& B, int ColId, TFltV& Result) const {
28  Assert(B.GetRows() >= RowN && Result.Len() >= ColN);
29  int i, j, len; TFlt *ResV = Result.BegI();
30  for (j = 0; j < ColN; j++) {
31  const TIntFltKdV& ColV = ColSpVV[j];
32  len = ColV.Len(); ResV[j] = 0.0;
33  for (i = 0; i < len; i++) {
34  ResV[j] += ColV[i].Dat * B(ColV[i].Key, ColId);
35  }
36  }
37 }
38 
39 void TSparseColMatrix::PMultiplyT(const TFltV& Vec, TFltV& Result) const {
40  Assert(Vec.Len() >= RowN && Result.Len() >= ColN);
41  int i, j, len; TFlt *VecV = Vec.BegI(), *ResV = Result.BegI();
42  for (j = 0; j < ColN; j++) {
43  const TIntFltKdV& ColV = ColSpVV[j];
44  len = ColV.Len(); ResV[j] = 0.0;
45  for (i = 0; i < len; i++) {
46  ResV[j] += ColV[i].Dat * VecV[ColV[i].Key];
47  }
48  }
49 }
50 
52 // Sparse-Row-Matrix
53 TSparseRowMatrix::TSparseRowMatrix(const TStr& MatlabMatrixFNm) {
54  FILE *F = fopen(MatlabMatrixFNm.CStr(), "rt"); IAssert(F != NULL);
56  RowN = 0; ColN = 0;
57  while (! feof(F)) {
58  int row=-1, col=-1; float val;
59  if (fscanf(F, "%d %d %f\n", &row, &col, &val) == 3) {
60  IAssert(row > 0 && col > 0);
61  MtxV.Add(TTriple<TInt, TInt, TSFlt>(row, col, val));
62  RowN = TMath::Mx(RowN, row);
63  ColN = TMath::Mx(ColN, col);
64  }
65  }
66  fclose(F);
67  // create matrix
68  MtxV.Sort();
69  RowSpVV.Gen(RowN);
70  int cnt = 0;
71  for (int row = 1; row <= RowN; row++) {
72  while (cnt < MtxV.Len() && MtxV[cnt].Val1 == row) {
73  RowSpVV[row-1].Add(TIntFltKd(MtxV[cnt].Val2-1, MtxV[cnt].Val3()));
74  cnt++;
75  }
76  }
77 }
78 
79 void TSparseRowMatrix::PMultiplyT(const TFltVV& B, int ColId, TFltV& Result) const {
80  Assert(B.GetRows() >= RowN && Result.Len() >= ColN);
81  for (int i = 0; i < ColN; i++) Result[i] = 0.0;
82  for (int j = 0; j < RowN; j++) {
83  const TIntFltKdV& RowV = RowSpVV[j]; int len = RowV.Len();
84  for (int i = 0; i < len; i++) {
85  Result[RowV[i].Key] += RowV[i].Dat * B(j,ColId);
86  }
87  }
88 }
89 
90 void TSparseRowMatrix::PMultiplyT(const TFltV& Vec, TFltV& Result) const {
91  Assert(Vec.Len() >= RowN && Result.Len() >= ColN);
92  for (int i = 0; i < ColN; i++) Result[i] = 0.0;
93  for (int j = 0; j < RowN; j++) {
94  const TIntFltKdV& RowV = RowSpVV[j]; int len = RowV.Len();
95  for (int i = 0; i < len; i++) {
96  Result[RowV[i].Key] += RowV[i].Dat * Vec[j];
97  }
98  }
99 }
100 
101 void TSparseRowMatrix::PMultiply(const TFltVV& B, int ColId, TFltV& Result) const {
102  Assert(B.GetRows() >= ColN && Result.Len() >= RowN);
103  for (int j = 0; j < RowN; j++) {
104  const TIntFltKdV& RowV = RowSpVV[j];
105  int len = RowV.Len(); Result[j] = 0.0;
106  for (int i = 0; i < len; i++) {
107  Result[j] += RowV[i].Dat * B(RowV[i].Key, ColId);
108  }
109  }
110 }
111 
112 void TSparseRowMatrix::PMultiply(const TFltV& Vec, TFltV& Result) const {
113  Assert(Vec.Len() >= ColN && Result.Len() >= RowN);
114  for (int j = 0; j < RowN; j++) {
115  const TIntFltKdV& RowV = RowSpVV[j];
116  int len = RowV.Len(); Result[j] = 0.0;
117  for (int i = 0; i < len; i++) {
118  Result[j] += RowV[i].Dat * Vec[RowV[i].Key];
119  }
120  }
121 }
122 
124 // Full-Col-Matrix
125 TFullColMatrix::TFullColMatrix(const TStr& MatlabMatrixFNm): TMatrix() {
126  TLAMisc::LoadMatlabTFltVV(MatlabMatrixFNm, ColV);
127  RowN=ColV[0].Len(); ColN=ColV.Len();
128  for (int i = 0; i < ColN; i++) {
129  IAssertR(ColV[i].Len() == RowN, TStr::Fmt("%d != %d", ColV[i].Len(), RowN));
130  }
131 }
132 
133 void TFullColMatrix::PMultiplyT(const TFltVV& B, int ColId, TFltV& Result) const {
134  Assert(B.GetRows() >= RowN && Result.Len() >= ColN);
135  for (int i = 0; i < ColN; i++) {
136  Result[i] = TLinAlg::DotProduct(B, ColId, ColV[i]);
137  }
138 }
139 
140 void TFullColMatrix::PMultiplyT(const TFltV& Vec, TFltV& Result) const {
141  Assert(Vec.Len() >= RowN && Result.Len() >= ColN);
142  for (int i = 0; i < ColN; i++) {
143  Result[i] = TLinAlg::DotProduct(Vec, ColV[i]);
144  }
145 }
146 
147 void TFullColMatrix::PMultiply(const TFltVV& B, int ColId, TFltV& Result) const {
148  Assert(B.GetRows() >= ColN && Result.Len() >= RowN);
149  for (int i = 0; i < RowN; i++) { Result[i] = 0.0; }
150  for (int i = 0; i < ColN; i++) {
151  TLinAlg::AddVec(B(i, ColId), ColV[i], Result, Result);
152  }
153 }
154 
155 void TFullColMatrix::PMultiply(const TFltV& Vec, TFltV& Result) const {
156  Assert(Vec.Len() >= ColN && Result.Len() >= RowN);
157  for (int i = 0; i < RowN; i++) { Result[i] = 0.0; }
158  for (int i = 0; i < ColN; i++) {
159  TLinAlg::AddVec(Vec[i], ColV[i], Result, Result);
160  }
161 }
162 
164 // Basic Linear Algebra Operations
165 double TLinAlg::DotProduct(const TFltV& x, const TFltV& y) {
166  IAssertR(x.Len() == y.Len(), TStr::Fmt("%d != %d", x.Len(), y.Len()));
167  double result = 0.0; int Len = x.Len();
168  for (int i = 0; i < Len; i++)
169  result += x[i] * y[i];
170  return result;
171 }
172 
173 double TLinAlg::DotProduct(const TFltVV& X, int ColIdX, const TFltVV& Y, int ColIdY) {
174  Assert(X.GetRows() == Y.GetRows());
175  double result = 0.0, len = X.GetRows();
176  for (int i = 0; i < len; i++)
177  result = result + X(i,ColIdX) * Y(i,ColIdY);
178  return result;
179 }
180 
181 double TLinAlg::DotProduct(const TFltVV& X, int ColId, const TFltV& Vec) {
182  Assert(X.GetRows() == Vec.Len());
183  double result = 0.0; int Len = X.GetRows();
184  for (int i = 0; i < Len; i++)
185  result += X(i,ColId) * Vec[i];
186  return result;
187 }
188 
189 double TLinAlg::DotProduct(const TIntFltKdV& x, const TIntFltKdV& y) {
190  const int xLen = x.Len(), yLen = y.Len();
191  double Res = 0.0; int i1 = 0, i2 = 0;
192  while (i1 < xLen && i2 < yLen) {
193  if (x[i1].Key < y[i2].Key) i1++;
194  else if (x[i1].Key > y[i2].Key) i2++;
195  else { Res += x[i1].Dat * y[i2].Dat; i1++; i2++; }
196  }
197  return Res;
198 }
199 
200 double TLinAlg::DotProduct(const TFltV& x, const TIntFltKdV& y) {
201  double Res = 0.0; const int xLen = x.Len(), yLen = y.Len();
202  for (int i = 0; i < yLen; i++) {
203  const int key = y[i].Key;
204  if (key < xLen) Res += y[i].Dat * x[key];
205  }
206  return Res;
207 }
208 
209 double TLinAlg::DotProduct(const TFltVV& X, int ColId, const TIntFltKdV& y) {
210  double Res = 0.0; const int n = X.GetRows(), yLen = y.Len();
211  for (int i = 0; i < yLen; i++) {
212  const int key = y[i].Key;
213  if (key < n) Res += y[i].Dat * X(key,ColId);
214  }
215  return Res;
216 }
217 
218 void TLinAlg::LinComb(const double& p, const TFltV& x,
219  const double& q, const TFltV& y, TFltV& z) {
220 
221  Assert(x.Len() == y.Len() && y.Len() == z.Len());
222  const int Len = x.Len();
223  for (int i = 0; i < Len; i++) {
224  z[i] = p * x[i] + q * y[i]; }
225 }
226 
227 void TLinAlg::ConvexComb(const double& p, const TFltV& x, const TFltV& y, TFltV& z) {
228  AssertR(0.0 <= p && p <= 1.0, TFlt::GetStr(p));
229  LinComb(p, x, 1.0 - p, y, z);
230 }
231 
232 void TLinAlg::AddVec(const double& k, const TFltV& x, const TFltV& y, TFltV& z) {
233  LinComb(k, x, 1.0, y, z);
234 }
235 
236 void TLinAlg::AddVec(const double& k, const TIntFltKdV& x, const TFltV& y, TFltV& z) {
237  Assert(y.Len() == z.Len());
238  z = y; // first we set z to be y
239  // and than we add x to z (==y)
240  const int xLen = x.Len(), yLen = y.Len();
241  for (int i = 0; i < xLen; i++) {
242  const int ii = x[i].Key;
243  if (ii < yLen) {
244  z[ii] = k * x[i].Dat + y[ii];
245  }
246  }
247 }
248 
249 void TLinAlg::AddVec(const double& k, const TIntFltKdV& x, TFltV& y) {
250  const int xLen = x.Len(), yLen = y.Len();
251  for (int i = 0; i < xLen; i++) {
252  const int ii = x[i].Key;
253  if (ii < yLen) {
254  y[ii] += k * x[i].Dat;
255  }
256  }
257 }
258 
259 void TLinAlg::AddVec(double k, const TFltVV& X, int ColIdX, TFltVV& Y, int ColIdY){
260  Assert(X.GetRows() == Y.GetRows());
261  const int len = Y.GetRows();
262  for (int i = 0; i < len; i++) {
263  Y(i,ColIdY) = Y(i,ColIdY) + k * X(i, ColIdX);
264  }
265 }
266 
267 void TLinAlg::AddVec(double k, const TFltVV& X, int ColId, TFltV& Result){
268  Assert(X.GetRows() == Result.Len());
269  const int len = Result.Len();
270  for (int i = 0; i < len; i++) {
271  Result[i] = Result[i] + k * X(i, ColId);
272  }
273 }
274 
275 void TLinAlg::AddVec(const TIntFltKdV& x, const TIntFltKdV& y, TIntFltKdV& z) {
277 }
278 
279 double TLinAlg::SumVec(const TFltV& x) {
280  const int len = x.Len();
281  double Res = 0.0;
282  for (int i = 0; i < len; i++) {
283  Res += x[i];
284  }
285  return Res;
286 }
287 
288 double TLinAlg::SumVec(double k, const TFltV& x, const TFltV& y) {
289  Assert(x.Len() == y.Len());
290  const int len = x.Len();
291  double Res = 0.0;
292  for (int i = 0; i < len; i++) {
293  Res += k * x[i] + y[i];
294  }
295  return Res;
296 }
297 
298 double TLinAlg::EuclDist2(const TFltV& x, const TFltV& y) {
299  Assert(x.Len() == y.Len());
300  const int len = x.Len();
301  double Res = 0.0;
302  for (int i = 0; i < len; i++) {
303  Res += TMath::Sqr(x[i] - y[i]);
304  }
305  return Res;
306 }
307 
308 double TLinAlg::EuclDist2(const TFltPr& x, const TFltPr& y) {
309  return TMath::Sqr(x.Val1 - y.Val1) + TMath::Sqr(x.Val2 - y.Val2);
310 }
311 
312 double TLinAlg::EuclDist(const TFltV& x, const TFltV& y) {
313  return sqrt(EuclDist2(x, y));
314 }
315 
316 double TLinAlg::EuclDist(const TFltPr& x, const TFltPr& y) {
317  return sqrt(EuclDist2(x, y));
318 }
319 
320 double TLinAlg::Norm2(const TFltV& x) {
321  return DotProduct(x, x);
322 }
323 
324 double TLinAlg::Norm(const TFltV& x) {
325  return sqrt(Norm2(x));
326 }
327 
329  const double xNorm = Norm(x);
330  if (xNorm > 0.0) { MultiplyScalar(1/xNorm, x, x); }
331 }
332 
333 double TLinAlg::Norm2(const TIntFltKdV& x) {
334  double Result = 0;
335  for (int i = 0; i < x.Len(); i++) {
336  Result += TMath::Sqr(x[i].Dat);
337  }
338  return Result;
339 }
340 
341 double TLinAlg::Norm(const TIntFltKdV& x) {
342  return sqrt(Norm2(x));
343 }
344 
346  MultiplyScalar(1/Norm(x), x, x);
347 }
348 
349 double TLinAlg::Norm2(const TFltVV& X, int ColId) {
350  return DotProduct(X, ColId, X, ColId);
351 }
352 
353 double TLinAlg::Norm(const TFltVV& X, int ColId) {
354  return sqrt(Norm2(X, ColId));
355 }
356 
357 double TLinAlg::NormL1(const TFltV& x) {
358  double norm = 0.0; const int Len = x.Len();
359  for (int i = 0; i < Len; i++)
360  norm += TFlt::Abs(x[i]);
361  return norm;
362 }
363 
364 double TLinAlg::NormL1(double k, const TFltV& x, const TFltV& y) {
365  Assert(x.Len() == y.Len());
366  double norm = 0.0; const int len = x.Len();
367  for (int i = 0; i < len; i++) {
368  norm += TFlt::Abs(k * x[i] + y[i]);
369  }
370  return norm;
371 }
372 
373 double TLinAlg::NormL1(const TIntFltKdV& x) {
374  double norm = 0.0; const int Len = x.Len();
375  for (int i = 0; i < Len; i++)
376  norm += TFlt::Abs(x[i].Dat);
377  return norm;
378 }
379 
381  const double xNorm = NormL1(x);
382  if (xNorm > 0.0) { MultiplyScalar(1/xNorm, x, x); }
383 }
384 
386  const double xNorm = NormL1(x);
387  if (xNorm > 0.0) { MultiplyScalar(1/xNorm, x, x); }
388 }
389 
390 double TLinAlg::NormLinf(const TFltV& x) {
391  double norm = 0.0; const int Len = x.Len();
392  for (int i = 0; i < Len; i++)
393  norm = TFlt::GetMx(TFlt::Abs(x[i]), norm);
394  return norm;
395 }
396 
397 double TLinAlg::NormLinf(const TIntFltKdV& x) {
398  double norm = 0.0; const int Len = x.Len();
399  for (int i = 0; i < Len; i++)
400  norm = TFlt::GetMx(TFlt::Abs(x[i].Dat), norm);
401  return norm;
402 }
403 
405  const double xNormLinf = NormLinf(x);
406  if (xNormLinf > 0.0) { MultiplyScalar(1.0/xNormLinf, x, x); }
407 }
408 
410  const double xNormLInf = NormLinf(x);
411  if (xNormLInf> 0.0) { MultiplyScalar(1.0/xNormLInf, x, x); }
412 }
413 
414 void TLinAlg::MultiplyScalar(const double& k, const TFltV& x, TFltV& y) {
415  Assert(x.Len() == y.Len());
416  int Len = x.Len();
417  for (int i = 0; i < Len; i++)
418  y[i] = k * x[i];
419 }
420 
421 void TLinAlg::MultiplyScalar(const double& k, const TIntFltKdV& x, TIntFltKdV& y) {
422  Assert(x.Len() == y.Len());
423  int Len = x.Len();
424  for (int i = 0; i < Len; i++)
425  y[i].Dat = k * x[i].Dat;
426 }
427 
428 void TLinAlg::Multiply(const TFltVV& A, const TFltV& x, TFltV& y) {
429  Assert(A.GetCols() == x.Len() && A.GetRows() == y.Len());
430  int n = A.GetRows(), m = A.GetCols();
431  for (int i = 0; i < n; i++) {
432  y[i] = 0.0;
433  for (int j = 0; j < m; j++)
434  y[i] += A(i,j) * x[j];
435  }
436 }
437 
438 void TLinAlg::Multiply(const TFltVV& A, const TFltV& x, TFltVV& C, int ColId) {
439  Assert(A.GetCols() == x.Len() && A.GetRows() == C.GetRows());
440  int n = A.GetRows(), m = A.GetCols();
441  for (int i = 0; i < n; i++) {
442  C(i,ColId) = 0.0;
443  for (int j = 0; j < m; j++)
444  C(i,ColId) += A(i,j) * x[j];
445  }
446 }
447 
448 void TLinAlg::Multiply(const TFltVV& A, const TFltVV& B, int ColId, TFltV& y) {
449  Assert(A.GetCols() == B.GetRows() && A.GetRows() == y.Len());
450  int n = A.GetRows(), m = A.GetCols();
451  for (int i = 0; i < n; i++) {
452  y[i] = 0.0;
453  for (int j = 0; j < m; j++)
454  y[i] += A(i,j) * B(j,ColId);
455  }
456 }
457 
458 void TLinAlg::Multiply(const TFltVV& A, const TFltVV& B, int ColIdB, TFltVV& C, int ColIdC){
459  Assert(A.GetCols() == B.GetRows() && A.GetRows() == C.GetRows());
460  int n = A.GetRows(), m = A.GetCols();
461  for (int i = 0; i < n; i++) {
462  C(i,ColIdC) = 0.0;
463  for (int j = 0; j < m; j++)
464  C(i,ColIdC) += A(i,j) * B(j,ColIdB);
465  }
466 }
467 
468 void TLinAlg::MultiplyT(const TFltVV& A, const TFltV& x, TFltV& y) {
469  Assert(A.GetRows() == x.Len() && A.GetCols() == y.Len());
470  int n = A.GetCols(), m = A.GetRows();
471  for (int i = 0; i < n; i++) {
472  y[i] = 0.0;
473  for (int j = 0; j < m; j++)
474  y[i] += A(j,i) * x[j];
475  }
476 }
477 
478 void TLinAlg::Multiply(const TFltVV& A, const TFltVV& B, TFltVV& C) {
479  Assert(A.GetRows() == C.GetRows() && B.GetCols() == C.GetCols() && A.GetCols() == B.GetRows());
480 
481  int n = C.GetRows(), m = C.GetCols(), l = A.GetCols();
482  for (int i = 0; i < n; i++) {
483  for (int j = 0; j < m; j++) {
484  double sum = 0.0;
485  for (int k = 0; k < l; k++)
486  sum += A(i,k)*B(k,j);
487  C(i,j) = sum;
488  }
489  }
490 }
491 
492 // general matrix multiplication (GEMM)
493 void TLinAlg::Gemm(const double& Alpha, const TFltVV& A, const TFltVV& B, const double& Beta,
494  const TFltVV& C, TFltVV& D, const int& TransposeFlags) {
495 
496  bool tA = (TransposeFlags & GEMM_A_T) == GEMM_A_T;
497  bool tB = (TransposeFlags & GEMM_B_T) == GEMM_B_T;
498  bool tC = (TransposeFlags & GEMM_C_T) == GEMM_C_T;
499 
500  // setting dimensions
501  int a_i = tA ? A.GetRows() : A.GetCols();
502  int a_j = tA ? A.GetCols() : A.GetRows();
503 
504  int b_i = tB ? A.GetRows() : A.GetCols();
505  int b_j = tB ? A.GetCols() : A.GetRows();
506 
507  int c_i = tC ? A.GetRows() : A.GetCols();
508  int c_j = tC ? A.GetCols() : A.GetRows();
509 
510  int d_i = D.GetCols();
511  int d_j = D.GetRows();
512 
513  // assertions for dimensions
514  bool Cnd = a_i == c_j && b_i == c_i && a_i == b_j && c_i == d_i && c_j == d_j;
515  if (!Cnd) {
516  Assert(Cnd);
517  }
518 
519  double Aij, Bij, Cij;
520 
521  // rows
522  for (int j = 0; j < a_j; j++) {
523  // cols
524  for (int i = 0; i < a_i; i++) {
525  // not optimized for speed (!)
526  double sum = 0;
527  for (int k = 0; k < a_i; k++) {
528  Aij = tA ? A.At(j, k) : A.At(k, j);
529  Bij = tB ? B.At(k, i) : B.At(i, k);
530  sum += Alpha * Aij * Bij;
531  }
532  Cij = tC ? C.At(i, j) : C.At(j, i);
533  sum += Beta * Cij;
534  D.At(i, j) = sum;
535  }
536  }
537 }
538 
539 void TLinAlg::Transpose(const TFltVV& A, TFltVV& B) {
540  Assert(B.GetRows() == A.GetCols() && B.GetCols() == A.GetRows());
541  for (int i = 0; i < A.GetCols(); i++) {
542  for (int j = 0; j < A.GetRows(); j++) {
543  B.At(i, j) = A.At(j, i);
544  }
545  }
546 }
547 
548 void TLinAlg::Inverse(const TFltVV& A, TFltVV& B, const TLinAlgInverseType& DecompType) {
549  switch (DecompType) {
550  case DECOMP_SVD:
551  InverseSVD(A, B);
552  }
553 }
554 
555 void TLinAlg::InverseSVD(const TFltVV& M, TFltVV& B) {
556  // create temp matrices
557  TFltVV U, V;
558  TFltV E;
559  TSvd SVD;
560 
561  U.Gen(M.GetRows(), M.GetRows());
562  V.Gen(M.GetCols(), M.GetCols());
563 
564  // do the SVD decompostion
565  SVD.Svd(M, U, E, V);
566 
567  // calculate reciprocal values for diagonal matrix = inverse diagonal
568  for (int i = 0; i < E.Len(); i++) {
569  E[i] = 1 / E[i];
570  }
571 
572  // calculate pseudoinverse: M^(-1) = V * E^(-1) * U'
573  for (int i = 0; i < U.GetCols(); i++) {
574  for (int j = 0; j < V.GetRows(); j++) {
575  double sum = 0;
576  for (int k = 0; k < U.GetCols(); k++) {
577  sum += E[k] * V.At(i, k) * U.At(j, k);
578  }
579  B.At(i, j) = sum;
580  }
581  }
582 }
583 
585  IAssert(Q.Len() > 0);
586  int m = Q.Len(); // int n = Q[0].Len();
587  for (int i = 0; i < m; i++) {
588  printf("%d\r",i);
589  for (int j = 0; j < i; j++) {
590  double r = TLinAlg::DotProduct(Q[i], Q[j]);
591  TLinAlg::AddVec(-r,Q[j],Q[i],Q[i]);
592  }
593  TLinAlg::Normalize(Q[i]);
594  }
595  printf("\n");
596 }
597 
598 void TLinAlg::GS(TFltVV& Q) {
599  int m = Q.GetCols(), n = Q.GetRows();
600  for (int i = 0; i < m; i++) {
601  for (int j = 0; j < i; j++) {
602  double r = TLinAlg::DotProduct(Q, i, Q, j);
603  TLinAlg::AddVec(-r,Q,j,Q,i);
604  }
605  double nr = TLinAlg::Norm(Q,i);
606  for (int k = 0; k < n; k++)
607  Q(k,i) = Q(k,i) / nr;
608  }
609 }
610 
611 void TLinAlg::Rotate(const double& OldX, const double& OldY, const double& Angle, double& NewX, double& NewY) {
612  NewX = OldX*cos(Angle) - OldY*sin(Angle);
613  NewY = OldX*sin(Angle) + OldY*cos(Angle);
614 }
615 
616 void TLinAlg::AssertOrtogonality(const TVec<TFltV>& Vecs, const double& Threshold) {
617  int m = Vecs.Len();
618  for (int i = 0; i < m; i++) {
619  for (int j = 0; j < i; j++) {
620  double res = DotProduct(Vecs[i], Vecs[j]);
621  if (TFlt::Abs(res) > Threshold)
622  printf("<%d,%d> = %.5f", i,j,res);
623  }
624  double norm = Norm2(Vecs[i]);
625  if (TFlt::Abs(norm-1) > Threshold)
626  printf("||%d|| = %.5f", i, norm);
627  }
628 }
629 
630 void TLinAlg::AssertOrtogonality(const TFltVV& Vecs, const double& Threshold) {
631  int m = Vecs.GetCols();
632  for (int i = 0; i < m; i++) {
633  for (int j = 0; j < i; j++) {
634  double res = DotProduct(Vecs, i, Vecs, j);
635  if (TFlt::Abs(res) > Threshold)
636  printf("<%d,%d> = %.5f", i, j, res);
637  }
638  double norm = Norm2(Vecs, i);
639  if (TFlt::Abs(norm-1) > Threshold)
640  printf("||%d|| = %.5f", i, norm);
641  }
642  printf("\n");
643 }
644 
646 // Numerical Linear Algebra
647 double TNumericalStuff::sqr(double a) {
648  return a == 0.0 ? 0.0 : a*a;
649 }
650 
651 double TNumericalStuff::sign(double a, double b) {
652  return b >= 0.0 ? fabs(a) : -fabs(a);
653 }
654 
655 void TNumericalStuff::nrerror(const TStr& error_text) {
656  printf("NR_ERROR: %s", error_text.CStr());
657  throw new TNSException(error_text);
658 }
659 
660 double TNumericalStuff::pythag(double a, double b) {
661  double absa = fabs(a), absb = fabs(b);
662  if (absa > absb)
663  return absa*sqrt(1.0+sqr(absb/absa));
664  else
665  return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+sqr(absa/absb)));
666 }
667 
669  int l,k,j,i;
670  double scale,hh,h,g,f;
671  for (i=n;i>=2;i--) {
672  l=i-1;
673  h=scale=0.0;
674  if (l > 1) {
675  for (k=1;k<=l;k++)
676  scale += fabs(a(i-1,k-1).Val);
677  if (scale == 0.0) //Skip transformation.
678  e[i]=a(i-1,l-1);
679  else {
680  for (k=1;k<=l;k++) {
681  a(i-1,k-1) /= scale; //Use scaled a's for transformation.
682  h += a(i-1,k-1)*a(i-1,k-1);
683  }
684  f=a(i-1,l-1);
685  g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
686  IAssertR(_isnan(g) == 0, TFlt::GetStr(h));
687  e[i]=scale*g;
688  h -= f*g; //Now h is equation (11.2.4).
689  a(i-1,l-1)=f-g; //Store u in the ith row of a.
690  f=0.0;
691  for (j=1;j<=l;j++) {
692  // Next statement can be omitted if eigenvectors not wanted
693  a(j-1,i-1)=a(i-1,j-1)/h; //Store u=H in ith column of a.
694  g=0.0; //Form an element of A  u in g.
695  for (k=1;k<=j;k++)
696  g += a(j-1,k-1)*a(i-1,k-1);
697  for (k=j+1;k<=l;k++)
698  g += a(k-1,j-1)*a(i-1,k-1);
699  e[j]=g/h; //Form element of p in temporarily unused element of e.
700  f += e[j]*a(i-1,j-1);
701  }
702  hh=f/(h+h); //Form K, equation (11.2.11).
703  for (j=1;j<=l;j++) { //Form q and store in e overwriting p.
704  f=a(i-1,j-1);
705  e[j]=g=e[j]-hh*f;
706  for (k=1;k<=j;k++) { //Reduce a, equation (11.2.13).
707  a(j-1,k-1) -= (f*e[k]+g*a(i-1,k-1));
708  Assert(_isnan(static_cast<double>(a(j-1,k-1))) == 0);
709  }
710  }
711  }
712  } else
713  e[i]=a(i-1,l-1);
714  d[i]=h;
715  }
716  // Next statement can be omitted if eigenvectors not wanted
717  d[1]=0.0;
718  e[1]=0.0;
719  // Contents of this loop can be omitted if eigenvectors not
720  // wanted except for statement d[i]=a[i][i];
721  for (i=1;i<=n;i++) { //Begin accumulation of transformationmatrices.
722  l=i-1;
723  if (d[i]) { //This block skipped when i=1.
724  for (j=1;j<=l;j++) {
725  g=0.0;
726  for (k=1;k<=l;k++) //Use u and u=H stored in a to form PQ.
727  g += a(i-1,k-1)*a(k-1,j-1);
728  for (k=1;k<=l;k++) {
729  a(k-1,j-1) -= g*a(k-1,i-1);
730  Assert(_isnan(static_cast<double>(a(k-1,j-1))) == 0);
731  }
732  }
733  }
734  d[i]=a(i-1,i-1); //This statement remains.
735  a(i-1,i-1)=1.0; //Reset row and column of a to identity matrix for next iteration.
736  for (j=1;j<=l;j++) a(j-1,i-1)=a(i-1,j-1)=0.0;
737  }
738 }
739 
741  int m,l,iter,i,k; // N = n+1;
742  double s,r,p,g,f,dd,c,b;
743  // Convenient to renumber the elements of e
744  for (i=2;i<=n;i++) e[i-1]=e[i];
745  e[n]=0.0;
746  for (l=1;l<=n;l++) {
747  iter=0;
748  do {
749  // Look for a single small subdiagonal element to split the matrix.
750  for (m=l;m<=n-1;m++) {
751  dd=TFlt::Abs(d[m])+TFlt::Abs(d[m+1]);
752  if ((double)(TFlt::Abs(e[m])+dd) == dd) break;
753  }
754  if (m != l) {
755  if (iter++ == 60) nrerror("Too many iterations in EigSymmetricTridiag");
756  //Form shift.
757  g=(d[l+1]-d[l])/(2.0*e[l]);
758  r=pythag(g,1.0);
759  //This is dm - ks.
760  g=d[m]-d[l]+e[l]/(g+sign(r,g));
761  s=c=1.0;
762  p=0.0;
763  // A plane rotation as in the original QL, followed by
764  // Givens rotations to restore tridiagonal form
765  for (i=m-1;i>=l;i--) {
766  f=s*e[i];
767  b=c*e[i];
768  e[i+1]=(r=pythag(f,g));
769  // Recover from underflow.
770  if (r == 0.0) {
771  d[i+1] -= p;
772  e[m]=0.0;
773  break;
774  }
775  s=f/r;
776  c=g/r;
777  g=d[i+1]-p;
778  r=(d[i]-g)*s+2.0*c*b;
779  d[i+1]=g+(p=s*r);
780  g=c*r-b;
781  // Next loop can be omitted if eigenvectors not wanted
782  for (k=0;k<n;k++) {
783  f=z(k,i);
784  z(k,i)=s*z(k,i-1)+c*f;
785  z(k,i-1)=c*z(k,i-1)-s*f;
786  }
787  }
788  if (r == 0.0 && i >= l) continue;
789  d[l] -= p;
790  e[l]=g;
791  e[m]=0.0;
792  }
793  } while (m != l);
794  }
795 }
796 
798  Assert(A.GetRows() == A.GetCols());
799  int n = A.GetRows(); p.Reserve(n,n);
800 
801  int i,j,k;
802  double sum;
803  for (i=1;i<=n;i++) {
804  for (j=i;j<=n;j++) {
805  for (sum=A(i-1,j-1),k=i-1;k>=1;k--) sum -= A(i-1,k-1)*A(j-1,k-1);
806  if (i == j) {
807  if (sum <= 0.0)
808  nrerror("choldc failed");
809  p[i-1]=sqrt(sum);
810  } else A(j-1,i-1)=sum/p[i-1];
811  }
812  }
813 }
814 
815 void TNumericalStuff::CholeskySolve(const TFltVV& A, const TFltV& p, const TFltV& b, TFltV& x) {
816  IAssert(A.GetRows() == A.GetCols());
817  int n = A.GetRows(); x.Reserve(n,n);
818 
819  int i,k;
820  double sum;
821 
822  // Solve L * y = b, storing y in x
823  for (i=1;i<=n;i++) {
824  for (sum=b[i-1],k=i-1;k>=1;k--)
825  sum -= A(i-1,k-1)*x[k-1];
826  x[i-1]=sum/p[i-1];
827  }
828 
829  // Solve L^T * x = y
830  for (i=n;i>=1;i--) {
831  for (sum=x[i-1],k=i+1;k<=n;k++)
832  sum -= A(k-1,i-1)*x[k-1];
833  x[i-1]=sum/p[i-1];
834  }
835 }
836 
838  IAssert(A.GetRows() == A.GetCols());
839  TFltV p; CholeskyDecomposition(A, p);
840  CholeskySolve(A, p, b, x);
841 }
842 
844  IAssert(A.GetRows() == A.GetCols());
845  int n = A.GetRows(); TFltV x(n);
846 
847  int i, j, k; double sum;
848  for (i = 0; i < n; i++) {
849  // solve L * y = e_i, store in x
850  // elements from 0 to i-1 are 0.0
851  for (j = 0; j < i; j++) x[j] = 0.0;
852  // solve l_ii * y_i = 1 => y_i = 1/l_ii
853  x[i] = 1/p[i];
854  // solve y_j for j > i
855  for (j = i+1; j < n; j++) {
856  for (sum = 0.0, k = i; k < j; k++)
857  sum -= A(j,k) * x[k];
858  x[j] = sum / p[j];
859  }
860 
861  // solve L'* x = y, store in upper triangule of A
862  for (j = n-1; j >= i; j--) {
863  for (sum = x[j], k = j+1; k < n; k++)
864  sum -= A(k,j)*x[k];
865  x[j] = sum/p[j];
866  }
867  for (int j = i; j < n; j++) A(i,j) = x[j];
868  }
869 
870 }
871 
873  IAssert(A.GetRows() == A.GetCols());
874  TFltV p;
875  // first we calculate cholesky decomposition of A
876  CholeskyDecomposition(A, p);
877  // than we solve system A x_i = e_i for i = 1..n
878  InverseSubstitute(A, p);
879 }
880 
882  IAssert(A.GetRows() == A.GetCols());
883  int n = A.GetRows(); TFltV x(n), p(n);
884 
885  int i, j, k; double sum;
886  // copy upper triangle to lower one as we'll overwrite upper one
887  for (i = 0; i < n; i++) {
888  p[i] = A(i,i);
889  for (j = i+1; j < n; j++)
890  A(j,i) = A(i,j);
891  }
892  // solve
893  for (i = 0; i < n; i++) {
894  // solve R * x = e_i, store in x
895  // elements from 0 to i-1 are 0.0
896  for (j = n-1; j > i; j--) x[j] = 0.0;
897  // solve l_ii * y_i = 1 => y_i = 1/l_ii
898  x[i] = 1/p[i];
899  // solve y_j for j > i
900  for (j = i-1; j >= 0; j--) {
901  for (sum = 0.0, k = i; k > j; k--)
902  sum -= A(k,j) * x[k];
903  x[j] = sum / p[j];
904  }
905  for (int j = 0; j <= i; j++) A(j,i) = x[j];
906  }
907 }
908 
909 void TNumericalStuff::LUDecomposition(TFltVV& A, TIntV& indx, double& d) {
910  Assert(A.GetRows() == A.GetCols());
911  int n = A.GetRows(); indx.Reserve(n,n);
912 
913  int i=0,imax=0,j=0,k=0;
914  double big,dum,sum,temp;
915  TFltV vv(n); // vv stores the implicit scaling of each row.
916  d=1.0; // No row interchanges yet.
917 
918  // Loop over rows to get the implicit scaling information.
919  for (i=1;i<=n;i++) {
920  big=0.0;
921  for (j=1;j<=n;j++)
922  if ((temp=TFlt::Abs(A(i-1,j-1))) > big) big=temp;
923  if (big == 0.0) nrerror("Singular matrix in routine LUDecomposition");
924  vv[i-1]=1.0/big;
925  }
926 
927  for (j=1;j<=n;j++) {
928  for (i=1;i<j;i++) {
929  sum=A(i-1,j-1);
930  for (k=1;k<i;k++) sum -= A(i-1,k-1)*A(k-1,j-1);
931  A(i-1,j-1)=sum;
932  }
933  big=0.0; //Initialize for the search for largest pivot element.
934  for (i=j;i<=n;i++) {
935  sum=A(i-1,j-1);
936  for (k=1;k<j;k++)
937  sum -= A(i-1,k-1)*A(k-1,j-1);
938  A(i-1,j-1)=sum;
939 
940  //Is the figure of merit for the pivot better than the best so far?
941  if ((dum=vv[i-1] * TFlt::Abs(sum)) >= big) {
942  big=dum;
943  imax=i;
944  }
945  }
946 
947  //Do we need to interchange rows?
948  if (j != imax) {
949  //Yes, do so...
950  for (k=1;k<=n;k++) {
951  dum=A(imax-1,k-1);
952  A(imax-1,k-1)=A(j-1,k-1); // Tadej: imax-1,k looks wrong
953  A(j-1,k-1)=dum;
954  }
955  //...and change the parity of d.
956  d = -d;
957  vv[imax-1]=vv[j-1]; //Also interchange the scale factor.
958  }
959  indx[j-1]=imax;
960 
961  //If the pivot element is zero the matrix is singular (at least to the precision of the
962  //algorithm). For some applications on singular matrices, it is desirable to substitute
963  //TINY for zero.
964  if (A(j-1,j-1) == 0.0) A(j-1,j-1)=1e-20;
965 
966  //Now, finally, divide by the pivot element.
967  if (j != n) {
968  dum=1.0/(A(j-1,j-1));
969  for (i=j+1;i<=n;i++) A(i-1,j-1) *= dum;
970  }
971  } //Go back for the next column in the reduction.
972 }
973 
974 void TNumericalStuff::LUSolve(const TFltVV& A, const TIntV& indx, TFltV& b) {
975  Assert(A.GetRows() == A.GetCols());
976  int n = A.GetRows();
977  int i,ii=0,ip,j;
978  double sum;
979  for (i=1;i<=n;i++) {
980  ip=indx[i-1];
981  sum=b[ip-1];
982  b[ip-1]=b[i-1];
983  if (ii)
984  for (j=ii;j<=i-1;j++) sum -= A(i-1,j-1)*b[j-1];
985  else if (sum) ii=i;b[i-1]=sum;
986  }
987  for (i=n;i>=1;i--) {
988  sum=b[i-1];
989  for (j=i+1;j<=n;j++) sum -= A(i-1,j-1)*b[j-1];
990  b[i-1]=sum/A(i-1,i-1);
991  }
992 }
993 
995  TIntV indx; double d;
996  LUDecomposition(A, indx, d);
997  x = b;
998  LUSolve(A, indx, x);
999 }
1000 
1002 // Sparse-SVD
1004  const TFltVV& Vec, int ColId, TFltV& Result) {
1005  TFltV tmp(Matrix.GetRows());
1006  // tmp = A * Vec(:,ColId)
1007  Matrix.Multiply(Vec, ColId, tmp);
1008  // Vec = A' * tmp
1009  Matrix.MultiplyT(tmp, Result);
1010 }
1011 
1013  const TFltV& Vec, TFltV& Result) {
1014  TFltV tmp(Matrix.GetRows());
1015  // tmp = A * Vec
1016  Matrix.Multiply(Vec, tmp);
1017  // Vec = A' * tmp
1018  Matrix.MultiplyT(tmp, Result);
1019 }
1020 
1022  int NumSV, int IterN, TFltV& SgnValV) {
1023 
1024  int i, j, k;
1025  int N = Matrix.GetCols(), M = NumSV;
1026  TFltVV Q(N, M);
1027 
1028  // Q = rand(N,M)
1029  TRnd rnd;
1030  for (i = 0; i < N; i++) {
1031  for (j = 0; j < M; j++)
1032  Q(i,j) = rnd.GetUniDev();
1033  }
1034 
1035  TFltV tmp(N);
1036  for (int IterC = 0; IterC < IterN; IterC++) {
1037  printf("%d..", IterC);
1038  // Gram-Schmidt
1039  TLinAlg::GS(Q);
1040  // Q = A'*A*Q
1041  for (int ColId = 0; ColId < M; ColId++) {
1042  MultiplyATA(Matrix, Q, ColId, tmp);
1043  for (k = 0; k < N; k++) Q(k,ColId) = tmp[k];
1044  }
1045  }
1046 
1047  SgnValV.Reserve(NumSV,0);
1048  for (i = 0; i < NumSV; i++)
1049  SgnValV.Add(sqrt(TLinAlg::Norm(Q,i)));
1050  TLinAlg::GS(Q);
1051 }
1052 
1054  const int& NumEig, TFltV& EigValV,
1055  const bool& DoLocalReortoP, const bool& SvdMatrixProductP) {
1056 
1057  if (SvdMatrixProductP) {
1058  // if this fails, use transposed matrix
1059  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1060  } else {
1061  IAssert(Matrix.GetRows() == Matrix.GetCols());
1062  }
1063 
1064  const int N = Matrix.GetCols(); // size of matrix
1065  TFltV r(N), v0(N), v1(N); // current vector and 2 previous ones
1066  TFltV alpha(NumEig, 0), beta(NumEig, 0); // diagonal and subdiagonal of T
1067 
1068  printf("Calculating %d eigen-values of %d x %d matrix\n", NumEig, N, N);
1069 
1070  // set starting vector
1071  //TRnd Rnd(0);
1072  for (int i = 0; i < N; i++) {
1073  r[i] = 1/sqrt((double)N); // Rnd.GetNrmDev();
1074  v0[i] = v1[i] = 0.0;
1075  }
1076  beta.Add(TLinAlg::Norm(r));
1077 
1078  for (int j = 0; j < NumEig; j++) {
1079  printf("%d\r", j+1);
1080  // v_j -> v_(j-1)
1081  v0 = v1;
1082  // v_j = (1/beta_(j-1)) * r
1083  TLinAlg::MultiplyScalar(1/beta[j], r, v1);
1084  // r = A*v_j
1085  if (SvdMatrixProductP) {
1086  // A = Matrix'*Matrix
1087  MultiplyATA(Matrix, v1, r);
1088  } else {
1089  // A = Matrix
1090  Matrix.Multiply(v1, r);
1091  }
1092  // r = r - beta_(j-1) * v_(j-1)
1093  TLinAlg::AddVec(-beta[j], v0, r, r);
1094  // alpha_j = vj'*r
1095  alpha.Add(TLinAlg::DotProduct(v1, r));
1096  // r = r - v_j * alpha_j
1097  TLinAlg::AddVec(-alpha[j], v1, r, r);
1098  // reortogonalization if neessary
1099  if (DoLocalReortoP) { } //TODO
1100  // beta_j = ||r||_2
1101  beta.Add(TLinAlg::Norm(r));
1102  // compoute approximatie eigenvalues T_j
1103  // test bounds for convergence
1104  }
1105  printf("\n");
1106 
1107  // prepare matrix T
1108  TFltV d(NumEig + 1), e(NumEig + 1);
1109  d[1] = alpha[0]; d[0] = e[0] = e[1] = 0.0;
1110  for (int i = 1; i < NumEig; i++) {
1111  d[i+1] = alpha[i]; e[i+1] = beta[i]; }
1112  // solve eigne problem for tridiagonal matrix with diag d and subdiag e
1113  TFltVV S(NumEig+1,NumEig+1); // eigen-vectors
1114  TLAMisc::FillIdentity(S); // make it identity
1115  TNumericalStuff::EigSymmetricTridiag(d, e, NumEig, S); // solve
1116  //TLAMisc::PrintTFltV(d, "AllEigV");
1117 
1118  // check convergence
1119  TFltKdV AllEigValV(NumEig, 0);
1120  for (int i = 1; i <= NumEig; i++) {
1121  const double ResidualNorm = TFlt::Abs(S(i-1, NumEig-1) * beta.Last());
1122  if (ResidualNorm < 1e-5)
1123  AllEigValV.Add(TFltKd(TFlt::Abs(d[i]), d[i]));
1124  }
1125 
1126  // prepare results
1127  AllEigValV.Sort(false); EigValV.Gen(NumEig, 0);
1128  for (int i = 0; i < AllEigValV.Len(); i++) {
1129  if (i == 0 || (TFlt::Abs(AllEigValV[i].Dat/AllEigValV[i-1].Dat) < 0.9999))
1130  EigValV.Add(AllEigValV[i].Dat);
1131  }
1132 }
1133 
1134 void TSparseSVD::Lanczos(const TMatrix& Matrix, int NumEig,
1135  int Iters, const TSpSVDReOrtoType& ReOrtoType,
1136  TFltV& EigValV, TFltVV& EigVecVV, const bool& SvdMatrixProductP) {
1137 
1138  if (SvdMatrixProductP) {
1139  // if this fails, use transposed matrix
1140  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1141  } else {
1142  IAssert(Matrix.GetRows() == Matrix.GetCols());
1143  }
1144  IAssertR(NumEig <= Iters, TStr::Fmt("%d <= %d", NumEig, Iters));
1145 
1146  //if (ReOrtoType == ssotFull) printf("Full reortogonalization\n");
1147  int i, N = Matrix.GetCols(), K = 0; // K - current dimension of T
1148  double t = 0.0, eps = 1e-6; // t - 1-norm of T
1149 
1150  //sequence of Ritz's vectors
1151  TFltVV Q(N, Iters);
1152  double tmp = 1/sqrt((double)N);
1153  for (i = 0; i < N; i++) {
1154  Q(i,0) = tmp;
1155  }
1156  //converget Ritz's vectors
1157  TVec<TFltV> ConvgQV(Iters);
1158  TIntV CountConvgV(Iters);
1159  for (i = 0; i < Iters; i++) CountConvgV[i] = 0;
1160  // const int ConvgTreshold = 50;
1161 
1162  //diagonal and subdiagonal of T
1163  TFltV d(Iters+1), e(Iters+1);
1164  //eigenvectors of T
1165  //TFltVV V;
1166  TFltVV V(Iters, Iters);
1167 
1168  // z - current Lanczos's vector
1169  TFltV z(N), bb(Iters), aa(Iters), y(N);
1170  //printf("svd(%d,%d)...\n", NumEig, Iters);
1171 
1172  if (SvdMatrixProductP) {
1173  // A = Matrix'*Matrix
1174  MultiplyATA(Matrix, Q, 0, z);
1175  } else {
1176  // A = Matrix
1177  Matrix.Multiply(Q, 0, z);
1178  }
1179 
1180  for (int j = 0; j < (Iters-1); j++) {
1181  //printf("%d..\r",j+2);
1182 
1183  //calculates (j+1)-th Lanczos's vector
1184  // aa[j] = <Q(:,j), z>
1185  aa[j] = TLinAlg::DotProduct(Q, j, z);
1186  //printf(" %g -- ", aa[j].Val); //HACK
1187 
1188  TLinAlg::AddVec(-aa[j], Q, j, z);
1189  if (j > 0) {
1190  // z := -aa[j] * Q(:,j) + z
1191  TLinAlg::AddVec(-bb[j-1], Q, j-1, z);
1192 
1193  //reortogonalization
1194  if (ReOrtoType == ssotSelective || ReOrtoType == ssotFull) {
1195  for (i = 0; i <= j; i++) {
1196  // if i-tj vector converget, than we have to ortogonalize against it
1197  if ((ReOrtoType == ssotFull) ||
1198  (bb[j-1] * TFlt::Abs(V(K-1, i)) < eps * t)) {
1199 
1200  ConvgQV[i].Reserve(N,N); CountConvgV[i]++;
1201  TFltV& vec = ConvgQV[i];
1202  //vec = Q * V(:,i)
1203  for (int k = 0; k < N; k++) {
1204  vec[k] = 0.0;
1205  for (int l = 0; l < K; l++) {
1206  vec[k] += Q(k,l) * V(l,i);
1207  }
1208  }
1209  TLinAlg::AddVec(-TLinAlg::DotProduct(ConvgQV[i], z), ConvgQV[i], z ,z);
1210  }
1211  }
1212  }
1213  }
1214 
1215  //adds (j+1)-th Lanczos's vector to Q
1216  bb[j] = TLinAlg::Norm(z);
1217  if (!(bb[j] > 1e-10)) {
1218  printf("Rank of matrix is only %d\n", j+2);
1219  printf("Last singular value is %g\n", bb[j].Val);
1220  break;
1221  }
1222 
1223  for (i = 0; i < N; i++) {
1224  Q(i, j+1) = z[i] / bb[j];
1225  }
1226 
1227  //next Lanzcos vector
1228  if (SvdMatrixProductP) {
1229  // A = Matrix'*Matrix
1230  MultiplyATA(Matrix, Q, j+1, z);
1231  } else {
1232  // A = Matrix
1233  Matrix.Multiply(Q, j+1, z);
1234  }
1235 
1236  //calculate T (K x K matrix)
1237  K = j + 2;
1238  // calculate diagonal
1239  for (i = 1; i < K; i++) d[i] = aa[i-1];
1240  d[K] = TLinAlg::DotProduct(Q, K-1, z);
1241  // calculate subdiagonal
1242  e[1] = 0.0;
1243  for (i = 2; i <= K; i++) e[i] = bb[i-2];
1244 
1245  //calculate 1-norm of T
1246  t = TFlt::GetMx(TFlt::Abs(d[1]) + TFlt::Abs(e[2]), TFlt::Abs(e[K]) + TFlt::Abs(d[K]));
1247  for (i = 2; i < K; i++) {
1248  t = TFlt::GetMx(t, TFlt::Abs(e[i]) + TFlt::Abs(d[i]) + TFlt::Abs(e[i+1]));
1249  }
1250 
1251  //set V to identity matrix
1252  //V.Gen(K,K);
1253  for (i = 0; i < K; i++) {
1254  for (int k = 0; k < K; k++) {
1255  V(i,k) = 0.0;
1256  }
1257  V(i,i) = 1.0;
1258  }
1259 
1260  //eigenvectors of T
1262  }//for
1263 
1264  //printf("\n");
1265 
1266  // Finds NumEig largest eigen values
1267  TFltIntKdV sv(K);
1268  for (i = 0; i < K; i++) {
1269  sv[i].Key = TFlt::Abs(d[i+1]);
1270  sv[i].Dat = i;
1271  }
1272  sv.Sort(false);
1273 
1274  TFltV uu(Matrix.GetRows());
1275  const int FinalNumEig = TInt::GetMn(NumEig, K);
1276  EigValV.Reserve(FinalNumEig,0);
1277  EigVecVV.Gen(Matrix.GetCols(), FinalNumEig);
1278  for (i = 0; i < FinalNumEig; i++) {
1279  //printf("s[%d] = %20.15f\r", i, sv[i].Key.Val);
1280  int ii = sv[i].Dat;
1281  double sigma = d[ii+1].Val;
1282  // calculate singular value
1283  EigValV.Add(sigma);
1284  // calculate i-th right singular vector ( V := Q * W )
1285  TLinAlg::Multiply(Q, V, ii, EigVecVV, i);
1286  }
1287  //printf("done \n");
1288 }
1289 
1290 void TSparseSVD::Lanczos2(const TMatrix& Matrix, int MaxNumEig,
1291  int MaxSecs, const TSpSVDReOrtoType& ReOrtoType,
1292  TFltV& EigValV, TFltVV& EigVecVV, const bool& SvdMatrixProductP) {
1293 
1294  if (SvdMatrixProductP) {
1295  // if this fails, use transposed matrix
1296  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1297  } else {
1298  IAssert(Matrix.GetRows() == Matrix.GetCols());
1299  }
1300  //IAssertR(NumEig <= Iters, TStr::Fmt("%d <= %d", NumEig, Iters));
1301 
1302  //if (ReOrtoType == ssotFull) printf("Full reortogonalization\n");
1303  int i, N = Matrix.GetCols(), K = 0; // K - current dimension of T
1304  double t = 0.0, eps = 1e-6; // t - 1-norm of T
1305 
1306  //sequence of Ritz's vectors
1307  TFltVV Q(N, MaxNumEig);
1308  double tmp = 1/sqrt((double)N);
1309  for (i = 0; i < N; i++)
1310  Q(i,0) = tmp;
1311  //converget Ritz's vectors
1312  TVec<TFltV> ConvgQV(MaxNumEig);
1313  TIntV CountConvgV(MaxNumEig);
1314  for (i = 0; i < MaxNumEig; i++) CountConvgV[i] = 0;
1315  // const int ConvgTreshold = 50;
1316 
1317  //diagonal and subdiagonal of T
1318  TFltV d(MaxNumEig+1), e(MaxNumEig+1);
1319  //eigenvectors of T
1320  //TFltVV V;
1321  TFltVV V(MaxNumEig, MaxNumEig);
1322 
1323  // z - current Lanczos's vector
1324  TFltV z(N), bb(MaxNumEig), aa(MaxNumEig), y(N);
1325  //printf("svd(%d,%d)...\n", NumEig, Iters);
1326 
1327  if (SvdMatrixProductP) {
1328  // A = Matrix'*Matrix
1329  MultiplyATA(Matrix, Q, 0, z);
1330  } else {
1331  // A = Matrix
1332  Matrix.Multiply(Q, 0, z);
1333  }
1334  TExeTm ExeTm;
1335  for (int j = 0; j < (MaxNumEig-1); j++) {
1336  printf("%d [%s]..\r",j+2, ExeTm.GetStr());
1337  if (ExeTm.GetSecs() > MaxSecs) { break; }
1338 
1339  //calculates (j+1)-th Lanczos's vector
1340  // aa[j] = <Q(:,j), z>
1341  aa[j] = TLinAlg::DotProduct(Q, j, z);
1342  //printf(" %g -- ", aa[j].Val); //HACK
1343 
1344  TLinAlg::AddVec(-aa[j], Q, j, z);
1345  if (j > 0) {
1346  // z := -aa[j] * Q(:,j) + z
1347  TLinAlg::AddVec(-bb[j-1], Q, j-1, z);
1348 
1349  //reortogonalization
1350  if (ReOrtoType == ssotSelective || ReOrtoType == ssotFull) {
1351  for (i = 0; i <= j; i++) {
1352  // if i-tj vector converget, than we have to ortogonalize against it
1353  if ((ReOrtoType == ssotFull) ||
1354  (bb[j-1] * TFlt::Abs(V(K-1, i)) < eps * t)) {
1355 
1356  ConvgQV[i].Reserve(N,N); CountConvgV[i]++;
1357  TFltV& vec = ConvgQV[i];
1358  //vec = Q * V(:,i)
1359  for (int k = 0; k < N; k++) {
1360  vec[k] = 0.0;
1361  for (int l = 0; l < K; l++)
1362  vec[k] += Q(k,l) * V(l,i);
1363  }
1364  TLinAlg::AddVec(-TLinAlg::DotProduct(ConvgQV[i], z), ConvgQV[i], z ,z);
1365  }
1366  }
1367  }
1368  }
1369 
1370  //adds (j+1)-th Lanczos's vector to Q
1371  bb[j] = TLinAlg::Norm(z);
1372  if (!(bb[j] > 1e-10)) {
1373  printf("Rank of matrix is only %d\n", j+2);
1374  printf("Last singular value is %g\n", bb[j].Val);
1375  break;
1376  }
1377  for (i = 0; i < N; i++)
1378  Q(i, j+1) = z[i] / bb[j];
1379 
1380  //next Lanzcos vector
1381  if (SvdMatrixProductP) {
1382  // A = Matrix'*Matrix
1383  MultiplyATA(Matrix, Q, j+1, z);
1384  } else {
1385  // A = Matrix
1386  Matrix.Multiply(Q, j+1, z);
1387  }
1388 
1389  //calculate T (K x K matrix)
1390  K = j + 2;
1391  // calculate diagonal
1392  for (i = 1; i < K; i++) d[i] = aa[i-1];
1393  d[K] = TLinAlg::DotProduct(Q, K-1, z);
1394  // calculate subdiagonal
1395  e[1] = 0.0;
1396  for (i = 2; i <= K; i++) e[i] = bb[i-2];
1397 
1398  //calculate 1-norm of T
1399  t = TFlt::GetMx(TFlt::Abs(d[1]) + TFlt::Abs(e[2]), TFlt::Abs(e[K]) + TFlt::Abs(d[K]));
1400  for (i = 2; i < K; i++)
1401  t = TFlt::GetMx(t, TFlt::Abs(e[i]) + TFlt::Abs(d[i]) + TFlt::Abs(e[i+1]));
1402 
1403  //set V to identity matrix
1404  //V.Gen(K,K);
1405  for (i = 0; i < K; i++) {
1406  for (int k = 0; k < K; k++)
1407  V(i,k) = 0.0;
1408  V(i,i) = 1.0;
1409  }
1410 
1411  //eigenvectors of T
1413  }//for
1414  printf("... calc %d.", K);
1415  // Finds NumEig largest eigen values
1416  TFltIntKdV sv(K);
1417  for (i = 0; i < K; i++) {
1418  sv[i].Key = TFlt::Abs(d[i+1]);
1419  sv[i].Dat = i;
1420  }
1421  sv.Sort(false);
1422 
1423  TFltV uu(Matrix.GetRows());
1424  const int FinalNumEig = K; //TInt::GetMn(NumEig, K);
1425  EigValV.Reserve(FinalNumEig,0);
1426  EigVecVV.Gen(Matrix.GetCols(), FinalNumEig);
1427  for (i = 0; i < FinalNumEig; i++) {
1428  //printf("s[%d] = %20.15f\r", i, sv[i].Key.Val);
1429  int ii = sv[i].Dat;
1430  double sigma = d[ii+1].Val;
1431  // calculate singular value
1432  EigValV.Add(sigma);
1433  // calculate i-th right singular vector ( V := Q * W )
1434  TLinAlg::Multiply(Q, V, ii, EigVecVV, i);
1435  }
1436  printf(" done\n");
1437 }
1438 
1439 
1441  const int& CalcSV, TFltV& SngValV, const bool& DoLocalReorto) {
1442 
1443  SimpleLanczos(Matrix, CalcSV, SngValV, DoLocalReorto, true);
1444  for (int SngValN = 0; SngValN < SngValV.Len(); SngValN++) {
1445  //IAssert(SngValV[SngValN] >= 0.0);
1446  if (SngValV[SngValN] < 0.0) {
1447  printf("bad sng val: %d %g\n", SngValN, SngValV[SngValN]());
1448  SngValV[SngValN] = 0;
1449  }
1450  SngValV[SngValN] = sqrt(SngValV[SngValN].Val);
1451  }
1452 }
1453 
1454 void TSparseSVD::LanczosSVD(const TMatrix& Matrix, int NumSV,
1455  int Iters, const TSpSVDReOrtoType& ReOrtoType,
1456  TFltV& SgnValV, TFltVV& LeftSgnVecVV, TFltVV& RightSgnVecVV) {
1457 
1458  // solve eigen problem for Matrix'*Matrix
1459  Lanczos(Matrix, NumSV, Iters, ReOrtoType, SgnValV, RightSgnVecVV, true);
1460  // calculate left singular vectors and sqrt singular values
1461  const int FinalNumSV = SgnValV.Len();
1462  LeftSgnVecVV.Gen(Matrix.GetRows(), FinalNumSV);
1463  TFltV LeftSgnVecV(Matrix.GetRows());
1464  for (int i = 0; i < FinalNumSV; i++) {
1465  if (SgnValV[i].Val < 0.0) { SgnValV[i] = 0.0; }
1466  const double SgnVal = sqrt(SgnValV[i]);
1467  SgnValV[i] = SgnVal;
1468  // calculate i-th left singular vector ( U := A * V * S^(-1) )
1469  Matrix.Multiply(RightSgnVecVV, i, LeftSgnVecV);
1470  for (int j = 0; j < LeftSgnVecV.Len(); j++) {
1471  LeftSgnVecVV(j,i) = LeftSgnVecV[j] / SgnVal; }
1472  }
1473  //printf("done \n");
1474 }
1475 
1476 void TSparseSVD::Project(const TIntFltKdV& Vec, const TFltVV& U, TFltV& ProjVec) {
1477  const int m = U.GetCols(); // number of columns
1478 
1479  ProjVec.Gen(m, 0);
1480  for (int j = 0; j < m; j++) {
1481  double x = 0.0;
1482  for (int i = 0; i < Vec.Len(); i++)
1483  x += U(Vec[i].Key, j) * Vec[i].Dat;
1484  ProjVec.Add(x);
1485  }
1486 }
1487 
1489 // Sigmoid
1490 double TSigmoid::EvaluateFit(const TFltIntKdV& data, const double A, const double B)
1491 {
1492  double J = 0.0;
1493  for (int i = 0; i < data.Len(); i++)
1494  {
1495  double zi = data[i].Key; int yi = data[i].Dat;
1496  double e = exp(-A * zi + B);
1497  double denum = 1.0 + e;
1498  double prob = (yi > 0) ? (1.0 / denum) : (e / denum);
1499  J -= log(prob < 1e-20 ? 1e-20 : prob);
1500  }
1501  return J;
1502 }
1503 
1504 void TSigmoid::EvaluateFit(const TFltIntKdV& data, const double A, const double B, double& J, double& JA, double& JB)
1505 {
1506  // J(A, B) = \sum_{i : y_i = 1} ln [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} [ln [1 + e^{-Az_i + B}] - {-Az_i + B}]
1507  // = \sum_i ln [1 + e^{-Az_i + B}] - \sum_{i : y_i = -1} {-Az_i + B}.
1508  // partial J / partial A = \sum_i (-z_i) e^{-Az_i + B} / [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} Az_i.
1509  // partial J / partial B = \sum_i e^{-Az_i + B} / [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} (-1).
1510  J = 0.0; double sum_all_PyNeg = 0.0, sum_all_ziPyNeg = 0.0, sum_yNeg_zi = 0.0, sum_yNeg_1 = 0.0;
1511  for (int i = 0; i < data.Len(); i++)
1512  {
1513  double zi = data[i].Key; int yi = data[i].Dat;
1514  double e = exp(-A * zi + B);
1515  double denum = 1.0 + e;
1516  double prob = (yi > 0) ? (1.0 / denum) : (e / denum);
1517  J -= log(prob < 1e-20 ? 1e-20 : prob);
1518  sum_all_PyNeg += e / denum;
1519  sum_all_ziPyNeg += zi * e / denum;
1520  if (yi < 0) { sum_yNeg_zi += zi; sum_yNeg_1 += 1; }
1521  }
1522  JA = -sum_all_ziPyNeg + sum_yNeg_zi;
1523  JB = sum_all_PyNeg - sum_yNeg_1;
1524 }
1525 
1526 void TSigmoid::EvaluateFit(const TFltIntKdV& data, const double A, const double B, const double U,
1527  const double V, const double lambda, double& J, double& JJ, double& JJJ)
1528 {
1529  // Let E_i = e^{-(A + lambda U) z_i + (B + lambda V)}. Then we have
1530  // J(lambda) = \sum_i ln [1 + E_i] - \sum_{i : y_i = -1} {-(A + lambda U)z_i + (B + lambda V)}.
1531  // J'(lambda) = \sum_i (V - U z_i) E_i / [1 + E_i] - \sum_{i : y_i = -1} {V - U z_i).
1532  // = \sum_i (V - U z_i) [1 - 1 / [1 + E_i]] - \sum_{i : y_i = -1} {V - U z_i).
1533  // J"(lambda) = \sum_i (V - U z_i)^2 E_i / [1 + E_i]^2.
1534  J = 0.0; JJ = 0.0; JJJ = 0.0;
1535  for (int i = 0; i < data.Len(); i++)
1536  {
1537  double zi = data[i].Key; int yi = data[i].Dat;
1538  double e = exp(-A * zi + B);
1539  double denum = 1.0 + e;
1540  double prob = (yi > 0) ? (1.0 / denum) : (e / denum);
1541  J -= log(prob < 1e-20 ? 1e-20 : prob);
1542  double VU = V - U * zi;
1543  JJ += VU * (e / denum); if (yi < 0) JJ -= VU;
1544  JJJ += VU * VU * e / denum / denum;
1545  }
1546 }
1547 
1549  // Let z_i be the projection of the i'th training example, and y_i \in {-1, +1} be its class label.
1550  // Our sigmoid is: P(Y = y | Z = z) = 1 / [1 + e^{-Az + B}]
1551  // and we want to maximize \prod_i P(Y = y_i | Z = z_i)
1552  // = \prod_{i : y_i = 1} 1 / [1 + e^{-Az_i + B}] \prod_{i : y_i = -1} e^{-Az_i + B} / [1 + e^{-Az_i + B}]
1553  // or minimize its negative logarithm,
1554  // J(A, B) = \sum_{i : y_i = 1} ln [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} [ln [1 + e^{-Az_i + B}] - {-Az_i + B}]
1555  // = \sum_i ln [1 + e^{-Az_i + B}] - \sum_{i : y_i = -1} {-Az_i + B}.
1556  // partial J / partial A = \sum_i (-z_i) e^{-Az_i + B} / [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} Az_i.
1557  // partial J / partial B = \sum_i e^{-Az_i + B} / [1 + e^{-Az_i + B}] + \sum_{i : y_i = -1} (-1).
1558  double minProj = data[0].Key, maxProj = data[0].Key;
1559  {for (int i = 1; i < data.Len(); i++) {
1560  double zi = data[i].Key; if (zi < minProj) minProj = zi; if (zi > maxProj) maxProj = zi; }}
1561  //const bool dump = false;
1562  A = 1.0; B = 0.5 * (minProj + maxProj);
1563  double bestJ = 0.0, bestA = 0.0, bestB = 0.0, lambda = 1.0;
1564  for (int nIter = 0; nIter < 50; nIter++)
1565  {
1566  double J, JA, JB; TSigmoid::EvaluateFit(data, A, B, J, JA, JB);
1567  if (nIter == 0 || J < bestJ) { bestJ = J; bestA = A; bestB = B; }
1568  // How far should we move?
1569  //if (dump) printf("Iter %2d: A = %.5f, B = %.5f, J = %.5f, partial = (%.5f, %.5f)\n", nIter, A, B, J, JA, JB);
1570  double norm = TMath::Sqr(JA) + TMath::Sqr(JB);
1571  if (norm < 1e-10) break;
1572  const int cl = -1; // should be -1
1573 
1574  double Jc = TSigmoid::EvaluateFit(data, A + cl * lambda * JA / norm, B + cl * lambda * JB / norm);
1575  //if (dump) printf(" At lambda = %.5f, Jc = %.5f\n", lambda, Jc);
1576  if (Jc > J) {
1577  while (lambda > 1e-5) {
1578  lambda = 0.5 * lambda;
1579  Jc = TSigmoid::EvaluateFit(data, A + cl * lambda * JA / norm, B + cl * lambda * JB / norm);
1580  //if (dump) printf(" At lambda = %.5f, Jc = %.5f\n", lambda, Jc);
1581  } }
1582  else if (Jc < J) {
1583  while (lambda < 1e5) {
1584  double lambda2 = 2 * lambda;
1585  double Jc2 = TSigmoid::EvaluateFit(data, A + cl * lambda2 * JA / norm, B + cl * lambda2 * JB / norm);
1586  //if (dump) printf(" At lambda = %.5f, Jc = %.5f\n", lambda2, Jc2);
1587  if (Jc2 > Jc) break;
1588  lambda = lambda2; Jc = Jc2; } }
1589  if (Jc >= J) break;
1590  A += cl * lambda * JA / norm; B += cl * lambda * JB / norm;
1591  //if (dump) printf(" Lambda = %.5f, new A = %.5f, new B = %.5f, new J = %.5f\n", lambda, A, B, Jc);
1592  }
1593  A = bestA; B = bestB;
1594 }
1595 
1597 // Useful stuff (hopefuly)
1598 void TLAMisc::SaveCsvTFltV(const TFltV& Vec, TSOut& SOut) {
1599  for (int ValN = 0; ValN < Vec.Len(); ValN++) {
1600  SOut.PutFlt(Vec[ValN]); SOut.PutCh(',');
1601  }
1602  SOut.PutLn();
1603 }
1604 
1605 void TLAMisc::SaveMatlabTFltIntKdV(const TIntFltKdV& SpV, const int& ColN, TSOut& SOut) {
1606  const int Len = SpV.Len();
1607  for (int ValN = 0; ValN < Len; ValN++) {
1608  SOut.PutStrLn(TStr::Fmt("%d %d %g", SpV[ValN].Key+1, ColN+1, SpV[ValN].Dat()));
1609  }
1610 }
1611 
1612 void TLAMisc::SaveMatlabTFltV(const TFltV& m, const TStr& FName) {
1613  PSOut out = TFOut::New(FName);
1614  const int RowN = m.Len();
1615  for (int RowId = 0; RowId < RowN; RowId++) {
1616  out->PutStr(TFlt::GetStr(m[RowId], 20, 18));
1617  out->PutCh('\n');
1618  }
1619  out->Flush();
1620 }
1621 
1622 void TLAMisc::SaveMatlabTIntV(const TIntV& m, const TStr& FName) {
1623  PSOut out = TFOut::New(FName);
1624  const int RowN = m.Len();
1625  for (int RowId = 0; RowId < RowN; RowId++) {
1626  out->PutInt(m[RowId]);
1627  out->PutCh('\n');
1628  }
1629  out->Flush();
1630 }
1631 
1632 void TLAMisc::SaveMatlabTFltVVCol(const TFltVV& m, int ColId, const TStr& FName) {
1633  PSOut out = TFOut::New(FName);
1634  const int RowN = m.GetRows();
1635  for (int RowId = 0; RowId < RowN; RowId++) {
1636  out->PutStr(TFlt::GetStr(m(RowId,ColId), 20, 18));
1637  out->PutCh('\n');
1638  }
1639  out->Flush();
1640 }
1641 
1642 
1643 void TLAMisc::SaveMatlabTFltVV(const TFltVV& m, const TStr& FName) {
1644  PSOut out = TFOut::New(FName);
1645  const int RowN = m.GetRows();
1646  const int ColN = m.GetCols();
1647  for (int RowId = 0; RowId < RowN; RowId++) {
1648  for (int ColId = 0; ColId < ColN; ColId++) {
1649  out->PutStr(TFlt::GetStr(m(RowId,ColId), 20, 18));
1650  out->PutCh(' ');
1651  }
1652  out->PutCh('\n');
1653  }
1654  out->Flush();
1655 }
1656 
1658  int RowN, int ColN, const TStr& FName) {
1659 
1660  PSOut out = TFOut::New(FName);
1661  for (int RowId = 0; RowId < RowN; RowId++) {
1662  for (int ColId = 0; ColId < ColN; ColId++) {
1663  out->PutStr(TFlt::GetStr(m(RowId,ColId), 20, 18)); out->PutCh(' ');
1664  }
1665  out->PutCh('\n');
1666  }
1667  out->Flush();
1668 }
1669 
1670 void TLAMisc::LoadMatlabTFltVV(const TStr& FNm, TVec<TFltV>& ColV) {
1671  PSIn SIn = TFIn::New(FNm);
1673  int Row = 0, Col = 0; ColV.Clr();
1674  Lx.GetSym(syFlt, syEof, syEoln);
1675  //printf("%d x %d\r", Row, ColV.Len());
1676  while (Lx.Sym != syEof) {
1677  if (Lx.Sym == syFlt) {
1678  if (ColV.Len() > Col) {
1679  IAssert(ColV[Col].Len() == Row);
1680  ColV[Col].Add(Lx.Flt);
1681  } else {
1682  IAssert(Row == 0);
1683  ColV.Add(TFltV::GetV(Lx.Flt));
1684  }
1685  Col++;
1686  } else if (Lx.Sym == syEoln) {
1687  IAssert(Col == ColV.Len());
1688  Col = 0; Row++;
1689  if (Row%100 == 0) {
1690  //printf("%d x %d\r", Row, ColV.Len());
1691  }
1692  } else {
1693  Fail;
1694  }
1695  Lx.GetSym(syFlt, syEof, syEoln);
1696  }
1697  //printf("\n");
1698  IAssert(Col == ColV.Len() || Col == 0);
1699 }
1700 
1701 void TLAMisc::LoadMatlabTFltVV(const TStr& FNm, TFltVV& MatrixVV) {
1702  TVec<TFltV> ColV; LoadMatlabTFltVV(FNm, ColV);
1703  if (ColV.Empty()) { MatrixVV.Clr(); return; }
1704  const int Rows = ColV[0].Len(), Cols = ColV.Len();
1705  MatrixVV.Gen(Rows, Cols);
1706  for (int RowN = 0; RowN < Rows; RowN++) {
1707  for (int ColN = 0; ColN < Cols; ColN++) {
1708  MatrixVV(RowN, ColN) = ColV[ColN][RowN];
1709  }
1710  }
1711 }
1712 
1713 
1714 void TLAMisc::PrintTFltV(const TFltV& Vec, const TStr& VecNm) {
1715  printf("%s = [", VecNm.CStr());
1716  for (int i = 0; i < Vec.Len(); i++) {
1717  printf("%.5f", Vec[i]());
1718  if (i < Vec.Len() - 1) { printf(", "); }
1719  }
1720  printf("]\n");
1721 }
1722 
1723 
1724 void TLAMisc::PrintTFltVV(const TFltVV& A, const TStr& MatrixNm) {
1725  printf("%s = [\n", MatrixNm.CStr());
1726  for (int j = 0; j < A.GetRows(); j++) {
1727  for (int i = 0; i < A.GetCols(); i++) {
1728  printf("%f\t", A.At(i, j).Val);
1729  }
1730  printf("\n");
1731  }
1732  printf("]\n");
1733 }
1734 
1735 void TLAMisc::PrintTIntV(const TIntV& Vec, const TStr& VecNm) {
1736  printf("%s = [", VecNm.CStr());
1737  for (int i = 0; i < Vec.Len(); i++) {
1738  printf("%d", Vec[i]());
1739  if (i < Vec.Len() - 1) printf(", ");
1740  }
1741  printf("]\n");
1742 }
1743 
1744 
1745 void TLAMisc::FillRnd(TFltV& Vec, TRnd& Rnd) {
1746  int Len = Vec.Len();
1747  for (int i = 0; i < Len; i++)
1748  Vec[i] = Rnd.GetNrmDev();
1749 }
1750 
1752  IAssert(M.GetRows() == M.GetCols());
1753  int Len = M.GetRows();
1754  for (int i = 0; i < Len; i++) {
1755  for (int j = 0; j < Len; j++) M(i,j) = 0.0;
1756  M(i,i) = 1.0;
1757  }
1758 }
1759 
1760 void TLAMisc::FillIdentity(TFltVV& M, const double& Elt) {
1761  IAssert(M.GetRows() == M.GetCols());
1762  int Len = M.GetRows();
1763  for (int i = 0; i < Len; i++) {
1764  for (int j = 0; j < Len; j++) M(i,j) = 0.0;
1765  M(i,i) = Elt;
1766  }
1767 }
1768 
1769 int TLAMisc::SumVec(const TIntV& Vec) {
1770  const int Len = Vec.Len();
1771  int res = 0;
1772  for (int i = 0; i < Len; i++)
1773  res += Vec[i];
1774  return res;
1775 }
1776 
1777 double TLAMisc::SumVec(const TFltV& Vec) {
1778  const int Len = Vec.Len();
1779  double res = 0.0;
1780  for (int i = 0; i < Len; i++)
1781  res += Vec[i];
1782  return res;
1783 }
1784 
1785 void TLAMisc::ToSpVec(const TFltV& Vec, TIntFltKdV& SpVec,
1786  const double& CutSumPrc) {
1787 
1788  // determine minimal element value
1789  IAssert(0.0 <= CutSumPrc && CutSumPrc <= 1.0);
1790  const int Elts = Vec.Len();
1791  double EltSum = 0.0;
1792  for (int EltN = 0; EltN < Elts; EltN++) {
1793  EltSum += TFlt::Abs(Vec[EltN]); }
1794  const double MnEltVal = CutSumPrc * EltSum;
1795  // create sparse vector
1796  SpVec.Clr();
1797  for (int EltN = 0; EltN < Elts; EltN++) {
1798  if (TFlt::Abs(Vec[EltN]) > MnEltVal) {
1799  SpVec.Add(TIntFltKd(EltN, Vec[EltN]));
1800  }
1801  }
1802  SpVec.Pack();
1803 }
1804 
1805 void TLAMisc::ToVec(const TIntFltKdV& SpVec, TFltV& Vec, const int& VecLen) {
1806  Vec.Gen(VecLen); Vec.PutAll(0.0);
1807  int Elts = SpVec.Len();
1808  for (int EltN = 0; EltN < Elts; EltN++) {
1809  if (SpVec[EltN].Key < VecLen) {
1810  Vec[SpVec[EltN].Key] = SpVec[EltN].Dat;
1811  }
1812  }
1813 }
#define IAssert(Cond)
Definition: bd.h:262
static double EuclDist(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:312
Definition: lx.h:127
static void Transpose(const TFltVV &A, TFltVV &B)
Definition: linalg.cpp:539
TSpSVDReOrtoType
Definition: linalg.h:406
void Clr()
Definition: ds.h:2246
static double Norm(const TFltV &x)
Definition: linalg.cpp:324
#define IAssertR(Cond, Reason)
Definition: bd.h:265
static PSOut New(const TStr &FNm, const bool &Append=false)
Definition: fl.cpp:442
static void AssertOrtogonality(const TVec< TFltV > &Vecs, const double &Threshold)
Definition: linalg.cpp:616
static void SaveMatlabTFltIntKdV(const TIntFltKdV &SpV, const int &ColN, TSOut &SOut)
Definition: linalg.cpp:1605
Definition: tm.h:355
TLinAlgInverseType
Definition: linalg.h:268
virtual int PutCh(const char &Ch)=0
Definition: dt.h:11
static double SumVec(const TFltV &x)
Definition: linalg.cpp:279
Definition: ds.h:130
static double sqr(double a)
Definition: linalg.cpp:647
static double EvaluateFit(const TFltIntKdV &data, const double A, const double B)
Definition: linalg.cpp:1490
static const T & Mx(const T &LVal, const T &RVal)
Definition: xmath.h:32
virtual void PMultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:27
static void GS(TVec< TFltV > &Q)
Definition: linalg.cpp:584
static void FillRnd(TFltV &Vec)
Definition: linalg.h:521
static double GetMx(const double &Flt1, const double &Flt2)
Definition: dt.h:1444
void MultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:34
static void Lanczos(const TMatrix &Matrix, int NumEig, int Iters, const TSpSVDReOrtoType &ReOrtoType, TFltV &EigValV, TFltVV &EigVecVV, const bool &SvdMatrixProductP=false)
Definition: linalg.cpp:1134
static void ToVec(const TIntFltKdV &SpVec, TFltV &Vec, const int &VecLen)
Definition: linalg.cpp:1805
#define Fail
Definition: bd.h:238
static void SimpleLanczosSVD(const TMatrix &Matrix, const int &CalcSV, TFltV &SngValV, const bool &DoLocalReortoP=false)
Definition: linalg.cpp:1440
TVec< TIntFltKdV > RowSpVV
Definition: linalg.h:95
static void Lanczos2(const TMatrix &Matrix, int MaxNumEig, int MaxSecs, const TSpSVDReOrtoType &ReOrtoType, TFltV &EigValV, TFltVV &EigVecVV, const bool &SvdMatrixProductP=false)
Definition: linalg.cpp:1290
static void SaveMatlabTFltVV(const TFltVV &m, const TStr &FName)
Definition: linalg.cpp:1643
double Val
Definition: dt.h:1388
Definition: bits.h:119
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
TKeyDat< TInt, TFlt > TIntFltKd
Definition: ds.h:381
Definition: lx.h:51
TFlt A
Definition: linalg.h:456
static void SolveLinearSystem(TFltVV &A, const TFltV &b, TFltV &x)
Definition: linalg.cpp:994
static void SaveMatlabTFltVVCol(const TFltVV &m, int ColId, const TStr &FName)
Definition: linalg.cpp:1632
TKeyDat< TFlt, TFlt > TFltKd
Definition: ds.h:396
static double Sqr(const double &x)
Definition: xmath.h:12
static void LanczosSVD(const TMatrix &Matrix, int NumSV, int Iters, const TSpSVDReOrtoType &ReOrtoType, TFltV &SgnValV, TFltVV &LeftSgnVecVV, TFltVV &RightSgnVecVV)
Definition: linalg.cpp:1454
int PutLn(const int &Lns=1)
Definition: fl.cpp:158
static void MultiplyT(const TFltVV &A, const TFltV &x, TFltV &y)
Definition: linalg.cpp:468
static void InverseSVD(const TFltVV &A, TFltVV &B)
Definition: linalg.cpp:555
TLxSym Sym
Definition: lx.h:149
static void InverseTriagonal(TFltVV &A)
Definition: linalg.cpp:881
Definition: dt.h:1386
bool Empty() const
Tests whether the vector is empty.
Definition: ds.h:570
static void Normalize(TFltV &x)
Definition: linalg.cpp:328
static void NormalizeLinf(TFltV &x)
Definition: linalg.cpp:404
Definition: lx.h:126
static void Project(const TIntFltKdV &Vec, const TFltVV &U, TFltV &ProjVec)
Definition: linalg.cpp:1476
static void SaveMatlabTFltVVMjrSubMtrx(const TFltVV &m, int rowN, int colN, const TStr &FName)
Definition: linalg.cpp:1657
static void SparseMerge(const TKeyDatV &SrcV1, const TKeyDatV &SrcV2, TKeyDatV &DstV)
Definition: linalg.h:548
TFullColMatrix()
Definition: linalg.h:146
static PSIn New(const TStr &FNm)
Definition: fl.cpp:290
static int SumVec(const TIntV &Vec)
Definition: linalg.cpp:1769
Definition: lx.h:126
static void Rotate(const double &OldX, const double &OldY, const double &Angle, double &NewX, double &NewY)
Definition: linalg.cpp:611
static void MultiplyScalar(const double &k, const TFltV &x, TFltV &y)
Definition: linalg.cpp:414
static void Svd(const TFltVV &InMtx, TFltVV &LSingV, TFltV &SingValV, TFltVV &RSingV)
Definition: xmath.cpp:1232
void Clr(const bool &DoDel=true, const TSizeTy &NoDelLim=-1)
Clears the contents of the vector.
Definition: ds.h:1022
static void NormalizeL1(TFltV &x)
Definition: linalg.cpp:380
static int GetMn(const int &Int1, const int &Int2)
Definition: dt.h:1183
void Sort(const bool &Asc=true)
Sorts the elements of the vector.
Definition: ds.h:1318
static void LUSolve(const TFltVV &A, const TIntV &indx, TFltV &b)
Definition: linalg.cpp:974
static double Norm2(const TFltV &x)
Definition: linalg.cpp:320
void PutAll(const TVal &Val)
Sets all elements of the vector to value Val.
Definition: ds.h:1229
int GetRows() const
Definition: linalg.h:45
static void CholeskyDecomposition(TFltVV &A, TFltV &p)
Definition: linalg.cpp:797
static double NormLinf(const TFltV &x)
Definition: linalg.cpp:390
static void AddVec(const double &k, const TFltV &x, const TFltV &y, TFltV &z)
Definition: linalg.cpp:232
Definition: lx.h:51
virtual void PMultiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:101
TSigmoid()
Definition: linalg.h:474
Definition: lx.h:129
Definition: xmath.h:397
static void PrintTIntV(const TIntV &Vec, const TStr &VecNm)
Definition: linalg.cpp:1735
#define Assert(Cond)
Definition: bd.h:251
int GetCols() const
Definition: linalg.h:47
const TVal & Last() const
Returns a reference to the last element of the vector.
Definition: ds.h:579
static void nrerror(const TStr &error_text)
Definition: linalg.cpp:655
int PutInt(const int &Int)
Definition: fl.cpp:93
static void EigSymmetricTridiag(TFltV &d, TFltV &e, int n, TFltVV &z)
Definition: linalg.cpp:740
TVec< TIntFltKdV > ColSpVV
Definition: linalg.h:60
static void PrintTFltVV(const TFltVV &A, const TStr &MatrixNm)
Definition: linalg.cpp:1724
static void SolveSymetricSystem(TFltVV &A, const TFltV &b, TFltV &x)
Definition: linalg.cpp:837
TSizeTy GetRows() const
Definition: ds.h:2252
Definition: fl.h:128
static void SaveMatlabTIntV(const TIntV &m, const TStr &FName)
Definition: linalg.cpp:1622
virtual void PMultiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:3
Definition: lx.h:45
int PutFlt(const double &Flt)
Definition: fl.cpp:109
Definition: linalg.h:8
virtual void Flush()=0
static void Multiply(const TFltVV &A, const TFltV &x, TFltV &y)
Definition: linalg.cpp:428
static double Abs(const double &Flt)
Definition: dt.h:1430
static void SymetricToTridiag(TFltVV &a, int n, TFltV &d, TFltV &e)
Definition: linalg.cpp:668
static void CholeskySolve(const TFltVV &A, const TFltV &p, const TFltV &b, TFltV &x)
Definition: linalg.cpp:815
static void SimpleLanczos(const TMatrix &Matrix, const int &NumEig, TFltV &EigValV, const bool &DoLocalReortoP=false, const bool &SvdMatrixProductP=false)
Definition: linalg.cpp:1053
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
double GetNrmDev()
Definition: dt.cpp:63
TVec< TFltV > ColV
Definition: linalg.h:131
static double NormL1(const TFltV &x)
Definition: linalg.cpp:357
static void LoadMatlabTFltVV(const TStr &FNm, TVec< TFltV > &ColV)
Definition: linalg.cpp:1670
virtual void PMultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:79
static void SaveCsvTFltV(const TFltV &Vec, TSOut &SOut)
Definition: linalg.cpp:1598
static double pythag(double a, double b)
Definition: linalg.cpp:660
Definition: dt.h:412
TIter BegI() const
Returns an iterator pointing to the first element in the vector.
Definition: ds.h:593
static TStr Fmt(const char *FmtStr,...)
Definition: dt.cpp:1599
TFlt B
Definition: linalg.h:457
void Pack()
Reduces vector capacity (frees memory) to match its size.
Definition: ds.h:1057
static void SaveMatlabTFltV(const TFltV &m, const TStr &FName)
Definition: linalg.cpp:1612
int PutStr(const char *CStr)
Definition: fl.cpp:117
static void FillIdentity(TFltVV &M)
Definition: linalg.cpp:1751
double GetSecs() const
Definition: tm.h:366
double GetUniDev()
Definition: dt.h:30
TVal1 Val1
Definition: ds.h:34
static void LUDecomposition(TFltVV &A, TIntV &indx, double &d)
Definition: linalg.cpp:909
TVal2 Val2
Definition: ds.h:35
virtual void PMultiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:147
static double DotProduct(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:165
#define AssertR(Cond, Reason)
Definition: bd.h:258
TStr GetStr() const
Definition: dt.h:1462
void Gen(const TSizeTy &_Vals)
Constructs a vector (an array) of _Vals elements.
Definition: ds.h:523
static void InverseSubstitute(TFltVV &A, const TFltV &p)
Definition: linalg.cpp:843
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TLxSym GetSym(const TFSet &Expect)
Definition: lx.cpp:315
static void Inverse(const TFltVV &A, TFltVV &B, const TLinAlgInverseType &DecompType)
Definition: linalg.cpp:548
static TVec< TFlt, int > GetV(const TFlt &Val1)
Returns a vector on element Val1.
Definition: ds.h:848
static void ToSpVec(const TFltV &Vec, TIntFltKdV &SpVec, const double &CutWordWgtSumPrc=0.0)
Definition: linalg.cpp:1785
static void OrtoIterSVD(const TMatrix &Matrix, int NumSV, int IterN, TFltV &SgnValV)
Definition: linalg.cpp:1021
const char * GetStr() const
Definition: tm.h:368
char * CStr()
Definition: dt.h:479
int PutStrLn(const TStr &Str, const bool &ForceInLn=false)
Definition: fl.h:161
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
static void Gemm(const double &Alpha, const TFltVV &A, const TFltVV &B, const double &Beta, const TFltVV &C, TFltVV &D, const int &TransposeFlags)
Definition: linalg.cpp:493
static void LinComb(const double &p, const TFltV &x, const double &q, const TFltV &y, TFltV &z)
Definition: linalg.cpp:218
static void ConvexComb(const double &p, const TFltV &x, const TFltV &y, TFltV &z)
Definition: linalg.cpp:227
TSizeTy GetCols() const
Definition: ds.h:2253
static void InverseSymetric(TFltVV &A)
Definition: linalg.cpp:872
double Flt
Definition: lx.h:151
void Gen(const TSizeTy &_XDim, const TSizeTy &_YDim)
Definition: ds.h:2247
static void PrintTFltV(const TFltV &Vec, const TStr &VecNm)
Definition: linalg.cpp:1714
virtual void PMultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.cpp:133
static void MultiplyATA(const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
Definition: linalg.cpp:1003
static double EuclDist2(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:298
const TVal & At(const TSizeTy &X, const TSizeTy &Y) const
Definition: ds.h:2256
static double sign(double a, double b)
Definition: linalg.cpp:651