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

3
#include <iostream>
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
  
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

124 125 126 127
  ////////////////////////////////////////////////////////////
  // output of sparse matrix
  ////////////////////////////////////////////////////////////
  size_t nnz=0;
128
  uint32_t *rs = (uint32_t*)malloc(N_nodes*4);
129
  for(int i=0; i<N_nodes; i++) {
130
    uint32_t n=0;
131
    for(int j=0; j<N_nodes; j++) {
132
      n += gsl_matrix_float_get(W,i,j)>1.0;
133
    }
134 135
    nnz+=n;
    rs[i]=n;
136 137
  }
  strcat(FN,".sparse");
138
  fprintf(stdout, "Write to File %s\n",FN);
139
  fhnd = fopen(FN, "w");
140 141 142 143
  if(fhnd==NULL) {
      cerr << "Couldn't open outfile file " << FN << "\n";
      exit(-1);
  }
144
#if defined(SPARSE_COO)
145
  fwrite("SPARSE\x00\x00",8,1,fhnd);
146 147 148 149 150
  size_t nr,nc;
  nr=nc=N_nodes;
  fwrite(&nr,8,1,fhnd);
  fwrite(&nc,8,1,fhnd);
  fwrite(&nnz,8,1,fhnd);
151 152 153 154 155
  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) {
156 157 158 159 160
        fwrite(&i, 4, 1, fhnd);
        fwrite(&j, 4, 1, fhnd);
      }
    }
  }
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
#elif defined(SPARSE_MM) 	// MatrixMarket format
  fprintf(fhnd,"%%MatrixMarket matrix coordinate binary general\n");
  fprintf(fhnd,"%%\n%%   see also http://math.nist.gov/MatrixMarket/formats.html#MMformat\n%%\n");
  size_t nr,nc;
  nr=nc=N_nodes;
  fprintf(fhnd,"%lu\t%lu\t%lu\n",nr,nc,nnz);
  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) {
        fprintf(fhnd,"%lu\t%lu\n",i,j);
      }
    }
  }
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
#elif 1 //defined(SPARSE_CSR)
  fwrite("SPARSE\x00\x02",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);
  fwrite(rs, 4, N_nodes, fhnd);
  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) {
        fwrite(&j, 4, 1, fhnd);
      }
    }
  }
#endif
194
  fclose(fhnd);
195
  if (rs != NULL) free(rs);
196 197 198 199

  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__);
200
  fflush(stdout);
201 202 203 204 205 206

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

207 208
if (N_nodes > 35000) return 0;

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
  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";
  
241 242 243
  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__);
244
}