%function [normdata, cms, calls] = norm_tetra(data_filename, positions_filename, num_sets);
%
% Normalize Polony sequencing data and display tetrahedra.
%
% Written by Greg Porreca, 01-31-2006; Church Lab, Harvard University
% Developed and originally written by Jay Shendure
%
%
function [normdata, cms, calls, coords_3d] = norm_tetra(raw_filtered_data, num_sets);

     sub_percentile = 1; %percentile value to subtract from all data (i.e. 1 means 1st percentile is subtracted)

     %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];

     num_cycles = (num_sets)/4;
     num_beads = size(raw_filtered_data, 1);

     %initialize return variables
     cms = zeros(num_beads, num_cycles + 7);
     calls = zeros(num_beads, num_cycles + 7);
     normdata = zeros(num_beads, (num_cycles*4));

     %build up matrix of identifying info for each bead; this will be
     %appended to normdata, cms, and calls
     id_matrix = raw_filtered_data(:,1:7);

     %remove id columns from raw_data for now; we will never reorder it, so
     %we can add id_matrix back later
     raw_filtered_data = raw_filtered_data(:,8:end);

     disp('    normalizing data...');

     %sort data so we can easily access percentile values
     raw_sorted_data = zeros(num_beads, num_sets);
     for(i=1:num_sets)
         raw_sorted_data(:,i) = sortrows(raw_filtered_data(:,i));
     end

     %compute a vector of values at the sub_percentile percentile
     sub_values = raw_sorted_data(floor(size(raw_filtered_data,1)/100)*sub_percentile,:);
     %sub_values = zeros(1, num_sets);
     %subtract the sub_percentile value from each row in a column, and start
     %creating the normdata matrix
     for(i=1:num_sets)
         normdata(:,i) = raw_filtered_data(:,i) - sub_values(i);
     end

     %turn all non-negative elements to 0, and report the changes
     normdata = normdata .* (normdata > 0);

     %now sum each column (set), divide by 300,000,000 and store this vector
     %set_sums = sum(raw_filtered_data) / 300000000;
     set_sums = sum(normdata) / 300000000;

     %now divide each value in each column by the appropriate sum just
     %computed
     for(i=1:num_sets)
         normdata(:,i) = (normdata(:,i) / set_sums(i));
     end

     %compute norms
     norms = zeros(num_beads, num_sets);
     for(i=1:4:num_sets)
         curr_norm = sqrt(sum( (normdata(:,i:i+3).*normdata(:,i:i+3))'))';
         norms(:,i:i+3) = curr_norm(:,ones(1,4)); %take the sum of each consecutive 4 columns, and put it into 4 corresponding positions in norms
     end
     normdata = normdata ./ (norms + eps);


     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %DATA IS NOW NORMALIZED; THE NEXT STEP IS TO MAKE BASECALLS, COMPUTE
     %CLUSTER MEDIANS, AND COMPUTE CMS SCORES FROM THESE MEDIANS
     %
     %
     %make initial basecalls by picking the highest of the 4 normalized vals for each cycle
     disp('    making initial basecalls...');
     calls(:,1:7) = id_matrix;
     for(i=1:num_cycles)
         [vals indx] = max(normdata(:,((i-1)*4)+1:((i-1)*4)+4)');
         calls(:,(7+i)) = indx';
     end

     %find cluster medians and calculate cms scores for each bead
     disp('    computing quality scores...');
     cms(:,1:7) = id_matrix;
     cluster_medians = zeros(num_cycles*4, 4);

     for(i=1:num_cycles)
         for(j=1:4) %base-identity
             cluster_medians(((i-1)*4)+j, :) = median(normdata(find(calls(:,i+7)==j),((i-1)*4)+1:((i-1)*4)+4)); %cluster_medians is NxM, N=#of sets, M=#of identities (4)
         end
     end

     %reassign calls based on minimum distance to cluster medians
     %
     %first, for each cycle, calculate the 4 distances to the 4 cluster
     %medians
     for(i=1:num_cycles) %iterate over the number of cycles (bases)
         for(j=1:4) %iterate over number of identities
             curr_cluster_med = cluster_medians(((i-1)*4)+j,:); %median value for the current cycle:identity is 4-D
             curr_cluster_med_matrix = curr_cluster_med(ones(size(normdata,1),1),:); %make a Nx4 matrix of these 4 vals (N=#of beads) for matrix subtract
             cm_distances(:,((i-1)*4)+j) = sqrt(sum(((normdata(:,(((i-1)*4)+1):(((i-1)*4)+4)) - curr_cluster_med_matrix) .* (normdata(:,(((i-1)*4)+1):(((i-1)*4)+4)) - curr_cluster_med_matrix))'))';
         end
     end
     %cm_distances now holds the distance of each bead, at each cycle, to
     %each of the 4 cluster medians (so, for example, cols 1-4 are the
     %distances of each bead to identity 1, 2, 3, and 4 for base 1; base 2
     %is cols 5-8, etc)

     %now we need to re-assign basecalls based on the identity closest to
     %its cluster median
     calls2(:,1:7) = id_matrix;
     for(i=1:num_cycles)
         [vals indx] = min(cm_distances(:,((i-1)*4)+1:((i-1)*4)+4)');
         cms(:,i+7) = vals'.*1000;
         calls2(:,(7+i)) = indx';
     end

     normdata = [id_matrix normdata];
     calls = calls2;

	 %generate list of 3-d coords for display
	 coords_3d = zeros(num_beads, num_cycles * 3);
	 for(i=1:num_cycles)
		 xyz = zeros(num_beads, 3);
		 for(j=1:4)
			 A = normdata(:,((i-1)*4)+j+7);
			 B = tetcoor(j,:);
			 xyz = xyz + B(ones(num_beads, 1),:) .* A(:, ones(3, 1));
		 end
		 coords_3d(:,((i-1)*3)+1:((i-1)*3)+3) = xyz;
	 end
	 
