function shake(ipr0)
	% shake(ipr)
	% 
	% Apply shake step in shake the box
	% Based on an initial guess of the particle distribution, iteratively
	% refine particle position and intenisty to match the image
	%
	% The Levenberg-Marquardt algorithm is implemented to iteratively solve
	% for the optimum particle distribution. 
	%
    % Input:
    %   camera      1 x n_c     camera model structure
    %   otf         1 x n_c     otf model structure
    %   ipr.img_meas     1 x n_c     cell array of particle images
    %
    %   ipr, reconstruction structure containing following fields
    %   re_part     3 x n_p     particle object space positions
    %   dx_part     3 x n_p     particle displacement relative to frame 2
    %   A_part      1 x n_p     particle brightness
    %
    %   opts        IPR options structure containing the following fields
    %   dx_max      maximum displacement (3 x 1)
    %   

	%% load calibration
    n_cameras		= ipr0.setup.n_cameras;
	camera			= ipr0.setup.camera;
    otf             = ipr0.setup.otf;
    re_lim          = ipr0.setup.re_bb_lim;
    re_fov          = abs(re_lim(2,:) - re_lim(1,:));
    
    opts            = ipr0.opts;
    n_iter			= opts.n_shake;
	
	% support
	support			= otf(1).support;
	S				= size(support);
	n_support_px	= S(1) * S(2);
    
	% dimensionless time of each step (normalsied by delta t)
	t_hat			= -1:1;

    %% define sorting plane
    % normal distance from this plane defines ordering of particles
    % we define this as the longest dimension of the reconstruction volume
    [~,i_long]      = max(re_fov);
    e_sort          = zeros(3,1);
    e_sort(i_long)  = 1;
	
	%% image space support
	% to increase sparsity, we only explicitly optimise the fit of image
	% values which are non-zero. This avoids optimising on regions with no
	% particles
	img_mask		= cell(1,n_cameras);
	for ic = 1 : n_cameras
		img_mask{ic}= ipr0.img_meas{ic} > 0;
	end
	
	%% refine position over n_iter iterations
	n_part			= ipr0.count();
	ipr0.lambda	    = opts.lambda0 * ones(1,n_part, 'double');
    % set bounds on displacement of particles during shaking procedure
	ipr0.re_part_lb	= bsxfun(@minus, ipr0.re_part, opts.dx_max_step(:));
	ipr0.re_part_ub	= bsxfun(@plus , ipr0.re_part, opts.dx_max_step(:));
		
	for n = 1 : n_iter
		%% report
		if opts.b_verbose
			fprintf('I%02d S%02d projecting and measuring residual\n', ipr0.riter, n);
		end

        %% count particles
        n_part          = ipr0.count(); 

		%% sort particles to optimise sparsity structure
		% particles which are likely to interact with each other will tend
		% to lie on the same epipolar plane. by taking the dot-product
		% of position with the sorting vector, we measure the height of each 
		% particle from the reference plane
		% then, by re-sorting according to this height, we place particles
		% that are likely to interact with each other at similar indices 
		% in the Jacobian matrix. this results in a banded structure for
		% J'J, which is better suited for solution with a sparse solver.
		
        ipr0.sort(e_sort);
        
        %% project particles onto images
        pimg0                   = ipr0.project();
        img_res_before          = ipr0.pimg_res;
        
		%% iterate to obtain jacobian, gramian and residual
        % jacobian J(ii,jj) represents change in projected image with respect to
        % the model variables re_part, dx_part
        % 
        % it is of dimension (n_px x 6*n_part)
        %
        % each particle has 6 columns in the Jacobian matrix
        % model parameters are plaid like
        % [x1 y1 z1 dx1 dy1 dz1, x2 y2 z2 dx2 dy2 dz2, ...]
        % 
				
        % Gramian moment matrix
		JtJ						= sparse([], [], [], 6*n_part, 6*n_part);
        % residual
		JtDy					= zeros(6*n_part,1);
        
		% add gpu capability
		JtJ						= sparse(JtJ);
		JtDy					= sparse(JtDy);
		n_px					= zeros(1, n_cameras);

		for ic = 1 : n_cameras
			%% size of this image
			n_i					= size(ipr0.img_meas{ic},1);
			n_j					= size(ipr0.img_meas{ic},2);
            n_px(ic)            = n_i*n_j;
            
            %% accumulate jacobian over time
			% jacobian is of size n_px x 7*n_part
			%
			% expected number of elements is:
			%	#visible particles x #support px x 6 
			n_visible_sum		= size(pimg0(1,ic).ind,2) ...
								+ size(pimg0(2,ic).ind,2) ...
								+ size(pimg0(3,ic).ind,2);
			big_jac             = sparse([], [], [], n_px(ic), 6*n_part, n_support_px*n_visible_sum*6);

			for rgb = 1 : 3
                %% map jacobian back to particle images
                % indices of particles 
                ind                 = pimg0(rgb,ic).ind;
                % row and column coordinates of particle images
                row_map				= pimg0(rgb,ic).map(:,:,1);
                col_map				= pimg0(rgb,ic).map(:,:,2);
                % take care of out of bounds elements in map
                % dep_map is always in bounds by construction
                map_ok				= row_map >= 1 & row_map <= n_i ...
                                    & col_map >= 1 & col_map <= n_j;
                % create map of indices where map is valid
				n_visible			= size(ind, 2);
                ind_map				= zeros([n_support_px n_visible]);
                ind_map(map_ok)     = sub2ind([n_i n_j], row_map(map_ok), col_map(map_ok));

                %% make big-ass Jacobian matrix
                jac_row_idx			= ind_map;
                jac_col_idx			= repmat(ind, n_support_px, 1);
                n_ok				= nnz(map_ok);
                pimg_jac_x			= double(pimg0(rgb,ic).jac(:,:,1));
                pimg_jac_y			= double(pimg0(rgb,ic).jac(:,:,2));
                pimg_jac_z			= double(pimg0(rgb,ic).jac(:,:,3));

				% jacobian with respect to model parameters, encoded
				% x, y, z, dx, dy, dz
                % jacobian with respect to particle brightness
				rows				= jac_row_idx(map_ok);
				cols				= 6*(jac_col_idx(map_ok)-1);
                % jacobian with respect to particle position at time t0
                big_jac             = big_jac + sparse(rows, cols+1, pimg_jac_x(map_ok), n_px(ic), n_part*6);
                big_jac             = big_jac + sparse(rows, cols+2, pimg_jac_y(map_ok), n_px(ic), n_part*6);
                big_jac             = big_jac + sparse(rows, cols+3, pimg_jac_z(map_ok), n_px(ic), n_part*6);
                % jacobian with respect to velocity
                if rgb == 1 || rgb == 3
                    big_jac         = big_jac + sparse(rows, cols+4, t_hat(rgb)*pimg_jac_x(map_ok), n_px(ic), n_part*6);
                    big_jac         = big_jac + sparse(rows, cols+5, t_hat(rgb)*pimg_jac_y(map_ok), n_px(ic), n_part*6);
                    big_jac         = big_jac + sparse(rows, cols+6, t_hat(rgb)*pimg_jac_z(map_ok), n_px(ic), n_part*6);
				end
            end
            
            %% compute moment matrix, J'J
			big_JtJ				= big_jac.' * big_jac;
            big_JtJ				= gather(big_JtJ);
			% accumulate system jacobian 
            JtJ					= JtJ  + big_JtJ;

            %% compute J' * (y - y^)
            % this is the term on the RHS of the Levenberg-Marquard update
            big_JtDy			= big_jac.' * double(ipr0.img_res{ic}(:));
			big_JtDy			= gather(big_JtDy);
            % accumulate residual
            JtDy				= JtDy + big_JtDy(:);
			clear big_JtJ big_JtDy;
			clear big_jac;
		end
		JtJ_bw					= bandwidth(JtJ);
		Jtj_bandden				= JtJ_bw / size(JtJ,1);
		fprintf('I%02d S%02d JtJ sparsity = %0.4f%% (%d elements)\n', ipr0.riter, n, 100*nnz(JtJ)/numel(JtJ), nnz(JtJ));
		fprintf('I%02d S%02d JtJ bandwidth= %5d (den=%0.3f%%)\n', ipr0.riter, n, JtJ_bw, 100*Jtj_bandden);
		
		%% update model
		% given residual and Jacobian, we calculate an update to the position
		% velocity and intensity (diameter^2) of each particle
		if opts.b_verbose
			fprintf('I%02d S%02d shaking the box!\n', ipr0.riter, n);
		end

		% add Levenberg-Marquardt parameter
		lambda_rep				= repelem(ipr0.lambda(:), 6);
		JtJ_diag				= diag(diag(JtJ).*lambda_rep);
		JtJ_plus_JtJ_lambda     = JtJ + JtJ_diag;
		% solve for change in model parameter
		delta_beta				= JtJ_plus_JtJ_lambda \ JtDy;
		% extract model parameter
        delta_pos				= [delta_beta(1:6:end,1) delta_beta(2:6:end,1) delta_beta(3:6:end,1)]';
		delta_dx				= [delta_beta(4:6:end,1) delta_beta(5:6:end,1) delta_beta(6:6:end,1)]';
		
		% limit maximum change in position to dx_max_step
		delta_pos_len			= sqrt(sum(delta_pos.^2,1));
		delta_pos_hat			= bsxfun(@rdivide, delta_pos, delta_pos_len);
		delta_pos				= bsxfun(@times, delta_pos_hat, min(delta_pos_len, opts.dx_max_step));
		
		% report update
		if opts.b_verbose
			fprintf('I%02d S%02d rms update pos %0.3f %0.3f %0.3f mm\n', ipr0.riter, n, rms(delta_pos,2));
			fprintf('I%02d S%02d max update pos %0.3f %0.3f %0.3f mm\n', ipr0.riter, n, max(abs(delta_pos),[],2));
			fprintf('I%02d S%02d rms update vel %0.3f %0.3f %0.3f mm/frame\n', ipr0.riter, n, rms(delta_dx,2));
			fprintf('I%02d S%02d max update vel %0.3f %0.3f %0.3f mm/frame\n', ipr0.riter, n, max(abs(delta_dx),[],2));
		end
		
        % create update
        ipr1                    = copy(ipr0);
        ipr1.re_part            = ipr1.re_part + delta_pos;
        ipr1.dx_part            = ipr1.dx_part + delta_dx;
		
		% constrain shift
		% we do not move position beyond bounds set at the beginning
        ipr1.re_part            = max(ipr1.re_part, ipr1.re_part_lb);
        ipr1.re_part            = min(ipr1.re_part, ipr1.re_part_ub);
        % displacement is constrained by opts.dx_lim
		ipr1.dx_part            = max(ipr1.dx_part, opts.dx_lim(:,1)); % lower bound 
        ipr1.dx_part            = min(ipr1.dx_part, opts.dx_lim(:,2)); % upper bound
        
		%% evaluate residual of new system
		fprintf('I%02d S%02d checking residual of new system...\n', ipr0.riter, n);
		pimg1                   = ipr1.project();

		%% apply intensity correction
		fprintf('I%02d S%02d applying intensity correction...\n', ipr0.riter, n);
		I_rel_mean				= ones(1, n_part, ipr1.opts.precision);
		for ic = 1 : n_cameras
			%% size of this image
			n_i					= size(ipr1.img_meas{ic},1);
			n_j					= size(ipr1.img_meas{ic},2);

            %% accumulate jacobian over time
			% jacobian is of size n_px x n_part
			%
			% expected number of elements is:
			%	#visible particles x #support px 

			% indices of particles 
			ind                 = pimg1(1,ic).ind;
			assert(isequal(pimg1(1,ic).ind, pimg1(2,ic).ind));
			assert(isequal(pimg1(1,ic).ind, pimg1(3,ic).ind));
			n_visible			= size(ind, 2);
			jac_local			= sparse(n_i*n_j,n_visible);
			for rgb = 1 : 3
                %% map jacobian back to particle images
                % row and column coordinates of particle images
                row_map				= pimg1(rgb,ic).map(:,:,1);
                col_map				= pimg1(rgb,ic).map(:,:,2);
                % take care of out of bounds elements in map
                % dep_map is always in bounds by construction
                map_ok				= row_map >= 1 & row_map <= n_i ...
                                    & col_map >= 1 & col_map <= n_j;
                % create map of indices where map is valid
                ind_px				= zeros([n_support_px n_visible]);
                ind_px(map_ok)      = sub2ind([n_i n_j], row_map(map_ok), col_map(map_ok));

                %% make Jacobian matrix
				rows				= ind_px;
				cols     			= repmat(1:n_visible, n_support_px, 1);
                vals				= double(pimg1(rgb,ic).img);
				jac_local			= jac_local  + sparse(rows(map_ok), cols(map_ok), vals(map_ok), n_i*n_j, n_visible);
			end
			%% solve linear least squares problem
			% this identifies the ratio by which we should make each
			% particle brighter/darker in order to 
			I_rel				= lsqlin(jac_local, double(ipr1.img_meas{ic}(:))).';
			% limit relative change in brightness to 1/e - e
			I_rel				= min(exp(+1), max(exp(-1), I_rel));
			I_rel_mean(1,ind)	= I_rel_mean(1,ind) .* I_rel;
		end
		% geometric mean across cameras
		I_rel_mean			= I_rel_mean.^(1./ipr1.n_views(1,:)); 
		% multiplicative update with relaxation factor
 		ipr1.d_part			= ipr1.d_part .* sqrt(I_rel_mean);

% 		%% apply intensity correction
% 		% following each shake, a multiplicative intensity correction is
% 		% applied
% 		% I' = I * sqrt(sum(res) / sum(proj))
% 		% see section 2.2.3 of Schanz et al. 2016
% 		fprintf('I%02d S%02d applying intensity correction...\n', ipr0.riter, n);
% 		pimg_res_sum			= zeros(1, n_part, ipr0.opts.precision);
% 		pimg_ref_sum			= zeros(1, n_part, ipr0.opts.precision);
% 		pimg_sum				= zeros(1, n_part, ipr0.opts.precision);
% 		I_rel_mean				= ones(1, n_part);
% 		for ic = 1 : n_cameras
% 			% index of particle
% 			ind					= pimg1(1,ic).ind;
% 			% reference image and projection in each red/green/blue exposure
% 			pimg_proj			= cat(1, pimg1(1,ic).img, pimg1(2,ic).img, pimg1(3,ic).img);
% 			pimg_meas			= cat(1, pimg1(1,ic).ref, pimg1(2,ic).ref, pimg1(3,ic).ref);
% 			pimg_res			= pimg_meas - pimg_proj;
% 			% ideal multiplicative solution, find
% 			%	pimg_meas = I_rel * pimg_proj
% 			% where pimg_proj is the projected image
% 			%       pimg_meas is the measured image
% 			% amp is determined in the least squares sense
% 			XtX					= dot(pimg_proj, pimg_proj, 1);
% 			Xty					= dot(pimg_proj, pimg_meas, 1);
% 			I_rel				= Xty ./ XtX;
% 			I_rel				= max(0, I_rel);
% 			I_rel_mean(1,ind)	= I_rel_mean(1,ind) .* I_rel;		
% 		end
% 		% geometric average across views
% 		I_rel_mean				= I_rel_mean.^(1./ipr0.n_views(1,:));
% 		% multiplicative update with relaxation factor
% 		ipr0.d_part				= ipr0.d_part .* sqrt(I_rel_mean);

		%% apply update of estimate
		% we only accept changes that reduce the residual in the particle
		% image. if the residual is improved above a certain threshold,
		% then the change is accepted. Otherwise, the change is rejected, and
		% the lambda parameter will be adjusted 

		% determine which shifts have improved fit
		b_accept					= ipr1.pimg_res <= ipr0.pimg_res;
		n_accept					= nnz(b_accept);

		% apply shifts that improve fit and recalculate lambda
		ipr1.select( b_accept);
		ipr0.select(~b_accept);
		ipr1.lambda = max(opts.lambda_min, ipr1.lambda / opts.Ldown);
		ipr0.lambda = min(opts.lambda_max, ipr0.lambda * opts.Lup  );
		ipr0.append(ipr1);

		if opts.b_verbose
			fprintf('I%02d S%02d accepted %d of %d\n', ipr0.riter, n, n_accept, n_part);
		end
		
		%% report norm of residual
		if opts.b_verbose
			fprintf('I%02d S%02d RMS    residual %0.3f counts per voxel\n', ipr0.riter, n, sqrt(mean(ipr0.pimg_res)));
			fprintf('I%02d S%02d median residual %0.3f counts per voxel\n', ipr0.riter, n, median(sqrt(ipr0.pimg_res)));
		end
	end
end

