Merge branch 'master' of https://git.iiens.net/trachino2017/sgi-choubi-nuja

matlab.h

+#pragma once
+
+#include <random>
+#include <iostream>
+#include <fstream>
+#include "../types.h"
+#include <complex>
+/*
+ * Code translated from Matlab base functions
+ * matlab translation guide : https://eigen.tuxfamily.org/dox/AsciiQuickReference.txt
+ */
+
+//-------------------------------------------------------------------------------------------------
+// VARIOUS
+
+const number pi = M_PI;
+
+// squares a number
+template<typename T> inline T squared(const T& x)
+{
+    return x*x;
+}
+
+//-------------------------------------------------------------------------------------------------
+// ALGEBRA
+
+// build a vector with values in the given interval
+Vector make_vector(const int fromValue, const int toValue, const int by = 1)
+{
+    int size = 1 + (toValue - fromValue) / by;
+    Vector result(size);
+
+    int currentValue = fromValue;
+    for(int i = 0; i < size; i++)
+    {
+        result(i) = currentValue;
+        currentValue += by;
+    }
+
+    return result;
+}
+
+// repeats a column vector n times in order to build a matrix
+Matrix repmat(const Vector& v, const int colNumber)
+{
+    Matrix result(v.size(), colNumber);
+
+    for(int col = 0; col < colNumber; col++)
+    {
+        result.col(col) = v;
+    }
+
+    return result;
+}
+
+// solves m = m1 \ m2
+inline ComplexMatrix mldivide(const ComplexMatrix& m1, const ComplexMatrix& m2)
+{
+    //return m1.colPivHouseholderQr().solve(m2);
+    return m1.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(m2); // you can use JacobiSVD for improved precision
+}
+
+// returns an hermitian matrix where v is the first row and conj(v) the first column
+Matrix toeplitz(const Vector& v)
+{
+    auto size = v.size();
+    Matrix result(size, size);
+
+    for(int row=0; row<size; row++)
+    {
+        for(int col=0; col<size; col++)
+        {
+            // complex
+            /*number value = v(std::abs(row - col));
+            if(row > col)
+            {
+                value = Sstd::conj(value);
+            }
+            result(row,col) = value;*/
+
+            // real
+            result(row,col) = v(std::abs(row - col));
+        }
+    }
+
+    return result;
+}
+
+// flip the column and rows
+inline Matrix flip(Matrix mat)
+{
+    Matrix matflipr = mat.rowwise().reverse(); // fliplr
+    Matrix matflipc = matflipr.colwise().reverse(); // flipud
+    return matflipc;
+}
+
+//-------------------------------------------------------------------------------------------------
+// POLYNOMIAL
+
+/*
+ * evaluates a polynomial on the given data
+ * use Horner's scheme
+ * NOTE: one could use a more numerically stable algorithm
+ */
+Vector polyval(const Polynomial& pol, const Vector& x)
+{
+    Vector result = Vector::Zero(x.size());
+
+    for(const auto& coef : pol)
+    {
+        result = result.binaryExpr(x, [coef](const number& result, const number& x){return result*x + coef;});
+    }
+
+    return result;
+}
+
+//-------------------------------------------------------------------------------------------------
+
+// RANDOM HELPING FUNCTION 
+
+const unsigned int n = 65535;//0xFFFF;
+const unsigned int m = 16807;
+const unsigned int b = 0;
+static long  seed = 43690;//0xAAAA;
+// generate an uniformaly distributed random number sampled from [0;0x7FFFFFFF]
+long uniform() {
+   long unsigned int hi = m * (seed & n);
+   long unsigned int lo = m * (seed >> 16);
+
+   lo+= (hi & 32767) << 16;
+   lo+= hi >> 15;
+
+   if (lo > 2147483647) //0x7FFFFFFF)
+      lo -= 2147483647; //0x7FFFFFFF;
+   seed=(long) lo;
+   return seed ;
+}
+
+
+//--------------------------------------------------------------------------------------------------
+
+// RANDOM
+
+// generate an uniformaly distributed random number sampled from [-1;1]
+inline number randu()
+{
+  //on le veut entre -1 et 1 / pour l'instant on l'a entre 0 et 1 
+    number result = uniform()/2147483647;  //(double)0x7FFFFFFF;
+    // number result = (uniform() - (number)3xFFFFFFF.8)/(number)3xFFFFFFF.8);
+     return result;
+  //  throw std::logic_error("Function not yet implemented.");
+}
+
+// generates a vector of the given size, the entries of the vector are all complex number
+// their real and imaginary parts are uniformaly distributed random number sampled from [-1;1]
+ComplexVector randu(int size)
+{
+    ComplexVector complexVector(size);
+    std::complex<number> mycomplex;   
+    for(int i = 0; i < size; i++)
+    {
+        std::complex<number> mycomplex(randu(), randu()); 
+        complexVector(i) = mycomplex;
+    }
+    return complexVector;
+
+    //throw std::logic_error("Function not yet implemented.");
+}
+
+// generate a normally distributed random number
+// sampled from a distribution with mean 0 and std 1
+inline number randn()
+{
+
+  number res;
+  number U0 = uniform()/2147483647;//(double)0x7FFFFFFF;
+  number U1 = uniform()/2147483647;//(double)0x7FFFFFFF;
+  res = sqrt(-2 * log(U0)) * cos(2 * M_PI * U1);
+  
+  return res;
+
+   // throw std::logic_error("Function not yet implemented.");
+}
+
+// generates a vector of the given size, the entries of the vector are all complex number
+// their real and imaginary parts are sampled from a normal distribution with mean 0 and std 1
+ComplexVector randn(int size)
+{
+    // NOTE: this can be built on the previous function
+    ComplexVector complexVector(size);
+    std::complex<number> mycomplex;   
+    for(int i = 0; i < size; i++)
+    {
+        std::complex<number> mycomplex(randn(), randn()); 
+        complexVector(i) = mycomplex;
+    }
+    return complexVector;
+    //throw std::logic_error("Function not yet implemented.");
+}
+
+
+
+//-------------------------------------------------------------------------------------------------
+// IMPORT
+
+/// import a vector of reals from a file
+Vector loadRealVector(const std::string& file)
+{
+    std::vector<number> result;
+    std::ifstream input(file);
+
+    number x;
+    while (input >> x)
+    {
+        result.push_back(x);
+    }
+
+    Vector output = Eigen::Map<Vector, Eigen::Unaligned>(result.data(), result.size());
+    return output;
+}
+
+/// import a vector of complex from a file
+ComplexVector loadComplexVector(const std::string& file)
+{
+    std::vector<complexNumber> result;
+    std::ifstream input(file);
+
+    number r;
+    number i;
+    while (input >> r)
+    {
+        input >> i;
+        complexNumber x(r,i);
+        result.push_back(x);
+    }
+
+    ComplexVector output = Eigen::Map<ComplexVector, Eigen::Unaligned>(result.data(), result.size());
+    return output;
+}
+
+//-------------------------------------------------------------------------------------------------