source: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/RadialBasisFunctions/MatrixUtilities.cs @ 14386

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab 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 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab 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 HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Data;

namespace HeuristicLab.Algorithms.DataAnalysis {
  internal static class MatrixUtilities {

    public static T[,] To2D<T>(T[][] source) {
      try {
        var firstDim = source.Length;
        var secondDim = source.GroupBy(row => row.Length).Single().Key; // throws InvalidOperationException if source is not rectangular
        var result = new T[firstDim, secondDim];
        for (var i = 0; i < firstDim; ++i)
          for (var j = 0; j < secondDim; ++j)
            result[i, j] = source[i][j];
        return result;
      }
      catch (InvalidOperationException) {
        throw new InvalidOperationException("The given jagged array is not rectangular.");
      }
    }

    public static DoubleMatrix SortRows(this DoubleMatrix x, int c1) {
      var sortable = new List<double[]>();
      for (var i = 0; i < x.Rows; i++)
        sortable.Add(SelectRow(x, i).ToArray());
      return new DoubleMatrix(To2D(sortable.OrderBy(p => p[c1]).ToArray()));
    }
    public static DoubleMatrix Ones(int rows, int cols) {
      var val = new double[rows, cols];
      for (var i = 0; i < rows; i++)
        for (var j = 0; j < cols; j++)
          val[i, j] = 1;
      return new DoubleMatrix(val);
    }
    public static DoubleMatrix Identity(int rows, int cols) {
      var val = new double[rows, cols];
      for (var i = 0; i < Math.Min(rows, cols); i++)
        val[i, i] = 1;
      return new DoubleMatrix(val);
    }
    public static DoubleMatrix Elementwise(this DoubleMatrix m, Func<double, double> func) {
      var val = new double[m.Rows, m.Columns];
      for (var i = 0; i < m.Rows; i++)
        for (var j = 0; j < m.Columns; j++)
          val[i, j] = func.Invoke(m[i, j]);
      return new DoubleMatrix(val);
    }
    public static DoubleMatrix Elementwise(this DoubleMatrix m, DoubleMatrix m1, Func<double, double, double> func) {
      if (m.Rows != m1.Rows || m.Columns != m1.Columns) throw new ArgumentException("Matrices are of different sizes");
      var val = new double[m.Rows, m.Columns];
      for (var i = 0; i < m.Rows; i++)
        for (var j = 0; j < m.Columns; j++)
          val[i, j] = func.Invoke(m[i, j], m1[i, j]);
      return new DoubleMatrix(val);
    }
    public static DoubleMatrix Elementwise(this DoubleMatrix m, DoubleMatrix m1, DoubleMatrix m2, Func<double, double, double, double> func) {
      if (m.Rows != m1.Rows || m.Columns != m1.Columns || m1.Rows != m2.Rows || m1.Columns != m2.Columns) throw new ArgumentException("Matrices are of different sizes");
      var val = new double[m.Rows, m.Columns];
      for (var i = 0; i < m.Rows; i++)
        for (var j = 0; j < m.Columns; j++)
          val[i, j] = func.Invoke(m[i, j], m1[i, j], m2[i, j]);
      return new DoubleMatrix(val);
    }
    public static DoubleMatrix SelectRow(this DoubleMatrix m, int row) {
      var m1 = new DoubleMatrix(1, m.Columns);
      for (var i = 0; i < m.Columns; i++) m1[0, i] = m[row, i];
      return m1;
    }
    public static DoubleMatrix SelectColumn(this DoubleMatrix m, int col) {
      var m1 = new DoubleMatrix(m.Rows, 1);
      for (var i = 0; i < m.Rows; i++) m1[i, 0] = m[i, col];
      return m1;
    }
    public static DoubleMatrix RowBind(this DoubleMatrix m1, DoubleMatrix m2) {
      if (m1.Columns != m2.Columns) throw new ArgumentException("Matrices need to have the same number of columns for RowBind");
      var m = new DoubleMatrix(m1.Rows + m2.Rows, m1.Columns);
      for (var i = 0; i < m1.Rows; i++)
        for (var j = 0; j < m.Columns; j++)
          m[i, j] = m1[i, j];
      for (var i = 0; i < m2.Rows; i++)
        for (var j = 0; j < m.Columns; j++)
          m[i + m1.Rows, j] = m2[i, j];
      return m;

    }
    public static DoubleMatrix ColumnBind(this DoubleMatrix m1, DoubleMatrix m2) {
      if (m1.Rows != m2.Rows) throw new ArgumentException("Matrices need to have the same number of rows for ColumnBind");
      var m = new DoubleMatrix(m1.Rows, m1.Columns + m1.Columns);
      for (var i = 0; i < m1.Rows; i++) {
        for (var j = 0; j < m1.Columns; j++)
          m[i, j] = m1[i, j];
        for (var j = 0; j < m2.Columns; j++)
          m[i, j + m1.Columns] = m2[i, j];
      }
      return m;
    }
    public static DoubleMatrix Select(this DoubleMatrix m, int[] rows, int[] cols) {
      var res = new DoubleMatrix(rows.Length, cols.Length);
      for (var i = 0; i < res.Rows; i++) {
        for (var j = 0; j < res.Columns; j++)
          res[i, j] = m[rows[i], cols[j]];
      }
      return res;
    }
    public static DoubleMatrix Select(this DoubleMatrix m, int[] rows) {
      return Select(m, rows, GetIndexVector(m.Columns));
    }
    public static int[] GetIndexVector(int size) {
      var res = new int[size];
      for (var i = 0; i < size; i++)
        res[i] = i;
      return res;
    }
    public static DoubleMatrix Transpose(this DoubleMatrix m) {
      var res = new DoubleMatrix(m.Columns, m.Rows);
      for (var i = 0; i < res.Rows; i++)
        for (var j = 0; j < res.Columns; j++)
          res[i, j] = m[j, i];
      return res;
    }
    public static DoubleMatrix ToRow(IEnumerable<double> values) {
      var arr = values.ToArray();
      var res = new DoubleMatrix(1, arr.Length);
      for (var i = 0; i < arr.Length; i++) res[0, i] = arr[i];
      return res;
    }
    public static double[] Min(this DoubleMatrix m, int dim) {
      if (dim != 0 && dim != 1) throw new ArgumentException("Can only find the minimum along rows or colums");
      var res = new double[dim == 0 ? m.Rows : m.Columns];
      for (var i = 0; i < res.Length; i++)
        res[i] = (dim == 0 ? SelectRow(m, i) : SelectColumn(m, i)).Min();
      return res;
    }
    public static double[] Max(this DoubleMatrix m, int dim) {
      if (dim != 0 && dim != 1) throw new ArgumentException("Can only find the minimum along rows or colums");
      var res = new double[dim == 0 ? m.Rows : m.Columns];
      for (var i = 0; i < res.Length; i++)
        res[i] = (dim == 0 ? SelectRow(m, i) : SelectColumn(m, i)).Max();
      return res;
    }
    public static double[] GetElementVector(int size) {
      var res = new double[size];
      for (var d = 1; d <= size; d++) {
        res[d - 1] = d;
      }
      return res;
    }
    public static double[] GetRange(double min, double max, int length) {
      if (length < 2 || min > max) throw new ArgumentException("can not create range of length: " + length + " from " + min + " to " + max);
      var res = new double[length];
      var step = (max - min) / (length - 1);
      res[0] = min;
      for (var i = 1; i < length; i++)
        res[i] = step + res[i - 1];
      return res;
    }
    public static int[] GetRange(int min, int max) {
      if (min > max) throw new ArgumentException("can not create range from " + min + " to " + max);
      var res = new int[max - min + 1];
      var i = 0;
      for (var v = min; v <= max; v++)
        res[i++] = v;
      return res;
    }
    public static bool IsEqualTo(this DoubleMatrix m1, DoubleMatrix m2) {
      if (m1.Rows != m2.Rows) return false;
      if (m1.Columns != m2.Columns) return false;
      return m1.SequenceEqual(m2);
    }
    public static DoubleMatrix ToColumnVector(this IEnumerable<double> values) {
      var arr = values.ToArray();
      var res = new DoubleMatrix(arr.Length, 1);
      for (var i = 0; i < arr.Length; i++)
        res[i, 0] = arr[i];
      return res;
    }
    public static DoubleMatrix ToMatrix<T>(this IEnumerable<T> rows) where T : IEnumerable<double> {
      return new DoubleMatrix(To2D(rows.Select(x => x.ToArray()).ToArray())); //not the fastest way to do it

    }

    /// <summary>
    /// m need to be symmetric positive definit
    /// </summary>
    /// <param name="m"></param>
    /// <returns></returns>
    public static DoubleMatrix Invert(this DoubleMatrix m) {
      int info;
      alglib.matinvreport report;
      var upper = m.CloneAsMatrix();
      alglib.rmatrixinverse(ref upper, out info, out report);
      if (info != 1) throw new ArgumentException("Could not invert matrix. Is ist quadratic symmetric positive definite?");
      return new DoubleMatrix(upper);
    }
    public static DoubleMatrix Mul(this DoubleMatrix left, DoubleMatrix right) {
      if (left.Columns != right.Rows) throw new ArgumentException("Dims dont match");
      var res = new double[left.Rows, right.Columns];
      alglib.rmatrixgemm(left.Rows, right.Columns, left.Columns, 1, left.CloneAsMatrix(), 0, 0, 0, right.CloneAsMatrix(), 0, 0, 0, 0, ref res, 0, 0);
      return new DoubleMatrix(res);
    }
    public static DoubleMatrix Solve(this DoubleMatrix a, DoubleMatrix b) {
      if (b.Columns != 1 || a.Rows != b.Rows) throw new ArgumentException("This system can not be solved due to dimensions");
      var b1 = b.ToArray();
      double[] x;
      int info;
      alglib.densesolverlsreport denseSolveRep;//what is this used for? In HeuristicLab.Algorithms.DataAnalysis.GaussianProcessSurrogate is a use-case similar to this one. What do you normally do with the rep?
      alglib.rmatrixsolvels(a.CloneAsMatrix(), a.Rows, a.Columns, b1, 0.0, out info, out denseSolveRep, out x);
      if (info != 1) {
        throw new ArgumentException("This system can not be solved due to content or numerics");
      }

      return x.ToColumnVector();
    }

  }
}