function [breakpoint,smallest_summed_deviance,confidence]=change_point_analysis_deviance_method(input_data,bootstrap_runs,do_plotting)
%this is modeled on o'connor, wittenberg and wang, 2005
%possible future additions include the consideration of a second candidate change point, to be identified at the same time as the first one, and bootstrap to obtain confidence intervals

smallest_summed_deviance=sum((input_data-mean(input_data)).^2);
breakpoint=0;

for candidate_breakpoint=1:length(input_data)-1
    summed_deviance=sum((input_data(1:candidate_breakpoint)-mean(input_data(1:candidate_breakpoint))).^2)...
        +sum((input_data((candidate_breakpoint+1):end)-mean(input_data((candidate_breakpoint+1):end))).^2);
    if summed_deviance<smallest_summed_deviance
        smallest_summed_deviance=summed_deviance;
        breakpoint=candidate_breakpoint;
    end
end

if do_plotting
	plot(input_data, 'ko')
	hold on
	plot([1 breakpoint],[mean(input_data(1:breakpoint)) mean(input_data(1:breakpoint))],'r')
	plot([breakpoint+1 length(input_data)],[mean(input_data((breakpoint+1):end)) mean(input_data((breakpoint+1):end))],'r')
end

if bootstrap_runs
	original_total_deviance=sum((input_data-mean(input_data)).^2);
	smallest_shuffled_summed_deviances=[];
	for i=1:bootstrap_runs
        smallest_shuffled_summed_deviance=original_total_deviance;
        shuffled_data=shuffle_randomization(input_data);
        %shuffled_breakpoint=0;
		for candidate_breakpoint=1:length(shuffled_data)-1
            summed_deviance=sum((shuffled_data(1:candidate_breakpoint)-mean(shuffled_data(1:candidate_breakpoint))).^2)...
                +sum((shuffled_data((candidate_breakpoint+1):end)-mean(shuffled_data((candidate_breakpoint+1):end))).^2);
            if summed_deviance<smallest_shuffled_summed_deviance
                smallest_shuffled_summed_deviance=summed_deviance;
                %shuffled_breakpoint=candidate_breakpoint;
            end
		end
        smallest_shuffled_summed_deviances=[smallest_shuffled_summed_deviances smallest_shuffled_summed_deviance];
	end
	confidence=100*(length(find(smallest_shuffled_summed_deviances>smallest_summed_deviance))/bootstrap_runs);
else
    confidence=nan;
end