<commit_before>#include <iostream>
#include <cstdint>
#include <cmath>
#include <vector>
#include <initializer_list>
#include <iostream>
#include <map>
#include <fstream>
#include <memory>
#include "Eigen/Dense"

#include "model.h"


using namespace std;

double DirichletMultinomialLogProbability(double alphas[4], ReadData data) {
	// TODO: Cache most of the math here
	// TODO: Does not include the multinomail coefficient
	int read_count = data.reads[0]+data.reads[1]+data.reads[2]+data.reads[3];
	double alpha_total = alphas[0]+alphas[1]+alphas[2]+alphas[3];
	double result = 0.0;
	for(int i : {0,1,2,3}) {
		for(int x = 0; x < data.reads[i]; ++x) {
			result += log(alphas[i]+x);
		}
	}
	for(int x = 0; x < read_count; ++x)
		result -= log(alpha_total+x);
	return result;
}

DiploidProbs DiploidPopulation(const ModelParams &params, int ref_allele) {
	ReadData d;
	DiploidProbs result;
	double alphas[4];
	for(int i : {0,1,2,3})
		alphas[i] = params.theta*params.nuc_freq[i];
	for(int i : {0,1,2,3}) {
		for(int j=0;j<i;++j) {
			d.key = 0;
			d.reads[ref_allele] = 1;
			d.reads[i] += 1;
			d.reads[j] += 1;
			result[i+j*4] = DirichletMultinomialLogProbability(alphas, d);
			result[j+i*4] = result[i+j*4];
		}
		d.key = 0;
		d.reads[ref_allele] = 1;
		d.reads[i] += 2;
		result[i+i*4] = DirichletMultinomialLogProbability(alphas, d);
	}
	return result.exp();
}

TransitionMatrix F81(const ModelParams &params){
	double beta = 1.0;
	for(auto d : params.nuc_freq)
		beta -= d*d;
	beta = 1.0/beta;
	beta = exp(-beta*params.mutation_rate);
	Eigen::Matrix4d m;
	for(int i : {0,1,2,3}) {
		for(int j : {0,1,2,3}) {
			m(i,j) = params.nuc_freq[i]*(1.0-beta);
		}
		m(i,i) += beta;
    }
    return m;
}



MutationMatrix MutationAccumulation_(const ModelParams &params, bool and_mut){
	TransitionMatrix m = F81(params);
    if(params.ploidy_ancestor == 1){//haploid->haploid design
        if(and_mut){
            return m;
        }
        Eigen::ArrayXXd result = MutationMatrix_HH::Zero();
        //identity matrix initialzer is for Matrices only, these are (despite teh
        //name) arrays so buid it up:
        for(int i : {0,1,2,3}){
            for (int j : {0,1,2,3}){
                if(i == j){
                    result(i,j) = m(i,j);
                }
            }
        }
        return result;
    }
    if(params.ploidy_descendant == 1){//diploid -> haploid design
        result = MutationMarix(16,4)
        for(int i : {0,1,2,3}) {
		    for(int j : {0,1,2,3}) {
			    for(int k : {0,1,2,3}) {
				    result(i*4+j,k) = 0.0;
        				if(!and_mut || i != k)
					result(i*4+j,k) += 0.5*m(i,k);
		        		if(!and_mut || j != k)
					result(i*4+j,k) += 0.5*m(j,k);
                }
            }
	    }
	    return result;
    }
    //only diploid->diploid left 
    result = MutationMarix(16,16);
	for(int i : {0,1,2,3}) {
		for(int j : {0,1,2,3}) {
			for(int k : {0,1,2,3}) {
    			for(int l : {0,1,2,3}) {                
    				result(i*4+j,k*4+l) = 0.0;
    				if(!and_mut || i != k || j != l ){//TODO: check this is right transition prob
					    result(i*4+j,k*4+l) += m(i,k) * m(j,l);
                    }
                }
            }
        }
    }
    return result;
}


DiploidProbs DiploidSequencing(const ModelParams &params, int ref_allele, ReadData data) {
	DiploidProbs result;
	double alphas_total = (1.0-params.phi_diploid)/params.phi_diploid;
	for(int i : {0,1,2,3}) {
		for(int j=0;j<i;++j) {
			double alphas[4];
			for(int k : {0,1,2,3}) {
				if(k == i || k == j)
					alphas[k] = (0.5-params.error_prob/3.0)*alphas_total;
				else
					alphas[k] = (params.error_prob/3.0)*alphas_total;
			}
			result[i*4+j] = DirichletMultinomialLogProbability(alphas, data);
			result[j*4+i] = result[i*4+j];
		}
		double alphas[4];
		for(int k : {0,1,2,3}) {
			if(k == i)
				alphas[k] = (1.0-params.error_prob)*alphas_total;
			else
				alphas[k] = params.error_prob/3.0*alphas_total;
		}
		result[i*4+i] = DirichletMultinomialLogProbability(alphas, data);
	}
	double scale = result.maxCoeff();
	return (result - scale).exp();
}

HaploidProbs HaploidSequencing(const ModelParams &params, int ref_allele, ReadData data) {
	HaploidProbs result;
	double alphas_total = (1.0-params.phi_haploid)/params.phi_haploid;
	for(int i : {0,1,2,3}) {
		double alphas[4];
		for(int k : {0,1,2,3}) {
			if(k == i)
				alphas[k] = (1.0-params.error_prob)*alphas_total;
			else
				alphas[k] = params.error_prob/3.0*alphas_total;
		}
		result[i] = DirichletMultinomialLogProbability(alphas, data);
	}
	double scale = result.maxCoeff();
	return (result - scale).exp();
}

double TetMAProbability(const ModelParams &params, const ModelInput site_data) {
	MutationMatrix_DD m = MutationAccumulation_DD(params, false);
	MutationMatrix_DD mt = MutationAccumulation_DD(params, true);
	MutationMatrix_DD mn = m-mt;	
	DiploidProbs pop_genotypes = DiploidPopulation(params, site_data.reference);
		
	auto it = site_data.all_reads.begin();
	DiploidProbs anc_genotypes = DiploidSequencing(params, site_data.reference, *it);
	anc_genotypes *= pop_genotypes;
	DiploidProbs num_genotypes = anc_genotypes;
	for(++it; it != site_data.all_reads.end(); ++it) {
		DiploidProbs p = DiploidSequencing(params, site_data.reference, *it);
		anc_genotypes *= (m.matrix()*p.matrix()).array();
		num_genotypes *= (mn.matrix()*p.matrix()).array();

	}

//    cerr << "\n" << anc_genotypes.sum() << '\t' << num_genotypes.sum() << endl;
	return 1.0 - num_genotypes.sum()/anc_genotypes.sum();
}

double TetMAProbOneMutation(const ModelParams &params, const ModelInput site_data) {
	MutationMatrix_DD m = MutationAccumulation_DD(params, false);
	MutationMatrix_DD mt = MutationAccumulation_DD(params, true);
	MutationMatrix_DD mn = m-mt;	
	DiploidProbs pop_genotypes = DiploidPopulation(params, site_data.reference);	
	auto it = site_data.all_reads.begin();
	DiploidProbs anc_genotypes = DiploidSequencing(params, site_data.reference, *it);
	anc_genotypes *= pop_genotypes;

  	DiploidProbs denom = anc_genotypes;   //product of p(Ri|A)
    
    DiploidProbs nomut_genotypes = anc_genotypes; //Product of p(Ri & noMutatoin|A)
    DiploidProbs mut_genotypes = DiploidProbs::Zero();      //Sum of p(Ri&Mutation|A=x)
	for(++it; it != site_data.all_reads.end(); ++it) {
        DiploidProbs p = DiploidSequencing(params, site_data.reference, *it);
        DiploidProbs dgen =  (mn.matrix()*p.matrix()).array();
        DiploidProbs agen = (m.matrix()*p.matrix()).array();
        nomut_genotypes *= dgen;
        mut_genotypes += (agen/dgen - 1); //(agen+dgen)/agen
        denom *= agen;
    }
    double result = (nomut_genotypes * mut_genotypes).sum() / denom.sum();
    return(result);
}


//int main(){
//      ModelParams p = { 
//        0.0001, 
//        {0.38, 0.12, 0.12, 0.38}, 
//        1e-2,
//        0.01,
//        0.01,
//        0.05,
//    };
//   MutationMatrix_DH m = MutationAccumulation_DH(p, true);
//   MutationMatrix_DH mn = MutationAccumulation_DH(p, false);
////   cerr << m << endl;
//   cout << mn- m << endl;
//   return 0;
//}
// Uncommon and compile with this:
// clang++ -std=c++11 -Ithird-party/bamtools/src/ -Lboost_progam_options model.cc
//
// to play around with / debug results.
//int main(){
//    ModelParams p = { 
//        0.0001, 
//        {0.38, 0.12, 0.12, 0.38}, 
//        1e-8,
//        0.01,
//        0.01,
//        0.05,
//    };
//    ModelInput two_vars = { 2, 
//        {
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0,  0,  0, 30},
//        { 0, 0,   0, 30},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//
//        }
//    };
//
//
//    ModelInput  one_vars = { 1,
//        {
//        { 0, 30,  0,  0},
//        { 0,  0,  0, 30},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0, 0}}
//    };    
//
//   ModelInput no_vars = {1,
//      {{0,   0,  0, 5},       
//       {3,   0,  0, 4},       
//       {8,   0,  0, 0},
//       {3,   0,  0, 3},     
//       {0,  0,  0, 2},       
//       {2,  0,  0, 3},
//       {0,  0,  0, 0},  
//       {0,  0,  0, 5},       
//       {0,  0,  0, 4},
//       {1,  0,  0, 3},
//       {0,  0,  0, 28},
//       {0 ,  0,  0, 3}}
//   };
//    
//    cout << "___With the no-variant data___" << endl;
////    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,no_vars)<< endl;
////    cout << "P(any|data)= " << TetMAProbability(p,no_vars) << endl;
//
//    cout << TetMAProbability(p,no_vars) << endl;
//    
//    cout << "___With the one-variant data___" << endl;  
////    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,one_vars)<< endl;
////    cout << "P(any|data)= " << TetMAProbability(p,one_vars) << endl;
//      cout << TetMAProbability(p,one_vars) << endl;
//    
//    cout << "___With the two-variant data___" << endl;
//    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,two_vars)<< endl;
//    cout << "P(any|data)= " << TetMAProbability(p,two_vars) << endl;   
//
//    
//    cout << "calculating the same number once: " << TetMAProbOneMutation(p,two_vars) << endl;
//    cout << "then another time: " << TetMAProbOneMutation(p,two_vars) << endl;
//    TetMAProbOneMutation(p,no_vars);
//    cout << "And once more after calling from the the no-vars data: " << TetMAProbOneMutation(p,two_vars) << endl;
//    return 0;
//}
////






<commit_msg>Haploid 'population' probs fxn for #5<commit_after>#include <iostream>
#include <cstdint>
#include <cmath>
#include <vector>
#include <initializer_list>
#include <iostream>
#include <map>
#include <fstream>
#include <memory>
#include "Eigen/Dense"

#include "model.h"


using namespace std;

double DirichletMultinomialLogProbability(double alphas[4], ReadData data) {
	// TODO: Cache most of the math here
	// TODO: Does not include the multinomail coefficient
	int read_count = data.reads[0]+data.reads[1]+data.reads[2]+data.reads[3];
	double alpha_total = alphas[0]+alphas[1]+alphas[2]+alphas[3];
	double result = 0.0;
	for(int i : {0,1,2,3}) {
		for(int x = 0; x < data.reads[i]; ++x) {
			result += log(alphas[i]+x);
		}
	}
	for(int x = 0; x < read_count; ++x)
		result -= log(alpha_total+x);
	return result;
}

DiploidProbs DiploidPopulation(const ModelParams &params, int ref_allele) {
	ReadData d;
	DiploidProbs result;
	double alphas[4];
	for(int i : {0,1,2,3})
		alphas[i] = params.theta*params.nuc_freq[i];
	for(int i : {0,1,2,3}) {
		for(int j=0;j<i;++j) {
			d.key = 0;
			d.reads[ref_allele] = 1;
			d.reads[i] += 1;
			d.reads[j] += 1;
			result[i+j*4] = DirichletMultinomialLogProbability(alphas, d);
			result[j+i*4] = result[i+j*4];
		}
		d.key = 0;
		d.reads[ref_allele] = 1;
		d.reads[i] += 2;
		result[i+i*4] = DirichletMultinomialLogProbability(alphas, d);
	}
	return result.exp();
}

HaploidProbs HaploidPopulation(const ModelParams &params, int ref_allele) {
	ReadData d;
	HaploidProbs result;
	double alphas[4];
	for(int i : {0,1,2,3})
		alphas[i] = params.theta*params.nuc_freq[i];
    for(int i : {0,1,2,3}){
        d.key = 0;
        d.reads[i] = 1;
        d.reads[ref_allele] = 1;
        result[i] = DirichletMultinomialLogProbability(alphas, d);
    }
    return result;
}
        

}



TransitionMatrix F81(const ModelParams &params){
	double beta = 1.0;
	for(auto d : params.nuc_freq)
		beta -= d*d;
	beta = 1.0/beta;
	beta = exp(-beta*params.mutation_rate);
	Eigen::Matrix4d m;
	for(int i : {0,1,2,3}) {
		for(int j : {0,1,2,3}) {
			m(i,j) = params.nuc_freq[i]*(1.0-beta);
		}
		m(i,i) += beta;
    }
    return m;
}



MutationMatrix MutationAccumulation_(const ModelParams &params, bool and_mut){
	TransitionMatrix m = F81(params);
    if(params.ploidy_ancestor == 1){//haploid->haploid design
        if(and_mut){
            return m;
        }
        Eigen::ArrayXXd result = MutationMatrix_HH::Zero();
        //identity matrix initialzer is for Matrices only, these are (despite teh
        //name) arrays so buid it up:
        for(int i : {0,1,2,3}){
            for (int j : {0,1,2,3}){
                if(i == j){
                    result(i,j) = m(i,j);
                }
            }
        }
        return result;
    }
    if(params.ploidy_descendant == 1){//diploid -> haploid design
        result = MutationMarix(16,4)
        for(int i : {0,1,2,3}) {
		    for(int j : {0,1,2,3}) {
			    for(int k : {0,1,2,3}) {
				    result(i*4+j,k) = 0.0;
        				if(!and_mut || i != k)
					result(i*4+j,k) += 0.5*m(i,k);
		        		if(!and_mut || j != k)
					result(i*4+j,k) += 0.5*m(j,k);
                }
            }
	    }
	    return result;
    }
    //only diploid->diploid left 
    result = MutationMarix(16,16);
	for(int i : {0,1,2,3}) {
		for(int j : {0,1,2,3}) {
			for(int k : {0,1,2,3}) {
    			for(int l : {0,1,2,3}) {                
    				result(i*4+j,k*4+l) = 0.0;
    				if(!and_mut || i != k || j != l ){//TODO: check this is right transition prob
					    result(i*4+j,k*4+l) += m(i,k) * m(j,l);
                    }
                }
            }
        }
    }
    return result;
}


DiploidProbs DiploidSequencing(const ModelParams &params, int ref_allele, ReadData data) {
	DiploidProbs result;
	double alphas_total = (1.0-params.phi_diploid)/params.phi_diploid;
	for(int i : {0,1,2,3}) {
		for(int j=0;j<i;++j) {
			double alphas[4];
			for(int k : {0,1,2,3}) {
				if(k == i || k == j)
					alphas[k] = (0.5-params.error_prob/3.0)*alphas_total;
				else
					alphas[k] = (params.error_prob/3.0)*alphas_total;
			}
			result[i*4+j] = DirichletMultinomialLogProbability(alphas, data);
			result[j*4+i] = result[i*4+j];
		}
		double alphas[4];
		for(int k : {0,1,2,3}) {
			if(k == i)
				alphas[k] = (1.0-params.error_prob)*alphas_total;
			else
				alphas[k] = params.error_prob/3.0*alphas_total;
		}
		result[i*4+i] = DirichletMultinomialLogProbability(alphas, data);
	}
	double scale = result.maxCoeff();
	return (result - scale).exp();
}

HaploidProbs HaploidSequencing(const ModelParams &params, int ref_allele, ReadData data) {
	HaploidProbs result;
	double alphas_total = (1.0-params.phi_haploid)/params.phi_haploid;
	for(int i : {0,1,2,3}) {
		double alphas[4];
		for(int k : {0,1,2,3}) {
			if(k == i)
				alphas[k] = (1.0-params.error_prob)*alphas_total;
			else
				alphas[k] = params.error_prob/3.0*alphas_total;
		}
		result[i] = DirichletMultinomialLogProbability(alphas, data);
	}
	double scale = result.maxCoeff();
	return (result - scale).exp();
}

double TetMAProbability(const ModelParams &params, const ModelInput site_data) {
	MutationMatrix_DD m = MutationAccumulation_DD(params, false);
	MutationMatrix_DD mt = MutationAccumulation_DD(params, true);
	MutationMatrix_DD mn = m-mt;	
	DiploidProbs pop_genotypes = DiploidPopulation(params, site_data.reference);
		
	auto it = site_data.all_reads.begin();
	DiploidProbs anc_genotypes = DiploidSequencing(params, site_data.reference, *it);
	anc_genotypes *= pop_genotypes;
	DiploidProbs num_genotypes = anc_genotypes;
	for(++it; it != site_data.all_reads.end(); ++it) {
		DiploidProbs p = DiploidSequencing(params, site_data.reference, *it);
		anc_genotypes *= (m.matrix()*p.matrix()).array();
		num_genotypes *= (mn.matrix()*p.matrix()).array();

	}

//    cerr << "\n" << anc_genotypes.sum() << '\t' << num_genotypes.sum() << endl;
	return 1.0 - num_genotypes.sum()/anc_genotypes.sum();
}

double TetMAProbOneMutation(const ModelParams &params, const ModelInput site_data) {
	MutationMatrix m = MutationAccumulation_DD(params, false);
	MutationMatrix mt = MutationAccumulation_DD(params, true);
	MutationMatrix mn = m-mt;	
	DiploidProbs pop_genotypes = DiploidPopulation(params, site_data.reference);	
	auto it = site_data.all_reads.begin();
	DiploidProbs anc_genotypes = DiploidSequencing(params, site_data.reference, *it);
	anc_genotypes *= pop_genotypes;

  	DiploidProbs denom = anc_genotypes;   //product of p(Ri|A)
    
    DiploidProbs nomut_genotypes = anc_genotypes; //Product of p(Ri & noMutatoin|A)
    DiploidProbs mut_genotypes = DiploidProbs::Zero();      //Sum of p(Ri&Mutation|A=x)
	for(++it; it != site_data.all_reads.end(); ++it) {
        DiploidProbs p = DiploidSequencing(params, site_data.reference, *it);
        DiploidProbs dgen =  (mn.matrix()*p.matrix()).array();
        DiploidProbs agen = (m.matrix()*p.matrix()).array();
        nomut_genotypes *= dgen;
        mut_genotypes += (agen/dgen - 1); //(agen+dgen)/agen
        denom *= agen;
    }
    double result = (nomut_genotypes * mut_genotypes).sum() / denom.sum();
    return(result);
}


//int main(){
//      ModelParams p = { 
//        0.0001, 
//        {0.38, 0.12, 0.12, 0.38}, 
//        1e-2,
//        0.01,
//        0.01,
//        0.05,
//    };
//   MutationMatrix_DH m = MutationAccumulation_DH(p, true);
//   MutationMatrix_DH mn = MutationAccumulation_DH(p, false);
////   cerr << m << endl;
//   cout << mn- m << endl;
//   return 0;
//}
// Uncommon and compile with this:
// clang++ -std=c++11 -Ithird-party/bamtools/src/ -Lboost_progam_options model.cc
//
// to play around with / debug results.
//int main(){
//    ModelParams p = { 
//        0.0001, 
//        {0.38, 0.12, 0.12, 0.38}, 
//        1e-8,
//        0.01,
//        0.01,
//        0.05,
//    };
//    ModelInput two_vars = { 2, 
//        {
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0,  0,  0, 30},
//        { 0, 0,   0, 30},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//
//        }
//    };
//
//
//    ModelInput  one_vars = { 1,
//        {
//        { 0, 30,  0,  0},
//        { 0,  0,  0, 30},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0,  0},
//        { 0, 30,  0, 0}}
//    };    
//
//   ModelInput no_vars = {1,
//      {{0,   0,  0, 5},       
//       {3,   0,  0, 4},       
//       {8,   0,  0, 0},
//       {3,   0,  0, 3},     
//       {0,  0,  0, 2},       
//       {2,  0,  0, 3},
//       {0,  0,  0, 0},  
//       {0,  0,  0, 5},       
//       {0,  0,  0, 4},
//       {1,  0,  0, 3},
//       {0,  0,  0, 28},
//       {0 ,  0,  0, 3}}
//   };
//    
//    cout << "___With the no-variant data___" << endl;
////    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,no_vars)<< endl;
////    cout << "P(any|data)= " << TetMAProbability(p,no_vars) << endl;
//
//    cout << TetMAProbability(p,no_vars) << endl;
//    
//    cout << "___With the one-variant data___" << endl;  
////    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,one_vars)<< endl;
////    cout << "P(any|data)= " << TetMAProbability(p,one_vars) << endl;
//      cout << TetMAProbability(p,one_vars) << endl;
//    
//    cout << "___With the two-variant data___" << endl;
//    cout << "P(one|data)= "<<  TetMAProbOneMutation(p,two_vars)<< endl;
//    cout << "P(any|data)= " << TetMAProbability(p,two_vars) << endl;   
//
//    
//    cout << "calculating the same number once: " << TetMAProbOneMutation(p,two_vars) << endl;
//    cout << "then another time: " << TetMAProbOneMutation(p,two_vars) << endl;
//    TetMAProbOneMutation(p,no_vars);
//    cout << "And once more after calling from the the no-vars data: " << TetMAProbOneMutation(p,two_vars) << endl;
//    return 0;
//}
////







{'fraction_non_alphanumeric': 0.18567109106749347, 'fraction_numerical': 0.03838466663250141, 'mean_word_length': 2.462384669978708, 'pattern_counts': {'":': 0, '<': 213, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5650598', 'n_tokens_mistral': 8708, 'n_tokens_neox': 7914, 'n_words': 1558}