#include #include #include #include #include #include "calc_stats_2p.hpp" using namespace std; // calculate second order statistics of W containing two populations int calc_N_sec(gsl_matrix *W, int N_nodes[2], double (*N_edge)[2], double (*N_recip)[2], double (*N_conv)[2][2], double (*N_div)[2][2], double (*N_chain)[2][2]) { // create the four submatrices gsl_matrix_view W00 = gsl_matrix_submatrix (W, 0, 0, N_nodes[0], N_nodes[0]); gsl_matrix_view W01 = gsl_matrix_submatrix (W, 0, N_nodes[0], N_nodes[0], N_nodes[1]); gsl_matrix_view W10 = gsl_matrix_submatrix (W, N_nodes[0], 0, N_nodes[1], N_nodes[0]); gsl_matrix_view W11 = gsl_matrix_submatrix (W, N_nodes[0], N_nodes[0], N_nodes[1], N_nodes[1]); // create a matrix of pointers to the submatrices gsl_matrix *Wsub[2][2]; Wsub[0][0] = &W00.matrix; Wsub[0][1] = &W01.matrix; Wsub[1][0] = &W10.matrix; Wsub[1][1] = &W11.matrix; // calculate number of edge entries for(int ind1=0; ind1<2; ind1++) for(int ind2=0; ind2<2; ind2++) { N_edge[ind1][ind2] = 0.0; for(int k=0; k= ind1) { // N_recip is trace gsl_vector_view W2diag = gsl_matrix_diagonal(Wmult[ind1][ind1]); N_recip[ind1][ind2] = gsl_blas_dasum(&W2diag.vector); // need to subtract off N_chain N_chain[ind1][ind2][ind1] -= N_recip[ind1][ind2]; } else { // this can only occur for ind1=2, ind2=1 // so N_recip[1][2] has already been calculated // need to subtract it off N_chain N_chain[ind1][ind2][ind1] -= N_recip[ind2][ind1]; } for(int ind3=0; ind3<2; ind3++) { // finish N_chain for ind3 != ind1 if(ind3 != ind1) { // compute Wab*Wbc, store in Wmult gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, Wsub[ind1][ind2], Wsub[ind2][ind3], 0.0, Wmult[ind1][ind3]); // N_chain^abc is sum of Wab*Wbc // don't have to subtract N_recip here N_chain[ind1][ind2][ind3] = 0.0; for(int i=0; i= ind2) { // now use Wmult for Wab^T*Wac gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, Wsub[ind1][ind2], Wsub[ind1][ind3], 0.0, Wmult[ind2][ind3]); // first part of N_conv^abc is sum of Wab^T*Wac N_conv[ind1][ind2][ind3]=0; for(int i=0; i= ind1) alphahat_recip[ind1][ind2] = N_recip[ind1][ind2] / (denom*phat[ind1][ind2]*phat[ind2][ind1]) - 1.0; for(int ind3=0; ind3<2; ind3++) { denom = N_nodes[ind1]*(double) N_nodes[ind2]*(double)N_nodes[ind3]; if(ind1 == ind2) { if(ind1 == ind3) denom += -3*gsl_pow_2(N_nodes[ind1]) + 2*N_nodes[ind1]; else denom -= N_nodes[ind1]*N_nodes[ind3]; } else if(ind1 == ind3) denom -= N_nodes[ind1]*N_nodes[ind2]; else if(ind2 == ind3) denom -= N_nodes[ind1]*N_nodes[ind2]; if(ind3 >= ind2) { alphahat_conv[ind1][ind2][ind3] = N_conv[ind1][ind2][ind3] / (denom*phat[ind1][ind2]*phat[ind1][ind3]) - 1.0; alphahat_div[ind2][ind3][ind1] = N_div[ind2][ind3][ind1] / (denom*phat[ind2][ind1]*phat[ind3][ind1]) - 1.0; } alphahat_chain[ind1][ind2][ind3] = N_chain[ind1][ind2][ind3] / (denom*phat[ind1][ind2]*phat[ind2][ind3]) - 1.0; } } return 0; } int calc_gaus_covs(gsl_matrix *W, int N_nodes[2], double (*sigma)[2], double (*cov_recip)[2], double (*cov_conv)[2][2], double (*cov_div)[2][2], double (*cov_chain)[2][2], double &cov_other) { int N_shift[2]; N_shift[0]=0; N_shift[1]=N_nodes[0]; int N_nodes_tot= N_nodes[0]+N_nodes[1]; // calc covariances of W (assume everything mean zero) cov_other=0.0; for(int i=0; i<2; i++) for(int j=0; j<2; j++) { sigma[i][j]=0.0; cov_recip[i][j]=0.0; for(int k=0; k<2; k++) { cov_conv[i][j][k]=0.0; cov_div[i][j][k]=0.0; cov_chain[i][j][k]=0.0; } } for(int ntype1=0; ntype1<2; ntype1++) for(int i=N_shift[ntype1]; i