function [mu, sig, dmu, dsig] = SOFM_PS(fun,y,x,CovVec,CovVal,step1,epsHat)
%Written by Jan Christoph Krueger, Hamburg University of Technology
%Implementation of the principal sensitivity second-order
%fourth-moment method. Theoretical details are given in

%OUTPUT
% mu: mean
% dmu: gradient of mean with respect to design variables
% sig: standard deviation
% dsig: gradient of standard deviation with respect to design variables
%INPUT
% fun: objective function such that [f, dfdy, dfdx]=fun(y,x);
% y: Current design variables
% x: Current random variables
%CovVec: Matrix including the eigenvectors of the Covariance
%matrix
%CovVal: Column-Vector of eigenvalues of the Covariance matrix
%step1: Step Size deltax
%epsHat: normalized step size of second order principal
%sensitivity directions

    %memory allocation
    lx = size(CovVec, 2); 
    ly = length(y);
    dfxx = zeros(lx, lx); 
    dfyx = zeros(ly, lx);
    dfydata = zeros(ly, 2*lx); 
    dfxdata = zeros(lx, 2*lx);
    fdata = zeros(2*lx,1);

    % Step 1:Function evalutaion at unperturbed design
    [f, dfy0, dfx0] = fun(y,x); 
    dfx0 = CovVec' * dfx0;
    dmu = dfy0;
    mu=f;  

    %Step 2: Function evaluation at perturbed design
    dist = [-step1 * CovVec, step1 * CovVec];
    for i = 1 : 2*lx
        [fdata(i), dfydata(:,i), dfx1] = fun(y, x+dist(:,i));
        dfxdata(:,i) = CovVec' * dfx1;
    end

    %Step 3: Process function values to compute mean (gradient) and
    %standard deviation
    for i = 1 : lx
        f1 = fdata(lx + i); 
        f2 = fdata(i);
        dfy1 = dfydata(:, lx+i); 
        dfx1 = dfxdata(:, lx+i);
        dfy2 = dfydata(:, i); 
        dfx2 = dfxdata(:, i);
        dfxx(:, i) = (dfx1 - dfx2) / (2 * step1); 
        dfyx(:, i)=(dfy1 - dfy2) / (2 * step1);
        dmu = dmu + 0.5 * CovVal(i) * (dfy1 - 2*dfy0 + dfy2) / (step1^2);
        mu = mu + 0.5 * CovVal(i) * (f1 - 2*f + f2) / (step1^2);
    end %for
    sig = 0.5 * trace((CovVal.*dfxx) * (CovVal.*dfxx)) + dfx0'*(CovVal .* dfx0);
    sig=sqrt(sig);  

    %Compute derivative of standard deviation only if output is requested
    if nargout>2
        %Step 4: Compute and scale second order principal sensitivity directions
        deltax = CovVal .* dfxx .* CovVal';
        epsi = epsHat ./ sqrt(sum(deltax.^2));
        deltax = deltax .* epsi;
        
        %Step 4: Function evaluations at perturbed designs
        dfyx2 = zeros(ly,lx);
        dist = [-step1 * CovVec + CovVec * deltax, step1 * CovVec + CovVec*deltax];
        for i=1:2*lx
            [~, dfydata(:,i)] = fun(y, x+dist(:,i));
        end
        
        %Step 5: Process function values to compute the gradient of the
        %variance
        for i = 1 : lx
            dfy2 = dfydata(:,lx+i); 
            dfy3=dfydata(:,i);
            dfyx2(:,i) = (dfy2 - dfy3) / (2*step1); 
        end
        dsig = sum((dfyx2 - dfyx) ./ epsi, 2);
        dsig = dsig + 2 * dfyx * (CovVal .* dfx0); 
        dsig = dsig / (2 * sig);
    end%nargout
end%SOFM