(*Step 1: Generate the function to be retrieved and the "measurement matrix"*)
dim = 200; (*Number of points in the original data (input)*)
nm = 30; (*Number of measurements done*)
hf := (E^(-(x^2/2)) HermiteH)/Sqrt\Pi; (*Hermite polynomials*)
nb = 6; (*number of nonzero coefficients*)
nzc = Sort@RandomSampleRange40, nb (*non zero coefficients*)
c = RandomReal{-0.8, 0.8}, nb; (*coeffieicents*)
g = Sumcj hfnzcj, z, {j, 1, nb}; (*function we would like to retrieve (sparse in the Hermite polynomial basis by construction)*)
xideal = Normal@SparseArraynzc -> c, {dim}; (*ideal x vector we would like to retrieve*)
df = Tableg // N, {z, -7, 7, 14/(dim - 1)}; 
m = Sort@RandomSampleRangedim, nm; (*Points where we make the measurements (in the canonical basis)*)
Atmp = SparseArrayTable{j, mj}, {j, 1, nm} -> Table1, {j, 1, nm} , {nm, dim}; (*measurement matrix*)
\CapitalPhi = TableN@hfn, x, {n, 1, dim}, {x, -7, 7, 14/(dim - 1)};(*matrix of change of basis from Hermite to canonical (rows are the hermite polynomials)*)
A = Atmp.Transpose\CapitalPhi;
y = A.xideal;

(*Step 2: A working (but not optimized) implementation of the Orthogonal Matching Pursuit algorithm*)
(*initialization*)
\Epsilon0 = 10^-10; (*error threshold*)
x = Table0, {dim}; (*first guess for the coefficients is that they are all zero*)
r = y; (*so the residuals are identical to the measurements*)
nS = Rangedim; (*the complement of the support is everything*)

\Epsilon = TableNormAAll, j.r/NormAAll, j ^2 AAll, j - r^2, {j, 1, dim}nS; (*calculate the error for all columns of A that are not already in the support*)
j = Position\Epsilon, Min\Epsilon1, 1; (*select the one with the biggest impact on the error*)
\Epsilon = \Epsilonj; (*that one is the new error*)
nS = DropnS, {j}; (*update the complement of the support*)
S = DeleteCasesRangedim, Alternatives @@ nS; (*and thus update the support*)
AS = AAll, S; 
xS = InverseTransposeAS.AS.TransposeAS.y; (*find the best fit for the new estimate of x (least square fit)*)
r = y - A.x; (*update the residuals*)
tmp = x;

evo = ReapWhile\Epsilon > \Epsilon0, (*repeat until the error is small enough*)
      \Epsilon = TableNormAAll, j.r/NormAAll, j ^2 AAll, j - r^2, {j, 1, dim}nS;
      j = Position\Epsilon, Min\Epsilon1, 1;
      \Epsilon = \Epsilonj;
      nS = DropnS, {j};
      S = DeleteCasesRangedim, Alternatives @@ nS;
      AS = AAll, S; 
      xS = InverseTransposeAS.AS.TransposeAS.y;
      r = y - A.x;
      Sowx;
      ;2, 1;
evo = Prependevo, tmp;

(*Step 3: Generate the animation*)
p0 = TableGraphicsRow{
     Show
      ListPlotdf, Joined -> True, PlotStyle -> {Thick, Gray}, Axes -> False, PlotRange -> {-1, 1}, Epilog -> {PointSize0.02, Point({m, y} // Transpose)1 ;; k } 
     ,
     ListPlotTable0, {50}, PlotRange -> {{0, 50}, {-1, 1}}, Filling -> Axis, PlotStyle -> {Purple, PointSize0.02}, Axes -> False, Frame -> True, FrameLabel -> {"Element of the basis", "Coefficient"}, 
      LabelStyle -> {Bold, Black}
     }, ImageSize -> Large
   , {k, 1, nm};
p1 = Table
   GraphicsRow{
     Show
      ListPlotdf, Joined -> True, PlotStyle -> {Thick, Gray}, PlotRange -> {-1, 1}, ListPlot(1 - \Tau)*0 + \Tau Transpose\CapitalPhi.evo1, Joined -> True, PlotStyle -> {Thick, Orange}
      , Axes -> False, PlotRange -> {-1, 1}, Epilog -> {PointSize0.02, Point{m, y} // Transpose}
      
     ,
     ListPlot\Tau evo1, PlotRange -> {{0, 50}, {-1, 1}}, Filling -> Axis, PlotStyle -> {Purple, PointSize0.02}, Axes -> False, Frame -> True, FrameLabel -> {"Element of the basis", "Coefficient"}, 
      LabelStyle -> {Bold, Black}
     }, ImageSize -> Large
   , {\Tau, 0, 1, 0.1};
p2 = TableTable
    GraphicsRow{
      Show
       ListPlotdf, Joined -> True, PlotStyle -> {Thick, Gray}, PlotRange -> {-1, 1}, ListPlot(1 - \Tau)*Transpose\CapitalPhi.evok - 1 + \Tau Transpose\CapitalPhi.evok, 
        Joined -> True, PlotStyle -> {Thick, Orange}
       , Axes -> False, PlotRange -> {-1, 1}, Epilog -> {PointSize0.02, Point{m, y} // Transpose}
       
      ,
      ListPlot(1 - \Tau) evok - 1 + \Tau evok, PlotRange -> {{0, 50}, {-1, 1}}, Filling -> Axis, PlotStyle -> {Purple, PointSize0.02}, Axes -> False, Frame -> True, 
       FrameLabel -> {"Element of the basis", "Coefficient"}, LabelStyle -> {Bold, Black}
      }, ImageSize -> Large
    , {\Tau, 0, 1, 0.1}, {k, 2, DimensionsS1 - 1};
ListAnimateFlattenJoinp0, p1, p2