function [Surf,Scores,ScoreLogMultipliers] = MethodII(f,nusteps,sigmasteps,ySteps,Report)

% [Surf,Scores,ScoreLogMultipliers] = MethodII(f,nusteps,sigmasteps,ySteps,Report)
%
% Five input arguments, the last two being optional
%
% f = the observed data, corresponding to a tree of observed fluorescence
%     values; each row is a cell, and each column is an observation; NaN
%     indicates values not observed
% nuSteps = values of nu to test (pass a single value or a vector)
% sigmaSteps = values of sigma to test (pass a single value or a vector)
% ySteps = increments for numerical integration over the y's
%          should be of the form 0:X:Z for some X>0 and Z>0
%          defaults to 0:25:20000;
% Report = 0 - returns P(f|nu,sigma)
%        = 1 - returns log of P(f|nu,sigma)
%        = 2 - returns P(f|nu,sigma) normalized to sum to 1
%        defaults to 1
%
% Three output arguments, though the last two are more for debugging
%
% Surf = The requested result, depending on Report
%
% For more details, see: 
% Nitzan Rosenfeld, Theodore J. Perkins, Uri Alon, Michael B. Elowitz, and Peter S. Swain
% A Fluctuation Method to Quantify In Vivo Fluorescence Data
% Biophys. J.  2006 91:  759-766; Published ahead of print on April 28, 2006 as doi:10.1529/biophysj.105.073098.

if nargin<4
    DeltaY = 25;
    ySteps = 0:DeltaY:20000;
else
    DeltaY = ySteps(2)-ySteps(1);
end

if nargin<5
    Report = 1;
end

Scores = zeros(length(nusteps),length(sigmasteps));
ScoreLogMultipliers = zeros(length(nusteps),length(sigmasteps));
[NCells,NObs] = size(f);
tic;
for nui = 1:length(nusteps)
    for sigmai = 1:length(sigmasteps)
        if nui>1 | sigmai>1
            TimeSpent = toc;
            NDone = (nui-1)*length(sigmasteps)+(sigmai-1);
            TotalTimeEst = TimeSpent * length(sigmasteps) * length(nusteps) / NDone;
            EstTimeLeft = TotalTimeEst-TimeSpent;
            ETLHours = floor(EstTimeLeft/3600);
            ETLMins = floor((EstTimeLeft-ETLHours*3600)/60);
            ETLSecs = EstTimeLeft-ETLHours*3600-ETLMins*60;
            if sigmai==length(sigmasteps)
                disp(['nu ',num2str(nui),'/',num2str(length(nusteps)),' sigma ',num2str(sigmai),'/',num2str(length(sigmasteps)),' etl ',num2str(ETLHours),' hr ',num2str(ETLMins),' min ',num2str(ETLSecs), 'sec']);
            end
            %else
            %disp(['nu ',num2str(nui),'/',num2str(length(nusteps)),' sigma ',num2str(sigmai),'/',num2str(length(sigmasteps)),' etl (min) ?']);
        end
        Z = cell(2*NCells+1,1);
        Multipliers = ones(NCells,1);
        for c=NCells:-1:1
            Z{c} = MethodII_ProcessNode_V2(f,nusteps(nui),sigmasteps(sigmai),ySteps,c,Z{2*c},Z{2*c+1});
            %Z{c} = Z{c}*10^3;
            if ~isempty(Z{c})
                if max(Z{c}>0)
                    Multipliers(c) = 1/max(Z{c});
                    Z{c} = Z{c}*Multipliers(c);
                end
            end
        end
        Scores(nui,sigmai) = DeltaY*sum(Z{1});
        I = find(Multipliers>0);
        ScoreLogMultipliers(nui,sigmai) = -sum(log(Multipliers(I)));
    end
end

if Report==0
    Surf = Scores.*exp(ScoreLogMultipliers);
end
if Report==1
    Surf = log(Scores)+ScoreLogMultipliers;
end
if Report==2
    Surf = Scores.*exp(ScoreLogMultipliers-max(max(ScoreLogMultipliers)));
    Norm = sum(sum(Surf));
    if Norm>0
        Surf = Surf/Norm;
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Z = MethodII_ProcessNode_V2(f,nu,sigma,ySteps,Node,LChildZ,RChildZ)

%%% WARNING! THIS CODE ASSUMES YSTEPS STARTS AT 0 AND THAT IT GOES IN EVEN INCRMENTS OF SOME REAL X>0. 
%%% I.E., IT ASSUMES YSTEPS = [0 X 2X 3X 4X ... AX] FOR SOME REAL X>0 AND INTEGER A>0.

L = length(ySteps);
[R,C] = size(f);
DeltaY = ySteps(2);
if all(isnan(f(Node,:)))
    Z = [];
else
    Z = ones(1,L);
    
    % Account for noisy observations of y
    for c=1:C
        if ~isnan(f(Node,c))
            Z = Z*(1/sigma).*exp(-(ySteps-f(Node,c)).^2/(2*sigma^2));
        end
    end
    
    % Account for fluctuations in number of molecules of children
    if ~isempty(LChildZ) || ~isempty(RChildZ)
        % If there's nothing at the current node, there must be nothing at
        % the children too.
        %Z(1) = 1;
        if ~isempty(LChildZ)
            Z(1) = Z(1)*LChildZ(1);
        end
        if ~isempty(RChildZ)
            Z(1) = Z(1)*RChildZ(1);
        end
        % Otherwise, the amount at the children must sum to be the amount
        % at the current node.
        Temp = [0 DeltaY*(1/sqrt(nu))*(1./sqrt(ySteps(2:end)))];
        for i=2:L
            if Z(i)>0 & Temp(i)>0 
                ZMult = exp(-(2*ySteps(1:i)-ySteps(i)).^2/(2*nu*ySteps(i)));
                if ~isempty(LChildZ)
                    ZMult = ZMult.*LChildZ(1:i);
                end
                if ~isempty(RChildZ)
                    ZMult = ZMult.*RChildZ(i:-1:1);
                end
                Z(i) = Z(i)*sum(ZMult)*Temp(i);
            end
        end
    end
end



