function obj = constructObjectiveParameterKineticNormal(m, kbar, Vkbar, UseParams, normalized)

%TODO: implement PValue

% Clean up inputs
if nargin < 5
    normalized = [];
    if nargin < 4
        UseParams = [];
        if nargin < 3
            V = [];
            if nargin < 2
                kbar = [];
                if nargin < 1
                    m = [];
                end
            end
        end
    end
end

% Special case: return empty structure array if inputs are numeric
if isnumeric(m)
    temp = struct('update', [], 'discreteTimes', [], 'continuous', [], 'complex', [], 'G', [], 'dGdp', [], 'd2Gdp2', [], 'F', [], 'Fn', [], 'p', []);
    if isempty(m)
        m = 1;
    end
    if ~isscalar(m)
        obj(prod(m),1) = temp;
        obj = reshape(obj, m);
        return
    end
    if isempty(kbar)
        kbar = 1;
    end
    obj(m, kbar) = temp;
    return
end

% Check inputs
n = numel(kbar);
assert(all(size(Vkbar) == [n,n]), 'KroneckerBio:constructObjectiveParameterNormal:Vsize', 'Input "Vkbar" must be a square matrix of numel(kbar)')

% Find unique times
discreteTimes = 0;

% Constants
nx = m.nX;
nk = m.nP;

% Fix UseParams
[UseParams, nTk] = fixUseParams(UseParams, nk);

% Pre-select kinetic parameters
k = m.p;
Tk = k(UseParams);
Tkbar = kbar(UseParams);
VTbark = Vkbar(UseParams,UseParams);

% Normalization
if normalized
    Tk = log(Tk);
    VTbark = diag(Tkbar.^(-1)) * VTbark * diag(Tkbar.^(-1));
    Tkbar = log(Tkbar);
end

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

% Objective structure
obj.update = @update;

obj.discreteTimes = discreteTimes;
obj.continuous    = false;
obj.complex       = false;

obj.G      = @G;
obj.dGdp   = @dGdk;
obj.d2Gdp2 = @d2Gdk2;

obj.F  = @F;
obj.Fn = @Fn;
obj.p  = @p;

%% Update
    function objNew = update(mNew)
        objNew = constructObjectiveParameterKineticNormal(mNew, kbar, Vkbar, UseParams, normalized);
    end

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

    function [val stopTimes] = G(t,sol,u,nSim)
        stopTimes = discreteTimes;
        diff = Tk - Tkbar;
        val = diff.' * (VTbark \ diff);
    end

    function val = dGdk(t,sol,u,nSim)
        val = zeros(nk,1);
        if normalized
            val(UseParams) = 2 * diag(k(UseParams).^(-1)) * (VTbark \ (Tk - Tkbar));
        else
            val(UseParams) = 2 * (VTbark \ (Tk - kbar));
        end
    end

    function val = d2Gdk2(t,sol,u,nSim)
        val = zeros(nk,nk);
        if normalized
            val(UseParams,UseParams) = 2 * diag(k(UseParams).^(-1)) * (VTbark \ diag(k(UseParams).^(-1))) - 2 * diag(k(UseParams).^(-2)) * diag(VTbark \ (Tk - kbar));
        else
            val(UseParams,UseParams) = 2 * inv(VTbark);
        end
    end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Information theory %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function p = p(sol, u)
        diff = Tk - Tkbar;
        p = (2*pi).^(-nTk/2) * det(VTbark).^(-1/2) * exp(-1/2 * diff.' * (VTbark \ diff));
    end

    function val = F(dxdTSol)
        % This objective only provides information on the first nTk parameters
        nT = (size(dxdTSol.y, 1) - nx) / nx;
        val = zeros(nT,nT);
        
        if ~normalized
            val(1:nTk,1:nTk) = inv(VTbark);
        else
            % Unnormalize
            TkbarTemp = exp(Tkbar);
            val(1:nTk,1:nTk) = diag(TkbarTemp) * inv(VTbark) * diag(TkbarTemp);
        end
    end

    function val = Fn(dxdTSol, T)
        % This objective only provides information on the first nTk parameters
        nT = (size(dxdTSol.y, 1) - nx) / nx;
        val = zeros(nT,nT);
        
        if ~normalized
            % Normalize
            val(1:nTk,1:nTk) = diag(Tkbar.^(-1)) * inv(VTbark) * diag(Tkbar.^(-1));
        else
            val(1:nTk,1:nTk) = inv(VTbark);
        end
    end
end