%%  wiebeFit.m
%
%   [FitCoeffs,CAparams,Wiebes] = wiebeFitBF(NumWiebes,Breakdown,CAs,HRR,CAparams,Guesses,Ranges,VecSize,a)
%
%   Fits a Wiebe function to a heat release rate curve via Brute Force
%
%   NumWiebes = # of sumltaneously brute-forced Wiebe functions, 1, 2, or 3
%   CAs: Vector of crank angles
%   HRR: Vector of heat release rates
%   CAparams: 3-element vector, 1st element = start of CA window, 2nd
%             element = end of CA window, 3rd element = CA interval
%   Guesses: NumWiebes-by-4 matrix of starting guesses for the wiebe
%            function parameters, 1st col = theta_0, 2nd col = delta theta,
%            3rd col = Q, 4th col = M
%   VecSize = Number of possibilities to try for each parameter
%   EfficiencyFactor = Efficiency factor "a" in the Wiebe functions


function [FitCoeffs,CAparams,Wiebe1s] = wiebeFit(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor)
    options = optimoptions('fmincon','FiniteDifferenceStepSize', 0.1);
    [BestCoeffs, fval] = fmincon(@(CoeffArray) OBJ(CoeffArray,NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor), [Guesses(:,1), Guesses(:,2), Guesses(:,3), Guesses(:,4)], options);
    % BestCoeffs(1,:) = cell2mat(CoeffArray(1,:)); %converts a cell array into an ordinary array (array in matrix)

    FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction
    FitCoeffs = zeros(size(BestCoeffs));
    for WiebeNum = 1:NumWiebes
        for CoeffNum = 1:4
            FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum));
        end
    end
    [Err, Wiebels] = OBJ(BestCoeffs, NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor);
end

function [Errs, Wiebels] = OBJ(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor)
    NewCAs = CAparams(1):CAparams(3):CAparams(2);
    NewHRR = zeros(1,length(NewCAs));
    for CAnum = 1:length(NewCAs)
        CA = NewCAs(CAnum);
        [~,CAind] = min(abs(CAs-CA)); %locate the index
        Err = CAs(CAind)-CA;
        if Err == 0
            NewHRR(CAnum) = HRR(CAind);
        else
            if CAind < 2
                Delta = HRR(CAind+1)-HRR(CAind);
                InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg
            else
                Delta = HRR(CAind)-HRR(CAind-1);
                InterpFac = Err/(CAs(CAind)-CAs(CAind-1));
            end
            NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac;
        end
    end
    NumCAs = length(NewCAs);

    NumCoeffs = NumWiebes*4;
    Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11    
    Coeffs = [];
    Coeffs = ones(NumWiebes,4,VecSize);
    for SweepNum = 1:VecSize
        Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges;
    end
    Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor);

%     NumPoss = size(CoeffArray,1);
    NumPoss = 1;
    
    Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension
    Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array

    Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension
    Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array

    if NumWiebes > 1;         
        for Poss2 = 1:NumPoss
            Coeff2 = CoeffArray(Poss2,:);
            Wiebe2 = Wiebes(2,:,Coeff2{:}); 
            NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2;
            for Poss3 = 1:NumPoss
                Coeff3 = CoeffArray(Poss3,:);
                Wiebe3 = Wiebes(3,:,Coeff3{:});
                NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3;
                Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR);
            end
        end
    else
        Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row
        [~,NewInds] = min(Errs);
    end  
end

% function [FitCoeffs,CAparams,Wiebe1s] = wiebeFit(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor)
%     options = optimoptions('fmincon','FiniteDifferenceStepSize', 0.1);
%     [BestCoeffs, fval] = fmincon(@(CoeffArray) OBJ(CoeffArray,NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor), [Guesses(:,1), Guesses(:,2), Guesses(:,3), Guesses(:,4)], options);
%     % BestCoeffs(1,:) = cell2mat(CoeffArray(1,:)); %converts a cell array into an ordinary array (array in matrix)
% 
%     FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction
%     FitCoeffs = zeros(size(BestCoeffs));
%     for WiebeNum = 1:NumWiebes
%         for CoeffNum = 1:4
%             FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum));
%         end
%     end
%     [Err, Wiebels] = OBJ(BestCoeffs, NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor);
% end
% 
% function [Errs, Wiebels] = OBJ(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor)
%     NewCAs = CAparams(1):CAparams(3):CAparams(2);
%     NewHRR = zeros(1,length(NewCAs));
%     for CAnum = 1:length(NewCAs)
%         CA = NewCAs(CAnum);
%         [~,CAind] = min(abs(CAs-CA)); %locate the index
%         Err = CAs(CAind)-CA;
%         if Err == 0
%             NewHRR(CAnum) = HRR(CAind);
%         else
%             if CAind < 2
%                 Delta = HRR(CAind+1)-HRR(CAind);
%                 InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg
%             else
%                 Delta = HRR(CAind)-HRR(CAind-1);
%                 InterpFac = Err/(CAs(CAind)-CAs(CAind-1));
%             end
%             NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac;
%         end
%     end
%     NumCAs = length(NewCAs);
% 
%     NumCoeffs = NumWiebes*4;
%     Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11    
%     Coeffs = [];
%     Coeffs = ones(NumWiebes,4,VecSize);
%     for SweepNum = 1:VecSize
%         Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges;
%     end
%     Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor);
% 
% %     NumPoss = size(CoeffArray,1);
%     NumPoss = 1;
%     
%     Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension
%     Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array
% 
%     Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension
%     Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array
% 
%     if NumWiebes > 1;         
%         for Poss2 = 1:NumPoss
%             Coeff2 = CoeffArray(Poss2,:);
%             Wiebe2 = Wiebes(2,:,Coeff2{:}); 
%             NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2;
%             for Poss3 = 1:NumPoss
%                 Coeff3 = CoeffArray(Poss3,:);
%                 Wiebe3 = Wiebes(3,:,Coeff3{:});
%                 NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3;
%                 Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR);
%             end
%         end
%     else
%         Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row
%         [~,NewInds] = min(Errs);
%     end  
% end

% function [FitCoeffs,CAparams,Wiebe1s] = TwiebeFitBruteForce(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor)
% 
% NewCAs = CAparams(1):CAparams(3):CAparams(2);
% NewHRR = zeros(1,length(NewCAs));
% for CAnum = 1:length(NewCAs)
%     CA = NewCAs(CAnum);
%     [~,CAind] = min(abs(CAs-CA)); %locate the index
%     Err = CAs(CAind)-CA;
%     if Err == 0
%         NewHRR(CAnum) = HRR(CAind);
%     else
%         if CAind < 2
%             Delta = HRR(CAind+1)-HRR(CAind);
%             InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg
%         else
%             Delta = HRR(CAind)-HRR(CAind-1);
%             InterpFac = Err/(CAs(CAind)-CAs(CAind-1));
%         end
%         NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac;
%     end
% end
% NumCAs = length(NewCAs);
% 
% NumCoeffs = NumWiebes*4;
% Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11    
% Coeffs = [];
% Coeffs = ones(NumWiebes,4,VecSize);
% for SweepNum = 1:VecSize
%     Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges;
% end
% Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor);
% 
% CoeffArray = num2cell(permn(1:VecSize,4)); %converts permn... into an cell array 4 (11^4 different combinations)
% NumPoss = size(CoeffArray,1);
% 
% Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension
% Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array
% 
% BestError = inf;
% BestCoeffs = zeros(NumWiebes,4);
% 
% if NumWiebes > 1;         
%     for Poss2 = 1:NumPoss
%         Coeff2 = CoeffArray(Poss2,:);
%         Wiebe2 = Wiebes(2,:,Coeff2{:}); 
%         NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2;
%         for Poss3 = 1:NumPoss
%             Coeff3 = CoeffArray(Poss3,:);
%             Wiebe3 = Wiebes(3,:,Coeff3{:});
%             NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3;
%             Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR);
%             [NewError,NewInds] = min(Errs);
%             if NewError < BestError
%                 BestError = NewError;
%                 BestCoeffs(1,:) = cell2mat(CoeffArray(NewInds,:)); %converts a cell array into an ordinary array
%                 BestCoeffs(2,:) = cell2mat(CoeffArray(Poss2,:));
%                 BestCoeffs(3,:) = cell2mat(CoeffArray(Poss3,:));
%             end
%         end
%     end
% else
%     Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row
%     [~,NewInds] = min(Errs);
%     BestCoeffs(1,:) = cell2mat(CoeffArray(NewInds,:)); %converts a cell array into an ordinary array (array in matrix)
% end
% 
% FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction
% FitCoeffs = zeros(size(BestCoeffs));
% for WiebeNum = 1:NumWiebes
%     for CoeffNum = 1:4
%         FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum));
%     end
% end
% 
% end