classdef Mus_V_DJ_cuts_insertions_dinuc_ntbias_model < VDJ_base_model
    % Model class that inherits from base model
    % Contains independent distributions for insertions,
    % joint distribution for palindromes and deletions,
    % and joint distribution for V,D,J genes.

    properties (SetAccess = public)
        %% Model settings
        max_assignments; % Max no. of assignments to explore during Expectation step of EM algorithm.
        max_insertions; % Maximum no. of insertions
        max_V_deletions; % Maximum no. of deletions
        max_D_deletions; % Maximum no. of deletions
        max_J_deletions; % Maximum no. of deletions
        max_excess_V_deletions; % Maximum excess deletions to consider
        max_excess_D_deletions; % Maximum excess deletions to consider
        max_excess_J_deletions; % Maximum excess deletions to consider

        max_palindrome; % Maximum half-palindrome length
        min_V_cut;
        max_V_cut;
        min_D_cut;
        max_D_cut;
        min_J_cut;
        max_J_cut;
        

        % Various entropies of the features from the model
        S_total;
        S_gene;
        S_V;
        S_DJ;
        S_D;
        S_J;
        S_insVD;
        S_insDJ;
        S_insVD_length;
        S_insDJ_length;
        S_insVD_nt;
        S_insDJ_nt;
        S_delV;
        S_delD;
        S_delJ;
        S_delDJ;



        negative_excess_deletions_max; % Maximum 'negative' deletions to consider
        % When considering more deletions than implied by alignment,
        % what are the minimum V,D and J lengths to maintain
        min_J_align_length;
        min_J_assign_length;
        min_V_length;
        % Region on left of sequence that is discarded
        % because of high position dependent error rate.
        high_error_region;
        % Whether sequences with too short D alignments are even
        % considered.
        use_no_D_match_seqs;
        % Read length in data set. Usually set by main model fitting
        % script.
        read_length;

        %% Model parameters
        PinsVD; % P(insertions)
        PinsDJ;
        RnucleotideVD_per_nucleotideVD_5prime;
        RnucleotideDJ_per_nucleotideDJ_3prime;

        PcutV_given_V;
        PcutJ_given_J;
        PcutDlcutDr_given_D;

        PV;
        PDJ; % Joint P( V, D, J gene choices)
        PVallele_given_gene % Probabilities of alleles given gene, for each gene
        PDallele_given_gene
        Rerror_per_sequenced_nucleotide % Error rate

    end
    %%
    methods
        function model = Mus_V_DJ_cuts_insertions_dinuc_ntbias_model(genV, genD, genJ)

            % These are default values for these settings. They can be
            % passed in at run time by calling the function settings()
            model.max_assignments = 5000;
            model.max_insertions = 30;

            model.max_V_deletions = 16;
            model.max_D_deletions = 16;
            model.max_J_deletions = 22;
            model.max_excess_V_deletions = 12;
            model.max_excess_J_deletions = 10;
            model.max_excess_D_deletions = 15;

            model.max_palindrome = 6;
            model.min_J_align_length = 2;
            model.min_J_assign_length = 1;

            model.min_V_length = 15;

            model.negative_excess_deletions_max = 0;
            model.high_error_region = 15;
            model.use_no_D_match_seqs = 1;
            model.read_length = 87;

            model.min_V_cut = -model.max_palindrome;
            model.max_V_cut = model.max_V_deletions;

            model.min_D_cut = -model.max_palindrome;
            model.max_D_cut = model.max_D_deletions;

            model.min_J_cut = -model.max_palindrome;
            model.max_J_cut = model.max_J_deletions;


            % Every variable value is indexed from 1.
            model.PinsVD = ones(model.max_insertions+1,1);
            model.PinsDJ = ones(model.max_insertions+1,1);

            model.PV=ones(max([genV.gene_index]),1);
            model.PDJ = ones(max([genD.gene_index]), max([genJ.gene_index])); % flat distrib on all V genomics segments

            model.PcutV_given_V = ones( model.max_V_cut - model.min_V_cut + 1, max([genV.gene_index]));
            model.PcutJ_given_J = ones( model.max_J_cut - model.min_J_cut + 1, max([genJ.gene_index]));
            model.PcutDlcutDr_given_D = ones( model.max_D_cut - model.min_D_cut + 1, model.max_D_cut - model.min_D_cut + 1, max([genD.gene_index]));
            for d=1:max([genD.gene_index])
                lseq = length(genD(d).seq);
                for cutDl=0:model.max_D_cut
                    for cutDr=0:model.max_D_cut;
                        if cutDl + cutDr > lseq
                            model.PcutDlcutDr_given_D( 1 - model.min_D_cut + cutDl, 1 - model.min_D_cut + cutDr , d) = 0;
                        end
                    end
                end
            end

            model.RnucleotideVD_per_nucleotideVD_5prime = ones(4,4)./4;
            model.RnucleotideDJ_per_nucleotideDJ_3prime = ones(4,4)./4;


            % Initialize probabilities of alleles given genes, to 1/n_alleles for each gene.
            V_genes = max([genV.gene_index]);
            max_V_alleles = max([genV.n_alleles])+1;
            model.PVallele_given_gene = zeros(max_V_alleles , V_genes);
            for a=1:length(genV)
                model.PVallele_given_gene(genV(a).allele,genV(a).gene_index) = 1/genV(a).n_alleles;
            end

            D_genes = max([genD.gene_index]);
            max_D_alleles = max([genD.n_alleles])+1;
            model.PDallele_given_gene = zeros(max_D_alleles , D_genes);
            for a=1:length(genD)
                model.PDallele_given_gene(genD(a).allele,genD(a).gene_index) = 1/genD(a).n_alleles;
            end

            model.Rerror_per_sequenced_nucleotide = 10^(-7); % error rate
            % Normalizes all distributions
            model = normalize_model(model);
            model = calculate_assignment_entropies(model);
        end

        % Fuction that explores different assignments of hidden variables
        function assigns = calculate_assignment_probabilities(model, seq, V, D, J, genV, genD, genJ, initial_assigns, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION)
            % if do_smoothing is false, then you get only 1 assignment per sequence.
            assigns = Mus_V_DJ_cuts_insertions_dinuc_ntbias_model_assignments(struct(model), seq, V, D, J, genV, genD, genJ, initial_assigns, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION);
            %disp(['Total assignments : ' num2str(in-1) ])
        end

        function model = calculate_assignment_entropies(model)
            model.S_V = entropy(model.PV);
            model.S_DJ = entropy(model.PDJ);
            model.S_J = entropy(sum(model.PDJ,1));
            model.S_D = model.S_DJ - model.S_J;

            model.S_gene = model.S_V + model.S_DJ;
            model.S_delV = entropy(model.PcutV_given_V*diag(model.PV)) - entropy(model.PV);

            PcutJJD = bsxfun(@times, model.PcutJ_given_J, shiftdim(model.PDJ',-1));
            PcutDlcutDrcutJJD = bsxfun(@times, model.PcutDlcutDr_given_D , shiftdim(shiftdim(PcutJJD,2),-2));

            model.S_delDJ = entropy(PcutDlcutDrcutJJD) - entropy(model.PDJ);

            model.S_delJ = entropy(PcutJJD) - entropy(model.PDJ);
            model.S_delD = model.S_delDJ - model.S_delJ;

            max_insVD = length(model.PinsVD)-1;
            max_insDJ = length(model.PinsDJ)-1;

            [p_monont_VD,E] =  eig(model.RnucleotideVD_per_nucleotideVD_5prime);
            eig_one_i = find( abs(diag(E)-1)<0.00001);
            p_monont_VD = p_monont_VD(:,eig_one_i)./sum(p_monont_VD(:,eig_one_i));

            [p_monont_DJ,~] =  eig(model.RnucleotideDJ_per_nucleotideDJ_3prime);
            eig_one_i = find( abs(diag(E)-1)<0.00001);
            p_monont_DJ = p_monont_DJ(:,eig_one_i)./sum(p_monont_DJ(:,eig_one_i));

            S_monont_VD = entropy(p_monont_VD);
            S_monont_DJ = entropy(p_monont_DJ);

            S_conditional_nt_VD = entropy(model.RnucleotideVD_per_nucleotideVD_5prime*diag(p_monont_VD)) - S_monont_VD;
            S_conditional_nt_DJ = entropy(model.RnucleotideDJ_per_nucleotideDJ_3prime*diag(p_monont_DJ)) - S_monont_DJ;

            S_nt_VD_vs_Lins_VD = zeros(1 + max_insVD,1);
            S_nt_VD_vs_Lins_VD(1) = 0;
            S_nt_VD_vs_Lins_VD(2:1+max_insVD) = S_monont_VD + [0:(max_insVD-1)]*S_conditional_nt_VD;

            S_nt_DJ_vs_Lins_DJ = zeros(1 + max_insDJ,1);
            S_nt_DJ_vs_Lins_DJ(1) = 0;
            S_nt_DJ_vs_Lins_DJ(2:1+max_insDJ) = S_monont_DJ + [0:(max_insDJ-1)]*S_conditional_nt_DJ;

            model.S_insVD_length = entropy(model.PinsVD);
            model.S_insDJ_length = entropy(model.PinsDJ);

            model.S_insVD_nt = dot(model.PinsVD, S_nt_VD_vs_Lins_VD);
            model.S_insDJ_nt = dot(model.PinsDJ, S_nt_DJ_vs_Lins_DJ);

            model.S_insVD = model.S_insVD_length + model.S_insVD_nt;
            model.S_insDJ = model.S_insDJ_length + model.S_insDJ_nt;

            model.S_total = model.S_gene + model.S_insVD + model.S_insDJ + model.S_delV + model.S_delDJ;

        end


        % counter adds up counts for the distributions. This function
        % initializes it.
        function counter = initialize_counter(model)

            % Call initialize function defined in base model
            counter = initialize_counter@VDJ_base_model(model);

            % Add any other fields you might want to keep track of, but are
            % not in model, below. Begin the name with an 'nP'
            % or 'nM', if it's a distribution or mean value, respectively.

            
            counter.nPv_alignment_index = zeros(size(model.PV,1),1);
            counter.nPj_alignment_index = zeros(size(model.PDJ,2),1);
            counter.nPd_alignment_index = zeros(201,size(model.PDJ,1));
            
            counter.nMmononucleotideVD = zeros(4,1);
            counter.nMmononucleotideDJ = zeros(4,1);
            counter.nMinsertionVD = zeros(4,1);
            counter.nMinsertionDJ = zeros(4,1);


            counter.nPVD_left_edge_dinucleotide = zeros(4,4);
            counter.nPVD_right_edge_dinucleotide = zeros(4,4);
            counter.nPDJ_left_edge_dinucleotide = zeros(4,4);
            counter.nPDJ_right_edge_dinucleotide = zeros(4,4);



            counter.nMtrinucleotideVD = zeros(4,4,4);
            counter.nMtrinucleotideDJ = zeros(4,4,4);

            counter.nPpVmax_delV_V =  zeros(model.max_palindrome + 1, model.max_V_deletions + 1, size(model.PV,1));
            counter.nPpJmax_delJ_J = zeros(model.max_palindrome + 1, model.max_J_deletions + 1, size(model.PDJ,2));

            counter.nPpDlmax_delDl_D = zeros(model.max_palindrome + 1, model.max_D_deletions + 1, size(model.PDJ,1));
            counter.nPpDrmax_delDr_D = zeros(model.max_palindrome + 1, model.max_D_deletions + 1, size(model.PDJ,1));

            counter.nMzeroD = 0;
            counter.nPVDJ = zeros(size(model.PV,1), size(model.PDJ,1), size(model.PDJ,2));
            counter.nPpVdelV = zeros( model.max_palindrome + 1, model.max_V_deletions + 1);
            counter.nPpDldelDl = zeros( model.max_palindrome + 1, model.max_D_deletions + 1);
            counter.nPpDrdelDr = zeros( model.max_palindrome + 1, model.max_D_deletions + 1);
            counter.nPpJdelJ = zeros( model.max_palindrome + 1, model.max_J_deletions + 1);


            counter.nPVcutV = zeros(size(model.PV,1),model.max_palindrome + model.max_V_deletions + 1);
            counter.nPDcutDl = zeros(size(model.PDJ,1),model.max_palindrome + model.max_D_deletions + 1);
            counter.nPDcutDr = zeros(size(model.PDJ,1),model.max_palindrome + model.max_D_deletions + 1);
            counter.nPJcutJ = zeros(size(model.PDJ,2),model.max_palindrome + model.max_J_deletions + 1);

            counter.nPDcutV = zeros(size(model.PDJ,1),model.max_palindrome + model.max_V_deletions + 1);
            counter.nPDcutJ = zeros(size(model.PDJ,1),model.max_palindrome + model.max_J_deletions + 1);

            counter.nPVcutJ = zeros(size(model.PV,1),model.max_palindrome + model.max_J_deletions + 1);
            counter.nPVcutDl = zeros(size(model.PV,1),model.max_palindrome + model.max_D_deletions + 1);
            counter.nPVcutDr = zeros(size(model.PV,1),model.max_palindrome + model.max_D_deletions + 1);

            counter.nPJcutDl = zeros(size(model.PDJ,2),model.max_palindrome + model.max_D_deletions + 1);
            counter.nPJcutDr = zeros(size(model.PDJ,2),model.max_palindrome + model.max_D_deletions + 1);
            counter.nPJcutV = zeros(size(model.PDJ,2),model.max_palindrome + model.max_V_deletions + 1);

            counter.nPinsVDcutV = zeros(model.max_insertions + 1,model.max_palindrome + model.max_V_deletions + 1);
            counter.nPinsDJcutV = zeros(model.max_insertions + 1,model.max_palindrome + model.max_V_deletions + 1);

            counter.nPinsVDcutDl = zeros(model.max_insertions + 1,model.max_palindrome + model.max_D_deletions + 1);
            counter.nPinsDJcutDl = zeros(model.max_insertions + 1,model.max_palindrome + model.max_D_deletions + 1);

            counter.nPinsVDcutDr = zeros(model.max_insertions + 1,model.max_palindrome + model.max_D_deletions + 1);
            counter.nPinsDJcutDr = zeros(model.max_insertions + 1,model.max_palindrome + model.max_D_deletions + 1);

            counter.nPinsVDcutJ = zeros(model.max_insertions + 1,model.max_palindrome + model.max_J_deletions + 1);
            counter.nPinsDJcutJ = zeros(model.max_insertions + 1,model.max_palindrome + model.max_J_deletions + 1);


            max_V_length = model.read_length;
            max_D_length = 16;
            max_J_length = 22;

            counter.nPV_align_length = zeros(max_V_length + 1,1);
            counter.nPD_align_length = zeros(max_D_length + 1,1);
            counter.nPJ_align_length = zeros(max_J_length + 1,1);

            counter.nMVV_err_pos = zeros(size(model.PV,1), max_V_length + model.max_V_deletions);
            counter.nMJJ_err_pos = zeros(size(model.PDJ,2), max_J_length);

            counter.nPJJ_align_length = zeros(size(model.PDJ,2), 1 + max_J_length);
            counter.nPVV_align_length = zeros(size(model.PV,1), 1 + max_V_length);


            counter.nPinsDJ_D_align_length = zeros(model.max_insertions +1, max_D_length + 1);
            counter.nPinsVD_D_align_length = zeros(model.max_insertions +1, max_D_length + 1);

            counter.nPinsDJ_J_align_length = zeros(model.max_insertions +1, max_J_length + 1);
            counter.nPinsVD_J_align_length = zeros(model.max_insertions +1, max_J_length + 1);

            counter.nPinsDJ_V_align_length = zeros(model.max_insertions +1, max_V_length + 1);
            counter.nPinsVD_V_align_length = zeros(model.max_insertions +1, max_V_length + 1);

            counter.nPDallele_D_align_length = zeros(3, max_D_length + 1);


            counter.nPdelVinsVD = zeros(model.max_V_deletions +1, model.max_insertions +1);
            counter.nPdelVinsDJ = zeros(model.max_V_deletions +1, model.max_insertions +1);
            counter.nPdelVdelDl = zeros(model.max_V_deletions +1, model.max_D_deletions +1);
            counter.nPdelVdelDr = zeros(model.max_V_deletions +1, model.max_D_deletions +1);
            counter.nPdelVdelJ = zeros(model.max_V_deletions +1, model.max_J_deletions +1);

            counter.nPdelJinsVD = zeros(model.max_J_deletions +1, model.max_insertions +1);
            counter.nPdelJinsDJ = zeros(model.max_J_deletions +1, model.max_insertions +1);
            counter.nPdelJdelDl = zeros(model.max_J_deletions +1, model.max_D_deletions +1);
            counter.nPdelJdelDr = zeros(model.max_J_deletions +1, model.max_D_deletions +1);

            counter.nPdelDlinsVD = zeros(model.max_D_deletions +1, model.max_insertions +1);
            counter.nPdelDlinsDJ = zeros(model.max_D_deletions +1, model.max_insertions +1);
            counter.nPdelDldelDr = zeros(model.max_D_deletions +1, model.max_D_deletions +1);

            counter.nPdelDrinsVD = zeros(model.max_D_deletions +1, model.max_insertions +1);
            counter.nPdelDrinsDJ = zeros(model.max_D_deletions +1, model.max_insertions +1);

            counter.nPinsVDinsDJ = zeros(model.max_insertions +1, model.max_insertions +1);

            counter.nPVdelV = zeros(size(model.PV,1), model.max_V_deletions+1 );
            counter.nPDdelDl = zeros(size(model.PDJ,1), model.max_D_deletions+1 );
            counter.nPDdelDr = zeros(size(model.PDJ,1), model.max_D_deletions+1 );
            counter.nPJdelJ = zeros(size(model.PDJ,2), model.max_J_deletions+1 );

            counter.nPVinsVD = zeros(size(model.PV,1), model.max_insertions+1 );
            counter.nPDinsVD = zeros(size(model.PDJ,1), model.max_insertions+1 );
            counter.nPDinsDJ = zeros(size(model.PDJ,1), model.max_insertions+1 );
            counter.nPJinsDJ = zeros(size(model.PDJ,2), model.max_insertions+1 );

            counter.nPVdelDl = zeros(size(model.PV,1), model.max_D_deletions+1 );
            counter.nPVdelDr = zeros(size(model.PV,1), model.max_D_deletions+1 );
            counter.nPVdelJ = zeros(size(model.PV,1), model.max_J_deletions+1 );

            counter.nPJdelV = zeros(size(model.PDJ,2), model.max_V_deletions+1 );
            counter.nPJdelDl = zeros(size(model.PDJ,2), model.max_D_deletions+1 );
            counter.nPJdelDr = zeros(size(model.PDJ,2), model.max_D_deletions+1 );

            counter.nPDdelV = zeros(size(model.PDJ,1), model.max_V_deletions+1 );
            counter.nPDdelJ = zeros(size(model.PDJ,1), model.max_J_deletions+1 );

            counter.nPVinsDJ = zeros(size(model.PV,1), model.max_insertions+1 );
            counter.nPJinsVD = zeros(size(model.PDJ,2), model.max_insertions+1 );

            counter.nPpVinsVD = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpVinsDJ = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpVdelDl = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPpVdelDr = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPpVdelJ = zeros(model.max_palindrome +1, model.max_J_deletions +1);
            counter.nPVpV = zeros(size(model.PV,1), model.max_palindrome+1 );
            counter.nPJpV = zeros(size(model.PDJ,2), model.max_palindrome+1 );
            counter.nPDpV = zeros(size(model.PDJ,1), model.max_palindrome+1 );

            counter.nPpJinsVD = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpJinsDJ = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpJdelV = zeros(model.max_palindrome +1, model.max_V_deletions +1);
            counter.nPpJdelDl = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPpJdelDr = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPVpJ = zeros(size(model.PV,1), model.max_palindrome+1 );
            counter.nPJpJ = zeros(size(model.PDJ,2), model.max_palindrome+1 );
            counter.nPDpJ = zeros(size(model.PDJ,1), model.max_palindrome+1 );

            counter.nPpDlinsVD = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpDlinsDJ = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpDldelV = zeros(model.max_palindrome +1, model.max_V_deletions +1);
            counter.nPpDldelJ = zeros(model.max_palindrome +1, model.max_J_deletions +1);
            counter.nPpDldelDr = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPVpDl = zeros(size(model.PV,1), model.max_palindrome+1 );
            counter.nPJpDl = zeros(size(model.PDJ,2), model.max_palindrome+1 );
            counter.nPDpDl = zeros(size(model.PDJ,1), model.max_palindrome+1 );

            counter.nPpDrinsVD = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpDrinsDJ = zeros(model.max_palindrome +1, model.max_insertions +1);
            counter.nPpDrdelV = zeros(model.max_palindrome +1, model.max_V_deletions +1);
            counter.nPpDrdelJ = zeros(model.max_palindrome +1, model.max_J_deletions +1);
            counter.nPpDrdelDl = zeros(model.max_palindrome +1, model.max_D_deletions +1);
            counter.nPVpDr = zeros(size(model.PV,1), model.max_palindrome+1 );
            counter.nPJpDr = zeros(size(model.PDJ,2), model.max_palindrome+1 );
            counter.nPDpDr = zeros(size(model.PDJ,1), model.max_palindrome+1 );

            counter.nPpVpDl = zeros(model.max_palindrome + 1, model.max_palindrome + 1);
            counter.nPpVpDr = zeros(model.max_palindrome + 1, model.max_palindrome + 1);
            counter.nPpVpJ = zeros(model.max_palindrome + 1, model.max_palindrome + 1);
            counter.nPpDlpDr = zeros(model.max_palindrome + 1, model.max_palindrome + 1);
            counter.nPpDlpJ = zeros(model.max_palindrome + 1, model.max_palindrome + 1);
            counter.nPpDrpJ = zeros(model.max_palindrome + 1, model.max_palindrome + 1);


            %
            counter.nMerror_vs_position = zeros(model.read_length,1);
            counter.nMcoverage = zeros(model.read_length,1);




        end

        % Function that compiles the assignments script.
        function call_codegen(model, genV, genD, genJ)
            model_name = class(model);

            counter = model.initialize_counter();
            assigns = initialize_assigns(model, counter);
            ps = fieldnames(assigns);
            for p=1:length(ps)
                A.(ps{p}) = coder.typeof( assigns.(ps{p}) , [], 1 );
            end
            max_D_aligns = 200;
            v_max = length(genV);
            d_max = length(genD);
            j_max = length(genJ);
            vg_max = max([genV.gene_index]);
            dg_max = max([genD.gene_index]);
            jg_max = max([genJ.gene_index]);
            V_allowed_errors_max=20;
            J_allowed_errors_max=10;
            D_allowed_errors_max=4;

            negative_excess_deletions_max = 3;

            V_max_deletions_alignment = 20;
            J_max_deletions_alignment = 18;
            D_max_deletions_alignment = 16;
            J_max_length = 22;

            seq = repmat('A',1,150);
            V.align_length = zeros(v_max,1);
            V.align_position = zeros(2,v_max);
            V.min_deletions = zeros(v_max,1);
            V.n_errors = zeros(v_max,1);
            V.error_positions = zeros(v_max, V_allowed_errors_max); % to force second dimension to be variable sized.
            V.p_region_max_length = zeros(v_max,V_max_deletions_alignment+1);
            V.excess_error_positions = zeros(v_max,negative_excess_deletions_max);
            V.alleles_all = [1:v_max]';
            V.alleles_from_distinct_genes = [1:vg_max]';

            J.align_length = zeros(j_max,1);
            J.align_position = zeros(2,j_max);
            J.primer_start_from_right = zeros(j_max,1);
            J.min_deletions = zeros(j_max,1);
            J.n_errors = zeros(j_max,1);
            J.error_positions = zeros(j_max, J_allowed_errors_max);
            J.p_region_max_length = zeros(j_max,J_max_length+1);
            J.excess_error_positions = zeros(j_max,negative_excess_deletions_max);
            J.alleles_all = [1:j_max]';
            J.alleles_from_distinct_genes = [1:jg_max]';

            D.align_length = zeros(max_D_aligns, d_max);
            D.score = zeros(max_D_aligns, d_max);
            D.n_errors = zeros(max_D_aligns, d_max);
            D.error_positions = zeros(max_D_aligns,D_allowed_errors_max,d_max);
            D.excess_error_positions_left = zeros(max_D_aligns, negative_excess_deletions_max, d_max);
            D.excess_error_positions_right = zeros(max_D_aligns, negative_excess_deletions_max, d_max);
            D.align_position_left = zeros(max_D_aligns, d_max);
            D.align_position_right = zeros(max_D_aligns, d_max);
            D.deletions_left = zeros(max_D_aligns, d_max);
            D.deletions_right = zeros(max_D_aligns, d_max);
            D.p_region_max_length_left = zeros(max_D_aligns, D_max_deletions_alignment +1, d_max);
            D.p_region_max_length_right = zeros(max_D_aligns, D_max_deletions_alignment +1,d_max);
            D.n_alignments = zeros(d_max,1);
            D.allele_order = [1:d_max];

            do_smoothing = true;
            force_all_alleles = true;
            probability_threshold_factor = 0.0001;
            no_error = false;
            ignore_deep_error = false;
            READ_LENGTH_CORRECTION = 0;

            ps = fieldnames(V);
            for p=1:length(ps)
                Vc.(ps{p}) = coder.typeof( V.(ps{p}) , [], 1 );
            end
            ps = fieldnames(D);
            for p=1:length(ps)
                Dc.(ps{p}) = coder.typeof( D.(ps{p}) , [], 1 );
            end
            ps = fieldnames(J);
            for p=1:length(ps)
                Jc.(ps{p}) = coder.typeof( J.(ps{p}) , [], 1 );
            end

            cfg = coder.config('mex');
            cfg.DynamicMemoryAllocation = 'AllVariableSizeArrays';

            eval(['codegen ' model_name '_assignments -config cfg -o ' model_name '_assignments -args {struct(model), coder.typeof(seq,[1 inf]), Vc, Dc, Jc, genV, genD, genJ, A, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION }']);
            %codegen VDJ_pdel_insVD_insDJ_model_assignments -config:mex -args {struct(model), C}
        end
    end % methods

end % classdef
