run_secorder.cpp 6.07 KB
Newer Older
1 2
#include <sys/time.h>

3
#include <iostream>
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
4
#include <stdint.h>
5 6 7 8 9 10 11
#include <string.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <gsl/gsl_rng.h>
#include "secorder_rec_1p.hpp"
#include "calc_stats_1p.hpp"

12
struct timeval T0;
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

using namespace std;

int main(int argc, char *argv[]) {
  

  ///////////////////////////////////////////////////////////////
  // Read seven input parameters 
  // N_nodes      number of nodes in the network
  // p            the probability of any one given connection
  // alpha_recip  determines covariance of reciprocal connections
  // alpha_conv   determines covariance of convergent connections
  // alpha_div    determines covariance of divergent connections
  // alpha_chain  determines covariance of chain connections
  // seed         seed for the random number generator (optional)
28
  // algorithm    S(standard), T(triangle), B16, B32, B64 (Block with size)
29 30
  ///////////////////////////////////////////////////////////////

31 32
  struct timeval t0,dT;
  gettimeofday(&T0,NULL);
33
 
34 35
  if ((argc < 7) || (argc > 9)) {
    cerr << "Requires six, seven, or eight parameters: N_nodes p alpha_recip alpha_conv alpha_div alpha_chain [seed] [algorithm]\n";
36 37 38 39 40 41 42 43 44 45 46 47 48
    exit(-1);
  }
  int N_nodes = atoi(argv[1]);
  double p = atof(argv[2]);
  double alpha_recip = atof(argv[3]);
  double alpha_conv = atof(argv[4]);
  double alpha_div = atof(argv[5]);
  double alpha_chain = atof(argv[6]);

  gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937);

  int deterministic_seed=0;
  int rng_seed = 0;
49 50 51 52 53 54 55
  int count = 7;
  if(argc > 7) {
      if (~isalpha(argv[count][0])) {
        deterministic_seed=1;
        rng_seed = atoi(argv[count]);
	count++;
      }
56
  }
57
  const char *algorithm=argv[count];
58 59 60 61 62 63 64 65 66

  // set the seed
  // if deterministic_seed, use rng_seed for seed
  if(deterministic_seed)
    gsl_rng_set(rng,rng_seed); 
  // else use time in seconds for the seed
  else
    gsl_rng_set(rng, time(NULL));
  
67
  gsl_matrix_float *W;
68 69 70 71 72

  gettimeofday(&t0,NULL);
  timersub(&t0,&T0,&dT);
  fprintf(stdout,"t=%li.%06li\t%s init (%s line %i)\n",dT.tv_sec,dT.tv_usec,__func__,__FILE__,__LINE__);

73
  W=secorder_rec_1p(N_nodes,p,  alpha_recip, alpha_conv, alpha_div, alpha_chain,
74
		    rng, algorithm);
75
  
76 77 78
  gettimeofday(&t0,NULL);
  timersub(&t0,&T0,&dT);
  fprintf(stdout,"t=%li.%06li\t%s computed (%s line %i)\n",dT.tv_sec,dT.tv_usec,__func__,__FILE__,__LINE__);
79
  fflush(stdout);
80

81 82 83 84 85 86 87 88 89 90 91
    // if failed to generate a matrix, write error and quit
  if(!W) {
    cerr << "Failed to generate the matrix\n";
    return -1;
  }

  // matrix files will be stored in data directory
  // with filename determined by the six or seven input parameters
  mkdir("data",0755);  // make data directory, if it doesn't exist
  char FNbase[200];
  if(deterministic_seed)
92
    sprintf(FNbase,"_%i_%1.3f_%1.3f_%1.3f_%1.3f_%1.3f_%i%s",N_nodes,p,alpha_recip, alpha_conv, alpha_div, alpha_chain, rng_seed, algorithm);
93
  else
94
    sprintf(FNbase,"_%i_%1.3f_%1.3f_%1.3f_%1.3f_%1.3f%s",N_nodes,p,alpha_recip, alpha_conv, alpha_div, alpha_chain, algorithm);
95
  
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
96
  char FN[2000];
97 98 99 100
  FILE *fhnd;
  strcpy(FN, "data/w");
  strcat(FN, FNbase);
  strcat(FN, ".dat");
101

102
#if 0
103 104
  // Large matrices are stored only in sparse from
  // Thise are used only for debugging
105 106 107 108 109 110 111 112
  fhnd = fopen(FN, "w");
  if(fhnd==NULL) {
    cerr << "Couldn't open outfile file " << FN << "\n";
    exit(-1);
  }

  for(int i=0; i<N_nodes;i++) {
    for(int j=0; j<N_nodes; j++) {
113
      fprintf(fhnd, "%i ", gsl_matrix_float_get(W,i,j)>1.0);
114 115 116 117 118
    }
    fprintf(fhnd,"\n");
  }
  fclose(fhnd);

119 120 121
  gettimeofday(&t0,NULL);
  timersub(&t0,&T0,&dT);
  fprintf(stdout,"t=%li.%06li\t%s output_full (%s line %i)\n",dT.tv_sec,dT.tv_usec,__func__,__FILE__,__LINE__);
122
#endif
123

Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
124 125 126 127 128 129
  ////////////////////////////////////////////////////////////
  // output of sparse matrix
  ////////////////////////////////////////////////////////////
  size_t nnz=0;
  for(int i=0; i<N_nodes; i++) {
    for(int j=0; j<N_nodes; j++) {
130
      nnz += gsl_matrix_float_get(W,i,j)>1.0;
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
131 132 133
    }
  }
  strcat(FN,".sparse");
134
  fprintf(stdout, "Write to File %s\n",FN);
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
135
  fhnd = fopen(FN, "w");
136 137 138 139
  if(fhnd==NULL) {
      cerr << "Couldn't open outfile file " << FN << "\n";
      exit(-1);
  }
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
140 141 142 143 144 145
  fwrite("SPARSE\x0\0x1",8,1,fhnd);
  size_t nr,nc;
  nr=nc=N_nodes;
  fwrite(&nr,8,1,fhnd);
  fwrite(&nc,8,1,fhnd);
  fwrite(&nnz,8,1,fhnd);
146 147 148 149 150
  for (int i=0; i<N_nodes; i++) {
    for (int j=0; j<N_nodes; j++) {
      int t = gsl_matrix_float_get(W,i,j) > 1.0;
      gsl_matrix_float_set(W,i,j,(float)t);
      if (t) {
Alois SCHLOEGL's avatar
Alois SCHLOEGL committed
151 152 153 154 155 156
        fwrite(&i, 4, 1, fhnd);
        fwrite(&j, 4, 1, fhnd);
      }
    }
  }
  fclose(fhnd);
157 158 159 160

  gettimeofday(&t0,NULL);
  timersub(&t0,&T0,&dT);
  fprintf(stdout,"t=%li.%06li\t%s output_sparse (%s line %i)\n",dT.tv_sec,dT.tv_usec,__func__,__FILE__,__LINE__);
161 162 163 164 165 166 167
  

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

168 169
if (N_nodes > 35000) return 0;

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
  cout << "Testing statistics of W ...";
  cout.flush();
  double phat, alphahat_recip, alphahat_conv, alphahat_div,
    alphahat_chain, alphahat_other;
  
  calc_phat_alphahat_1p(W, N_nodes, phat, alphahat_recip, 
			alphahat_conv, alphahat_div,
			alphahat_chain, alphahat_other);

  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, "%e %e %e %e %e %e\n", phat, alphahat_recip, 
	  alphahat_conv, alphahat_div, alphahat_chain, alphahat_other);
  fclose(fhnd);


  cout << "done\n";
  cout << "\nActual statistics of matrix:\n";
  cout << "phat = " << phat << "\n";
  cout << "alphahats:\n";
  cout << "alpha_recip = " << alphahat_recip
       << ", alpha_conv = " << alphahat_conv
       << ", alpha_div = " << alphahat_div
       << ", alpha_chain = " << alphahat_chain
       << ", alpha_other = " << alphahat_other << "\n";
  
202 203 204
  gettimeofday(&t0,NULL);
  timersub(&t0,&T0,&dT);
  fprintf(stdout,"t=%li.%06li\t%s test_alphas (%s line %i)\n",dT.tv_sec,dT.tv_usec,__func__,__FILE__,__LINE__);
205
}