SNAP Library 4.0, Developer Reference  2017-07-27 13:18:06
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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(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(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:2245
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:1441
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:1385
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:1383
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:1180
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:2251
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:1427
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:1459
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, TSizeTy > 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:476
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:2252
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:2246
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:2255
static double sign(double a, double b)
Definition: linalg.cpp:651