genran.C

#include "random.h"
#include <cstdio>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <unistd.h>
#include <ctime>

/*
    GenRan: non-Gaussian random field generator.
    Copyright (C) 2006 Peter Grassl
    Distributed under the MIT licence.

    Command-line driver. Reads simulation parameters from an .ini-style
    input file, calls generateRandomField, writes the resulting field to
    the user-specified output file and a (mean, sdev/mean) summary to
    stat.dat in the current directory. See README.md for input file syntax.
*/

#define MYPI 3.14159265


int main(int argc, char* argv[]){
  printf("GenRan: Non-gaussian random field generator \n");
  printf("Copyright (C) 2006 Peter Grassl\n");
  printf("Version 1.0-Devel\n\n");
  
  if(argc != 3){
    printf("wrong number of arguments: input file and output file needed\n");
    std::exit(1);
  }
  
  int nDim;
  nDim = getint_ini(argv[1],"numberOfDimensions");
  
  int numberReal[nDim]; 
  if(nDim == 1){
    numberReal[0] = getint_ini(argv[1],"realNumber1");
  }
  else if(nDim ==2){
    numberReal[0] = getint_ini(argv[1],"realNumber1");
    numberReal[1] = getint_ini(argv[1],"realNumber2");
  }
  else if(nDim ==3){
    numberReal[0] = getint_ini(argv[1],"realNumber1");
    numberReal[1] = getint_ini(argv[1],"realNumber2");
    numberReal[2] = getint_ini(argv[1],"realNumber3");  
  }


  double origin[nDim];
  if(nDim == 1){
    origin[0] = getdouble_ini(argv[1],"xOrigin");
  }
  else if(nDim == 2){
    origin[0] = getdouble_ini(argv[1],"xOrigin");    
    origin[1] = getdouble_ini(argv[1],"yOrigin");    
  }
  else if(nDim == 3){
    origin[0] = getdouble_ini(argv[1],"xOrigin");    
    origin[1] = getdouble_ini(argv[1],"yOrigin");    
    origin[2] = getdouble_ini(argv[1],"zOrigin");    
  }
  
  // The radix-2 FFT requires every realNumber to be a positive power of two.
  for (int k = 0; k < nDim; ++k) {
    const int n = numberReal[k];
    if (n <= 0 || (n & (n - 1)) != 0) {
      printf("Error: realNumber%d = %d must be a positive power of two.\n", k + 1, n);
      std::exit(1);
    }
  }
  
  //autocorrelation length
  double autoLength[nDim];
  if(nDim==1){
    autoLength[0] = getdouble_ini(argv[1],"autoLength1");
  }
  else if(nDim==2){
    autoLength[0] = getdouble_ini(argv[1],"autoLength1");
    autoLength[1] = getdouble_ini(argv[1],"autoLength2");
  }
  else if(nDim==3){
    autoLength[0] = getdouble_ini(argv[1],"autoLength1");
    autoLength[1] = getdouble_ini(argv[1],"autoLength2");
    autoLength[2] = getdouble_ini(argv[1],"autoLength3");
  }
  
  //Read in converter which transforms the realNumbers to the dimensions of the random field
  double converter = getdouble_ini(argv[1], "converter");

  // Marginal CDF: 1 = Gaussian, 2 = Weibull, 3 = grafted Weibull-Gaussian.
  int typeOfCDF = getint_ini(argv[1],"typeOfCDF");

  // randomParameter slots:
  //   typeOfCDF = 1 (Gaussian):
  //     [0] mean, [1] cov
  //   typeOfCDF = 2 (Weibull):
  //     [0] scaling, [1] modulus, [2] derived mean, [3] derived variance
  //   typeOfCDF = 3 (grafted):
  //     [0] scaling, [1] modulus, [2] graft probability, [3] cov,
  //     [4] rescaling factor, [5] gaussian-core mean, [6] gaussian-core std
  double randomParameter[7];
  double gaussianMean = 0., gaussianSDV = 0., tempParameter = 0.;

  if (typeOfCDF == 1) {
    randomParameter[0] = getdouble_ini(argv[1], "mean");
    randomParameter[1] = getdouble_ini(argv[1], "cov");
  } else if (typeOfCDF == 2) {
    randomParameter[0] = getdouble_ini(argv[1], "scaling");
    randomParameter[1] = getdouble_ini(argv[1], "modulus");
    randomParameter[2] = randomParameter[0] * exp(gammaln(1. + 1. / randomParameter[1]));
    randomParameter[3] = pow(randomParameter[0], 2.)
                         * (exp(gammaln(1. + 2. / randomParameter[1]))
                            - pow(exp(gammaln(1. + 1. / randomParameter[1])), 2.));
  } else if (typeOfCDF == 3) {
    randomParameter[0] = getdouble_ini(argv[1], "scaling");
    randomParameter[1] = getdouble_ini(argv[1], "modulus");
    randomParameter[2] = getdouble_ini(argv[1], "graft");
    randomParameter[3] = getdouble_ini(argv[1], "cov");
    tempParameter = computeRescaling(randomParameter, gaussianMean, gaussianSDV);
    randomParameter[4] = tempParameter;
    randomParameter[5] = gaussianMean;
    randomParameter[6] = gaussianSDV;
  } else {
    std::fprintf(stderr, "Error: unknown typeOfCDF %d (expected 1, 2, or 3)\n", typeOfCDF);
    std::exit(1);
  }

  long randomInteger;
  int seedFromFile = 0;
  const bool seedProvided = tryGetint_ini(argv[1], "ranint", &seedFromFile);
  if (!seedProvided) {
    randomInteger = -time(NULL);
  } else if (seedFromFile >= 0) {
    std::cout << "Random integer must be smaller than zero. Use time based random integer instead\n";
    sleep(2);
    randomInteger = -time(NULL);
  } else {
    randomInteger = seedFromFile;
  }
  

  int numberTotal = numberReal[0];
  for(int i=1;i<nDim;i++){
    numberTotal = numberTotal* numberReal[i];
  }
  
  // The simulator works in scaled spectral coordinates with a fixed
  // upper-cutoff wavenumber kappaUpper = 1.8π. The corresponding spatial
  // step deltaX = π/kappaUpper is the same in every dimension. The user-
  // facing `converter` (samples per unit length) is rescaled to scaled
  // units so output coordinates map back to user units cleanly, and the
  // autocorrelation length is converted from the user's spatial units to
  // the target spectrum's b parameter (b = 2·la/√π).
  double kappaUpper = 1.8 * MYPI;
  double deltaX[nDim];
  for (int t = 0; t < nDim; ++t) deltaX[t] = MYPI / kappaUpper;

  converter = converter * deltaX[0];
  for (int i = 0; i < nDim; ++i) {
    autoLength[i] = autoLength[i] * (2. / sqrt(MYPI)) * converter;
  }
  
  // Optional padding: simulate on a `padding`x larger grid in each dimension
  // and crop the user's requested region. Suppresses boundary contamination
  // from the FFT's implicit periodicity. Must be a power of two so that the
  // padded grid stays a power of two (required by the radix-2 FFT).
  int padding = 1;
  int paddingFromFile = 0;
  if (tryGetint_ini(argv[1], "padding", &paddingFromFile)) {
    padding = paddingFromFile;
    if (padding < 1 || (padding & (padding - 1)) != 0) {
      printf("Error: padding must be a positive power of two (1, 2, 4, ...). Got %d.\n", padding);
      std::exit(1);
    }
  }

  // Optional SDB-2011 iterative determination of the underlying Gaussian
  // PSDF for non-Gaussian fields. When iterate == 0 (the default) the user's
  // target spectrum is used directly as the Gaussian spectrum, which is the
  // legacy behaviour and is correct for typeOfCDF == 1.
  int iterate = 0;
  int iterTmp = 0;
  if (tryGetint_ini(argv[1], "iterate", &iterTmp)) {
    iterate = iterTmp;
  }
  int iterMaxIter = 100;
  if (tryGetint_ini(argv[1], "iterMaxIter", &iterTmp)) iterMaxIter = iterTmp;
  int iterHermiteOrder = 30;
  if (tryGetint_ini(argv[1], "iterHermiteOrder", &iterTmp)) iterHermiteOrder = iterTmp;
  double iterTolerance = 0.01;
  tryGetdouble_ini(argv[1], "iterTolerance", &iterTolerance);
  double iterBeta = 1.4;
  tryGetdouble_ini(argv[1], "iterBeta", &iterBeta);

  double* field = new double[numberTotal];

  generateRandomField(field, numberReal, autoLength, randomParameter, nDim,
                      typeOfCDF, randomInteger, padding, iterate, iterMaxIter,
                      iterTolerance, iterBeta, iterHermiteOrder);

  printf("generation of random field is done\n");

  // Summary statistics → stat.dat: mean and coefficient of variation.
  // moment() uses the legacy 1-based indexing convention; offset the
  // pointer locally and restore before delete[].
  double ave, adev, sdev, var, skew, curt;
  field--;
  moment(field, numberTotal, &ave, &adev, &sdev, &var, &skew, &curt);
  field++;
  std::ofstream outputStat("stat.dat");
  outputStat << ave << "\t" << sdev / ave << "\n \n";
  outputStat.close();

  //write the field into the output file
  std::ofstream outputField(argv[2]); 
    
  if(nDim==1){
    for(int k=0;k<numberReal[0];k++){            
      outputField << origin[0] + k/converter*deltaX[0] << "\t" << field[k] << "\n";
    } 
  }
  else if(nDim==2){
    outputField << numberReal[0] << "\t" << numberReal[1] << "\t" << deltaX[0] << "\t" << deltaX[1] <<"\n";
    for(int k=0;k<numberReal[0];k++){       
      for(int l=0;l<numberReal[1];l++){
        outputField << origin[0] + k/converter*deltaX[0] << "\t" << origin[1] + l/converter*deltaX[1] << "\t" <<  field[k*numberReal[1]+l] << "\n";
      }
      outputField << "\n";
    }
  }
  else if(nDim ==3){
    outputField << numberReal[0] << "\t" << numberReal[1]  << "\t" << numberReal[2] << "\t" << deltaX[0] << "\t" << deltaX[1] << "\t" << deltaX[2] <<"\n";
    for(int k=0;k<numberReal[0];k++){  
      for(int l=0;l<numberReal[1];l++){
        for(int m=0;m<numberReal[2];m++){     
          outputField << origin[0] + k/converter*deltaX[0] << "\t" << origin[1] + l/converter*deltaX[1] << "\t" << origin[2] + m/converter*deltaX[2] << "\t" << field[k*numberReal[1]*numberReal[2]+l*numberReal[2]+m] << "\n";
        }
      }
    }
  }
  
  outputField.close();
  delete[] field;
  return 0;
}