test_utils.cc

Go to the documentation of this file.

/* Copyright (C) 2001-2012 Stefan Buehler <sbuehler@ltu.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#include "test_utils.h"
#include <cmath>
#include "arts.h"
#include "lin_alg.h"
#include "matpackII.h"

using std::abs;


void add_noise(VectorView v, Numeric scale) {
  Rand<Numeric> rand(0, scale);

  for (Index i = 0; i < v.nelem(); i++) {
    v[i] += rand();
  }
}


void random_fill_matrix(MatrixView A, Numeric range, bool positive) {
  Index m = A.nrows();
  Index n = A.ncols();

  Rand<Numeric> rand(positive ? 0 : -range, range);

  for (Index i = 0; i < m; i++) {
    for (Index j = 0; j < n; j++) {
      A(i, j) = (Numeric)rand();
    }
  }
}
void random_fill_matrix(ComplexMatrixView A, Numeric range, bool positive) {
  Index m = A.nrows();
  Index n = A.ncols();

  Rand<Numeric> rand(positive ? 0 : -range, range);

  for (Index i = 0; i < m; i++) {
    for (Index j = 0; j < n; j++) {
      A(i, j) = Complex((Numeric)rand(), (Numeric)rand());
    }
  }
}


void random_fill_matrix(Sparse& A, Numeric range, bool positive) {
  Index m = A.nrows();
  Index n = A.ncols();

  Rand<Numeric> random_number(positive ? 0 : -range, range);

  Index nelem = max(m, n);

  for (Index i = 0; i < nelem; i++) {
    Index m1, n1;
    m1 = rand() % m;
    n1 = rand() % n;

    A.rw(m1, n1) = random_number();
  }
}


void random_fill_matrix(Matrix& A, Sparse& B, Numeric range, bool positive) {
  Index m = A.nrows();
  Index n = A.ncols();

  assert(B.nrows() == m);
  assert(B.ncols() == n);

  Rand<Numeric> random_number(positive ? 0 : -range, range);

  Index nelem = max(m, n);

  for (Index i = 0; i < nelem; i++) {
    Index m1, n1;
    m1 = rand() % m;
    n1 = rand() % n;

    A(m1, n1) = random_number();
    B.rw(m1, n1) = A(m1, n1);
  }
}


void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive) {
  random_fill_matrix(A, range, positive);
  Matrix M(A);
  A += transpose(M);
}
void random_fill_matrix_symmetric(ComplexMatrixView A,
                                  Numeric range,
                                  bool positive) {
  random_fill_matrix(A, range, positive);
  ComplexMatrix M(A);
  A += transpose(M);
}


void random_fill_matrix_pos_def(MatrixView A, Numeric range, bool positive) {
  Index n = A.ncols();

  // Ensure that A is square.
  assert(A.ncols() == A.nrows());

  // Generate random, pos. semi-def. Matrix
  random_fill_matrix(A, range, positive);
  Matrix M(A);
  mult(A, M, transpose(M));

  // Add identity matrix.
  for (Index i = 0; i < n; i++) {
    A(i, i) += 1.0;
  }
}


void random_fill_matrix_pos_semi_def(MatrixView A,
                                     Numeric range,
                                     bool positive) {
  random_fill_matrix(A, range, positive);
  Matrix M(A);
  mult(A, M, transpose(M));
}


void random_fill_vector(VectorView v, Numeric range, bool positive) {
  Index n = v.nelem();

  Rand<Numeric> rand(positive ? 0 : -range, range);

  for (Index i = 0; i < n; i++) {
    v[i] = rand();
  }
}


MatrixView random_submatrix(MatrixView A, int m, int n) {
  Index m0(A.nrows()), n0(A.ncols());
  assert((0 <= m) && (m <= m0));
  assert((0 <= n) && (m <= n0));

  Rand<Index> rand_m(0, (m0 - m - 1)), rand_n(0, (n0 - n - 1));
  Index m1, n1;
  m1 = rand_m();
  n1 = rand_n();

  Range r1(m1, m), r2(n1, n);
  return A(r1, r2);
}


Range random_range(Index n) {
  Rand<Index> extent(1, n);
  Index e = extent();

  Index s = 0;
  if (0 <= (n - e - 1)) {
    Rand<Index> start(0, n - e - 1);
    s = start();
  }
  return Range(s, e);
}


Numeric get_maximum_error(ConstVectorView v1,
                          ConstVectorView v2,
                          bool relative) {
  Index n = min(v1.nelem(), v2.nelem());

  Numeric max = 0.0, err = 0.0;

  for (Index i = 0; i < n; i++) {
    err = 0.0;

    if (relative) {
      if (v2[i] != 0.0) {
        err = abs((v2[i] - v1[i]) / v2[i]);
      }

    } else {
      err = abs(v2[i] - v2[i]);
    }

    if (err > max) {
      max = err;
    }
  }
  return err;
}


Numeric get_maximum_error(ConstMatrixView A1,
                          ConstMatrixView A2,
                          bool relative) {
  Index m = min(A1.nrows(), A2.nrows());
  Index n = min(A1.ncols(), A2.ncols());

  Numeric max = 0.0, err = 0.0;

  for (Index i = 0; i < m; i++) {
    for (Index j = 0; j < n; j++) {
      err = 0.0;

      if (relative) {
        if (A2(i, j) != 0.0) {
          err = abs((A2(i, j) - A1(i, j)) / A2(i, j));
        }

      } else {
        err = A2(i, j) - A1(i, j);
      }

      if (err > max) {
        max = err;
      }
    }
  }

  return max;
}
Numeric get_maximum_error(ConstComplexMatrixView A1,
                          ConstComplexMatrixView A2,
                          bool relative) {
  Index m = min(A1.nrows(), A2.nrows());
  Index n = min(A1.ncols(), A2.ncols());

  Numeric max = 0.0, err = 0.0;

  for (Index i = 0; i < m; i++) {
    for (Index j = 0; j < n; j++) {
      err = 0.0;

      if (relative) {
        if (A2(i, j).real() != 0.0 && A2(i, j).imag() != 0.0) {
          err = abs((A2(i, j) - A1(i, j)) / A2(i, j));
        }

      } else {
        err = abs(A2(i, j) - A1(i, j));
      }

      if (err > max) {
        max = err;
      }
    }
  }

  return max;
}