% Poisson_EstimationScript

% The embryos
if 1
    Emb = {'BECac6', 'BECaf3', 'BECah11', 'BECah6', 'BECah8', 'BECah9', 'BECai1', 'BECaj2', 'BECaj8', 'BECaj9', 'BECak1', 'FEScb13', 'FEScc01', 'FEScc34', 'FEScc42', 'FEScd15', 'FESce12', 'FESce20', 'FESce21', 'FEScf04', 'FEScg20', 'FESch02', 'FESch04', 'FESch06', 'HETaa16', 'HETaa4', 'HETaa9', 'HETab11', 'HETab12', 'HETab13', 'HETac14', 'HETaf5', 'HETaf6', 'HETaf9', 'HETbb1', 'HETbc23', 'HETbe21', 'HETbe6', 'HETbe9', 'ba3', 'bd5', 'bk1', 'bk15', 'br12', 'ca1', 'ca14', 'dm14', 'ds1', 'em37', 'fa12', 'fq1', 'fq4', 'gq10', 'gq13', 'gq14', 'gq21', 'gq4', 'gr12', 'gr16', 'gr3', 'hne2', 'hne4', 'hne8', 'jr1', 'jr16', 'kd17', 'kd26', 'kd31', 'kd9', 'kf3', 'kf9', 'kw14', 'kw2', 'kw7', 'lesl15', 'lesl16', 'lesl20', 'lesl21', 'lesl26', 'mt1', 'mu15', 'mu20', 'mu7', 'nb1', 'nb11', 'nb5', 'ne11', 'ne12', 'nk11', 'nk14', 'nk5', 'pd4', 'pe2', 'ra6', 'ra7', 'rb2', 'rb8', 'rf11', 'rf6', 'rge7', 'rge9', 'sasha3', 'sasha4', 'se10', 'sk2', 'sk5', 'sl30', 'sp1', 'ta2', 'ta3', 'tb6', 'tf4', 'tg11', 'tg14', 'tg23', 'th9', 'tim11', 'tim12', 'tim15', 'tim18', 'tim2', 'tn2', 'tu1', 'tu3', 'we11', 'wn1', 'wn17', 'wn8', 'wp2', 'wp7'};
end


% The following computes the Anterior-Posterior coordinate transformation
% used by the local regression scheme
if 0
    CoordParams = zeros(length(Emb),9);
    for i=1:length(Emb);
        % Visual output, just to keep track of where we are.  The full
        % optimization may take some days!
        disp([num2str(i),' ',Emb{i}]);
        
        % Load the nuclear data
        Fname = ['Poisson_Embryos/',Emb{i},'.txt'];
        QD = LoadQuantData(Fname);

        % For each channel, check if it is a pair-rule gene.  If so,
        % optimize the transformation.
        for Chan=1:3
            GNum = GeneNameToNum(QD.Genes{Chan});
            % If pair-rule gene, optimize coords
            if GNum==4 | GNum==5 | GNum==7 | GNum==10 | GNum==11 | GNum==12 | GNum==13
                S_AP0_DV0 = Poisson_OptimizeCoords(QD.Data(:,2),QD.Data(:,3),QD.Data(:,3+Chan));
                CoordParams(i,3*(Chan-1)+1) = S_AP0_DV0(1);
                CoordParams(i,3*(Chan-1)+2) = S_AP0_DV0(2);
                CoordParams(i,3*(Chan-1)+3) = S_AP0_DV0(3);
            end
        end
    end
    save -ascii Poisson_OptimizedCoords.txt CoordParams
end


% Next, we estimate the expected intensity of each nucleus in each embryo.
if 1
    IntensityEst = cell(1,length(Emb));
    CoordParams = load('-ascii','Poisson_OptimizedCoords.txt');
    for i=1:length(Emb)
        % Visual display of progress
        disp([num2str(i),' ',Emb{i}]);

        % Load the nuclear data
        Fname = ['Poisson_Embryos/',Emb{i},'.txt'];
        QD = LoadQuantData(Fname);

        % For each channel, estimate, with method depending on whether
        % its a pair-rule gene or not.
        for Chan=1:3
            GNum = GeneNameToNum(QD.Genes{Chan});
            DV = QD.Data(:,3);
            Expr = QD.Data(:,3+Chan);
            if GNum==4 | GNum==5 | GNum==7 | GNum==10 | GNum==11 | GNum==12 | GNum==13
                AP = Poisson_Recoordinate(QD.Data(:,2),QD.Data(:,3),...
                    CoordParams(i,3*(Chan-1)+1), ...
                    CoordParams(i,3*(Chan-1)+2), ...
                    CoordParams(i,3*(Chan-1)+3));
                Est = Poisson_KNNQuadfit(AP,DV/3,Expr);
            else
                AP = QD.Data(:,2);
                Est = Poisson_KNNQuadfit(AP,DV,Expr);
            end
            IntensityEst{i}(1:QD.NNuc,Chan) = Est;
        end
    end
    save Poisson_EstimatedIntensities.mat IntensityEst
end


% Construct the nu estimates
% First, we have to go through all embryos and collect everything by gene
if 1
    ObsAndEsts = cell(1,14);
    load Poisson_EstimatedIntensities.mat;
    for i=1:length(Emb)
        % Visual display of progress
        %disp([num2str(i),' ',Emb{i}]);

        % Load the nuclear data
        Fname = ['Poisson_Embryos/',Emb{i},'.txt'];
        QD = LoadQuantData(Fname);

        % Go through each channel
        for Chan=1:3
            GNum = GeneNameToNum(QD.Genes{Chan});
            ObsAndEsts{GNum} = [ObsAndEsts{GNum}; ...
                QD.Data(:,3+Chan) IntensityEst{i}(:,Chan)];
        end
    end
end
% Next, we go through by genes, and make our estimates
if 1
    for g=1:14
        if ~isempty(ObsAndEsts{g})
            % We must discard points with non-positive predictions
            I = find(ObsAndEsts{g}(:,2)>0);
            % Compute mean squared error (for curiosity or x-val)
            MSQ = mean((ObsAndEsts{g}(I,1)-ObsAndEsts{g}(I,2)).^2);
            % Construct the pointwise estimates
            NuEsts = ((ObsAndEsts{g}(I,1)-ObsAndEsts{g}(I,2)).^2)./ ...
                (ObsAndEsts{g}(I,2));
            % Discard 10% largest, and average the rest
            NuEsts = sort(NuEsts);
            NuEsts = NuEsts(1:round(0.9*length(NuEsts)));
            NuEst = mean(NuEsts);
            % Print out report
            disp([GeneNumToName(g),' ',num2str(MSQ),' ',num2str(NuEst),' ',num2str(255/NuEst)]);
        end
    end
end