function [budget, cost, Lperror, maxerror] = NonlinearApprox(Wavelet_Coefficients,Dilations,Shifts,Budget,Kernel,Kernel_Parameter,TL_Smoothness,Approx_Norm,Time_End)

%%Input Descriptions:
%Budget - Total budget of centers; 
%Kernel - String: options 'Gauss' or 'MQ'
%Kernel_Parameter - Gaussian exponent alpha, or MQ exponent alpha
%TL_Smoothness - s value for Triebel-Lizorkin space (roughly # of
%derivatives)
%Approx_Norm - Lp norm for estimating difference of f and approximant
%Wavelet_Coefficients - f_I coefficients of f in Meyer wavelet expansion
%Dilations/Shifts - vectors containing dilations and shifts of unit cube
%% here f = sum_j Wavelet_Coefficients_j*Meyer(2^(Dilations_j)*(x-Shifts_j))%%

%%Initialize with Default Options:
if nargin < 3
    disp('Error: Must give a minimum of 3 arguments');
end
if nargin < 4 || isempty(Budget)
    Budget = 20;
    disp('Budget assigned to be 20 centres');
end
if nargin < 5 || isempty(Kernel)
    Kernel = 'Gauss';
    disp('Kernel chosen to be the Gaussian with Parameter .5');
end
%{
if nargin < 6 | (isempty(Kernel_Parameter) & isequal(strcmp(Kernel,'Gauss'),1))
    Kernel_Parameter = .5;
    disp('Kernel_Parameter chosen to be .5 for the Gaussian kernel');
else
    Kernel_Parameter=-3/2;
    disp('Kernel_Parameter chosen to be -3/2 for the Multiquadric kernel');
end
%}
if nargin <  7 || isempty(TL_Smoothness)
    TL_Smoothness = 0;
    disp('Smoothness order chosen to be 0 by default');
end
if nargin < 8 || isempty(Approx_Norm)
    Approx_Norm = 1;
    disp('Approximation norm chosen to be L_1 by default');
end
if nargin < 9 || isempty(Time_End)
    Time_End = 6;
    disp('Interval in Time domain chosen to be [-6,6]');
end


%%Preliminaries: Define Meyer mother wavelet in time (Meyer) and frequency
%(Meyer_hat)%%
f1 = @(x) 2*(sin(2*pi/3*(x+5/4))-sin(4*pi/3*(x+5/4)))./(pi*(4*x+5))+2*(sin(2*pi/3*(x-1/4))-sin(4*pi/3*(x-1/4)))./(pi*(4*x-1));
f2 = @(x) (cos(8*pi/3*(x+1/8))-cos(4*pi/3*(x+1/8)))./(2*pi*(x+1/8))+(cos(4*pi/3*(x+7/8))-cos(8*pi/3*(x+7/8)))./(2*pi*(x+7/8));
Meyer = @(x) f1(x)+f2(x);
Meyer_hat = @(t) ((2*pi/3<abs(t) & abs(t)<=4*pi/3).*sin(pi/2*(3*abs(t)/(2*pi)-1)).*exp(1i*t/2))+((4*pi/3<abs(t) & abs(t)<=8*pi/3).*cos(pi/2*(3*abs(t)/(4*pi)-1)).*exp(1i*t/2));

%%Set Kernel to use in T_N%%
if strcmp(Kernel,'Gauss')==1
    Kernel = @(x) exp(-x.^2/Kernel_Parameter^2);
    Kernel_hat = @(x) Kernel_Parameter*sqrt(pi)*exp(-Kernel_Parameter^2*x.^2/4);
    disp('Using Gaussian Kernel');
elseif strcmp(Kernel,'MQ')==1
    Kernel = @(x) (x.^2+1).^Kernel_Parameter;
    Kernel_hat = @(x) sqrt(2*pi)*2^(1+Kernel_Parameter)/gamma(-Kernel_Parameter)*abs(x).^(-Kernel_Parameter-1/2).*besselk(abs(Kernel_Parameter+1/2),abs(x));
    disp('Using Multiquadric Kernel');
end

%%Set Triebel-Lizorkin space parameters%%
%p is the approximation norm; s is the smoothness order, and tau is given
%by the formula 1/tau = 1/p + s
tau = 1/(1/Approx_Norm+TL_Smoothness);

%%Evaluate f at a lot of points on interval [-time_end,time_end], which is arbitrary%%
time = linspace(-Time_End,Time_End,1000);
for i = 1:length(time)
    for k = 1:length(Shifts)
        terms(k) = Wavelet_Coefficients(k)*Meyer(2^Dilations(k)*(time(i)-Shifts(k)));
    end
     f(i) = sum(terms);
end

%%Determine Cost Distribution%%
%Initialize
cost = zeros(1,length(Shifts));
measure_dyadic_cube = 2.^Dilations;

%Calculate Cost vector (uses assumption that cubes are disjoint)%
cost = round(abs(Wavelet_Coefficients).^tau.*abs(measure_dyadic_cube).^(tau/Approx_Norm)./(sum(abs(Wavelet_Coefficients).^tau.*abs(measure_dyadic_cube).^(tau/Approx_Norm)).^(1/tau))*Budget);

%%Create approximant T_N%%
for i=1:length(time)
    for l = 1:length(cost)
        j = -cost(l):cost(l);
        for k=1:length(j)
            psi_kernel(l,k) = (cost(l)~=0)*real(1/(2*pi)*integral(@(x) Meyer_hat(x)./Kernel_hat(x).*exp(1i*x*j(k)/sqrt(cost(l))),-10,10));
            Kernel_vec(l,k) = (cost(l)~=0)*Kernel(2^Dilations(l)*(time(i)-Shifts(l))-j(k)/sqrt(cost(l)));
        end
        T_N_flat_terms(l) = (cost(l)~=0)*cost(l)^(-1/2)*sum(psi_kernel(l,:).*Kernel_vec(l,:));
    end
    S_f_N(i) = sum(Wavelet_Coefficients.*T_N_flat_terms);
end
plot(time,f,time,S_f_N,'-r');

%%To display Budget and Error, remove the comments below%%
%{
disp('Budget:');
disp(Budget);
disp('Lp error:'); 
disp((sum(abs(f-S_f_N).^Approx_Norm)*2*Time_End/length(time))^(1/Approx_Norm));
%}
budget = Budget;
Lperror = (sum(abs(f-S_f_N).^Approx_Norm)*2*Time_End/length(time))^(1/Approx_Norm);
maxerror = max(abs(f-S_f_N));

end