function obj = constructObjectiveChiSquare(m, outputs, times, sd, nonNegMeasurements, measurements, perfect)
%obj = constructObjectiveChiSquare(m, outputs, times, sd, measurements)
%returns X2, chi-square, the generalized least squares statistic

% sd: @(t,yInd,yVal) returns standard error given t, index of output, and
% value of output

% Clean up inputs
if nargin < 7
    perfect = [];
    if nargin < 6
        measurements = [];
        if nargin < 5
            nonNegMeasurements = [];
            if nargin < 4
                sd = [];
                if nargin < 3
                    times = [];
                    if nargin < 2
                        outputs = [];
                        if nargin < 1
                            m = [];
                        end
                    end
                end
            end
        end
    end
end

% Special case: return empty structure array if inputs are numeric
if isnumeric(m)
%    temp = struct('update', [], 'G', [], 'dGdx', [], 'dGdp', [], 'd2Gdx2', [], 'd2Gdp2', [], 'd2Gdpdx', [], 'd2Gdxdp', [], 'n', [], 'discreteTimes', [], 'continuous', [], 'complex', [], 'H', [], 'Hn', [], 'EH', [], 'EHn', [], 'F',  [], 'Fn', [], 'AddData', [], 'AddExpectedData', [], 'p', [], 'PValue', [], 'QuantifyPrediction', []);
    temp = struct('update', [], 'G', [], 'dGdx', [], 'dGdp', [], 'd2Gdx2', [], 'd2Gdp2', [], 'd2Gdpdx', [], 'd2Gdxdp', [], 'n', [], 'discreteTimes', [], 'continuous', [], 'complex', [], 'H', [], 'Hn', [], 'EH', [], 'EHn', [], 'F',  [], 'Fn', [], 'AddData', [], 'AddExpectedData', [], 'PValue', [], 'QuantifyPrediction', []);
    if isempty(m)
        m = 1;
    end
    if ~isscalar(m)
        obj(prod(m),1) = temp;
        obj = reshape(obj, m);
        return
    end
    if isempty(outputs)
        outputs = 1;
    end
    obj(m, outputs) = temp;
    return
end

% Check inputs
assert(isvector(outputs), 'KroneckerBio:constructObjectiveChiSquare:outputs', 'Input "outputs" must be a vector.')
n = length(outputs);
assert(isvector(times) && length(times) == n, 'KroneckerBio:constructObjectiveChiSquare:times', 'Input "times" must be a vector length of "outputs".')
assert((isvector(measurements) && length(measurements) == n) || isempty(measurements), 'KroneckerBio:constructObjectiveChiSquare:measurements', 'Input "measurements" must be a vector length of "outputs" or empty.')

% Gut m
temp = m;
clear m
m.nX = temp.nX;
m.nP = temp.nP;
m.c = temp.c;
clear temp

% Constants
nX = m.nX;

% Find unique times
discreteTimes = vec(unique(times)).';

% Compute covariance matrix
if ~isempty(measurements)
    V = sparse(n,n);
    for j = 1:n
        V(j,j) = sd(times(j), outputs(j), measurements(j))^2;
    end
end

% Default is nonNegMeasurements
if isempty(nonNegMeasurements)
    nonNegMeasurements = true;
end

% Default is real data
if isempty(perfect)
    perfect = false;
end

obj.update = @update;
obj.G = @G;
obj.dGdx = @dGdx;
obj.dGdp = @dGdp;
obj.d2Gdx2 = @d2Gdx2;
obj.d2Gdp2 = @d2Gdp2;
obj.d2Gdpdx = @d2Gdpdx;
obj.d2Gdxdp = @d2Gdxdp;
obj.n = n;
obj.discreteTimes = discreteTimes;
obj.continuous = false;
obj.complex = false;
obj.H = @H;
obj.Hn = @Hn;
obj.EH = @EH;
obj.EHn = @EHn;
obj.F = @F;
obj.Fn = @Fn;
%obj.FnSimbio = @FnSimbio;
obj.AddData = @AddData;
obj.AddExpectedData = @AddExpectedData;
%obj.p = @p;
obj.PValue = @pvalue;
obj.QuantifyPrediction = @QuantifyPrediction;
    function a=PValue
    end

%% Update
    function objNew = update(mNew)
        objNew = constructObjectiveChiSquare(mNew, outputs, times, sd, true, measurements, perfect);
    end

%% AddData
    function objNew = AddData(sol, u)
        nX = size(sol.c,2);
        
        % Evaluate the ODE solver structure
        if isempty(sol.idata)
            % The solution is already discretized
            x = sol.y(1:nX, lookup(discreteTimes, sol.x));
        else
            % The complete solution is provided
            x = deval(sol, discreteTimes, 1:nX); % x_t
        end
        
        % New measurements
        newMeasurements = zeros(n,1);
        for i = 1:n
            newMeasurements(i) = sol.c(outputs(i),:) * x(:,discreteTimes == times(i));
            newMeasurements(i) = newMeasurements(i) + randn*sd(times(i),outputs(i),newMeasurements(i));
        end
        
        if nonNegMeasurements
            newMeasurements(newMeasurements < 0) = 0;
        end
        
        objNew = constructObjectiveChiSquare(m, outputs, times, sd, nonNegMeasurements, newMeasurements, false);
    end

%% AddExpectedData
    function objNew = AddExpectedData(sol, u)
        % Evaluate the ODE solver structure
        if isempty(sol.idata)
            % The solution is already discretized
            x = sol.y(1:nX, lookup(discreteTimes, sol.x));
        else
            % The complete solution is provided
            x = deval(sol, discreteTimes, 1:nX); % x_t
        end
        
        % New measurements
        newMeasurements = zeros(n,1);
        for i = 1:n
            newMeasurements(i) = m.c(outputs(i),:) * x(:,discreteTimes == times(i));
        end
        
        % Build new objective function with perfect data
        objNew = constructObjectiveChiSquare(m, outputs, times, sd, nonNegMeasurements, newMeasurements, true);
    end

%% Probability density function
    function p = p(sol, u)
        % This is a centralized probability density function. That is, when
        % chi2 = n, then p = 1. The probability has been scaled to
        % accomodate this and is necessary to store the probability as a
        % floating point number. The relative probabilities between
        % solutions using the same objective function (or even objective
        % functions with the same number of data points) remains correct.
        % p is a direct consequence of chi-square
        
        chi2 = G([], sol, [], 1);
        %p = (2*pi)^(-n/2) * det(V)^(-1/2) * exp(-1/2*chi2); % old non-centralized
        p = exp(-1/2 * (chi2-n)); % centralized
    end

%% P-Value
    function pval = pvalue(sol, con, obj)
        chi2 = 0;
        ntot = 0;
        
        % Sum chi-square and n across all data
        for i = 1:length(sol)
            chi2 = chi2 + obj(i).G([], sol(i), con(i).u, 1);
            ntot = ntot + obj(i).n;
        end
        
        pval = chi2pvalue(chi2, ntot);
    end

%% QuantifyPrediction
    function result = QuantifyPrediction(sol1, sol2, type)
        % Special case: return empty structure array if inputs are numeric
        if isnumeric(sol1)
            result = zeros(sol1);
            return
        end
        
        % Check type
        if isempty(type)
            type = 'diffscalemaxmean';
        end
        
        % Evaluate the ODE solver structure
        if isempty(sol1.idata)
            % The solution is already discretized
            if all(sol1.x == discreteTimes)
                x1 = sol1.y(1:nX,:); % x_t
                x2 = sol2.y(1:nX,:); % x_t
            else
                error('Need to implement multiple obj functions here')
            end
        else
            % The complete solution is provided
            x1 = deval(sol1, discreteTimes, 1:nX); % x_t
            x2 = deval(sol2, discreteTimes, 1:nX); % x_t
        end
        
        % Convert to outputs
        y1 = zeros(n,1);
        y2 = zeros(n,1);
        for i = 1:n
            y1(i,1) = m.c(outputs(i),:)*x1(:,discreteTimes == times(i));
            y2(i,1) = m.c(outputs(i),:)*x2(:,discreteTimes == times(i));
        end
        
        switch lower(type)
            case 'diffscalemaxmean'
                max1 = max(x1,[],2).^-1;
                max1 = spdiags(max1, 0, nX,nX);
                
                % Scale species by the maximum
                x1 = max1 * x1;
                x2 = max1 * x2;
                
                % Convert to outputs
                y1 = zeros(n,1);
                y2 = zeros(n,1);
                for i = 1:n
                    y1(i,1) = m.c(outputs(i),:)*x1(:,discreteTimes == times(i));
                    y2(i,1) = m.c(outputs(i),:)*x2(:,discreteTimes == times(i));
                end
                
                % Difference
                result = mean(abs(y1 - y2), 1);
                
            case 'diffscalemaxmax'
                max1 = max(x1,[],2).^-1;
                max1 = spdiags(max1, 0, nX,nX);
                
                % Scale species by the maximum
                x1 = max1 * x1;
                x2 = max1 * x2;
                
                % NaN should be zero
                x1(isnan(x1)) = 0;
                x2(isnan(x2)) = 0;
                
                % Convert to outputs
                y1 = zeros(n,1);
                y2 = zeros(n,1);
                for i = 1:n
                    y1(i,1) = m.c(outputs(i),:)*x1(:,discreteTimes == times(i));
                    y2(i,1) = m.c(outputs(i),:)*x2(:,discreteTimes == times(i));
                end
                
                % Difference
                result = max(abs(y1 - y2), [], 1);
                
            case 'diffscalemaxtruemax'
                max1 = max(x1,[],2).^-1;
                max1 = spdiags(max1, 0, nX,nX);
                
                % Scale species by the maximum
                x1 = max1 * x1;
                x2 = max1 * x2;
                
                % NaN should be zero
                x1(isnan(x1)) = 0;
                x2(isnan(x2)) = 0;
                
                % Convert to outputs
                y1 = zeros(n,1);
                y2 = zeros(n,1);
                for i = 1:n
                    y1(i,1) = m.c(outputs(i),:)*x1(:,discreteTimes == times(i));
                    y2(i,1) = m.c(outputs(i),:)*x2(:,discreteTimes == times(i));
                end
                
                % Difference
                result = max(abs(y1 - y2), [], 1);
                
        end % switch
    end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Parameter fitting functions %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function [val stopTimes] = G(t,xSol,u,nSim)
        if ~perfect
            % Evaluate the ODE solver structure
            if isempty(xSol.idata)
                % The solution is already discretized
                if numel(xSol(nSim).x) == numel(discreteTimes)
                    x = xSol(nSim).y(1:nX,:); % x_t
                else
                    ind = lookup(discreteTimes, xSol(nSim).x);
                    x = xSol(nSim).y(1:nX,ind);
                end
            else
                % The complete solution is provided
                x = deval(xSol(nSim), discreteTimes, 1:nX); % x_t
            end
            
            % Get e for every point evaluated
            e = zeros(n,1);
            %temp
            %v = zeros(n,1);
            %/temp
            for i = 1:n
                %temp
                %y = m.c(outputs(i),:)*x(:,discreteTimes == times(i));
                %/temp
                e(i,1) = m.c(outputs(i),:)*x(:,discreteTimes == times(i)) - measurements(i);
                %temp
                %v(i,1) = sd(times(i), outputs(i), y)^2;
                %/temp
            end
            
            % Sort by absolute value
            [unused I] = sort(abs(e.'*sqrt(V)));
            e = e(I);
            Vsort = diag(V); % Extract diagonal
            Vsort = spdiags(Vsort(I),0,n,n); % Replace diagonal with sorted version
            %Proper way:
            %Vsort = spdiags(v(I),0,n,n); % Replace diagonal with sorted version
            
            % Compute chi-square
            val = e'*(Vsort\e);
            
            if nargout >= 2
                stopTimes = discreteTimes;
            end
        else
            % Perfect data requires an expected goal value
            val = EG(t,xSol,u,nSim);
            stopTimes = discreteTimes;
        end
    end

    function val = dGdx(t,xSol,u,nSim)
        %Find all data points that have a time that matches t
        currentIndexes = find(t == times);
        tlength = length(currentIndexes);
        if tlength > 0
            % Evaluate the ODE solver structure
            if isempty(xSol.idata)
                % The solution is already discretized
                % Get the solution at the matching timepoint
                x = xSol(nSim).y(1:nX, xSol.x == t);
            else
                % The complete solution is provided
                x = deval(xSol(nSim), t, 1:nX);
            end
            
            % Compute e for each datapoint that has a matching time to t
            e = zeros(tlength,1);
            C = zeros(tlength,nX);
            %temp
            %v = zeros(tlength,1);
            %/temp
            for i = 1:tlength
                C(i,:) = m.c(outputs(currentIndexes(i)),:);
                %temp
                %y = C(i,:)*x;
                %/temp
                e(i,1) = C(i,:)*x - measurements(currentIndexes(i));
                %temp
                %v(i,1) = sd(times(currentIndexes(i)), outputs(currentIndexes(i)), y)^2;
                %/temp
            end
            
            % Sort by absolute value
            %temp
            %Vsort = diag(v);
            %/temp
            Vsort = V(currentIndexes,currentIndexes);
            [unused I] = sort(abs(e.'*sqrt(Vsort)));
            e = e(I);
            C = C(I,:);
            Vsort = diag(Vsort); % Extract diagonal
            Vsort = spdiags(Vsort(I),0,tlength,tlength); % Replace diagonal with sorted version
            
            % Compute dX2/dx = 2*e'*W*dy/dx = 2*e'*W*C
            val = 2*e'*(Vsort\C);
        else
            val = zeros(1, nX);
        end
    end

    function val = dGdp(t,xSol,u,nSim)
        val = sparse(1, m.nP);
    end

    function val = d2Gdx2(t,xSol,u,nSim)
        %Find all data points that have a time that matches t
        currentIndexes = find(t == times);
        tlength = length(currentIndexes);
        if tlength < 0
            %compute d2x2/dx2 = 2*C'*W*C
            C = zeros(tlength,nX);
            for i = 1:tlength
                C(i,:) = m.c(outputs(currentIndexes(i)),:);
            end
            val = 2*C*(V(currentIndexes,currentIndexes)\C);
        else
            val = zeros(nX, nX);
        end
    end

    function val = d2Gdp2(t,xSol,u,nSim)
        val = sparse(m.nP, m.nP);
    end

    function val = d2Gdxdp(t,xSol,u,nSim)
        val = sparse(m.nP, nX);
    end

    function val = d2Gdpdx(t,xSol,u,nSim)
        val = sparse(nX, m.nP);
    end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Experimental design functions %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%Hessian of the objective function
    function val = H(xSol, dxdpSol, dx2dp2Sol)
        nP = size(dxdpSol.y, 1) / nX;
        
        % Evaluate the ODE solver's xSol structure to get x(time)
        x = deval(xSol, times, 1:nX);
        dxdp = deval(dxdpSol, times, 1:nX*nP);
        d2xdp2 = deval(dx2dp2Sol, times, 1:nX*nP*nP);

        % Get values for every point
        e = zeros(n,1);
        dydp = zeros(n, nP);
        d2ydp2 = zeros(n, nP*nP);
        for i = 1:n
            iC = m.c(outputs(i),:);
            % Compute e
            e(i) = iC*x(:,i) - measurements(i);
            
            % Compute dy/dp
            dydp(i,:) = iC * reshape(dxdp(:,i), nX, nP);
            
            % Compute d2y/dp2
            d2ydp2(i,:) = iC * reshape(d2xdp2(:,i), nX, nP*nP);
        end
        
        % Compute hessian
        val =  e.' * (V \ d2ydp2);
        val = reshape(val, nP, nP);
        val = val + dydp.' * (V \ dydp);
        val = 2 * symmat(val);
    end

%% Normalized hessian of the objective function
    function val = Hn(xSol, dxdpSol, dx2dp2Sol, p)
        nP = size(dxdpSol.y, 1) / nX;
        
        % Evaluate the ODE solver's xSol structure to get x(time)
        x = deval(xSol, times, 1:nX);
        dxdp = deval(dxdpSol, times, 1:nX*nP);
        d2xdp2 = deval(dx2dp2Sol, times, 1:nX*nP*nP);
        
        % Get values for every point
        e = zeros(n,1);
        dydp = zeros(n, nP); % y_p
        d2ydp2 = zeros(n, nP*nP); % y_pp
        for i = 1:n
            iC = m.c(outputs(i),:);
            % Compute e
            e(i) = iC*x(:,i) - measurements(i);
            
            % Compute dy/dp
            dydp(i,:) = iC * reshape(dxdp(:,i), nX, nP);
            
            % Compute d2y/dp2
            d2ydp2(i,:) = iC * reshape(d2xdp2(:,i), nX, nP*nP); % auto-reshape invoked
        end
        d2ydp2 = reshape(d2ydp2, n*nP, nP); % yp_p
        
        % Construct normalization matrices
        dpdlnp = spdiags(p,0,sparse(nP,nP)); % p along the diagonal
        linInd = (0:(nP-1)).' * (nP*nP+nP+1) + 1; % linear indexes for putting p into d2pdlnp2
        d2pdlnp2 = sparse(nP, nP*nP);
        d2pdlnp2(linInd) = p; % p along the 3D diagonal, or 2D along the linear indexes
        
        % Compute normalized hessian
        dydlnp = dydp*dpdlnp; % n_p
        val = dpdlnp*dpdlnp; % p_p
        val = d2ydp2*val; % np_p
        val = val + reshape( dydp*d2pdlnp2, n*nP,nP); % np_p
        val = reshape(val, n,nP*nP); % n_pp
        val = e.' * (V \ val); % _pp
        val = reshape(val, nP,nP); % p_p
        val = val + dydlnp.' * (V \ dydlnp); % p_p
        val = 2 * symmat(val); % p_p
    end

%% Expected hessian of the objective function
    function val = EH(dxdpSol)
        nP = (size(dxdpSol.y, 1) - nX) / nX;
        
        % Evaluate the ODE solver structure
        joint = deval(dxdpSol, times, 1:(nX*(nP+1))); % j_t
        x = joint(1:nX,:); % x_t
        dxdp = joint((nX+1):(nX*(nP+1)),:); % xp_t
        
        % Get dydp for every point
        y = zeros(n,1);
        dydp = zeros(n, nP);
        EV = sparse(n,n);
        for i = 1:n
            iC = m.c(outputs(i),:);
            % Compute y
            y(i) = iC*x(:,i);
            
            % Compute dy/dp
            dydp(i,:) = iC * reshape(dxdp(:,i), nX, nP);
            
            % Compute expected V matrix
            EV(i,i) = sd(times(i), outputs(i), y(i))^2;
        end

        % Compute expected normalized hessian
        val = dydp.' * (EV \ dydp);
        val = 2 * symmat(val);
    end

%% Expected normalized hessian of the objective function
    function val = EHn(dxdpSol, p)
        nP = (size(dxdpSol.y, 1) - nX) / nX;
        
        % Evaluate the ODE solver structure
        joint = deval(dxdpSol, times, 1:(nX*(nP+1))); % j_t
        x = joint(1:nX,:); % x_t
        dxdp = joint((nX+1):(nX*(nP+1)),:); % xp_t
        
        % Get dydp for every point
        y = zeros(n,1);
        dydp = zeros(n, nP);
        EV = sparse(n,n);
        for i = 1:n
            iC = m.c(outputs(i),:);
            % Compute y
            y(i) = iC*x(:,i);
            
            % Compute dy/dp
            dydp(i,:) = iC * reshape(dxdp(:,i), nX, nP);
            
            % Compute expected V matrix
            EV(i,i) = sd(times(i), outputs(i), y(i))^2;
        end

        % Construct normalization matrix
        val = spdiags(p,0,sparse(nP,nP)); % p along the diagonal
        
        % Compute expected normalized hessian
        val = dydp*val;
        val = val.' * (EV \ val);
        val = 2 * symmat(val);
    end

%% Fisher information matrix
    function val = F(dxdpSol)
        nP = (size(dxdpSol.y, 1) - nX) / nX;
        
        % Evaluate the ODE solver structure
        if ~isfield(dxdpSol, 'idata')
            % The solution is already discretized
            if all(dxdpSol.x == discreteTimes)
                % Matching solution
                x = dxdpSol.y(1:nX,:);
                dxdp = dxdpSol.y((nX+1):(nX*(nP+1)),:); %xp_t
            else
                error('Need to implement multiple obj functions here')
            end
        else
            % The complete solution is provided
            joint = deval(dxdpSol, discreteTimes, 1:(nX*(nP+1))); % j_t
            x = joint(1:nX,:); % x_t
            dxdp = joint((nX+1):(nX*(nP+1)),:); % xp_t
        end
        
        % Get dydp for every point
        y = zeros(n,1);
        dydp = zeros(n, nP);
        EV = sparse(n,n);
        for i = 1:n
            ind = find(discreteTimes == times(i), 1, 'first');
            ix = x(:,ind);
            idxdp = reshape(dxdp(:,ind), nX,nP); %x_p
            iC = m.c(outputs(i),:);
            
            % Compute y
            y(i) = iC*ix;
            
            % Compute dy/dp
            dydp(i,:) = iC * idxdp;
            
            % Compute expected V matrix
            EV(i,i) = sd(times(i), outputs(i), y(i))^2;
        end

        % Compute expected normalized hessian
        val = dydp.' * (EV \ dydp);
        val = symmat(val);
    end

%% Fisher information matrix
    function val = Fn(dxdpSol, p)
        nP = (size(dxdpSol.y, 1) - nX) / nX;
        
        % Evaluate the ODE solver structure
        if ~isfield(dxdpSol, 'idata')
            % The solution is already discretized
            if all(dxdpSol.x == discreteTimes)
                % Matching solution
                x = dxdpSol.y(1:nX,:);
                dxdp = dxdpSol.y((nX+1):(nX*(nP+1)),:); %xp_t
            else
                error('Need to implement multiple obj functions here')
            end
        else
            % The complete solution is provided
            joint = deval(dxdpSol, discreteTimes, 1:(nX*(nP+1))); % j_t
            x = joint(1:nX,:); % x_t
            dxdp = joint((nX+1):(nX*(nP+1)),:); % xp_t
        end
        
        % Get dydp for every point
        y = zeros(n,1);
        dydp = zeros(n, nP);
        EV = sparse(n,n);
        for i = 1:n
            ind = find(discreteTimes == times(i), 1, 'first');
            ix = x(:,ind);
            idxdp = reshape(dxdp(:,ind), nX,nP); %x_p
            iC = m.c(outputs(i),:);
            
            % Compute y
            y(i) = iC*ix;
            
            % Compute dy/dp
            dydp(i,:) = iC * idxdp;
            
            % Compute expected V matrix
            EV(i,i) = sd(times(i), outputs(i), y(i))^2;
        end

        % Construct normalization matrix
        val = spdiags(p,0,sparse(nP,nP)); % p along the diagonal
        
        % Compute expected normalized hessian
        val = dydp*val;
        val = val.' * (EV \ val);
        val = symmat(val);
    end

%% Fisher information matrix - SimBiology version
%     function val = FnSimbio(simData, p)
%         % simData is stack of species and sensitivities
%         nP = (size(simData.Data, 2) - nX) / nX;
%         
%         % Use linear interpolation to evaluate the SimData
%         x = piecewiselinear(times, simData.Time, simData.Data(:,1:nX)).';
%         dxdp = piecewiselinear(times, simData.Time, simData.Data(:,(nX+1):end)).';
%         
%         % Get dydp for every point
%         y = zeros(n,1);
%         dydp = zeros(n, nP);
%         EV = sparse(n,n);
%         for i = 1:n
%             iC = m.c(outputs(i),:);
%             % Compute y
%             y(i) = iC*x(:,i);
%             
%             % Compute dy/dp
%             dydp(i,:) = iC * reshape(dxdp(:,i), nX, nP);
%             
%             % Compute expected V matrix
%             EV(i,i) = sd(times(i), outputs(i), y(i))^2;
%         end
% 
%         % Construct normalization matrix
%         val = spdiags(p,0,sparse(nP,nP)); % p along the diagonal
%         
%         % Compute expected normalized hessian
%         val = dydp*val;
%         val = val.' * (EV \ val);
%         val = symmat(val);
%     end

%% Expected goal function
% The expected value of normalized chi-square given the current
% measurements as the bias. If measurements is unset, there is assumed to
% be no bias.
    function val = EG(t,xSol,u,nSim)
        % Evaluate the ODE solver structure
        if isempty(xSol.idata)
            % The solution is already discretized
            if all(xSol(nSim).x == discreteTimes)
                x = xSol(nSim).y(1:nX,:); % x_t
            else
                error('Need to implement multiple obj functions here')
            end
        else
            % The complete solution is provided
            x = deval(xSol(nSim), discreteTimes, 1:nX); % x_t
        end
        
        % Get bias d for every point evaluated
        d = zeros(n,1);
        if ~isempty(measurements)
            for i = 1:n
                d(i,1) = m.c(outputs(i),:)*x(:,discreteTimes == times(i)) - measurements(i);
            end
        end
        
        % compute expected chi-square
        val = n + d.' * (V \ d);
    end

end