function [x_sol,L_curve_xpts,L_curve_ypts,percent_errors,residual_norms] = continuation_scheme_with_wavelets(A,b,x,tau_max,tau_min,num_taus,wav_tol)

	% x is the true solution, which would not be available to you
	% here just used for testing (you can test with x too)

	m = size(A,1);
	n = size(A,2);

	% initial guess
	w0 = zeros(n,1);

	w_bar = w0;

	% estimate the spectral norm of the matrix A*Winv 
	% where Winv is the inverse haar wavelet transform
	v = randn(n,1);
	v = v/norm(v);

	% 50 iterations of power scheme should be enough
	% to estimate the spectral norm of M = A*Winv, which 
	% we do not form explicitly
	for k=1:50
		%w = M'*(M*v) = W*A'*A*Winv*v;
		vinv = haar_1d_inverse(v);
		w = haar_1d(A'*(A*vinv)); 
		normval = norm(w)/norm(v);
		v = w/norm(w);
	end
	alpha = sqrt(normval)

	% set up arrays holding stuff
	solution_norms = zeros(num_taus,1); 
	residual_norms = zeros(num_taus,1); 
	L_curve_xpts = zeros(num_taus,1); 
	L_curve_ypts = zeros(num_taus,1); 
	residual_norms = zeros(num_taus,1); 
	percent_errors = zeros(num_taus,1); 
	taus = zeros(num_taus,1);

	% tau step in log space
	S = (log(tau_max) - log(tau_min))/(num_taus - 1);
	w_bar = w0;
	TOL = 0;
	maxiters = 100; %% max iters at each tau

	% this just needs to be computed once, Mt*b = Winvt*At*b = W*At*b
	Mtb = (1/alpha^2)*haar_1d(A'*b);

        load('data/wavelet_mats2.mat');
        A = A*Winv;
        nest = normest(A);
        alpha = nest;
        A = A/alpha;
        b = b/alpha;
        AtA = A'*A;

	fprintf('run continuation scheme\n');
	for i=1:(num_taus)
		tau = exp(log(tau_max) - S*(i-1));
		taus(i) = tau;

		%[w_bar,num_iters_taken,tElapsed] = thr_fista_with_wavelets_solve(A,b,Mtb,alpha,tau,w0,maxiters); 
                TOL = 0;
                l = 1.5; 
                p = 1;
                lambda = tau;
                maxiters = 5;
                %[w_bar,num_iters_taken,tElapsed,residuals,Js] = irls_sys_solve(A,AtA,b,x,l,p,lambda,w0,TOL,maxiters); 
                [w_bar,num_iters_taken,tElapsed,funct_vals,perrors_local]=conj_grad(A,b,x,p,lambda,w0,maxiters);
                %[w_bar,num_iters_taken,tElapsed,funct_vals,perrors_local] = thr_fista_solve(A,b,x,tau,w0,TOL,maxiters); 
		w0 = w_bar;

		% transform back
		%x_bar = haar_1d_inverse(w_bar);
		x_bar = Winv*w_bar;
		residual_norms(i) = norm(A*x_bar - b); 
		solution_norms(i) = norm(w_bar,1);
		percent_errors(i) = 100*norm(x_bar - x)/norm(x);
		L_curve_xpts(i) = log(residual_norms(i));
		L_curve_ypts(i) = log(solution_norms(i));
	end

	% at the end, remove small coefficients, 
	% then convert back to regular coordinates from wavelets
	% we follow Donoho wavelet denoising advice
	w_bar = softThreshold(w_bar,(sqrt(2*log(n))/sqrt(n))*wav_tol);
	%x_sol = haar_1d_inverse(w_bar);
	x_sol = Winv*w_bar;
end

