function [truemindex, trueT, mcpmy, linearpmy, AIC, AICc, SIC, mcpmy_data] = submitDrawMAPKTopologies(seed)
%% Options
% Global
opts.UseModelSeeds = false;
opts.UseModelInputs = false;
opts.RelTol = 1e-4;
opts.AbsTol = 1e-8;
opts.Normalized = true;
opts.Verbose = 1;

% Fitting
opts.LowerBound = 1e-10;
opts.UpperBound = 1e8;
opts.TolOptim = 1e-6;
opts.UseAdjoint = false;
opts.MaxIter = 2000;
opts.MaxFunEvals = 5000;
opts.Restart = 0;

% Sampling
opts.MaxSampleStep   = 1;     % Uncertainty vectors with 95% CI stetching more than this fold change are truncated
opts.AdaptMaxSampleStep = true; % Change the max sample step based on the acceptance ratio
opts.AdaptProposal   = true;  % Re-linearize the system every StepsPerCheck to get a new proposal distribution
opts.AdaptThinning   = false;  % Change the thinning based on the acceptance ratio
opts.StepsPerCheck   = 100;   % Number of steps between checks on the acceptance rate

% Topology probability
opts.NeedFit = false;
opts.TopologyTol = 0.05;
opts.VariableScheduleTarget = 0.2;
opts.PathMinimumSamples = 20;
opts.PathMaximumSamples = 20;

%% Fetch a bunch of models
mDKDP = LoadModel('Ferrell_MAPK_DKDP.txt');
mDKPP = LoadModel('Ferrell_MAPK_DKPP.txt');
mPKDP = LoadModel('Ferrell_MAPK_PKDP.txt');
mPKPP = LoadModel('Ferrell_MAPK_PKPP.txt');

m = [mDKDP; mDKPP; mPKDP; mPKPP];
nTop = numel(m);

%% Define priors
% Topologies prior is uniform
opts.PriorTopology = zeros(nTop,1) + 1/nTop;

% Parameter prior is a spherical Gaussian
CI = 100; % Fold change to 95% CI
variance = (log(CI)/1.96).^2; % Log space variance

mux = cell(nTop,1);
Vx = cell(nTop,1);
for i = 1:nTop
    % Mean
    mux{i} = zeros(m(i).nk,1) + 1e-1;
    
    % Variance (normalized with no covariance) edit: unnormalized here?
    Vx{i} = diag(mux{i} .* (zeros(m(i).nk,1) + variance) .* mux{i});
end

%% Create prior objectives
objPrior = Gzero(nTop);
for i = 1:nTop
    objPrior(i) = objectiveLogNormalPriorOnKineticParameters(mux{i}, Vx{i});
end

%% Define test experiment
tF = 100;
u = 1;
con = Experiment(m(1), tF, [], false, false, u, [], [], [], 'test');

outputs = [1;2;3];
ntimes = 8;
lintimes = linspace(tF/ntimes, tF, ntimes).'; % evenly spaced time points
sd = sdLinear(0.01,0.1);

%% Set random seed
rng(seed);

%% Draw a random model
% Topology
truemindex = intrnd(1,nTop);

% Parameter set
truenk = m(truemindex).nk;
if opts.Normalized
    vinv = diag(mux{truemindex}.^(-1));
    Vlnv = vinv * Vx{truemindex} * vinv;
    trueT = exp(vec(mvnbndrndgibbs(log(mux{truemindex}), Vlnv, log(zeros(truenk,1) + opts.LowerBound), log(zeros(truenk,1) + opts.UpperBound), [], 1, 50)));
else
    trueT = vec(mvnbndrndgibbs(mux{truemindex}, Vx{truemindex}, zeros(truenk,1) + opts.LowerBound, zeros(truenk,1) + opts.UpperBound));
end

% Create true model
mtrue = Update(m(truemindex), trueT);

%% Generate Fake Data
[obj, data] = generateFakeData(mtrue, con, outputs, lintimes, sd, true, true, [], opts);

%% Fit all models to data
mfit = m;
for i = 1:nTop
    % Reset parameters by setting to the mean
    mfit(i) = Update(mfit(i), zeros(mfit(i).nk,1) + 0.1);
    
    % FitObjective
    mfit(i) = FitObjective(mfit(i), con, [obj; objPrior(i)], opts);
end

%% Computer Monte Carlo p_m|y
opts.Verbose = 2;
opts.TopologyMethod = 'Path';
opts.PathProposalDist = 'Linear';
opts.PathSchedule = 'Variable';
opts.PathEvalMethod = 'Importance';
[mcpmy, mcpmy_data] = TopologyProbability(mfit, con, obj, objPrior, [], [], opts);

%% Compute linear p_m|y
opts.Verbose = 1;
opts.TopologyMethod = 'Linear';
linearpmy = TopologyProbability(mfit, con, obj, objPrior, [], [], opts);

%% Fit all models to data only
for i = 1:nTop
    % FitObjective
    mfit(i) = FitObjective(mfit(i), con, obj, opts);
end

%% Compute AIC, AICc, and SIC
n = numel(lintimes) * numel(outputs);

pyTm = zeros(nTop,1);
AIC = zeros(nTop,1);
AICc = zeros(nTop,1);
SIC = zeros(nTop,1);
for iTop = 1:nTop
    pyTm(iTop) = ObjectiveLogLikelihood(mfit(iTop), con, obj, opts);
    
    nk = mfit(iTop).nk;
    
    AIC(iTop) = 2 * nk - 2 * pyTm(iTop);
    AICc(iTop) = AIC(iTop) + 2 * nk * (nk+1) / (n - nk - 1);
    SIC(iTop) = nk * log(n) - 2 * pyTm(iTop);
end
