function do_mixing_model_fitting(num_iterations, model_A_path, model_A_name, model_B_path, model_B_name, alignments_file_base_name, model_type, varargin)

warning('off','MATLAB:structOnObject')
%% Parse arguments
n_varargs = length (varargin);
expect_arg = false;
for i = 1:n_varargs
    arg = varargin{i};
    if ~expect_arg && ischar (arg) && arg (1) == '-'
        switch arg
            case '-prefix',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    prefix = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -prefix missing.');
                end
                
            case '-out_file',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    out_file_name = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -prefix missing.');
                end

            otherwise, error ('do_probabilistic_model_fitting:unknownarg', 'Unknown argument %q.', arg);
        end
    else
        if ~expect_arg
            error ('do_probabilistic_model_fitting:unexparg', 'Unexpected argument at position %d in varargin.', i);
        end
        expect_arg = false;
    end
end

if ~exist('model_type','var')
    disp('Specify model type using -model_type!')
    return
end

if ~exist('n_chunk','var')
    n_chunk = 100;
end

%% load genomic segments
[genV, genD, genJ] = load_Mus_TCR_genomic_segments();



%% figure out alignments files

alignments_file_indices = 1:100;

for f=alignments_file_indices
    fname = [alignments_file_base_name '.alignments.' num2str(f) '.mat'];
    if exist(fname,'file')
        alignments_files{f} = fname;
    else
        break;
    end
end

alignments_file_base_name
alignments_files{:}

basename = alignments_file_base_name;


if ~exist('out_file','var')
    if ~exist('prefix','var')
        p_i = regexp(basename,'\.','split');
        prefix = p_i{1};
        s_i = regexp(prefix,'/','split');
        prefix = s_i{end};
    end
    prefix = [prefix '_MIX_' model_type];
    out_filename = ['results/', prefix, '.mat']
else
    out_filename = out_file;
end


model_A = load(model_A_path, 'model');
model_B = load(model_B_path, 'model');

model_A = model_A.model;
model_B = model_B.model;

model_A.validate_model();
model_B.validate_model();

assert(length(alignments_files) < 5, 'Too many alignment files found'); %Currently running model for small data sets---don't want to run out of memory.

%% Load Alignments

alignments_file = alignments_files{1}
load(alignments_file,'V', 'D', 'J','Sequence');


for f=2:length(alignments_files)
    alignments_file = alignments_files{f}
    current_alignment_file = load(alignments_file,'V', 'D', 'J','Sequence');
    V = [V ; current_alignment_file.V];
    D = [D ; current_alignment_file.D];
    J = [J ; current_alignment_file.J];
    Sequence = [Sequence ; current_alignment_file.Sequence];
end


    do_smoothing = (1 > 0);

    if exist('model_settings','var')
        model = model.settings(model_settings);
    end

    load(alignments_files{1},'alignment_parameters');
    other_model_settings.read_length = alignment_parameters.read_length;
    other_model_settings.max_assignments = 5000;
    other_model_settings.high_error_region = 0;
    other_model_settings.min_V_length = 15;
    force_all_alleles = false;
    probability_threshold_factor = 0.0001;
    no_error = false;
    ignore_deep_error = false;
    READ_LENGTH_CORRECTION = 0;
    

    eval(['model = ' model_type '(genV, genD, genJ);']);
    
    if exist('model_settings','var')
        model = model.settings(model_settings);
    end
        model = model.settings(other_model_settings);


    %% Calculate Pgens from modelA and modelB

    counter_A = model.initialize_counter();
    initial_assigns_A = model.initialize_assigns(counter_A);

    n_assignments_A = zeros(1, length(V));
    loglikelihoods_A = zeros(1, length(V));
    log_generation_probabilities_A = zeros(1, length(V));
    assignment_entropies_A = zeros(1, length(V));
    allele_entropies_A = zeros(1, length(V));
    event_entropies_A = zeros(1, length(V));
    assignment_entropies_no_errors_A = zeros(1, length(V));
    mean_errors_A = zeros(1, length(V));

    MaxLL_A.Vgene = zeros(1, length(V));
    MaxLL_A.Dgene = zeros(1, length(V));
    MaxLL_A.Jgene = zeros(1, length(V));
    MaxLL_A.insVD = zeros(1, length(V));
    MaxLL_A.insDJ = zeros(1, length(V));
    MaxLL_A.cutV = zeros(1, length(V));
    MaxLL_A.cutJ = zeros(1, length(V));
    MaxLL_A.cutDl = zeros(1, length(V));
    MaxLL_A.cutDr = zeros(1, length(V));
    MaxLL_A.V_length = zeros(1, length(V));
    MaxLL_A.J_length = zeros(1, length(V));
    MaxLL_A.D_length = zeros(1, length(V));
    MaxLL_A.errors = zeros(1, length(V));


    
    
    
    
    
    
    counter_B = model.initialize_counter();
    initial_assigns_B = model.initialize_assigns(counter_B);

    n_assignments_B = zeros(1, length(V));
    loglikelihoods_B = zeros(1, length(V));
    log_generation_probabilities_B = zeros(1, length(V));
    assignment_entropies_B = zeros(1, length(V));
    allele_entropies_B = zeros(1, length(V));
    event_entropies_B = zeros(1, length(V));
    assignment_entropies_no_errors_B = zeros(1, length(V));
    mean_errors_B = zeros(1, length(V));

    MaxLL_B.Vgene = zeros(1, length(V));
    MaxLL_B.Dgene = zeros(1, length(V));
    MaxLL_B.Jgene = zeros(1, length(V));
    MaxLL_B.insVD = zeros(1, length(V));
    MaxLL_B.insDJ = zeros(1, length(V));
    MaxLL_B.cutV = zeros(1, length(V));
    MaxLL_B.cutJ = zeros(1, length(V));
    MaxLL_B.cutDl = zeros(1, length(V));
    MaxLL_B.cutDr = zeros(1, length(V));
    MaxLL_B.V_length = zeros(1, length(V));
    MaxLL_B.J_length = zeros(1, length(V));
    MaxLL_B.D_length = zeros(1, length(V));
    MaxLL_B.errors = zeros(1, length(V));
    
    
    
    
    
    

disp('Starting Loop over sequences');
parfor k=1:length(V)
warning('off','MATLAB:structOnObject');
    if mod(k, 5) == 0, disp(['Finished seq #', num2str(k), ' out of ', num2str(length(V))]); end;
    assigns_A = calculate_assignment_probabilities(model_A, Sequence{k}, V(k), D(k), J(k), genV, genD, genJ, initial_assigns_A, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION);
    n_assignments_A(k) = assigns_A.n_assignments;
    loglikelihoods_A(k) = log(assigns_A.likelihood);
    logPgens_A(k) = log(assigns_A.generation_probability);

    if assigns_A.n_assignments > 0
        p_ass = assigns_A.proba(1:assigns_A.n_assignments)./assigns_A.likelihood;
        p_ass_nz = p_ass(p_ass>0);
        assignment_entropies_A(k) = -dot(p_ass_nz,log(p_ass_nz));

        %[~,max_p_ass_i] = max(p_ass);
        max_p_ass_i = assigns_A.max_proba_index;
        MaxLL_A.Vgene(k) = assigns_A.V(max_p_ass_i);
        MaxLL_A.Dgene(k) = assigns_A.DJ(max_p_ass_i,1);
        MaxLL_A.Jgene(k) = assigns_A.DJ(max_p_ass_i,2);
        MaxLL_A.insVD(k) = assigns_A.insVD(max_p_ass_i);
        MaxLL_A.insDJ(k) = assigns_A.insDJ(max_p_ass_i);
        MaxLL_A.cutV(k) = assigns_A.cutV_given_V(max_p_ass_i,1);
        MaxLL_A.cutJ(k) = assigns_A.cutJ_given_J(max_p_ass_i,1);
        MaxLL_A.cutDl(k) = assigns_A.cutDlcutDr_given_D(max_p_ass_i,1);
        MaxLL_A.cutDr(k) = assigns_A.cutDlcutDr_given_D(max_p_ass_i,2);
        MaxLL_A.V_length(k) = assigns_A.V_align_length(max_p_ass_i);
        MaxLL_A.J_length(k) = assigns_A.J_align_length(max_p_ass_i);
        MaxLL_A.D_length(k) = assigns_A.D_align_length(max_p_ass_i);
        MaxLL_A.errors(k) = assigns_A.error(max_p_ass_i);

        mean_errors_A(k) = dot(p_ass,assigns_A.error(1:assigns_A.n_assignments));


        if assigns_A.generation_probability > 0
            p_gen = assigns_A.event_probability(1:assigns_A.n_assignments)./assigns_A.generation_probability;
            p_gen = p_gen(p_gen > 0);
            assignment_entropies_no_errors_A(k) = -dot(p_gen,log(p_gen));
        else
            assignment_entropies_no_errors_A(k) = NaN;
        end

        if assigns_A.likelihood > 0
            [counter_A, delta_counter] = model.update_counter(assigns_A, counter_A);
            allele_entropies_A(k) = entropy(delta_counter.nPVallele_given_gene) - entropy(delta_counter.nPV) + entropy(delta_counter.nPDallele_given_gene) - entropy(sum(delta_counter.nPDJ,2));
        end

        event_entropies_A(k) = assignment_entropies_no_errors_A(k) - allele_entropies_A(k);

    end


    
    
    
    assigns_B = calculate_assignment_probabilities(model_B, Sequence{k}, V(k), D(k), J(k), genV, genD, genJ, initial_assigns_B, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION);
    n_assignments_B(k) = assigns_B.n_assignments;
    loglikelihoods_B(k) = log(assigns_B.likelihood);
    logPgens_B(k) = log(assigns_B.generation_probability);

    if assigns_B.n_assignments > 0
        p_ass = assigns_B.proba(1:assigns_B.n_assignments)./assigns_B.likelihood;
        p_ass_nz = p_ass(p_ass>0);
        assignment_entropies_B(k) = -dot(p_ass_nz,log(p_ass_nz));

        %[~,max_p_Bss_i] = max(p_Bss);
        max_p_ass_i = assigns_B.max_proba_index;
        MaxLL_B.Vgene(k) = assigns_B.V(max_p_ass_i);
        MaxLL_B.Dgene(k) = assigns_B.DJ(max_p_ass_i,1);
        MaxLL_B.Jgene(k) = assigns_B.DJ(max_p_ass_i,2);
        MaxLL_B.insVD(k) = assigns_B.insVD(max_p_ass_i);
        MaxLL_B.insDJ(k) = assigns_B.insDJ(max_p_ass_i);
        MaxLL_B.cutV(k) = assigns_B.cutV_given_V(max_p_ass_i,1);
        MaxLL_B.cutJ(k) = assigns_B.cutJ_given_J(max_p_ass_i,1);
        MaxLL_B.cutDl(k) = assigns_B.cutDlcutDr_given_D(max_p_ass_i,1);
        MaxLL_B.cutDr(k) = assigns_B.cutDlcutDr_given_D(max_p_ass_i,2);
        MaxLL_B.V_length(k) = assigns_B.V_align_length(max_p_ass_i);
        MaxLL_B.J_length(k) = assigns_B.J_align_length(max_p_ass_i);
        MaxLL_B.D_length(k) = assigns_B.D_align_length(max_p_ass_i);
        MaxLL_B.errors(k) = assigns_B.error(max_p_ass_i);

        mean_errors_B(k) = dot(p_ass,assigns_B.error(1:assigns_B.n_assignments));


        if assigns_B.generation_probability > 0
            p_gen = assigns_B.event_probability(1:assigns_B.n_assignments)./assigns_B.generation_probability;
            p_gen = p_gen(p_gen > 0);
            assignment_entropies_no_errors_B(k) = -dot(p_gen,log(p_gen));
        else
            assignment_entropies_no_errors_B(k) = NaN;
        end

        if assigns_B.likelihood > 0
            [counter_B, delta_counter] = model.update_counter(assigns_B, counter_B);
            allele_entropies_B(k) = entropy(delta_counter.nPVallele_given_gene) - entropy(delta_counter.nPV) + entropy(delta_counter.nPDallele_given_gene) - entropy(sum(delta_counter.nPDJ,2));
        end

        event_entropies_B(k) = assignment_entropies_no_errors_B(k) - allele_entropies_B(k);

    end
    
    
end


%% Mixing model iterations

mix_fracts = zeros(num_iterations + 1, 2);

mix_fracts(1, :) = [0.5, 0.5];



mix_counters = zeros(num_iterations+1, 2);

if num_iterations > 1
    for iteration = 1:num_iterations-1

        mix_counter = [0, 0];

        for seq_index = 1:length(V)
            if (loglikelihoods_A(seq_index) > -Inf) && (loglikelihoods_B(seq_index) > -Inf)
                prob_A = exp(loglikelihoods_A(seq_index)) * mix_fracts(iteration, 1); %note that mix_fract(iteration) refers to the mix fract inferred from the previous iteration;
                prob_B = exp(loglikelihoods_B(seq_index)) * mix_fracts(iteration, 2);
                mix_counter = mix_counter + [prob_A / (prob_A + prob_B), prob_B / (prob_A + prob_B)];
            end
        end

        mix_counters(iteration + 1, :) = mix_counter;

        mix_fracts(iteration + 1, :) = mix_counter/sum(mix_counter);

    end
end

iteration = num_iterations; %We want to collect extra info on the last iteration
loglikelihoods = zeros(1, length(V));
log_generation_probabilities = zeros(1, length(V));
MaxLL_orig_pop = zeros(1, length(V)); %1 for model_A, 2 for model_B, 0 for assignment

mix_counter = [0, 0];
    
for seq_index = 1:length(V)
    if (loglikelihoods_A(seq_index) > -Inf) && (loglikelihoods_B(seq_index) > -Inf)
        prob_A = exp(loglikelihoods_A(seq_index)) * mix_fracts(iteration, 1); %note that mix_fract(iteration) refers to the mix fract inferred from the previous iteration;
        prob_B = exp(loglikelihoods_B(seq_index)) * mix_fracts(iteration, 2);
        mix_counter = mix_counter + [prob_A / (prob_A + prob_B), prob_B / (prob_A + prob_B)];
        
        loglikelihoods(seq_index) = log(prob_A + prob_B);
        log_generation_probabilities(seq_index) = log(exp(log_generation_probabilities_A(seq_index)) * mix_fracts(iteration, 1) + exp(log_generation_probabilities_B(seq_index)) * mix_fracts(iteration, 2));
        if prob_A > prob_B
            MaxLL_orig_pop(k) = 1;
        else
            MaxLL_orig_pop(k) = 2;
        end
        
        
    else
        loglikelihoods(seq_index) = -Inf;
        log_generation_probabilities(seq_index) = -Inf;
        MaxLL_orig_pop(seq_index) = 0;
        
    end
end

mix_counters(iteration + 1, :) = mix_counter;

mix_fracts(iteration + 1, :) = mix_counter/sum(mix_counter);

model_A_info.name = model_A_name;
model_A_info.model = model_A;
model_A_info.counter = counter_A;
model_A_info.n_assignments = n_assignments_A;
model_A_info.MaxLL = MaxLL_A;
model_A_info.mean_errors = mean_errors_A;
model_A_info.loglikelihoods = loglikelihoods_A;
model_A_info.log_generation_probilities = log_generation_probabilities_A;
model_A_info.event_entropies = event_entropies_A;
model_A_info.assignment_entropies = assignment_entropies_A;
model_A_info.assignment_entropies_no_errors = assignment_entropies_no_errors_A;
model_A_info.allele_entropies = allele_entropies_A;

model_B_info.name = model_B_name;
model_B_info.model = model_B;
model_B_info.counter = counter_B;
model_B_info.n_assignments = n_assignments_B;
model_B_info.MaxLL = MaxLL_B;
model_B_info.mean_errors = mean_errors_B;
model_B_info.loglikelihoods = loglikelihoods_B;
model_B_info.log_generation_probilities = log_generation_probabilities_B;
model_B_info.event_entropies = event_entropies_B;
model_B_info.assignment_entropies = assignment_entropies_B;
model_B_info.assignment_entropies_no_errors = assignment_entropies_no_errors_B;
model_B_info.allele_entropies = allele_entropies_B;



save(out_filename, 'mix_fracts','loglikelihoods','log_generation_probabilities', 'MaxLL_orig_pop', 'model_A_info', 'model_B_info', 'genV', 'genD', 'genJ', 'Sequence', 'V', 'D', 'J');


end
