% function [basecalls, basequals] = polbasecall(data_filename,num_sets, cmsthresh, primer_thresh, size_min, size_max, each_cycle_plot, final_tetrahedron_plot, output_fn)
%
% Call bases from .dataset files by normalizing, applying cms threshold
% (supply a percentage of data you want to pass the cms thresh as parameter
% cmsthresh), and re-normalizing the remaining subset.  Tetrahedra of the
% final data are diaplayed if tetrahedron_plot==1.  num_sets is the number
% of image sets, not including the primer images (which are assumed to
% occur first).  So, for an experiment where 26 bases are produced, one
% would set num_sets=26.
%
% To write the data to a file, specify the filename as the last (optinal)
% argument
%
% Written by Greg Porreca, 01-31-2006; Church Lab, Harvard University
%
%
function [basecalls, basequals] = polbasecall(data_filename, num_sets, cmsthresh, primer_thresh, size_min, size_max, each_cycle_plot, final_tetrahedron_plot, output_fn)

%primer_thresh = 8000;
%size_min = 1;
%size_max = 5;

    %define coefficients for transformation to 3 dimensions
     tetcoor = [0 0 1; 0 0.9428 -0.3333; -0.8165 -0.4714 -0.3333; 0.8165 -0.4714 -0.3333];


%load input file; we need to load it in chunks, both because we don't want
%to normalize globally and also because we don't have enough memory for the
%whole thing
%
%we will shoot for approx 20 frames per subset, and we will assume there
%are 600 frames.  this means we will load 1/30th of the full file in each
%batch.  We conpute this by taking the full size of the file, and dividing
%it by the number of sets + 7 (first 7 cols are bead id info).  This is
%twice the number of beads (since each field is 2 bytes).
est_numframes = 600;
numframes_per_subset = 20;
%num_subsets = (est_numframes / numframes_per_subset)-1;
num_subsets = (est_numframes / numframes_per_subset);

fileinfo = dir(data_filename);
filesize = fileinfo.bytes;
total_num_beads = filesize / 2;
total_num_beads = total_num_beads / (num_sets+7);
%subset_size = floor(total_num_beads / (num_subsets+1));
subset_size = floor(total_num_beads / (num_subsets));

disp(['==>polbasecall running over a total of ' num2str(total_num_beads) ' beads in ' num2str(num_subsets) ' subsets (' num2str(subset_size) ' per subset)']);

fid = fopen(data_filename);

basecalls = [];
basequals = [];
tetrahedron_coords = [];

for(j=1:num_subsets)
%for(j=1:1)
    disp(['=======> Currently on subset ' num2str(j) ' of ' num2str(num_subsets) ' <=======']);
    if(j<num_subsets)
        raw_data = uint16(fread(fid, subset_size * (num_sets+7), 'uint16'));
    else
        raw_data = uint16(fread(fid, inf, 'uint16'));
    end    
    raw_data_temp = reshape(raw_data, num_sets+7, [])';
	raw_data = [raw_data_temp(:,1:5) raw_data_temp(:,end) raw_data_temp(:,6:end-1)];
	clear raw_data_temp;
	raw_data = double(raw_data);

	%make sure nothing is funny
	if (mod(num_sets, 4)~=0)
		disp('You entered an invalid number of sets (must be a multiple of 4)');
		return;
	end
	%raw_data = raw_data(:,1:end-4);
	%num_sets = num_sets - 4;
	num_cycles = (num_sets)/4;
	
	%threshold data and report
	raw_filtered_data = raw_data(find(raw_data(:,7) > primer_thresh),:);
	raw_filtered_data = raw_filtered_data(find(raw_filtered_data(:,3) >= size_min),:); %apply bead size_min thresh
	raw_filtered_data = raw_filtered_data(find(raw_filtered_data(:,3) <= size_max),:); %apply bead size_max thresh
	raw_filtered_data = raw_filtered_data(find(raw_filtered_data(:,6) == 0),:); %only use beads not on the border of frames
	num_beads = size(raw_filtered_data, 1);
	
	disp(['==>Started using ' data_filename ' with primer thresh of ' num2str(primer_thresh) ' and ' num2str(size_min) '<=size<=' num2str(size_max)]);
	disp(['==>Out of ' num2str(size(raw_data,1)) ' objects, ' num2str(num_beads) ' met the thresholds above']);
	
	max(raw_filtered_data(:,1));
	
	%call norm_tetra to generate initial basecalls and cms data so we can
	%select the best subset
	[normdata, cms, calls, coords3d] = norm_tetra(raw_filtered_data, num_sets);
	cms_sums = cms(:,8:end);
	if(size(cms_sums,2) > 1)
		cms_sums = sum(cms_sums')';
    end
    
    if(each_cycle_plot==1)
        figure, hist(cms_sums, 100);
    end
	
	%now throw out the worst frames and re-normalize only the good frames	
	cms_framemean = [];
	for(i=min(cms(:,1)):max(cms(:,1)))
		if(size(cms(:,8:end),2)==1)
			cms_framemean(i) = mean(cms(find(cms(:,1)==i),8:end));
		else
			cms_framemean(i) = mean(sum(cms(find(cms(:,1)==i),8:end)')');
		end
	end
	cms_framemean_s = sortrows(cms_framemean');
%	cms_frame_thresh = cms_framemean_s(end-2);
    cms_frame_thresh = 290 * (num_sets / 4);    
    
    if(each_cycle_plot==1)
        figure, hist(double(cms_framemean(min(cms(:,1)):max(cms(:,1)))));
    end
    positions = load('position_list.dat');
    positions(:,2) = positions(:,2) .* -1;
%	signal_vs_position(positions(1:size(cms_framemean,2),2:3), cms_framemean,1,0);
	
	disp(['==>Selecting beads in frames which meet CMS mean threshold of ' num2str(cms_frame_thresh) ' (best 90% of frames)...']);
	good_frames = [];
	good_cms = [];
	good_raw_filtered_data = [];
	for(i=min(cms(:,1)):max(cms(:,1)))
		if(cms_framemean(i) < cms_frame_thresh)
			good_frames = [good_frames i];
			good_cms = [good_cms; cms(find(cms(:,1)==i),:)];
			good_raw_filtered_data = [good_raw_filtered_data; raw_filtered_data(find(raw_filtered_data(:,1)==i),:)];
		end
	end
	num_beads_precms = num_beads;
	num_beads = size(good_raw_filtered_data,1);
	disp(['==>Of ' num2str(num_beads_precms) ' beads in all frames, ' num2str(num_beads) ' were in "good CMS" frames (' num2str(size(find(cms_framemean < cms_frame_thresh),2)) '/' num2str(size(cms_framemean,2)) ' frames)']);
	
	disp(['==>Selecting the best global ' num2str(cmsthresh) '% of beads on the basis of sum CMS..']);
	cms_sums = good_cms(:,8:end);
	if(size(cms_sums,2) > 1)
		cms_sums = sum(cms_sums')';
    end

    if(each_cycle_plot==1)
    	figure, hist(cms_sums, 100);
    end

    sortedcms = sortrows(cms_sums);
	cms_thresh = sortedcms(floor((num_beads / 100) * cmsthresh));

	%select the best subset
	orig_size = size(good_raw_filtered_data,1);
	cms_vals = good_cms(:,8:end);
	if(size(cms_vals,2) > 1)
		cms_filtered_data = good_raw_filtered_data(find(sum(cms_vals')' < cms_thresh),:);
	else
		cms_filtered_data = good_raw_filtered_data(find(cms_vals < cms_thresh),:);
	end
	num_beads_precms = num_beads;
	num_beads = size(cms_filtered_data, 1);
	disp(['==>Of the ' num2str(num_beads_precms) ' beads, ' num2str(num_beads) ' passed CMS threshold of ' num2str(cms_thresh) '.']);
	
    
	%re-run the normalization and basecalling w/ this subset, generating
	%tetrahedra if appropriate
	[normdata2, subset_quals, subset_calls, coords3d] = norm_tetra(cms_filtered_data, num_sets);
    if(each_cycle_plot==1)
    	for(i=1:4:num_sets)
        	curr_3d = ((((i-1)/4) + 1)*3) - 2;
            tetrahedron(coords3d(:,curr_3d:curr_3d+2), subset_calls(:, ((i-1)/4) + 8 ));
        end
    end
        
    
	tetrahedron_coords = [tetrahedron_coords; coords3d];
	cms_sums = subset_quals(:,8:end);
    if(size(cms_sums,2) > 1)
        cms_sums = sum(cms_sums')';
    end
    
    if(each_cycle_plot==1)
        figure, hist(cms_sums, 100);
        pause;
    end
    
	basequals = [basequals; subset_quals];
	basecalls = [basecalls; subset_calls];
	close all;
end

%plot tetrahedra?
if(final_tetrahedron_plot == 1)
    qualsum = zeros(size(positions, 1),num_sets/4);

	%OUTPUT BMP TO QC DIRECTORY
	ci = fopen(['../../EXEC/base_order.dat'],'r');
	cycle_names = textscan(ci, '%s');
	fclose(ci);
	ci = fopen(['../../IMAGES/RAW/cycle_id_table.txt'], 'r');
	cycle_nums = textscan(ci, '%s%s%s', 'delimiter', '\t', 'whitespace', '');
	fclose(ci);
	%
	
	for(i=1:4:num_sets)
        
		%OUTPUT BMP TO QC DIRECTORY
		curr = cell2mat(cycle_names{1}(int32(i/4)+1));
		for(j=1:4:num_sets)
			if(size(curr,2)==size(cell2mat(cycle_nums{2}(j)),2))
				if(curr == cell2mat(cycle_nums{2}(j)))
					fret_num = str2num(cell2mat(cycle_nums{1}(j)));
					cy5_num = fret_num+1;
					cy3_num = fret_num+2;
					txred_num = fret_num+3;
				end
			end
		end
		if(exist('fret_num') > 0)
			fret_name = sprintf('%03d', fret_num);
			cy5_name = sprintf('%03d', cy5_num);
			cy3_name = sprintf('%03d', cy3_num);
			txred_name = sprintf('%03d', txred_num);
		
			cycle_title = [curr '--' fret_name '-' txred_name];
			cycle_filename = ['../../QC/tetra_' cycle_title '.bmp'];
		end
		%
		
		curr_3d = ((((i-1)/4) + 1)*3) - 2;
		tetrahedron(tetrahedron_coords(:,curr_3d:curr_3d+2), basecalls(:, ((i-1)/4) + 8 ));

		if(exist('fret_num') > 0)
			title(cycle_title);
			saveas(gcf, cycle_filename, 'bmp');
		end
	end
	
    for(i=1:num_sets/4)
        for(j=min(basequals(:,1)):max(basequals(:,1)))
            qualsum(j,i) = sum(basequals(find(basequals(:,1)==j),i+7));
        end
        signal_vs_position(positions(:,2:3), qualsum(:,i), 1, 0);
    end
end

%save the results?
if(exist('output_fn','var')==1)
	disp(['==>Writing data to file ' output_fn]);
	output_data = [basecalls basequals(:,8:end)];
	save(output_fn, 'output_data', '-ASCII');
end
