Sélectionner une révision Git
matlab.h 6,30 Kio
#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;
}
//-------------------------------------------------------------------------------------------------