function find_insertion_profile_each_seq(model_path, 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 '-alignments_files',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    alignments_files_list = varargin{i + 1};
                    alignments_files = regexp(alignments_files_list,',','split');
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -alignments_files missing.');
                end
            case '-alignments_file_base_name',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    alignments_file_base_name = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -alignments_file_base_name missing.');
                end
            case '-alignments_file_indices',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    alignments_file_indices = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -alignments_file_indices missing.');
                end
            case '-model_type',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    model_type = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -model_type missing.');
                end
            case '-chunk_size',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    n_chunk = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -chunk_size missing.');
                end
            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 '-model_settings',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    model_settings = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -model_settings missing.');
                end
            case '-start_from_flat_prior',
                if i + 1 <= n_varargs
                    expect_arg = true;
                    start_from_flat_prior = varargin{i + 1};
                else
                    error ('do_probabilistic_model_fitting:noarg', 'Required argument for -start_from_flat_prior 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
if ~exist('alignments_files','var')
    if ~exist('alignments_file_indices','var')
        alignments_file_indices = 1:100;
    end

    if ~exist('alignments_file_base_name','var')
        disp('Need either list of comma separated alignments files or a base name and indices!')
        assert(1==2);
    end

    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
end

alignments_file_base_name
alignments_files{:}

if exist('alignments_file_base_name','var')
    basename = alignments_file_base_name;
else
    basename = alignments_files{1};
end

if ~exist('prefix','var')
    p_i = regexp(basename,'\.','split');
    prefix = p_i{1};
    s_i = regexp(prefix,'/','split');
    prefix = s_i{end};
end

%% Iterate to fit probabilistic model

    out_filename = ['results/', prefix, '_insertion_counters.mat'];


    %% Load Initial Model
    model = 0;
    
    disp(['Loading model from ', model_path ])
    load([model_path '/model.mat'])
    % loads model: PinsVD_specific, PinsDJ_specific, PdelV, PdelDl, PdelDr, PdelJ, PV, PD,
    % PJ, Perror
    model.validate_model(); % See if model is valid - whatever that means.
    
    do_smoothing = true;

    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;

    if ~do_smoothing
        other_model_settings.max_assignments = 1;
    else
        other_model_settings.max_assignments = 5000;
    end

    if other_model_settings.read_length==101
        other_model_settings.high_error_region = 20;
        other_model_settings.min_V_length = 30;
    else
        other_model_settings.high_error_region = 0;
        other_model_settings.min_V_length = 15;
    end

    force_all_alleles = false;
    probability_threshold_factor = 0.0001;
    other_model_settings.negative_excess_deletions_max = 0;



    model = model.settings(other_model_settings);
    no_error = false;
    ignore_deep_error = false;
    READ_LENGTH_CORRECTION = 0;

    counter = model.initialize_counter();
    initial_assigns = model.initialize_assigns(counter);

    n_assignments = [];
    loglikelihoods = [];
    log_generation_probabilities = [];
    assignment_entropies = [];
    allele_entropies = [];
    event_entropies = [];
    assignment_entropies_no_errors = [];
    mean_errors = [];
    insVDinsDJ_counter_per_seq = [];

    MaxLL.Vgene = [];
    MaxLL.Dgene = [];
    MaxLL.Jgene = [];
    MaxLL.insVD = [];
    MaxLL.insDJ = [];
    MaxLL.cutV = [];
    MaxLL.cutJ = [];
    MaxLL.cutDl = [];
    MaxLL.cutDr = [];
    MaxLL.V_length = [];
    MaxLL.J_length = [];
    MaxLL.D_length = [];
    MaxLL.errors = [];



    for f=1:length(alignments_files)
        %Load Alignments
        if f>1
            clear('V','D','J','Sequence');
        end
        alignments_file = alignments_files{f}
        load(alignments_file,'V', 'D', 'J','Sequence');

        N = size(V,1)
        N_p = ceil(N/n_chunk);

        dummy_counter = model.initialize_counter();
        counters = repmat(dummy_counter, N_p ,1);
        %all_counters = repmat(dummy_counter,N,1);

        myV = cell(N_p,1);
        myD = cell(N_p,1);
        myJ = cell(N_p,1);
        mySequence = cell(N_p,1);
        for k=1:N_p
            myV{k} = V((k-1)*n_chunk +1 : min(k*n_chunk, N));
            myD{k} = D((k-1)*n_chunk +1 : min(k*n_chunk, N));
            myJ{k} = J((k-1)*n_chunk +1 : min(k*n_chunk, N));
            mySequence{k} = Sequence((k-1)*n_chunk +1 : min(k*n_chunk, N));
        end

        cc = ceil(n_chunk/5);
        %h_assigns = zeros(N,1);
        my_n_assignments = cell(N_p,1);
        my_loglikelihoods = cell(N_p,1);
        my_logPgens = cell(N_p,1);
        my_assignment_entropies = cell(N_p,1);
        my_assignment_entropies_no_errors = cell(N_p,1);
        my_allele_entropies = cell(N_p,1);
        my_event_entropies = cell(N_p,1);
        my_mean_errors = cell(N_p,1);
        my_insVDinsDJ_counter_per_seq = cell(N_p, 1);

        my_Vgene = cell(N_p,1);
        my_Dgene = cell(N_p,1);
        my_Jgene = cell(N_p,1);
        my_insVD = cell(N_p,1);
        my_insDJ = cell(N_p,1);
        my_cutV = cell(N_p,1);
        my_cutJ = cell(N_p,1);
        my_cutDl = cell(N_p,1);
        my_cutDr = cell(N_p,1);
        my_V_length = cell(N_p,1);
        my_J_length = cell(N_p,1);
        my_D_length = cell(N_p,1);
        my_errors = cell(N_p,1);



        disp('Starting Loop over sequences');
        parfor k=1:N_p
            warning('off','MATLAB:structOnObject');
            my_n_assignments{k} = zeros(length(myV{k}),1);
            my_loglikelihoods{k} = zeros(length(myV{k}),1);
            my_logPgens{k} = zeros(length(myV{k}),1);
            my_assignment_entropies{k} = zeros(length(myV{k}),1);
            my_assignment_entropies_no_errors{k} = zeros(length(myV{k}),1);
            my_allele_entropies{k} = zeros(length(myV{k}),1);
            my_event_entropies{k} = zeros(length(myV{k}),1);

            my_mean_errors{k} = zeros(length(myV{k}),1);
            my_insVDinsDJ_counter_per_seq{k} = zeros(length(myV{k}), model.max_insertions + 1, model.max_insertions + 1);

            my_Vgene{k} = zeros(length(myV{k}),1);
            my_Dgene{k} = zeros(length(myV{k}),1);
            my_Jgene{k} = zeros(length(myV{k}),1);
            my_insVD{k} = zeros(length(myV{k}),1);
            my_insDJ{k} = zeros(length(myV{k}),1);
            my_cutV{k} = zeros(length(myV{k}),1);
            my_cutJ{k} = zeros(length(myV{k}),1);
            my_cutDl{k} = zeros(length(myV{k}),1);
            my_cutDr{k} = zeros(length(myV{k}),1);
            my_V_length{k} = zeros(length(myV{k}),1);
            my_J_length{k} = zeros(length(myV{k}),1);
            my_D_length{k} = zeros(length(myV{k}),1);
            my_errors{k} = zeros(length(myV{k}),1);


            for i=1:length(myV{k})
                if mod(i,cc)==0, disp([num2str(i + (k-1)*n_chunk) '*']); end;
                % now explore space of all possible (sufficiently probable)
                % assignments and update the expectation values.
                assigns = calculate_assignment_probabilities(model, mySequence{k}{i}, myV{k}(i),myD{k}(i),myJ{k}(i), genV, genD, genJ, initial_assigns, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION);
                my_n_assignments{k}(i) = assigns.n_assignments;
                my_loglikelihoods{k}(i) = log(assigns.likelihood);
                my_logPgens{k}(i) = log(assigns.generation_probability);

                if assigns.n_assignments > 0
                    p_ass = assigns.proba(1:assigns.n_assignments)./assigns.likelihood;
                    p_ass_nz = p_ass(p_ass>0);
                    my_assignment_entropies{k}(i) = -dot(p_ass_nz,log(p_ass_nz));

                    %[~,max_p_ass_i] = max(p_ass);
                    max_p_ass_i = assigns.max_proba_index;
                    my_Vgene{k}(i) = assigns.V(max_p_ass_i);
                    my_Dgene{k}(i) = assigns.DJ(max_p_ass_i,1);
                    my_Jgene{k}(i) = assigns.DJ(max_p_ass_i,2);
                    my_insVD{k}(i) = assigns.insVD(max_p_ass_i);
                    my_insDJ{k}(i) = assigns.insDJ(max_p_ass_i);
                    my_cutV{k}(i) = assigns.cutV_given_V(max_p_ass_i,1);
                    my_cutJ{k}(i) = assigns.cutJ_given_J(max_p_ass_i,1);
                    my_cutDl{k}(i) = assigns.cutDlcutDr_given_D(max_p_ass_i,1);
                    my_cutDr{k}(i) = assigns.cutDlcutDr_given_D(max_p_ass_i,2);
                    my_V_length{k}(i) = assigns.V_align_length(max_p_ass_i);
                    my_J_length{k}(i) = assigns.J_align_length(max_p_ass_i);
                    my_D_length{k}(i) = assigns.D_align_length(max_p_ass_i);
                    my_errors{k}(i) = assigns.error(max_p_ass_i);
                    

                    my_mean_errors{k}(i) = dot(p_ass,assigns.error(1:assigns.n_assignments));


                    if assigns.generation_probability > 0
                        p_gen = assigns.event_probability(1:assigns.n_assignments)./assigns.generation_probability;
                        p_gen = p_gen(p_gen > 0);
                        my_assignment_entropies_no_errors{k}(i) = -dot(p_gen,log(p_gen));
                    else
                        my_assignment_entropies_no_errors{k}(i) = NaN;
                    end

                    if assigns.likelihood > 0
                        [counters(k), delta_counter] = model.update_counter(assigns, counters(k));
                        my_allele_entropies{k}(i) = entropy(delta_counter.nPVallele_given_gene) - entropy(delta_counter.nPV) + entropy(delta_counter.nPDallele_given_gene) - entropy(sum(delta_counter.nPDJ,2));
                        my_insVDinsDJ_counter_per_seq{k}(i, :, :) = delta_counter.nPinsVDinsDJ;
                    end

                    my_event_entropies{k}(i) = my_assignment_entropies_no_errors{k}(i) - my_allele_entropies{k}(i);

                end


            end

        end
        counter = model.sum_counters([counter; counters]);
        for k=1:N_p
            n_assignments = [n_assignments, my_n_assignments{k}'];
            loglikelihoods = [loglikelihoods, my_loglikelihoods{k}'];
            log_generation_probabilities = [log_generation_probabilities, my_logPgens{k}'];
            assignment_entropies = [assignment_entropies, my_assignment_entropies{k}'];
            assignment_entropies_no_errors = [assignment_entropies_no_errors, my_assignment_entropies_no_errors{k}'];
            allele_entropies = [allele_entropies, my_allele_entropies{k}'];
            event_entropies = [event_entropies, my_event_entropies{k}'];

            mean_errors = [mean_errors, my_mean_errors{k}'];
            insVDinsDJ_counter_per_seq = [insVDinsDJ_counter_per_seq; my_insVDinsDJ_counter_per_seq{k}];

            MaxLL.Vgene = [MaxLL.Vgene, my_Vgene{k}'];
            MaxLL.Dgene = [MaxLL.Dgene, my_Dgene{k}'];
            MaxLL.Jgene = [MaxLL.Jgene, my_Jgene{k}'];
            MaxLL.insVD = [MaxLL.insVD, my_insVD{k}'];
            MaxLL.insDJ = [MaxLL.insDJ, my_insDJ{k}'];
            MaxLL.cutV = [MaxLL.cutV, my_cutV{k}'];
            MaxLL.cutJ = [MaxLL.cutJ, my_cutJ{k}'];
            MaxLL.cutDl = [MaxLL.cutDl, my_cutDl{k}'];
            MaxLL.cutDr = [MaxLL.cutDr, my_cutDr{k}'];
            MaxLL.V_length = [MaxLL.V_length, my_V_length{k}'];
            MaxLL.J_length = [MaxLL.J_length, my_J_length{k}'];
            MaxLL.D_length = [MaxLL.D_length, my_D_length{k}'];
            MaxLL.errors = [MaxLL.errors, my_errors{k}'];

        end
        %hist(h_assigns,30);
    end


    %%% Normalize the expectation values calculated and save to file.
        model = model.get_model_from_counter(counter);
        model = model.calculate_assignment_entropies();

    sequence_repertoire_entropy = model.S_total - mean(event_entropies(n_assignments>0));

    save(out_filename, 'model','counter','n_assignments','MaxLL','mean_errors','loglikelihoods','log_generation_probabilities','event_entropies','assignment_entropies','assignment_entropies_no_errors','allele_entropies','sequence_repertoire_entropy', 'insVDinsDJ_counter_per_seq');


end
