Xmipp  v3.23.11-Nereus
Classes | Functions
Matrix2D Matrices
Collaboration diagram for Matrix2D Matrices:

Classes

class  SparseMatrix2D
 
class  Matrix2D< T >
 

Functions

template<typename T >
void ludcmp (const Matrix2D< T > &A, Matrix2D< T > &LU, Matrix1D< int > &indx, T &d)
 
template<typename T >
void lubksb (const Matrix2D< T > &LU, Matrix1D< int > &indx, Matrix1D< T > &b)
 
void svbksb (Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v, Matrix1D< double > &b, Matrix1D< double > &x)
 
template<typename T >
void svdcmp (const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
 

Matrices speed up macros

typedef Matrix2D< double > DMatrix
 
typedef Matrix2D< int > IMatrix
 
template<typename T >
bool operator== (const Matrix2D< T > &op1, const Matrix2D< T > &op2)
 
#define MATRIX2D_ARRAY(m)   ((m).mdata)
 
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
 
#define MAT_ELEM(m, i, j)   ((m).mdata[((size_t)i)*(m).mdimx+((size_t)j)])
 
#define MAT_XSIZE(m)   ((m).mdimx)
 
#define MAT_YSIZE(m)   ((m).mdimy)
 
#define MAT_SIZE(m)   ((m).mdim)
 
#define dMij(m, i, j)   MAT_ELEM(m, i, j)
 
#define dMn(m, n)   ((m).mdata[(n)])
 
#define M3x3_BY_V3x1(a, M, b)
 
#define M3x3_BY_M3x3(A, B, C)
 
#define M2x2_BY_V2x1(a, M, b)
 
#define M2x2_BY_CT(M2, M1, k)
 
#define M3x3_BY_CT(M2, M1, k)
 
#define M4x4_BY_CT(M2, M1, k)
 
#define M2x2_INV(Ainv, A)
 
#define M3x3_INV(Ainv, A)
 
#define M4x4_INV(Ainv, A)
 

Matrix Related functions

These functions are not methods of Matrix2D

void generalizedEigs (const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix1D< double > &D, Matrix2D< double > &P)
 
void firstEigs (const Matrix2D< double > &A, size_t M, Matrix1D< double > &D, Matrix2D< double > &P, bool Pneeded=true)
 
void lastEigs (const Matrix2D< double > &A, size_t M, Matrix1D< double > &D, Matrix2D< double > &P)
 
void eigsBetween (const Matrix2D< double > &A, size_t I1, size_t I2, Matrix1D< double > &D, Matrix2D< double > &P)
 
void allEigs (const Matrix2D< double > &A, std::vector< std::complex< double > > &eigs)
 
void connectedComponentsOfUndirectedGraph (const Matrix2D< double > &G, Matrix1D< int > &component)
 
template<typename T1 , typename T2 >
void typeCast (const Matrix2D< T1 > &v1, Matrix2D< T2 > &v2)
 
template<typename T1 >
void typeCast (const Matrix2D< T1 > &v1, Matrix2D< T1 > &v2)
 
void orthogonalizeColumnsGramSchmidt (Matrix2D< double > &M)
 
void normalizeColumns (Matrix2D< double > &A)
 
void normalizeColumnsBetween0and1 (Matrix2D< double > &A)
 
void subtractColumnMeans (Matrix2D< double > &A)
 
void matrixOperation_AtA (const Matrix2D< double > &A, Matrix2D< double > &B)
 
void matrixOperation_AAt (const Matrix2D< double > &A, Matrix2D< double > &C)
 
void matrixOperation_AB (const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
 
void matrixOperation_Ax (const Matrix2D< double > &A, const Matrix1D< double > &x, Matrix1D< double > &y)
 
void matrixOperation_ABt (const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
 
void matrixOperation_AtB (const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
 
void matrixOperation_Atx (const Matrix2D< double > &A, const Matrix1D< double > &x, Matrix1D< double > &y)
 
void matrixOperation_AtBt (const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix2D< double > &C)
 
void matrixOperation_XtAX_symmetric (const Matrix2D< double > &X, const Matrix2D< double > &A, Matrix2D< double > &B)
 
void matrixOperation_IplusA (Matrix2D< double > &A)
 
void matrixOperation_IminusA (Matrix2D< double > &A)
 
void eraseFirstColumn (Matrix2D< double > &A)
 
void keepColumns (Matrix2D< double > &A, int j0, int jF)
 
#define VIA_BILIB
 

Detailed Description

Macro Definition Documentation

◆ dMij

#define dMij (   m,
  i,
  j 
)    MAT_ELEM(m, i, j)

Matrix element: Element access

This is just a redefinition of the function above

dMij(m, -2, 1) = 1;
val = dMij(m, -2, 1);

Definition at line 139 of file matrix2d.h.

◆ dMn

#define dMn (   m,
  n 
)    ((m).mdata[(n)])

Matrix element: Element access

This is just a redefinition of the function above

Definition at line 146 of file matrix2d.h.

◆ FOR_ALL_ELEMENTS_IN_MATRIX2D

#define FOR_ALL_ELEMENTS_IN_MATRIX2D (   m)
Value:
for (size_t i=0; i<(m).mdimy; i++) \
for (size_t j=0; j<(m).mdimx; j++)
#define i
#define j
int m

For all elements in the array

This macro is used to generate loops for the matrix in an easy way. It defines internal indexes 'i' and 'j' which ranges the matrix using its mathematical definition (ie, logical access).

{
std::cout << m(i, j) << " ";
}

Definition at line 104 of file matrix2d.h.

◆ M2x2_BY_CT

#define M2x2_BY_CT (   M2,
  M1,
  k 
)
Value:
{ \
dMn(M2, 0) = dMn(M1, 0) * k; \
dMn(M2, 1) = dMn(M1, 1) * k; \
dMn(M2, 2) = dMn(M1, 2) * k; \
dMn(M2, 3) = dMn(M1, 3) * k; }
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define M2
#define M1
#define dMn(m, n)
Definition: matrix2d.h:146

Matrix (2x2) by constant (M2=M1*k)

You must create the result matrix with the appropriate size. You can reuse the matrix M1 to store the results (that is, M2x2_BY_CT(M, M, k);, is allowed).

Definition at line 237 of file matrix2d.h.

◆ M2x2_BY_V2x1

#define M2x2_BY_V2x1 (   a,
  M,
  b 
)
Value:
{ \
spduptmp0 = dMn(M, 0) * XX(b) + dMn(M, 1) * YY(b); \
spduptmp1 = dMn(M, 2) * XX(b) + dMn(M, 3) * YY(b); \
XX(a) = spduptmp0; \
YY(a) = spduptmp1; }
doublereal * b
#define XX(v)
Definition: matrix1d.h:85
#define YY(v)
Definition: matrix1d.h:93
#define dMn(m, n)
Definition: matrix2d.h:146
doublereal * a

Matrix (2x2) by vector (2x1) (a=M*b)

You must "load" the temporary variables, and create the result vector with the appropriate size. You can reuse the vector b to store the results (that is, M2x2_BY_V2x1(b, M, b);, is allowed).

double example
{
M.init_random(0, 1);
b.init_random(0, 1);
return a.sum();
}

Definition at line 225 of file matrix2d.h.

◆ M2x2_INV

#define M2x2_INV (   Ainv,
 
)
Value:
{ \
spduptmp0 = dMn(A,0) * dMn(A,3) - dMn(A,1) * dMn(A,2);\
if (spduptmp0==0.0) \
REPORT_ERROR(ERR_NUMERICAL,"2x2 matrix is not invertible"); \
spduptmp0 = 1.0 / spduptmp0; \
dMn(Ainv, 0) = dMn(A,3); \
dMn(Ainv, 1) = -dMn(A,1); \
dMn(Ainv, 2) = -dMn(A,2); \
dMn(Ainv, 3) = dMn(A,0); \
M2x2_BY_CT(Ainv, Ainv, spduptmp0); }
Error related to numerical calculation.
Definition: xmipp_error.h:179
#define dMn(m, n)
Definition: matrix2d.h:146

Inverse of a matrix (2x2)

Input and output matrix cannot be the same one. The output is supposed to be already resized.

Definition at line 286 of file matrix2d.h.

◆ M3x3_BY_CT

#define M3x3_BY_CT (   M2,
  M1,
  k 
)
Value:
{ \
dMn(M2, 0) = dMn(M1, 0) * k; \
dMn(M2, 1) = dMn(M1, 1) * k; \
dMn(M2, 2) = dMn(M1, 2) * k; \
dMn(M2, 3) = dMn(M1, 3) * k; \
dMn(M2, 4) = dMn(M1, 4) * k; \
dMn(M2, 5) = dMn(M1, 5) * k; \
dMn(M2, 6) = dMn(M1, 6) * k; \
dMn(M2, 7) = dMn(M1, 7) * k; \
dMn(M2, 8) = dMn(M1, 8) * k; }
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define M2
#define M1
#define dMn(m, n)
Definition: matrix2d.h:146

Matrix (3x3) by constant (M2=M1*k)

You must create the result matrix with the appropriate size. You can reuse the matrix M1 to store the results (that is, M2x2_BY_CT(M, M, k);, is allowed).

Definition at line 248 of file matrix2d.h.

◆ M3x3_BY_M3x3

#define M3x3_BY_M3x3 (   A,
  B,
 
)
Value:
{ \
spduptmp0 = dMn(B,0) * dMn(C,0) + dMn(B,1) * dMn(C,3) + dMn(B,2) * dMn(C,6); \
spduptmp1 = dMn(B,0) * dMn(C,1) + dMn(B,1) * dMn(C,4) + dMn(B,2) * dMn(C,7); \
spduptmp2 = dMn(B,0) * dMn(C,2) + dMn(B,1) * dMn(C,5) + dMn(B,2) * dMn(C,8); \
spduptmp3 = dMn(B,3) * dMn(C,0) + dMn(B,4) * dMn(C,3) + dMn(B,5) * dMn(C,6); \
spduptmp4 = dMn(B,3) * dMn(C,1) + dMn(B,4) * dMn(C,4) + dMn(B,5) * dMn(C,7); \
spduptmp5 = dMn(B,3) * dMn(C,2) + dMn(B,4) * dMn(C,5) + dMn(B,5) * dMn(C,8); \
spduptmp6 = dMn(B,6) * dMn(C,0) + dMn(B,7) * dMn(C,3) + dMn(B,8) * dMn(C,6); \
spduptmp7 = dMn(B,6) * dMn(C,1) + dMn(B,7) * dMn(C,4) + dMn(B,8) * dMn(C,7); \
spduptmp8 = dMn(B,6) * dMn(C,2) + dMn(B,7) * dMn(C,5) + dMn(B,8) * dMn(C,8); \
dMn(A, 0) = spduptmp0; \
dMn(A, 1) = spduptmp1; \
dMn(A, 2) = spduptmp2; \
dMn(A, 3) = spduptmp3; \
dMn(A, 4) = spduptmp4; \
dMn(A, 5) = spduptmp5; \
dMn(A, 6) = spduptmp6; \
dMn(A, 7) = spduptmp7; \
dMn(A, 8) = spduptmp8; }
#define dMn(m, n)
Definition: matrix2d.h:146

Matrix (3x3) by Matrix (3x3) (A=B*C)

You must "load" the temporary variables, and create the result vector with the appropriate size. You can reuse any of the multiplicands to store the results (that is, M3x3_BY_M3x3(A, A, B);, is allowed).

Definition at line 182 of file matrix2d.h.

◆ M3x3_BY_V3x1

#define M3x3_BY_V3x1 (   a,
  M,
  b 
)
Value:
{ \
spduptmp0 = dMn(M, 0) * XX(b) + dMn(M, 1) * YY(b) + dMn(M, 2) * ZZ(b); \
spduptmp1 = dMn(M, 3) * XX(b) + dMn(M, 4) * YY(b) + dMn(M, 5) * ZZ(b); \
spduptmp2 = dMn(M, 6) * XX(b) + dMn(M, 7) * YY(b) + dMn(M, 8) * ZZ(b); \
XX(a) = spduptmp0; YY(a) = spduptmp1; ZZ(a) = spduptmp2; }
doublereal * b
#define XX(v)
Definition: matrix1d.h:85
#define YY(v)
Definition: matrix1d.h:93
#define dMn(m, n)
Definition: matrix2d.h:146
#define ZZ(v)
Definition: matrix1d.h:101
doublereal * a

Matrix (3x3) by vector (3x1) (a=M*b)

You must "load" the temporary variables, and create the result vector with the appropriate size. You can reuse the vector b to store the results (that is, M3x3_BY_V3x1(b, M, b);, is allowed).

double example
{
M.init_random(0, 1);
b.init_random(0, 1);
return a.sum();
}

Definition at line 170 of file matrix2d.h.

◆ M3x3_INV

#define M3x3_INV (   Ainv,
 
)
Value:
{ \
dMn(Ainv, 0) = dMn(A,8)*dMn(A,4)-dMn(A,7)*dMn(A,5); \
dMn(Ainv, 1) = -(dMn(A,8)*dMn(A,1)-dMn(A,7)*dMn(A,2)); \
dMn(Ainv, 2) = dMn(A,5)*dMn(A,1)-dMn(A,4)*dMn(A,2); \
dMn(Ainv, 3) = -(dMn(A,8)*dMn(A,3)-dMn(A,6)*dMn(A,5)); \
dMn(Ainv, 4) = dMn(A,8)*dMn(A,0)-dMn(A,6)*dMn(A,2); \
dMn(Ainv, 5) = -(dMn(A,5)*dMn(A,0)-dMn(A,3)*dMn(A,2)); \
dMn(Ainv, 6) = dMn(A,7)*dMn(A,3)-dMn(A,6)*dMn(A,4); \
dMn(Ainv, 7) = -(dMn(A,7)*dMn(A,0)-dMn(A,6)*dMn(A,1)); \
dMn(Ainv, 8) = dMn(A,4)*dMn(A,0)-dMn(A,3)*dMn(A,1); \
spduptmp0 = dMn(A,0)*dMn(Ainv,0)+dMn(A,3)*dMn(Ainv,1)+dMn(A,6)*dMn(Ainv,2); \
if (spduptmp0==0.0) \
REPORT_ERROR(ERR_NUMERICAL,"3x3 matrix is not invertible"); \
spduptmp0 = 1.0 / spduptmp0; \
M3x3_BY_CT(Ainv, Ainv, spduptmp0); }
Error related to numerical calculation.
Definition: xmipp_error.h:179
#define dMn(m, n)
Definition: matrix2d.h:146

Inverse of a matrix (3x3)

Input and output matrix cannot be the same one. The output is supposed to be already resized.

Definition at line 302 of file matrix2d.h.

◆ M4x4_BY_CT

#define M4x4_BY_CT (   M2,
  M1,
  k 
)
Value:
{ \
dMn(M2, 0) = dMn(M1, 0) * k; \
dMn(M2, 1) = dMn(M1, 1) * k; \
dMn(M2, 2) = dMn(M1, 2) * k; \
dMn(M2, 3) = dMn(M1, 3) * k; \
dMn(M2, 4) = dMn(M1, 4) * k; \
dMn(M2, 5) = dMn(M1, 5) * k; \
dMn(M2, 6) = dMn(M1, 6) * k; \
dMn(M2, 7) = dMn(M1, 7) * k; \
dMn(M2, 8) = dMn(M1, 8) * k; \
dMn(M2, 9) = dMn(M1, 9) * k; \
dMn(M2,10) = dMn(M1,10) * k; \
dMn(M2,11) = dMn(M1,11) * k; \
dMn(M2,12) = dMn(M1,12) * k; \
dMn(M2,13) = dMn(M1,13) * k; \
dMn(M2,14) = dMn(M1,14) * k; \
dMn(M2,15) = dMn(M1,15) * k;}
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define M2
#define M1
#define dMn(m, n)
Definition: matrix2d.h:146

Matrix (4x4) by constant (M2=M1*k)

You must create the result matrix with the appropriate size. You can reuse the matrix M1 to store the results (that is, M2x2_BY_CT(M, M, k);, is allowed).

Definition at line 263 of file matrix2d.h.

◆ M4x4_INV

#define M4x4_INV (   Ainv,
 
)

Inverse of a matrix (4x4)

Input and output matrix cannot be the same one. The output is supposed to be already resized.

Definition at line 323 of file matrix2d.h.

◆ MAT_ELEM

#define MAT_ELEM (   m,
  i,
  j 
)    ((m).mdata[((size_t)i)*(m).mdimx+((size_t)j)])

Access to a matrix element v is the array, i and j define the element v_ij.

MAT_ELEM(m, 0, 0) = 1;
val = MAT_ELEM(m, 0, 0);

Definition at line 116 of file matrix2d.h.

◆ MAT_SIZE

#define MAT_SIZE (   m)    ((m).mdim)

Total elements of the matrix

Definition at line 128 of file matrix2d.h.

◆ MAT_XSIZE

#define MAT_XSIZE (   m)    ((m).mdimx)

X dimension of the matrix

Definition at line 120 of file matrix2d.h.

◆ MAT_YSIZE

#define MAT_YSIZE (   m)    ((m).mdimy)

Y dimension of the matrix

Definition at line 124 of file matrix2d.h.

◆ MATRIX2D_ARRAY

#define MATRIX2D_ARRAY (   m)    ((m).mdata)

Array access.

This macro gives you access to the array (T)

Definition at line 89 of file matrix2d.h.

◆ VIA_BILIB

#define VIA_BILIB

Definition at line 1196 of file matrix2d.h.

Typedef Documentation

◆ DMatrix

typedef Matrix2D<double> DMatrix

Definition at line 1163 of file matrix2d.h.

◆ IMatrix

typedef Matrix2D<int> IMatrix

Definition at line 1164 of file matrix2d.h.

Function Documentation

◆ allEigs()

void allEigs ( const Matrix2D< double > &  A,
std::vector< std::complex< double > > &  eigs 
)

Compute all eigenvalues, even if they are complex

Definition at line 345 of file matrix2d.cpp.

346 {
347  int N=(int)MAT_YSIZE(A);
348  alglib::real_2d_array a, vl, vr;
349  a.setcontent(N,N,MATRIX2D_ARRAY(A));
350  alglib::real_1d_array wr, wi;
351 
352  bool ok=rmatrixevd(a, N, 0, wr, wi, vl, vr);
353  if (!ok)
354  REPORT_ERROR(ERR_NUMERICAL,"Could not perform eigenvector decomposition");
355  eigs.clear();
356  for (int n=0; n<N; ++n)
357  eigs.push_back(std::complex<double>(wr(n),wi(n)));
358 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
bool rmatrixevd(const real_2d_array &a, const ae_int_t n, const ae_int_t vneeded, real_1d_array &wr, real_1d_array &wi, real_2d_array &vl, real_2d_array &vr)
Definition: linalg.cpp:2468
Error related to numerical calculation.
Definition: xmipp_error.h:179
void setcontent(ae_int_t irows, ae_int_t icols, const double *pContent)
Definition: ap.cpp:6660
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
int * n
doublereal * a

◆ connectedComponentsOfUndirectedGraph()

void connectedComponentsOfUndirectedGraph ( const Matrix2D< double > &  G,
Matrix1D< int > &  component 
)

Find connected components of a graph. Assuming that the matrix G represents an undirected graph (of size NxN), this function returns a vector (of size N) that indicates for each element which is the number of its connected component.

Definition at line 360 of file matrix2d.cpp.

361 {
362  size_t N=MAT_XSIZE(G);
363  component.resizeNoCopy(N);
364  component.initConstant(-1);
365 
366  int nextComponent=0;
367  bool workDone=false;
368  std::queue<size_t> toExplore;
369  do
370  {
371  workDone=false;
372  // Find next unvisited element
373  bool found=false;
374  size_t seed=0;
376  if (VEC_ELEM(component,i)<0)
377  {
378  seed=i;
379  found=true;
380  break;
381  }
382 
383  // If found, get its connected component
384  if (found)
385  {
386  int currentComponent=nextComponent;
387  nextComponent++;
388 
389  VEC_ELEM(component,seed)=currentComponent;
390  toExplore.push(seed);
391  while (toExplore.size()>0)
392  {
393  seed=toExplore.front();
394  toExplore.pop();
395  for (size_t j=seed+1; j<N; ++j)
396  if (MAT_ELEM(G,seed,j)>0)
397  {
398  if (VEC_ELEM(component,j)<0)
399  {
400  VEC_ELEM(component,j)=currentComponent;
401  toExplore.push(j);
402  }
403 
404  }
405  }
406  workDone=true;
407  }
408  } while (workDone);
409 }
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void initConstant(T val)
Definition: matrix1d.cpp:83

◆ eigsBetween()

void eigsBetween ( const Matrix2D< double > &  A,
size_t  I1,
size_t  I2,
Matrix1D< double > &  D,
Matrix2D< double > &  P 
)

Compute eigenvectors between two indexes of a real, symmetric matrix. Solves the problem Av=dv. Only the eigenvectors of the smallest eigenvalues between indexes I1 and I2 are returned as columns of P. Indexes start at 0.

Definition at line 324 of file matrix2d.cpp.

325 {
326  size_t M = I2 - I1 + 1;
327  int N=(int)MAT_YSIZE(A);
329  a.setcontent(N,N,MATRIX2D_ARRAY(A));
331 
332  bool ok=smatrixevdi(a, N, true, false, I1, I2, d, z);
333  if (!ok)
334  REPORT_ERROR(ERR_NUMERICAL,"Could not perform eigenvector decomposition");
335 
336  D.resizeNoCopy(M);
338  VEC_ELEM(D,i)=d(M-1-i);
339  memcpy(&VEC_ELEM(D,0),d.getcontent(),M*sizeof(double));
340  P.resizeNoCopy(N,M);
342  MAT_ELEM(P,i,j)=z(i,j);
343 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
#define i
doublereal * d
bool smatrixevdi(const real_2d_array &a, const ae_int_t n, const ae_int_t zneeded, const bool isupper, const ae_int_t i1, const ae_int_t i2, real_1d_array &w, real_2d_array &z)
Definition: linalg.cpp:2009
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
Error related to numerical calculation.
Definition: xmipp_error.h:179
double z
#define j
double * getcontent()
Definition: ap.cpp:6207
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void setcontent(ae_int_t irows, ae_int_t icols, const double *pContent)
Definition: ap.cpp:6660
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
doublereal * a

◆ eraseFirstColumn()

void eraseFirstColumn ( Matrix2D< double > &  A)

Erase first column

Definition at line 542 of file matrix2d.cpp.

543 {
544  Matrix2D<double> Ap;
545  Ap.resize(MAT_YSIZE(A),MAT_XSIZE(A)-1);
546  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
547  memcpy(&MAT_ELEM(Ap,i,0),&MAT_ELEM(A,i,1),MAT_XSIZE(Ap)*sizeof(double));
548  A=Ap;
549 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void resize(size_t Ydim, size_t Xdim, bool noCopy=false)
Definition: matrix2d.cpp:1022

◆ firstEigs()

void firstEigs ( const Matrix2D< double > &  A,
size_t  M,
Matrix1D< double > &  D,
Matrix2D< double > &  P,
bool  Pneeded = true 
)

First eigenvectors of a real, symmetric matrix. Solves the problem Av=dv. Only the eigenvectors of the largest M eigenvalues are returned as columns of P

Definition at line 284 of file matrix2d.cpp.

285 {
286  int N=(int)MAT_YSIZE(A);
288  a.setcontent(N,N,MATRIX2D_ARRAY(A));
290  bool ok=smatrixevdi(a, N, Pneeded, false, N-M, N-1, d, z);
291  if (!ok)
292  REPORT_ERROR(ERR_NUMERICAL,"Could not perform eigenvector decomposition");
293 
294  D.resizeNoCopy(M);
296  VEC_ELEM(D,i)=d(M-1-i);
297  if (Pneeded)
298  {
299  P.resizeNoCopy(N,M);
301  MAT_ELEM(P,i,j)=z(i,M-1-j);
302  }
303 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
#define i
doublereal * d
bool smatrixevdi(const real_2d_array &a, const ae_int_t n, const ae_int_t zneeded, const bool isupper, const ae_int_t i1, const ae_int_t i2, real_1d_array &w, real_2d_array &z)
Definition: linalg.cpp:2009
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
Error related to numerical calculation.
Definition: xmipp_error.h:179
double z
#define j
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void setcontent(ae_int_t irows, ae_int_t icols, const double *pContent)
Definition: ap.cpp:6660
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
doublereal * a

◆ generalizedEigs()

void generalizedEigs ( const Matrix2D< double > &  A,
const Matrix2D< double > &  B,
Matrix1D< double > &  D,
Matrix2D< double > &  P 
)

Generalized eigenvector decomposition. Solves the problem Av=dBv. The decomposition is such that A=B P D P^-1. A and B must be square matrices of the same size.

Definition at line 267 of file matrix2d.cpp.

268 {
269  int N=(int)MAT_YSIZE(A);
271  a.setcontent(N,N,MATRIX2D_ARRAY(A));
272  b.setcontent(N,N,MATRIX2D_ARRAY(B));
274  bool ok=smatrixgevd(a, N, true, b, true, true, 1, d, z);
275  if (!ok)
276  REPORT_ERROR(ERR_NUMERICAL,"Could not perform eigenvector decomposition");
277  D.resizeNoCopy(N);
278  memcpy(&VEC_ELEM(D,0),d.getcontent(),N*sizeof(double));
279  P.resizeNoCopy(A);
281  MAT_ELEM(P,i,j)=z(i,j);
282 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
bool smatrixgevd(const real_2d_array &a, const ae_int_t n, const bool isuppera, const real_2d_array &b, const bool isupperb, const ae_int_t zneeded, const ae_int_t problemtype, real_1d_array &d, real_2d_array &z)
Definition: linalg.cpp:6778
#define i
doublereal * d
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
doublereal * b
Error related to numerical calculation.
Definition: xmipp_error.h:179
double z
#define j
double * getcontent()
Definition: ap.cpp:6207
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void setcontent(ae_int_t irows, ae_int_t icols, const double *pContent)
Definition: ap.cpp:6660
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
doublereal * a

◆ keepColumns()

void keepColumns ( Matrix2D< double > &  A,
int  j0,
int  jF 
)

Keep columns between j0 and jF

Definition at line 551 of file matrix2d.cpp.

552 {
553  Matrix2D<double> Ap;
554  Ap.resize(MAT_YSIZE(A),jF-j0+1);
555  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
556  memcpy(&MAT_ELEM(Ap,i,0),&MAT_ELEM(A,i,j0),MAT_XSIZE(Ap)*sizeof(double));
557  A=Ap;
558 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void resize(size_t Ydim, size_t Xdim, bool noCopy=false)
Definition: matrix2d.cpp:1022

◆ lastEigs()

void lastEigs ( const Matrix2D< double > &  A,
size_t  M,
Matrix1D< double > &  D,
Matrix2D< double > &  P 
)

Last eigenvectors of a real, symmetric matrix. Solves the problem Av=dv. Only the eigenvectors of the smallest M eigenvalues are returned as columns of P

Definition at line 305 of file matrix2d.cpp.

306 {
307  int N=(int)MAT_YSIZE(A);
309  a.setcontent(N,N,MATRIX2D_ARRAY(A));
311  bool ok=smatrixevdi(a, N, true, false, 0, M-1, d, z);
312  if (!ok)
313  REPORT_ERROR(ERR_NUMERICAL,"Could not perform eigenvector decomposition");
314 
315  D.resizeNoCopy(M);
317  VEC_ELEM(D,i)=d(M-1-i);
318  memcpy(&VEC_ELEM(D,0),d.getcontent(),M*sizeof(double));
319  P.resizeNoCopy(N,M);
321  MAT_ELEM(P,i,j)=z(i,j);
322 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
#define i
doublereal * d
bool smatrixevdi(const real_2d_array &a, const ae_int_t n, const ae_int_t zneeded, const bool isupper, const ae_int_t i1, const ae_int_t i2, real_1d_array &w, real_2d_array &z)
Definition: linalg.cpp:2009
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
Error related to numerical calculation.
Definition: xmipp_error.h:179
double z
#define j
double * getcontent()
Definition: ap.cpp:6207
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void setcontent(ae_int_t irows, ae_int_t icols, const double *pContent)
Definition: ap.cpp:6660
#define MATRIX2D_ARRAY(m)
Definition: matrix2d.h:89
doublereal * a

◆ lubksb()

template<typename T >
void lubksb ( const Matrix2D< T > &  LU,
Matrix1D< int > &  indx,
Matrix1D< T > &  b 
)

LU Backsubstitution

Definition at line 596 of file matrix2d.cpp.

597 {
601 }
size_t size() const
Definition: matrix1d.h:508
T * adaptForNumericalRecipes2() const
Definition: matrix2d.h:875
void lubksb(const Matrix2D< T > &LU, Matrix1D< int > &indx, Matrix1D< T > &b)
Definition: matrix2d.cpp:596
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844

◆ ludcmp()

template<typename T >
void ludcmp ( const Matrix2D< T > &  A,
Matrix2D< T > &  LU,
Matrix1D< int > &  indx,
T &  d 
)

LU Decomposition

Definition at line 586 of file matrix2d.cpp.

587 {
588  LU = A;
589  if (VEC_XSIZE(indx)!=A.mdimx)
590  indx.resizeNoCopy(A.mdimx);
592  indx.adaptForNumericalRecipes(), &d);
593 }
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
T * adaptForNumericalRecipes2() const
Definition: matrix2d.h:875
doublereal * d
void ludcmp(const Matrix2D< T > &A, Matrix2D< T > &LU, Matrix1D< int > &indx, T &d)
Definition: matrix2d.cpp:586
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
size_t mdimx
Definition: matrix2d.h:410
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844

◆ matrixOperation_AAt()

void matrixOperation_AAt ( const Matrix2D< double > &  A,
Matrix2D< double > &  C 
)

Matrix operation: C=A*A^t.

Definition at line 449 of file matrix2d.cpp.

450 {
451  C.initZeros(MAT_YSIZE(A), MAT_YSIZE(A));
452  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
453  for (size_t j = i; j < MAT_YSIZE(A); ++j)
454  {
455  double aux=0.;
456  for (size_t k = 0; k < MAT_XSIZE(A); ++k)
457  aux += MAT_ELEM(A, i, k) * MAT_ELEM(A, j, k);
458  MAT_ELEM(C, j, i)=MAT_ELEM(C, i, j)=aux;
459  }
460 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void initZeros()
Definition: matrix2d.h:626

◆ matrixOperation_AB()

void matrixOperation_AB ( const Matrix2D< double > &  A,
const Matrix2D< double > &  B,
Matrix2D< double > &  C 
)

Matrix operation: C=A*B.

Definition at line 411 of file matrix2d.cpp.

412 {
413  C.initZeros(MAT_YSIZE(A), MAT_XSIZE(B));
414  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
415  for (size_t j = 0; j < MAT_XSIZE(B); ++j)
416  {
417  double aux=0.;
418  for (size_t k = 0; k < MAT_XSIZE(A); ++k)
419  aux += MAT_ELEM(A, i, k) * MAT_ELEM(B, k, j);
420  MAT_ELEM(C, i, j)=aux;
421  }
422 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void initZeros()
Definition: matrix2d.h:626

◆ matrixOperation_ABt()

void matrixOperation_ABt ( const Matrix2D< double > &  A,
const Matrix2D< double > &  B,
Matrix2D< double > &  C 
)

Matrix operation: C=A*B^t.

Definition at line 462 of file matrix2d.cpp.

463 {
464  C.initZeros(MAT_YSIZE(A), MAT_YSIZE(B));
465  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
466  for (size_t j = 0; j < MAT_YSIZE(B); ++j)
467  {
468  double aux=0.;
469  for (size_t k = 0; k < MAT_XSIZE(A); ++k)
470  aux += MAT_ELEM(A, i, k) * MAT_ELEM(B, j, k);
471  MAT_ELEM(C, i, j)=aux;
472  }
473 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void initZeros()
Definition: matrix2d.h:626

◆ matrixOperation_AtA()

void matrixOperation_AtA ( const Matrix2D< double > &  A,
Matrix2D< double > &  B 
)

Matrix operation: B=A^t*A.

Definition at line 436 of file matrix2d.cpp.

437 {
438  B.resizeNoCopy(MAT_XSIZE(A), MAT_XSIZE(A));
439  for (size_t i = 0; i < MAT_XSIZE(A); ++i)
440  for (size_t j = i; j < MAT_XSIZE(A); ++j)
441  {
442  double aux=0.;
443  for (size_t k = 0; k < MAT_YSIZE(A); ++k)
444  aux += MAT_ELEM(A, k, i) * MAT_ELEM(A, k, j);
445  MAT_ELEM(B, j, i) = MAT_ELEM(B, i, j) = aux;
446  }
447 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ matrixOperation_AtB()

void matrixOperation_AtB ( const Matrix2D< double > &  A,
const Matrix2D< double > &  B,
Matrix2D< double > &  C 
)

Matrix operation: C=A^t*B.

Definition at line 475 of file matrix2d.cpp.

476 {
477  C.resizeNoCopy(MAT_XSIZE(A), MAT_XSIZE(B));
478  for (size_t i = 0; i < MAT_XSIZE(A); ++i)
479  for (size_t j = 0; j < MAT_XSIZE(B); ++j)
480  {
481  double aux=0.;
482  for (size_t k = 0; k < MAT_YSIZE(A); ++k)
483  aux += MAT_ELEM(A, k, i) * MAT_ELEM(B, k, j);
484  MAT_ELEM(C, i, j)=aux;
485  }
486 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ matrixOperation_AtBt()

void matrixOperation_AtBt ( const Matrix2D< double > &  A,
const Matrix2D< double > &  B,
Matrix2D< double > &  C 
)

Matrix operation: C=A^t*Bt.

Definition at line 500 of file matrix2d.cpp.

501 {
502  C.initZeros(MAT_XSIZE(A), MAT_YSIZE(B));
503  for (size_t i = 0; i < MAT_XSIZE(A); ++i)
504  for (size_t j = 0; j < MAT_YSIZE(B); ++j)
505  {
506  double aux=0.;
507  for (size_t k = 0; k < MAT_YSIZE(A); ++k)
508  aux += MAT_ELEM(A, k, i) * MAT_ELEM(B, j, k);
509  MAT_ELEM(C, i, j)=aux;
510  }
511 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void initZeros()
Definition: matrix2d.h:626

◆ matrixOperation_Atx()

void matrixOperation_Atx ( const Matrix2D< double > &  A,
const Matrix1D< double > &  x,
Matrix1D< double > &  y 
)

Matrix operation: y=A^t*x.

Definition at line 488 of file matrix2d.cpp.

489 {
490  y.resizeNoCopy(MAT_XSIZE(A));
491  for (size_t i = 0; i < MAT_XSIZE(A); ++i)
492  {
493  double aux=0.;
494  for (size_t k = 0; k < MAT_YSIZE(A); ++k)
495  aux += MAT_ELEM(A, k, i) * VEC_ELEM(x, k);
496  VEC_ELEM(y, i)=aux;
497  }
498 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458

◆ matrixOperation_Ax()

void matrixOperation_Ax ( const Matrix2D< double > &  A,
const Matrix1D< double > &  x,
Matrix1D< double > &  y 
)

Matrix operation: y=A*x.

Definition at line 424 of file matrix2d.cpp.

425 {
426  y.initZeros(MAT_YSIZE(A));
427  for (size_t i = 0; i < MAT_YSIZE(A); ++i)
428  {
429  double aux=0.;
430  for (size_t k = 0; k < MAT_XSIZE(A); ++k)
431  aux += MAT_ELEM(A, i, k) * VEC_ELEM(x, k);
432  VEC_ELEM(y, i)=aux;
433  }
434 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
void initZeros()
Definition: matrix1d.h:592
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ matrixOperation_IminusA()

void matrixOperation_IminusA ( Matrix2D< double > &  A)

Matrix operation: A=I-A

Definition at line 533 of file matrix2d.cpp.

534 {
536  if (i == j)
537  MAT_ELEM(A, i, j) = 1 - MAT_ELEM(A, i, j);
538  else
539  MAT_ELEM(A, i, j) = -MAT_ELEM(A, i, j);
540 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j

◆ matrixOperation_IplusA()

void matrixOperation_IplusA ( Matrix2D< double > &  A)

Matrix operation: A=I+A

Definition at line 527 of file matrix2d.cpp.

528 {
529  for (size_t i=0; i<MAT_YSIZE(A); ++i)
530  MAT_ELEM(A,i,i)+=1;
531 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116

◆ matrixOperation_XtAX_symmetric()

void matrixOperation_XtAX_symmetric ( const Matrix2D< double > &  X,
const Matrix2D< double > &  A,
Matrix2D< double > &  B 
)

Matrix operation: B=X^t*A*X. We know that the result B must be symmetric

Definition at line 513 of file matrix2d.cpp.

514 {
515  Matrix2D<double> AX=A*X;
516  B.resizeNoCopy(MAT_XSIZE(X), MAT_XSIZE(X));
517  for (size_t i = 0; i < MAT_XSIZE(X); ++i)
518  for (size_t j = i; j < MAT_XSIZE(X); ++j)
519  {
520  double aux=0.;
521  for (size_t k = 0; k < MAT_YSIZE(X); ++k)
522  aux += MAT_ELEM(X, k, i) * MAT_ELEM(AX, k, j);
523  MAT_ELEM(B, j, i) = MAT_ELEM(B, i, j) = aux;
524  }
525 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ normalizeColumns()

void normalizeColumns ( Matrix2D< double > &  A)

Normalize columns. So that they have zero mean and unit variance.

Definition at line 188 of file matrix2d.cpp.

189 {
190  if (MAT_YSIZE(A)<=1)
191  return;
192 
193  // Compute the mean and standard deviation of each column
194  Matrix1D<double> avg, stddev;
195  avg.initZeros(MAT_XSIZE(A));
196  stddev.initZeros(MAT_XSIZE(A));
197 
199  {
200  double x=MAT_ELEM(A,i,j);
201  VEC_ELEM(avg,j)+=x;
202  VEC_ELEM(stddev,j)+=x*x;
203  }
204 
205  double iN=1.0/MAT_YSIZE(A);
207  {
208  VEC_ELEM(avg,i)*=iN;
209  VEC_ELEM(stddev,i)=sqrt(fabs(VEC_ELEM(stddev,i)*iN - VEC_ELEM(avg,i)*VEC_ELEM(avg,i)));
210  if (VEC_ELEM(stddev,i)>XMIPP_EQUAL_ACCURACY)
211  VEC_ELEM(stddev,i)=1.0/VEC_ELEM(stddev,i);
212  else
213  VEC_ELEM(stddev,i)=0.0;
214  }
215 
216  // Now normalize
218  MAT_ELEM(A,i,j)=(MAT_ELEM(A,i,j)-VEC_ELEM(avg,j))*VEC_ELEM(stddev,j);
219 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
void sqrt(Image< double > &op)
doublereal * x
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
void initZeros()
Definition: matrix1d.h:592
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ normalizeColumnsBetween0and1()

void normalizeColumnsBetween0and1 ( Matrix2D< double > &  A)

Normalize columns. So that the minimum is 0 and the maximum is 1

Definition at line 221 of file matrix2d.cpp.

222 {
223  double maxValue,minValue;
224  A.computeMaxAndMin(maxValue,minValue);
225  double iMaxValue=1.0/(maxValue-minValue);
227  MAT_ELEM(A,i,j)=(MAT_ELEM(A,i,j)-minValue)*iMaxValue;
228 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
void computeMaxAndMin(T &maxValue, T &minValue) const
Definition: matrix2d.cpp:1262
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define j

◆ operator==()

template<typename T >
bool operator== ( const Matrix2D< T > &  op1,
const Matrix2D< T > &  op2 
)

Definition at line 1167 of file matrix2d.h.

1168 {
1169  return op1.equal(op2);
1170 }
bool equal(const Matrix2D< T > &op, double accuracy=XMIPP_EQUAL_ACCURACY) const
Definition: matrix2d.cpp:1202

◆ orthogonalizeColumnsGramSchmidt()

void orthogonalizeColumnsGramSchmidt ( Matrix2D< double > &  M)

Gram Schmidt orthogonalization by columns

Definition at line 560 of file matrix2d.cpp.

561 {
562  for(size_t j1=0; j1<MAT_XSIZE(M); j1++)
563  {
564  // Normalize column j1
565  double norm=0;
566  for (size_t i=0; i<MAT_YSIZE(M); i++)
567  norm+=MAT_ELEM(M,i,j1)*MAT_ELEM(M,i,j1);
568  double K=1.0/sqrt(norm);
569  for(size_t i = 0; i<MAT_YSIZE(M); i++)
570  MAT_ELEM(M,i,j1) *=K;
571 
572  // Update rest of columns
573  for(size_t j2=j1+1; j2<MAT_XSIZE(M); j2++)
574  {
575  // Compute the dot product
576  double K = 0;
577  for (size_t i=0; i<MAT_YSIZE(M); i++)
578  K+=MAT_ELEM(M,i,j1)*MAT_ELEM(M,i,j2);
579  for (size_t i=0; i<MAT_YSIZE(M); i++)
580  MAT_ELEM(M,i,j2) -= K*MAT_ELEM(M,i,j1);
581  }
582  }
583 }
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
void sqrt(Image< double > &op)
T norm(const std::vector< T > &v)
Definition: vector_ops.h:399
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define MAT_XSIZE(m)
Definition: matrix2d.h:120
constexpr int K

◆ subtractColumnMeans()

void subtractColumnMeans ( Matrix2D< double > &  A)

Subtract mean of columns. So that they have zero mean.

Definition at line 230 of file matrix2d.cpp.

231 {
232  if (MAT_YSIZE(A)<1)
233  return;
234 
235  // Compute the mean and standard deviation of each column
236  Matrix1D<double> avg;
237  avg.initZeros(MAT_XSIZE(A));
238 
240  VEC_ELEM(avg,j)+=MAT_ELEM(A,i,j);
241 
242  double iN=1.0/MAT_YSIZE(A);
244  VEC_ELEM(avg,i)*=iN;
245 
246  // Now normalize
248  MAT_ELEM(A,i,j)=MAT_ELEM(A,i,j)-VEC_ELEM(avg,j);
249 }
#define FOR_ALL_ELEMENTS_IN_MATRIX2D(m)
Definition: matrix2d.h:104
#define MAT_YSIZE(m)
Definition: matrix2d.h:124
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define i
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v)
Definition: matrix1d.h:72
void initZeros()
Definition: matrix1d.h:592
#define j
#define MAT_XSIZE(m)
Definition: matrix2d.h:120

◆ svbksb()

void svbksb ( Matrix2D< double > &  u,
Matrix1D< double > &  w,
Matrix2D< double > &  v,
Matrix1D< double > &  b,
Matrix1D< double > &  x 
)

SVD Backsubstitution

Definition at line 174 of file matrix2d.cpp.

176 {
177  // Call to the numerical recipes routine. Results will be stored in X
181  u.mdimy, u.mdimx,
184 }
T * adaptForNumericalRecipes2() const
Definition: matrix2d.h:875
size_t mdimy
Definition: matrix2d.h:413
size_t mdimx
Definition: matrix2d.h:410
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
void svbksb(Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v, Matrix1D< double > &b, Matrix1D< double > &x)
Definition: matrix2d.cpp:174

◆ svdcmp()

template<typename T >
void svdcmp ( const Matrix2D< T > &  a,
Matrix2D< double > &  u,
Matrix1D< double > &  w,
Matrix2D< double > &  v 
)

SVD Decomposition

Definition at line 125 of file matrix2d.cpp.

129 {
130  // svdcmp only works with double
131  typeCast(a, u);
132 
133  // Set size of matrices
134  w.initZeros(u.mdimx);
135  v.initZeros(u.mdimx, u.mdimx);
136 
137  // Call to the numerical recipes routine
138 #ifdef VIA_NR
139 
140  svdcmp(u.mdata,
141  u.mdimy, u.mdimx,
142  w.vdata,
143  v.mdata);
144 #endif
145 
146 #ifdef VIA_BILIB
147 
148  int status;
150  u.mdimy, u.mdimx,
151  w.vdata,
152  v.mdata,
153  5000, &status);
154 #endif
155 }
T * mdata
Definition: matrix2d.h:395
size_t mdimy
Definition: matrix2d.h:413
void initZeros()
Definition: matrix1d.h:592
void svdcmp(const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
Definition: matrix2d.cpp:125
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
int SingularValueDecomposition(double *U, long Lines, long Columns, double W[], double *V, long MaxIterations, int *Status)
void initZeros()
Definition: matrix2d.h:626
size_t mdimx
Definition: matrix2d.h:410
T * vdata
The array itself.
Definition: matrix1d.h:258

◆ typeCast() [1/2]

template<typename T1 , typename T2 >
void typeCast ( const Matrix2D< T1 > &  v1,
Matrix2D< T2 > &  v2 
)

Conversion from one type to another.

If we have an integer array and we need a double one, we can use this function. The conversion is done through a type casting of each element If n >= 0, only the nth volumes will be converted, otherwise all NSIZE volumes

Definition at line 1245 of file matrix2d.h.

1246 {
1247  if (v1.mdim == 0)
1248  {
1249  v2.clear();
1250  return;
1251  }
1252 
1253  if (v1.mdimx!=v2.mdimx || v1.mdimy!=v2.mdimy)
1254  v2.resizeNoCopy(v1);
1255  for (unsigned long int n = 0; n < v1.mdim; n++)
1256  v2.mdata[n] = static_cast< T2 > (v1.mdata[n]);
1257 }
T * mdata
Definition: matrix2d.h:395
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
void clear()
Definition: matrix2d.h:473
size_t mdimy
Definition: matrix2d.h:413
size_t mdim
Definition: matrix2d.h:416
size_t mdimx
Definition: matrix2d.h:410
int * n

◆ typeCast() [2/2]

template<typename T1 >
void typeCast ( const Matrix2D< T1 > &  v1,
Matrix2D< T1 > &  v2 
)

Conversion from one type to another. In some cases, the two types are the same. So a faster way is simply by assignment.

Definition at line 1263 of file matrix2d.h.

1264 {
1265  v2=v1;
1266 }
double v1