run_ei_balanced.cpp 8.63 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
#include <iostream>
#include <string.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_blas.h>
#include "secorder_rec_2p.hpp"
#include "calc_stats_2p.hpp"

using namespace std;



int main(int argc, char *argv[]) {
  
  // optional flags
  int deterministic_seed = 0; // if nonzero, use for random num gen seed

  int N_nodes[2] = {2000,500};
  double p[2][2] = {{0.01, 0.01},{0.01,0.01}};
  double alpha_recip[2][2] = {{0.0, 0.0},{99,0.0}};


  // form of large N covariance matrix centered around population x
  // conv_x00 conv_x01 chain_0x0 chain_1x0
  //          conv_x11 chain_0x1 chain_1x1
  //                     div_00x   div_01x
  //                               div_11x


  ///////////////////////////////////////////////
  // correlations centered around population zero
  // diagonal entries
  double alpha_conv_hom_00 = 0.5;  // zeros onto zero (1,1)
  double alpha_conv_hom_01 = 0.5;  // ones onto zero  (2,2)
  double alpha_div_hom_00  = 0.5;  // zero onto zeros (3,3)
  double alpha_div_hom_10  = 0.5;  // zero onto ones  (4,4)

  // cross entries for divergence and convergence
  // constrained by the diagonal entries
  double cc_conv_0_01 = 0.5;  // mix onto zero  (2,1)
  double cc_div_01_0  = 0.0;  // zero onto mix  (4,3)

  // remaining four cross entries are chains through population 0
  // constrained by diagonal entries
  double cc_chain_000 = 0.8;  // (3,1)
  double cc_chain_100 = 0.0;  // (4,1)
  double cc_chain_001 = 0.8;  // (3,2)
  double cc_chain_101 = 0.0;  // (4,2)

  ///////////////////////////////////////////////
  // correlations centered around population one
  // diagonal entries
  double alpha_conv_hom_10 = 0.5;  // zeros onto one (1,1)
  double alpha_conv_hom_11 = 0.5;  // ones onto one  (2,2)
  double alpha_div_hom_01  = 0.5;  // one onto zeros (3,3)
  double alpha_div_hom_11  = 0.5;  // one onto ones  (4,4)

  // cross entries for divergence and convergence
  // constrained by the diagonal entries
  double cc_conv_1_01 = 0.5;  // mix onto one  (2,1)
  double cc_div_01_1  = 0.0;  // one onto mix  (4,3)

  // remaining four cross entries are chains through population 1
  // constrained by diagonal entries
  double cc_chain_010 = 0.0;  // (3,1)
  double cc_chain_111 = 0.0;  // (4,1)
  double cc_chain_110 = 0.0;  // (3,2)
  double cc_chain_011 = 0.0;  // (4,2)


  // double alpha_conv[2][2][2] = {{{0.8,0.1},{99,0.8}},{{0.8,0.1},{99,0.8}}};
  // double alpha_div[2][2][2] = {{{0.8,0.8},{0.1,0.1}},{{99,99},{0.8,0.8}}};
  // //double alpha_chain[2][2][2] = {{{0.5,0.2},{0.3,0.4}},{{0.2,0.1},{0.4,0.5}}};
  // double alpha_chain[2][2][2] = {{{0.0,0.0},{0.0,0.0}},{{0.0,0.0},{0.5,0.0}}};

  double alpha_conv_hom[2][2] = {{alpha_conv_hom_00,alpha_conv_hom_01},{alpha_conv_hom_10,alpha_conv_hom_11}};
  double cc_conv_mixed[2][2][2] = {{{99,cc_conv_0_01},{99,99}},{{99,cc_conv_1_01},{99,99}}};
  double alpha_div_hom[2][2] = {{alpha_div_hom_00,alpha_div_hom_01},{alpha_div_hom_10,alpha_div_hom_11}};
  double cc_div_mixed[2][2][2] = {{{99,99},{cc_div_01_0,cc_div_01_1}},{{99,99},{99,99}}};
  double cc_chain[2][2][2] = {{{cc_chain_000,cc_chain_001},{cc_chain_010,cc_chain_011}},{{cc_chain_100,cc_chain_101},{cc_chain_110,cc_chain_111}}};


  //{{{000,001},{010,011}},{{100,101},{110,111}}}

  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);

  // set the seed
  // if deterministic_seed, use that for seed
  if(deterministic_seed)
    gsl_rng_set(r,deterministic_seed); 
  // else use time in seconds for the seed
  else
    gsl_rng_set(r, time(NULL));
  

  // matrix files will be stored in data directory
  // with filename determined by the six input parameters
  mkdir("data",0755);  // make data directory, if it doesn't exist
  char FNbase[200];

  sprintf(FNbase,"_%i_%i_%1.3f_%1.3f_%1.3f_%1.3f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f_%1.2f",
	  N_nodes[0], N_nodes[1],
	  p[0][0],p[0][1],p[1][0],p[1][1],
	  alpha_recip[0][0], alpha_recip[0][1], alpha_recip[1][1],
	  alpha_conv_hom[0][0], alpha_conv_hom[0][1],
	  alpha_conv_hom[1][0], alpha_conv_hom[1][1],
	  cc_conv_mixed[0][0][1], cc_conv_mixed[1][0][1], 
	  alpha_div_hom[0][0], alpha_div_hom[0][1], 
	  alpha_div_hom[1][0], alpha_div_hom[1][1],
	  cc_div_mixed[0][1][0], cc_div_mixed[0][1][1], 
	  cc_chain[0][0][0], cc_chain[0][0][1], 
	  cc_chain[0][1][0], cc_chain[0][1][1],
	  cc_chain[1][0][0], cc_chain[1][0][1],
	  cc_chain[1][1][0], cc_chain[1][1][1]);


  // generate second order matrix W
  gsl_matrix *W;
  W = secorder_rec_2p(N_nodes,p, alpha_recip, alpha_conv_hom, cc_conv_mixed,
		      alpha_div_hom, cc_div_mixed, cc_chain,
		      r);

  // if failed to generate a matrix, write error and quit
  if(!W) {
    cerr << "Failed to generate the matrix\n";
    return -1;
  }

  char FN[200];
  FILE *fhnd;


  ////////////////////////////////////////////////////////////
  // Calculate the covariance structure of the matrix 
  // This should approximately agree with the input alphas
  ////////////////////////////////////////////////////////////

  cout << "Testing statistics of matrix...";
  cout.flush();


  double phat[2][2];
  double alphahat_recip[2][2];
  double alphahat_conv[2][2][2];
  double alphahat_div[2][2][2];
  double alphahat_chain[2][2][2];

  calc_phat_alphahat(W, N_nodes, phat, alphahat_recip, 
		     alphahat_conv, alphahat_div, alphahat_chain);

  cout << "done\n";

  strcpy(FN, "data/stats");
  strcat(FN, FNbase);
  strcat(FN, ".dat");
  fhnd = fopen(FN, "w");
  if(fhnd==NULL) {
    cerr << "Couldn't open outfile file " << FN << "\n";
    exit(-1);
  }

  
  fprintf(fhnd, "%i %i ", N_nodes[1], N_nodes[2]);
  
  cout << "Actual statistics of matrix:\n";
  cout << "phat = ";
  for(int i=0; i<2; i++)
    for(int j=0; j<2; j++) {
      cout << phat[i][j] << " ";
      fprintf(fhnd, "%e ", phat[i][j]);
    }
  cout << "\n";
  cout << "alphahat_recip = ";
  for(int i=0; i<2; i++)
    for(int j=i; j<2; j++) {
      cout << alphahat_recip[i][j] << " ";
      fprintf(fhnd, "%e ", alphahat_recip[i][j]);
    }
  cout << "\n";
  cout << "alphahat_conv = "; 
  for(int i=0; i<2; i++)
    for(int j=0; j<2; j++)
      for(int k=j; k<2; k++) {
	cout << alphahat_conv[i][j][k] << " ";
	fprintf(fhnd, "%e ", alphahat_conv[i][j][k]);
      }
  cout << "\n";
  cout << "alphahat_div = ";
  for(int i=0; i<2; i++)
    for(int j=i; j<2; j++)
      for(int k=0; k<2; k++) {
	cout << alphahat_div[i][j][k] << " ";
	fprintf(fhnd, "%e ", alphahat_div[i][j][k]);
      }
  cout << "\n";
  cout << "alphahat_chain = ";
  for(int i=0; i<2; i++)
    for(int j=0; j<2; j++)
      for(int k=0; k<2; k++) {
	cout << alphahat_chain[i][j][k] << " ";
	fprintf(fhnd, "%e ", alphahat_chain[i][j][k]);
      }
  cout << "\n";
  cout.flush();

  fclose(fhnd);


  // balance excitation and inhibition
  int balance_average=1;

  if(balance_average) { 
    // balance so expected value of inputs is zero

    double ee_tot = p[0][0]*N_nodes[0];
    double ei_tot = p[0][1]*N_nodes[1];
    
    double ie_tot = p[1][0]*N_nodes[0];
    double ii_tot = p[1][1]*N_nodes[1];
    
    // calculate magnitudes to make the expected values balance
    //double Wee_mag = 1.0;
    double Wei_mag = -ee_tot/ei_tot;
    
    //double Wie_mag = 1.0;
    double Wii_mag = -ie_tot/ii_tot;
    
    gsl_matrix_view Wei = gsl_matrix_submatrix
      (W, 0, N_nodes[0], N_nodes[0], N_nodes[1]);
    gsl_matrix_scale(&Wei.matrix, Wei_mag);
    
    gsl_matrix_view Wii = gsl_matrix_submatrix
      (W, N_nodes[0], N_nodes[0], N_nodes[1], N_nodes[1]);
    gsl_matrix_scale(&Wii.matrix, Wii_mag);
  }
  else {
    // balance so each row sum is exactly zero

    for(int i=0; i<N_nodes[0]+N_nodes[1]; i++) {
      gsl_vector_view erow = gsl_matrix_subrow(W, i, 0, N_nodes[0]);
      gsl_vector_view irow = gsl_matrix_subrow(W, i, N_nodes[0],N_nodes[1]);
      
      // sum of excitatory inputs
      double esum = gsl_blas_dasum(&erow.vector);
      double isum = gsl_blas_dasum(&irow.vector);
      
      if(isum) {
	gsl_vector_scale(&irow.vector, -esum/isum);
      }
      else {
	gsl_vector_set_zero(&erow.vector);
	cout << "Set row " << i << " to zero\n";
	// gsl_vector_set(&irow.vector, gsl_rng_uniform_int(r,N_nodes[1]), -esum);
	// cout << "Set a element of row " << i << " to " << -esum << "\n";
      }
    }
  }

      


  // write matrix to a file
  strcpy(FN, "data/w");
  strcat(FN, FNbase);
  strcat(FN, ".dat");
  fhnd = fopen(FN, "w");
  if(fhnd==NULL) {
    cerr << "Couldn't open outfile file " << FN << "\n";
    exit(-1);
  }


  for(int i=0; i<N_nodes[0]+N_nodes[1];i++) {
    for(int j=0; j<N_nodes[0]+N_nodes[1]; j++) {
      fprintf(fhnd, "%f ", gsl_matrix_get(W,i,j));
    }
    fprintf(fhnd,"\n");
  }
  fclose(fhnd);

  gsl_matrix_free(W);


  return 0;

}