example08.cpp

#include <numcpp.h>
#include <numcpp/solvers.h>
using namespace numcpp;

int main()
{
  size_t n = 256;
  size_t N = n*n;
  size_t M = N / 4;

  auto a = array<cdouble>({1,2,3,4});
  auto idx = array<size_t>({0,1,2,3});

  SparseFFTArray<cdouble> F2(idx,4);

  print(dot(F2,a));
  print(solve(F2,a));

  //Array<cdouble> p = reshape(phantom(n),N);
  Array<cdouble> p = zeros(N);
  p(n*n/2 + n/4) = 1;

  auto tmp2 = eval(abs(reshape(p,n,n)));
  export_image( tmp2, "CSIm.png");

  //auto indices = eval(range(0,N-1));
  auto indices = randi<size_t>(N-1, M);

  SparseFFTArray<cdouble> F(indices,n,n);

  auto y = dot(F,p);
  auto noise = randn<double>(M);
  noise *= max(abs(y))*0.5;
  y += noise;

  {
    auto z = solve(F, y);
    export_image( eval( abs(reshape(z,n,n) ) ), "CSNaiv.png");
    std::cout << "error naiv: " << nrmsd(p, z) << " max " <<  max(abs(z)) << std::endl;
    std::cout << z.getMem() <<std::endl;
  }

  {
    auto z = SL0(F, y, 30, 3, 2);
    export_image( eval( abs(reshape(z,n,n) )), "CSSL0.png");
    std::cout << "error SL0: " << nrmsd(p, z) << " max " <<  max(abs(z)) << std::endl;
    std::cout << z.getMem() <<std::endl;
  }

  {
    auto z = fista(F, y, 30, 3, 1e-8);
    export_image( eval ( abs(reshape(z,n,n) )), "CSFista.png");
    std::cout << "error Fista: " << nrmsd(p, z) << " max " <<  max(abs(z)) << std::endl;
    std::cout << z.getMem() <<std::endl;
  }

}

/*!
@page example07 Compressed Sensing Example

@code
int main()
{
  size_t N = 256;

  auto p = phantom(N);
  export_image(p, "phantom.png");

  auto y = dwt(p, wavelets::Haar);
  export_image(y, "waveletHaar.png");

  auto y = dwt(p, wavelets::Daubechies4);
  export_image(y, "waveletDaubechies4.png");

  auto z = idwt(y, wavelets::Daubechies4);
  export_image(z, "reco.png");
}
@endcode
*/