%% calcCFFricGMM
%
% GMM Objective function for model with frictions affecting pre-period
% earnings levels. 
%
% Inputs:
% 1. Data structure with vector of bunching amounts, tax structure
% 2. Parameter vector, p(1) = elasticity, p(2) = adjustment cost.
% 3. 

function [obj, bHat] = calcCFFricGMM(ep, wMat, data, optBounds, optInner, nDraws,density)

    % Preliminaries
        e = ep(1);
        phi = ep(2);

    % Calculate bunching at first kink

        % Find delta zStar1
        deltaZStar1 = data.tax1.zStar.*( (1-data.tax1.t0)./(1-data.tax1.t1) ).^e;

        % calculate zn0Upper, zn0Lower, and zn1Upper, bounds on
        % counterfactual earnings, and marginal adjuster (given optimal
        % earnings z, set by zGrid). 
        b1s = zeros(length(data.tax1.zStar),1);
        for ii = 1:length(data.tax1.zStar)

            if ~strcmp(density, 'flexible')
                zGrid = linspace(data.tax1.zStar(ii), deltaZStar1(ii), nDraws);  
            else

                % Map uniforms into relevant range
                prL = interp1(data.income + data.tax1.zStar(ii), data.dist, data.tax1.zStar(ii));
                prU = min(interp1(data.income + data.tax1.zStar(ii), data.dist, deltaZStar1(ii)),1);
                zDrawsShift = prL + (prU-prL).*data.zDraws;

                % Find hInverse(uniforms)
                zGrid = min(interp1(data.dist, data.income + data.tax1.zStar(ii), zDrawsShift),...
                    max(data.income+data.tax1.zStar(ii)));

            end
            assert (sum(isnan(zGrid))==0);
            [z0L, z0U] = calcBounds(zGrid, e, phi, data.tax0, optBounds);
            
            tempTax.zStar = data.tax1.zStar(ii);
            tempTax.t0 = data.tax1.t0(ii);
            tempTax.t1 = data.tax1.t1(ii);
            [~, z1U] = calcBounds(zGrid, e, phi, tempTax, optBounds);
            if ~strcmp(density, 'flexible')
                % trapezoid rule to evaluate integral
                    % nb: integral should multiply by stepsize*/(b-a), but
                    % stepsize is bia/(nDraws-1), so it cancels out
                approxInt1 = trapz(min((z0U - z1U)./(z0U - z0L),1))/(nDraws-1); 
                bSaez = deltaZStar1(ii) - data.tax1.zStar(ii);
            else
                approxInt1 = mean( min((z0U - z1U)./(z0U - z0L),1) );
                bSaez = calcBunchFlex(deltaZStar1(ii), data.tax1.zStar(ii), data)./data.densAK;
            end
            
            b1s(ii) = bSaez.*approxInt1;
        end
        
        
    % Calculate bunching at second kink
    
        % Find highest buncher
        phiBC = max(phi, 0);
        deltaZStar2 = data.tax2.zStar.*( (1-data.tax2.t0)/(1-data.tax2.t1) )^e + phiBC;
        zStart = deltaZStar2;
        [upperZ2,~,flag] = fsolve(@(zOld) uOpt(zOld,e,data.tax2) - uAK(zOld, e, data.tax2) - phiBC, zStart, optInner);
        upperZ2 = min(upperZ2, max(deltaZStar1)); % must have <= bunching in period 2 than 1
        if (flag ~= 1)
            flag
            e
            phi
        end
        assert (flag==1)
        
        % calculate zn0Upper, zn0Lower, and zn1Upper 
            % note this is inefficient since we calcualted these above
            % But calculating them twice keeps the grid sizes the same 
        if ~strcmp(density, 'flexible')
            zGrid = linspace(data.tax1.zStar(end), upperZ2, nDraws);
        else
            
            % Map uniforms into relevant range
            prL = interp1(data.income + data.tax1.zStar(end), data.dist, data.tax1.zStar(end));
            prU = min(interp1(data.income + data.tax1.zStar(end), data.dist, upperZ2),1);
            zDrawsShift = prL + (prU-prL).*data.zDraws;

            % Find hInverse(uniforms)
            zGrid = min(interp1(data.dist, data.income + data.tax1.zStar(end), zDrawsShift),...
                max(data.income + data.tax1.zStar(end)));

        end
        
        [z0L, z0U] = calcBounds(zGrid, e, phi, data.tax0, optBounds);
        tempTax.zStar = data.tax1.zStar(end);
        tempTax.t0 = data.tax1.t0(end);
        tempTax.t1 = data.tax1.t1(end);
        [~, z1U] = calcBounds(zGrid, e, phi, tempTax, optBounds);
        
        if ~strcmp(density, 'flexible')
            approxInt2 = trapz(min((z0U - z1U)./(z0U - z0L),1))/(nDraws-1); 
            bNoFric = upperZ2 - data.tax1.zStar(end);
        else
            approxInt2 = mean( min((z0U - z1U)./(z0U - z0L),1) );
            bNoFric = calcBunchFlex(upperZ2, data.tax1.zStar(end), data)./data.densAK;
        end

        
        b2s = bNoFric.*approxInt2;
        
    % Return quadratic loss 
        bHat = [b1s ; b2s];
        obj = (bHat - data.bunch)'*wMat*(bHat - data.bunch);

end

