Commit fdd2036c authored by Alois SCHLOEGL's avatar Alois SCHLOEGL

initial import, assumed to be closest to the original code; obtained from...

initial import, assumed to be closest to the original code; obtained from /fs3/group/jonasgrp/Jose/SONET
parents
#include<iostream>
#include<stdio.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_math.h>
#include<gsl/gsl_roots.h>
#include<gsl/gsl_cdf.h>
#include<gsl/gsl_integration.h>
#include "calc_rhos.hpp"
using namespace std;
// function and structure declarations
double integrate_gaussian(double rho0, void *parameters);
struct Qfunction_params { double xth, rho; };
struct int_gauss_params{ double p1, p2, sec_mom; };
////////////////////////////////////////////////////
// calc_rho_given_alpha
// need to calculate value of rho, the correlation
// between two Gaussian random variables
// so that the covariance of the resulting Bernoulli
// variables will be alpha*p1*p2
// solve this 1D equation for rho numerically
// (The variance of each Gaussian is determined by
// the requirement that the Bernoulli obtained by
// thresholding at 1 has probability pj of being 1.)
/////////////////////////////////////////////////////
double calc_rho_given_alpha(double p1, double p2, double alpha,
int &status) {
// if correlation or a probability is zero, return rho=0
if(alpha==0 || p1==0 || p2==0) {
status = 0;
return 0.0;
}
if(alpha < -1) {
cerr << "alpha < -1, cannot calc rho" << endl;
status = -1;
return 0.0;
}
// desired second moment among bernoulli random variables
double b_sec_mom = p1*p2*(1+alpha);
// set up gsl function FF to point to function
// integrate_gaussian with
// parameters determined by arguments p and bcorr
struct int_gauss_params para;
para.p1 = p1;
para.p2 = p2;
para.sec_mom = b_sec_mom;
gsl_function FF;
FF.function = &integrate_gaussian;
FF.params = &para;
int iter=0, max_iter=1000;
double rho=0;
// we know rho has to be in interval [x_lo,x_hi]
double x_lo = -1, x_hi=1;
if (alpha > 0)
x_lo=0;
else x_hi=0;
// Initialize solver and iterate until converge to solution
const gsl_root_fsolver_type *T=gsl_root_fsolver_brent;
gsl_root_fsolver *s = gsl_root_fsolver_alloc(T);
gsl_root_fsolver_set(s, &FF, x_lo, x_hi);
do {
iter++;
status = gsl_root_fsolver_iterate (s);
if(status)
break;
rho = gsl_root_fsolver_root (s);
x_lo = gsl_root_fsolver_x_lower (s);
x_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (x_lo, x_hi, 0, 0.00001);
}
while (status == GSL_CONTINUE && iter < max_iter);
gsl_root_fsolver_free (s);
return rho;
}
// Q function, when integrated from xth to infinity,
// gives the second moment of the Bernoulli random variables
double Qfunction(double x, void *parameters){
struct Qfunction_params *para = (struct Qfunction_params *)parameters;
double rho = para->rho;
double xth = para->xth;
double Q=gsl_cdf_gaussian_Q(xth-rho*x, sqrt(1-rho*rho));
double g= 1/sqrt(2*M_PI)*exp(-0.5*x*x)*Q;
return g;
}
// Take the integral of the Gaussian that corresponds to the
// second moment of the Bernoulli random variables.
// Subtract off their required value so that function is
// zero when found correct value of rho
double integrate_gaussian(double rho0, void *parameters){
struct int_gauss_params * para = (struct int_gauss_params *) parameters;
double p1=para->p1;
double p2=para->p2;
double sec_mom=para->sec_mom;
double rho=rho0;
double xth1=gsl_cdf_ugaussian_Qinv(p1);
double xth2=gsl_cdf_ugaussian_Qinv(p2);
struct Qfunction_params qparams;
qparams.xth=xth1;
qparams.rho=rho;
gsl_function G;
G.function = &Qfunction;
G.params = &qparams;
gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
double res, err;
gsl_integration_qagiu(&G, xth2, 1e-7, 1e-7, 1000, ww, &res, &err);
gsl_integration_workspace_free (ww);
res -= sec_mom;
return res;
}
double calc_rho_given_alpha(double p1, double p2, double alpha,
int &status);
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include<gsl/gsl_eigen.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_cdf.h>
#include "calc_sqrtcov_gen.hpp"
using namespace std;
int calc_sqrtcov_given_rhos_large_N
(int N_pops, int *N_nodes, double *sigma, double *rho_recip, double *rho_conv,
double *rho_div, double *rho_chain,
double *sqrt_diag, double *sqrt_recip, double *sqrt_conv,
double *sqrt_div, double *sqrt_chain,
double *sqrt_diag_last, double *sqrt_recip_last, double *sqrt_conv_last,
double *sqrt_div_last, double *sqrt_chain_last,
double cov_mult=1);
int calc_sqrtcov_given_rhos_refine
(int N_pops, int *N_nodes, double *sigma, double *rho_recip, double *rho_conv,
double *rho_div, double *rho_chain,
double *sqrt_diag, double *sqrt_recip, double *sqrt_conv,
double *sqrt_div, double *sqrt_chain);
int calc_sqrtcov_given_rhos
(int N_pops, int *N_nodes, double *sigma, double *rho_recip, double *rho_conv,
double *rho_div, double *rho_chain,
double *sqrt_diag, double *sqrt_recip, double *sqrt_conv,
double *sqrt_div, double *sqrt_chain) {
double *sqrt_diag_last=new double[N_pops*N_pops];
double *sqrt_recip_last=new double[N_pops*N_pops];
double *sqrt_conv_last=new double[N_pops*N_pops*N_pops];
double *sqrt_div_last=new double[N_pops*N_pops*N_pops];
double *sqrt_chain_last=new double[N_pops*N_pops*N_pops];
// initialize all sqrt's to zero
for(int i=0; i<N_pops; i++)
for(int j=0; j<N_pops; j++) {
sqrt_diag[i*N_pops+j]=sqrt_recip[i*N_pops+j]=0.0;
sqrt_diag_last[i*N_pops+j]=sqrt_recip_last[i*N_pops+j]=0.0;
for(int k=0; k<N_pops; k++) {
sqrt_conv[(i*N_pops+j)*N_pops+k]=
sqrt_div[(i*N_pops+j)*N_pops+k]=
sqrt_chain[(i*N_pops+j)*N_pops+k]=0.0;
sqrt_conv_last[(i*N_pops+j)*N_pops+k]=
sqrt_div_last[(i*N_pops+j)*N_pops+k]=
sqrt_chain_last[(i*N_pops+j)*N_pops+k]=0.0;
}
}
int status;
status= calc_sqrtcov_given_rhos_large_N
(N_pops, N_nodes, sigma, rho_recip, rho_conv, rho_div, rho_chain,
sqrt_diag, sqrt_recip, sqrt_conv, sqrt_div, sqrt_chain,
sqrt_diag_last, sqrt_recip_last,
sqrt_conv_last, sqrt_div_last, sqrt_chain_last,
0);
if(status) {
delete [] sqrt_diag_last;
delete [] sqrt_recip_last;
delete [] sqrt_conv_last;
delete [] sqrt_div_last;
delete [] sqrt_chain_last;
return status;
}
cout << "Found large N square root\n";
cout.flush();
// status = calc_sqrtcov_given_rhos_refine
// (N_pops, N_nodes, sigma, rho_recip, rho_conv, rho_div, rho_chain,
// sqrt_diag, sqrt_recip, sqrt_conv, sqrt_div, sqrt_chain);
for(int i=0; i<N_pops; i++)
for(int j=0; j<N_pops; j++) {
sqrt_diag_last[i*N_pops+j]=sqrt_diag[i*N_pops+j];
sqrt_recip_last[i*N_pops+j]=sqrt_recip[i*N_pops+j];
for(int k=0; k<N_pops; k++) {
sqrt_conv_last[(i*N_pops+j)*N_pops+k]
=sqrt_conv[(i*N_pops+j)*N_pops+k];
sqrt_div_last[(i*N_pops+j)*N_pops+k]
=sqrt_div[(i*N_pops+j)*N_pops+k];
sqrt_chain_last[(i*N_pops+j)*N_pops+k]
=sqrt_chain[(i*N_pops+j)*N_pops+k];
}
}
double d_cov_mult=1.0;
double cov_mult=d_cov_mult;
for(int i=0; i<2000; i++) {
status= calc_sqrtcov_given_rhos_large_N
(N_pops, N_nodes, sigma, rho_recip, rho_conv, rho_div, rho_chain,
sqrt_diag, sqrt_recip, sqrt_conv, sqrt_div, sqrt_chain,
sqrt_diag_last, sqrt_recip_last,
sqrt_conv_last, sqrt_div_last, sqrt_chain_last,
cov_mult);
if(status) {
// cov_mult /= 2.0;
// d_cov_mult /= 2.0;
if(status>0)
cout << "Found " << status << " negative evals\n";
else
exit(-1);
//cout << "Reduced d_cov_mult to " << d_cov_mult << "\n";
}
//else {
double diff=0.0;
double sum=0.0;
// it doesn't hurt to calculate diff of non-existent entries
// since they are all zero
for(int i=0; i<N_pops; i++)
for(int j=0; j<N_pops; j++) {
diff += gsl_pow_2(sqrt_diag_last[i*N_pops+j]-sqrt_diag[i*N_pops+j]);
sum += gsl_pow_2(sqrt_diag[i*N_pops+j]);
sqrt_diag_last[i*N_pops+j]=sqrt_diag[i*N_pops+j];
diff += gsl_pow_2(sqrt_recip_last[i*N_pops+j]-sqrt_recip[i*N_pops+j]);
sum += gsl_pow_2(sqrt_recip[i*N_pops+j]);
sqrt_recip_last[i*N_pops+j]=sqrt_recip[i*N_pops+j];
for(int k=0; k<N_pops; k++) {
diff += gsl_pow_2(sqrt_conv_last[(i*N_pops+j)*N_pops+k]
-sqrt_conv[(i*N_pops+j)*N_pops+k]);
sum += gsl_pow_2(sqrt_conv[(i*N_pops+j)*N_pops+k]);
sqrt_conv_last[(i*N_pops+j)*N_pops+k]
=sqrt_conv[(i*N_pops+j)*N_pops+k];
diff += gsl_pow_2(sqrt_div_last[(i*N_pops+j)*N_pops+k]
-sqrt_div[(i*N_pops+j)*N_pops+k]);
sum += gsl_pow_2(sqrt_div[(i*N_pops+j)*N_pops+k]);
sqrt_div_last[(i*N_pops+j)*N_pops+k]
=sqrt_div[(i*N_pops+j)*N_pops+k];
diff += gsl_pow_2(sqrt_chain_last[(i*N_pops+j)*N_pops+k]
-sqrt_chain[(i*N_pops+j)*N_pops+k]);
sum += gsl_pow_2(sqrt_chain[(i*N_pops+j)*N_pops+k]);
sqrt_chain_last[(i*N_pops+j)*N_pops+k]
=sqrt_chain[(i*N_pops+j)*N_pops+k];
}
}
cout << "diff = " << diff << ", sum = " << sum << ", ratio = "
<< diff/sum << "\n";
if(cov_mult==1.0) {
if(diff/sum < 1E-5)
break;
}
else { //if(diff/sum < 1E-5) {
cov_mult += d_cov_mult;
if(cov_mult >1.0)
cov_mult=1.0;
cout << "cov_mult = " << cov_mult << "\n";
}
//}
}
// cout << "before refine\n";
// cout << "sqrt_diag = \n";
// for(int i=0; i<N_pops; i++) {
// for(int j=0; j<N_pops; j++)
// cout << sqrt_diag[i*N_pops+j] << " ";
// cout << "\n";
// }
// cout << "\n";
// cout << "sqrt_recip = \n";
// for(int i=0; i<N_pops; i++) {
// for(int j=i; j<N_pops; j++)
// cout << sqrt_recip[i*N_pops+j] << " ";
// cout << "\n";
// }
// cout << "\n";
// cout << "sqrt_conv = \n";
// for(int i=0; i<N_pops; i++) {
// for(int j=0; j<N_pops; j++) {
// for(int k=j; k<N_pops; k++)
// cout << sqrt_conv[(i*N_pops+j)*N_pops+k] << " ";
// cout << "\n";
// }
// cout << "\n";
// }
// cout << "sqrt_div = \n";
// for(int i=0; i<N_pops; i++) {
// for(int j=i; j<N_pops; j++) {
// for(int k=0; k<N_pops; k++)
// cout << sqrt_div[(i*N_pops+j)*N_pops+k] << " ";
// cout << "\n";
// }
// cout << "\n";
// }
// cout << "sqrt_chain = \n";
// for(int i=0; i<N_pops; i++) {
// for(int j=0; j<N_pops; j++) {
// for(int k=0; k<N_pops; k++)
// cout << sqrt_chain[(i*N_pops+j)*N_pops+k] << " ";
// cout << "\n";
// }
// cout << "\n";
// }
// status = calc_sqrtcov_given_rhos_refine
// (N_pops, N_nodes, sigma, rho_recip, rho_conv, rho_div, rho_chain,
// sqrt_diag, sqrt_recip, sqrt_conv, sqrt_div, sqrt_chain);
// double diff=0.0;
// double sum=0.0;
// // it doesn't hurt to calculate diff of non-existent entries
// // since they are all zero
// for(int i=0; i<N_pops; i++)
// for(int j=0; j<N_pops; j++) {
// diff += gsl_pow_2(sqrt_diag_last[i*N_pops+j]-sqrt_diag[i*N_pops+j]);
// sum += gsl_pow_2(sqrt_diag[i*N_pops+j]);
// sqrt_diag_last[i*N_pops+j]=sqrt_diag[i*N_pops+j];
// diff += gsl_pow_2(sqrt_recip_last[i*N_pops+j]-sqrt_recip[i*N_pops+j]);
// sum += gsl_pow_2(sqrt_recip[i*N_pops+j]);
// sqrt_recip_last[i*N_pops+j]=sqrt_recip[i*N_pops+j];
// for(int k=0; k<N_pops; k++) {
// diff += gsl_pow_2(sqrt_conv_last[(i*N_pops+j)*N_pops+k]
// -sqrt_conv[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_conv[(i*N_pops+j)*N_pops+k]);
// sqrt_conv_last[(i*N_pops+j)*N_pops+k]
// =sqrt_conv[(i*N_pops+j)*N_pops+k];
// diff += gsl_pow_2(sqrt_div_last[(i*N_pops+j)*N_pops+k]
// -sqrt_div[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_div[(i*N_pops+j)*N_pops+k]);
// sqrt_div_last[(i*N_pops+j)*N_pops+k]
// =sqrt_div[(i*N_pops+j)*N_pops+k];
// diff += gsl_pow_2(sqrt_chain_last[(i*N_pops+j)*N_pops+k]
// -sqrt_chain[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_chain[(i*N_pops+j)*N_pops+k]);
// sqrt_chain_last[(i*N_pops+j)*N_pops+k]
// =sqrt_chain[(i*N_pops+j)*N_pops+k];
// }
// }
// cout << "diff = " << diff << ", sum = " << sum << ", ratio = "
// << diff/sum << "\n";
// for(int i=0; i<10; i++) {
// status= calc_sqrtcov_given_rhos_large_N
// (N_pops, N_nodes, sigma, rho_recip, rho_conv, rho_div, rho_chain,
// sqrt_diag, sqrt_recip, sqrt_conv, sqrt_div, sqrt_chain,
// sqrt_diag_last, sqrt_recip_last,
// sqrt_conv_last, sqrt_div_last, sqrt_chain_last);
// if(0) {//status) {
// delete [] sqrt_diag_last;
// delete [] sqrt_recip_last;
// delete [] sqrt_conv_last;
// delete [] sqrt_div_last;
// delete [] sqrt_chain_last;
// return status;
// }
// double diff=0.0;
// double sum=0.0;
// // it doesn't hurt to calculate diff of non-existent entries
// // since they are all zero
// for(int i=0; i<N_pops; i++)
// for(int j=0; j<N_pops; j++) {
// diff += gsl_pow_2(sqrt_diag_last[i*N_pops+j]-sqrt_diag[i*N_pops+j]);
// sum += gsl_pow_2(sqrt_diag[i*N_pops+j]);
// sqrt_diag_last[i*N_pops+j]=sqrt_diag[i*N_pops+j];
// diff += gsl_pow_2(sqrt_recip_last[i*N_pops+j]-sqrt_recip[i*N_pops+j]);
// sum += gsl_pow_2(sqrt_recip[i*N_pops+j]);
// sqrt_recip_last[i*N_pops+j]=sqrt_recip[i*N_pops+j];
// for(int k=0; k<N_pops; k++) {
// diff += gsl_pow_2(sqrt_conv_last[(i*N_pops+j)*N_pops+k]
// -sqrt_conv[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_conv[(i*N_pops+j)*N_pops+k]);
// sqrt_conv_last[(i*N_pops+j)*N_pops+k]
// =sqrt_conv[(i*N_pops+j)*N_pops+k];
// diff += gsl_pow_2(sqrt_div_last[(i*N_pops+j)*N_pops+k]
// -sqrt_div[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_div[(i*N_pops+j)*N_pops+k]);
// sqrt_div_last[(i*N_pops+j)*N_pops+k]
// =sqrt_div[(i*N_pops+j)*N_pops+k];
// diff += gsl_pow_2(sqrt_chain_last[(i*N_pops+j)*N_pops+k]
// -sqrt_chain[(i*N_pops+j)*N_pops+k]);
// sum += gsl_pow_2(sqrt_chain[(i*N_pops+j)*N_pops+k]);
// sqrt_chain_last[(i*N_pops+j)*N_pops+k]
// =sqrt_chain[(i*N_pops+j)*N_pops+k];
// }
// }
// cout << "diff = " << diff << ", sum = " << sum << ", ratio = "
// << diff/sum << "\n";
// if(diff/sum < 1E-12)
// break;
// }
delete [] sqrt_diag_last;
delete [] sqrt_recip_last;
delete [] sqrt_conv_last;
delete [] sqrt_div_last;
delete [] sqrt_chain_last;
return 0;//status;
}
// calculate the components of the sqrt of the covariance matrix
// assuming that the number of neurons N_nodes are large
int calc_sqrtcov_given_rhos_large_N
(int N_pops, int *N_nodes, double *sigma, double *rho_recip, double *rho_conv,
double *rho_div, double *rho_chain,
double *sqrt_diag, double *sqrt_recip, double *sqrt_conv,
double *sqrt_div, double *sqrt_chain,
double *sqrt_diag_last, double *sqrt_recip_last, double *sqrt_conv_last,
double *sqrt_div_last, double *sqrt_chain_last,
double cov_mult) {
int found_neg_eval=0;
double largest_neg_eval=0.0;
int N_nodes_tot=0;
for(int i=0; i<N_pops; i++)
N_nodes_tot += N_nodes[i];
double *N_nodes_sqrt = new double[N_pops];
for(int i=0; i<N_pops; i++)
N_nodes_sqrt[i] = sqrt((double)N_nodes[i]);
// for large N, the equations for sqrt_conv, sqrt_div, and sqrt_chain
// decouple from the rest. More over, they decouple into N_pops groups
// of N_pops*(2*N_pops+1) equations based on the pop of
// the central neuron in the motif
// moreover, these N_pops*(2*N_pops+1) equations can be written as
// solving for the square root of a 2*N_pops x 2*N_pops covariance matrix
// Hence one can quickly determine if the equations have a real solution
// and find that real solution
const size_t nmat=2*N_pops;
gsl_eigen_symmv_workspace *work_eig=gsl_eigen_symmv_alloc(nmat);
gsl_matrix *A = gsl_matrix_alloc(nmat,nmat); // the 4x4 covariance matrix
gsl_matrix *sqrtA =gsl_matrix_alloc(nmat,nmat); // its square root
gsl_vector *evals=gsl_vector_alloc(nmat); // evals of A
gsl_matrix *evecs=gsl_matrix_alloc(nmat,nmat); // evects of A
for(int pop1=0; pop1 <N_pops; pop1++) {
// Need to scale the rho's appropriately so the equations
// can be written as finding the square root of covariance matrix A
for(int pop2=0; pop2<N_pops; pop2++) {
int first_pop12, second_pop12;
if(pop1 <= pop2) {
first_pop12=pop1;
second_pop12=pop2;
}
else {
first_pop12=pop2;
second_pop12=pop1;
}
for(int pop3=pop2; pop3<N_pops; pop3++ ) {
int first_pop13, second_pop13;
if(pop1 <= pop3) {
first_pop13=pop1;
second_pop13=pop3;
}
else {
first_pop13=pop3;
second_pop13=pop1;
}
double temp_conv = sigma[pop1*N_pops+pop2]*sigma[pop1*N_pops+pop3]
*rho_conv[(pop1*N_pops+pop2)*N_pops+pop3];
// subtract off additional components if sqrt's are nonzero
double finite_N_correct
= -(sqrt_diag_last[pop1*N_pops+pop2]+sqrt_diag_last[pop1*N_pops+pop3])
* sqrt_conv_last[(pop1*N_pops+pop2)*N_pops+pop3]
- sqrt_recip_last[first_pop12*N_pops+second_pop12]
* sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop3]
- sqrt_recip_last[first_pop13*N_pops+second_pop13]
* sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2]
- sqrt_div_last[(first_pop13*N_pops+second_pop13)*N_pops+pop2]
* sqrt_chain_last[(pop1*N_pops+pop3)*N_pops+pop2]
- sqrt_div_last[(first_pop12*N_pops+second_pop12)*N_pops+pop3]
* sqrt_chain_last[(pop1*N_pops+pop2)*N_pops+pop3];
finite_N_correct
+=sqrt_conv_last[(pop1*N_pops+first_pop12)*N_pops+second_pop12]
*sqrt_conv_last[(pop1*N_pops+first_pop13)*N_pops+second_pop13]
+sqrt_chain_last[(pop1*N_pops+pop1)*N_pops+pop2]
*sqrt_chain_last[(pop1*N_pops+pop1)*N_pops+pop3];
finite_N_correct
+=sqrt_conv_last[(pop1*N_pops+pop2)*N_pops+pop2]
*sqrt_conv_last[(pop1*N_pops+pop2)*N_pops+pop3]
+sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop2]
*sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop3];
finite_N_correct
+=sqrt_conv_last[(pop1*N_pops+pop2)*N_pops+pop3]
*sqrt_conv_last[(pop1*N_pops+pop3)*N_pops+pop3]
+sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop3];
temp_conv += finite_N_correct*cov_mult;
gsl_matrix_set(A,pop3,pop2,
N_nodes_sqrt[pop2]*N_nodes_sqrt[pop3]*temp_conv);
double temp_div =sigma[pop2*N_pops+pop1]*sigma[pop3*N_pops+pop1]
*rho_div[(pop2*N_pops+pop3)*N_pops+pop1];
// subtract off additional components if sqrt's are nonzero
finite_N_correct
=-(sqrt_diag_last[pop2*N_pops+pop1]+sqrt_diag_last[pop3*N_pops+pop1])
*sqrt_div_last[(pop2*N_pops+pop3)*N_pops+pop1]
-sqrt_recip_last[first_pop12*N_pops+second_pop12]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2]
-sqrt_recip_last[first_pop13*N_pops+second_pop13]
*sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop3]
-sqrt_conv_last[(pop2*N_pops+first_pop13)*N_pops+second_pop13]
*sqrt_chain_last[(pop2*N_pops+pop3)*N_pops+pop1]
-sqrt_conv_last[(pop3*N_pops+first_pop12)*N_pops+second_pop12]
*sqrt_chain_last[(pop3*N_pops+pop2)*N_pops+pop1];
finite_N_correct
+=sqrt_div_last[(first_pop12*N_pops+second_pop12)*N_pops+pop1]
*sqrt_div_last[(first_pop13*N_pops+second_pop13)*N_pops+pop1]
+sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop1]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop1];
finite_N_correct
+=sqrt_div_last[(pop2*N_pops+pop2)*N_pops+pop1]
*sqrt_div_last[(pop2*N_pops+pop3)*N_pops+pop1]
+sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop2]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2];
finite_N_correct
+=sqrt_div_last[(pop2*N_pops+pop3)*N_pops+pop1]
*sqrt_div_last[(pop3*N_pops+pop3)*N_pops+pop1]
+sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop3]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop3];
temp_div += finite_N_correct*cov_mult;
gsl_matrix_set(A,pop3+N_pops,pop2+N_pops,
N_nodes_sqrt[pop2]*N_nodes_sqrt[pop3]*temp_div);
}
for(int pop3=0; pop3<N_pops; pop3++ ) {
int first_pop13, second_pop13, first_pop23, second_pop23;
if(pop1 <= pop3) {
first_pop13=pop1;
second_pop13=pop3;
}
else {
first_pop13=pop3;
second_pop13=pop1;
}
if(pop2 <= pop3) {
first_pop23=pop2;
second_pop23=pop3;
}
else {
first_pop23=pop3;
second_pop23=pop2;
}
double temp_chain = sigma[pop3*N_pops+pop1]*sigma[pop1*N_pops+pop2]
*rho_chain[(pop3*N_pops+pop1)*N_pops+pop2];
// subtract off additional components if sqrt's are nonzero
double finite_N_correct
=- (sqrt_diag_last[pop3*N_pops+pop1]+sqrt_diag_last[pop1*N_pops+pop2])
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2]
-sqrt_recip_last[first_pop13*N_pops+second_pop13]
*sqrt_conv_last[(pop1*N_pops+first_pop23)*N_pops+second_pop23]
-sqrt_recip_last[first_pop12*N_pops+second_pop12]
*sqrt_div_last[(first_pop23*N_pops+second_pop23)*N_pops+pop1]
-sqrt_conv_last[(pop3*N_pops+first_pop12)*N_pops+second_pop12]
*sqrt_div_last[(first_pop13*N_pops+second_pop13)*N_pops+pop2]
-sqrt_chain_last[(pop2*N_pops+pop3)*N_pops+pop1]
*sqrt_chain_last[(pop1*N_pops+pop2)*N_pops+pop3];
// subtract off overcounting from when pop2 match pop3,pop1,pop2
finite_N_correct
+= sqrt_conv_last[(pop1*N_pops+first_pop23)*N_pops+second_pop23]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop3]
+sqrt_div_last[(pop3*N_pops+pop3)*N_pops+pop1]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2];
finite_N_correct
+= sqrt_conv_last[(pop1*N_pops+first_pop12)*N_pops+second_pop12]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop1]
+sqrt_div_last[(first_pop13*N_pops+second_pop13)*N_pops+pop1]
*sqrt_chain_last[(pop1*N_pops+pop1)*N_pops+pop2];
finite_N_correct
+= sqrt_conv_last[(pop1*N_pops+pop2)*N_pops+pop2]
*sqrt_chain_last[(pop3*N_pops+pop1)*N_pops+pop2]
+sqrt_div_last[(first_pop23*N_pops+second_pop23)*N_pops+pop1]
*sqrt_chain_last[(pop2*N_pops+pop1)*N_pops+pop2];
temp_chain += finite_N_correct*cov_mult;
gsl_matrix_set(A,pop3+N_pops,pop2,
N_nodes_sqrt[pop2]*N_nodes_sqrt[pop3]*temp_chain);
}
}
// to calculate square root of A
// 1. take it's eigen decomposition
// 2. take the square root of its eigenvalues
// 3. reconstruct with new eigenvalues to get a square root of A
gsl_eigen_symmv(A, evals, evecs, work_eig);
double max_eval=-1E100;
double min_eval=1E100;
for(size_t i=0; i<nmat; i++) {
double the_eval = gsl_vector_get(evals,i);
if(the_eval > max_eval)
max_eval=the_eval;
if(the_eval < min_eval)
min_eval=the_eval;
if(the_eval <= 0) {
if(the_eval > -1E-7) {
// allow eigenvalues to be slightly negative due to
// roundoff error
the_eval=0;
}
else if(the_eval > -1E-0) {
found_neg_eval++;
if(the_eval < largest_neg_eval)
largest_neg_eval = the_eval;
//cout << "found_neg_eval = " << found_neg_eval << "\n";
the_eval=0;
}
else {
// if have a negative eigenvalue, can't take square root
// system of equations does not have a real solution
// (at least in limit of large N)
cout << "Found a negative eval(" << i <<")=" << the_eval
<< " for pop1=" << pop1 << "\n";
delete [] N_nodes_sqrt;
gsl_eigen_symmv_free(work_eig);
gsl_matrix_free(A);
gsl_matrix_free(sqrtA);
gsl_matrix_free(evecs);
gsl_vector_free(evals);
return -1;
}
}
// scale eigenvector by fourth root of eval so
// reconstruction with be based on square root of eval