clear all, close all, clc
%% Generate Data
Beta = [10; 28; 8/3]; % Lorenz's parameters (chaotic)
n = 3;
x0=[-8; 8; 27];  % Initial condition
tspan=[.01:.01:50];
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,n));
[t,x]=ode45(@(t,x) lorenz(t,x,Beta),tspan,x0,options);

%% Compute Derivative
for i=1:length(x)
    dx(i,:) = lorenz(0,x(i,:),Beta);
end

%% Build library and compute sparse regression
polyorder = 3;
Theta = poolData(x,n,polyorder);      % up to third order polynomials


%lambda = 0.025;      % lambda is our sparsification knob.
%Xi = sparsifyDynamics(Theta,dx,lambda,n)
%poolDataLIST({'x','y','z'},Xi,n,polyorder);



%% Apply LASSO
snap = size(Theta,2);       % total number of terms 
Xi = zeros(snap,n);
for k=1:n                   % for each component in state variable X

    [XiAll FitInfo] = lasso(Theta,dx(:,k),'CV',10);

    XiSparse = XiAll(:,FitInfo.Index1SE);
    ind = find(abs(XiSparse)>0);

    XiDeBiased = pinv(Theta(:,ind))*dx(:,k);

    Xi(ind,k) = XiDeBiased;

end
poolDataLIST({'x','y','z'},Xi,n,polyorder);


