text:
/* #define _XOPEN_SOURCE 500 */
#ifndef _FILE_OFFSET_BITS
#define _FILE_OFFSET_BITS 64
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <inttypes.h>
#include <limits.h>

#ifndef SIZE_MAX
#define SIZE_MAX (~(size_t)0)
#endif

#ifdef _OPENMP
#include <omp.h>
#endif

/* for open/close and other file io*/
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

/*for sendfile (in-kernel file to file copy, supposedly very fast) */
#ifdef USE_SENDFILE
#include <sys/sendfile.h>
#endif

#if defined(USE_MMAP)
#include <sys/mman.h>
#elif defined(USE_SENDFILE)
#include <sys/sendfile.h>
#else
#endif

#ifdef USE_WRITEV
#ifndef USE_MMAP
#error USE_MMAP must be enabled with USE_WRITEV
#endif
#include <sys/uio.h>
#endif

#include <gsl/gsl_rng.h>

#include "macros.h"
#include "utils.h"
#include "progressbar.h"
#include "gadget_utils.h"

/* Copied straight from https://fossies.org/dox/gsl-2.2.1/shuffle_8c_source.html*/
/* Adapted to generate array indices. The returned random indices are in increasing order */
int gsl_ran_arr_index (const gsl_rng * r, size_t * dest, const size_t k, const size_t n)
{
  /* Choose k out of n items, return an array x[] of the k items.
      These items will preserve the relative order of the original
      input -- you can use shuffle() to randomize the output if you
      wish */

  if (k > n) {
	fprintf(stderr,"k=%zu is greater than n=%zu. Cannot sample more than n items\n",
			k, n);
	return EXIT_FAILURE;
  }
  if(k == n) {
	for(size_t i=0;i<n;i++) {
	  dest[i] = i;
	}
  } else {
	size_t j=0;
	for (size_t i = 0; i < n && j < k; i++) {
	  if ((n - i) * gsl_rng_uniform (r) < k - j) {
		dest[j] = i;
		j++ ;
	  }
	}
  }

  return EXIT_SUCCESS;
}

int write_random_subsample_of_field(int in_fd, int out_fd, off_t in_offset, off_t out_offset, const int dest_npart, const size_t itemsize, size_t *random_indices
#ifdef USE_MMAP
									,char *in_memblock
#endif
)
{
#ifdef USE_MMAP
  in_memblock += in_offset;//These are here so that I don't lose the original pointer from mmap
#endif

#ifndef USE_MMAP
  size_t buf[itemsize];
#endif

  int interrupted=0;

#ifndef _OPENMP
  init_my_progressbar(dest_npart,&interrupted);
#endif

#ifdef USE_WRITEV
  //vector writes. Loop has to be re-written in units of IOV_MAX
  for(int i=0;i<dest_npart;i+=IOV_MAX) {
#ifndef _OPENMP
      my_progressbar(i,&interrupted);//progress-bar will be jumpy
#endif
      const int nleft = ((dest_npart - i) > IOV_MAX) ? IOV_MAX:(dest_npart - i);
      struct iovec iov[IOV_MAX];
      for(int j=0;j<nleft;j++){
          const size_t ind = random_indices[j];
          const size_t offset_in_field = ind * itemsize;
          iov[j].iov_base = in_memblock + offset_in_field;
          iov[j].iov_len = itemsize;
      }
      random_indices += nleft;

      ssize_t bytes_written = writev(out_fd, iov, nleft);
      XRETURN(bytes_written == nleft*itemsize, EXIT_FAILURE, "Expected to write bytes = %zu but wrote %zd instead\n",nleft*itemsize, bytes_written);
  }
#else  //Not using WRITEV -> 3 options here i) MMAP + write ii) sendfile iii) pread + write
  for(int i=0;i<dest_npart;i++) {
#ifndef _OPENMP
	my_progressbar(i,&interrupted);
#endif
	const size_t ind = random_indices[i];
	size_t offset_in_field = ind * itemsize;
#if defined(USE_MMAP)
	ssize_t bytes_written = write(out_fd, in_memblock + offset_in_field, itemsize);
#elif defined(USE_SENDFILE)
	off_t input_offset = in_offset + offset_in_field;
	ssize_t bytes_written = sendfile(out_fd, in_fd, &input_offset, itemsize);
#else
	ssize_t bytes_read = pread(in_fd, &buf, itemsize, in_offset + offset_in_field);
	XRETURN(bytes_read == itemsize, EXIT_FAILURE, "Expected to read bytes = %zu but read %zd instead\n",itemsize, bytes_read);
	ssize_t bytes_written = write(out_fd, &buf, itemsize);
#endif//default case (pread + write)

	XRETURN(bytes_written == itemsize, EXIT_FAILURE, "Expected to write bytes = %zu but wrote %zd instead\n",itemsize, bytes_written);
	out_offset += itemsize;
  }
#endif //end of WRITEV

#ifndef _OPENMP
  finish_myprogressbar(&interrupted);
#endif

  return EXIT_SUCCESS;
}

int subsample_single_gadgetfile(const int dest_npart, const char *inputfile, const char *outputfile, const size_t id_bytes, const gsl_rng *rng, const double fraction, const int64_t nparttotal)
{
  if(dest_npart <= 0) {
	fprintf(stderr,"Error: Desired number of particles =%d in the subsampled file must be > 0\n",dest_npart);
	return EXIT_FAILURE;
  }
  struct timespec tstart, t0, t1;
  current_utc_time(&tstart);

  int in_fd=-1;
  int out_fd=-1;
  int status = EXIT_FAILURE;
  in_fd = open(inputfile, O_RDONLY);
  if(in_fd < 0) {
	fprintf(stderr,"Error (in function %s, line # %d) while opening input file = `%s'\n",__FUNCTION__,__LINE__,inputfile);
	perror(NULL);
	return EXIT_FAILURE;
  }

#ifdef USE_MMAP
  struct stat sb;
  if(fstat(in_fd, &sb) < 0) {
	perror(NULL);
	return EXIT_FAILURE;
  }
  char *in_memblock = mmap(NULL, sb.st_size, PROT_READ, MAP_SHARED, in_fd, 0);
#endif

  //Check that that the output file does not exist.
  {
	FILE *fp = fopen(outputfile,"r");
	if(fp != NULL) {
	  fclose(fp);
	  fprintf(stderr,"Warning: Output file = `%s' should not exist. "
			  "Aborting so as to avoid accidentally over-writing regular fles\n",outputfile);
	  return EXIT_FAILURE;
	}
  }

  out_fd = open(outputfile, O_CREAT | O_WRONLY | O_TRUNC, S_IRUSR | S_IWUSR);//set the mode since the file is being created
  if(out_fd < 0) {
	fprintf(stderr,"Error (in function %s, line # %d) while opening output file = `%s'\n",__FUNCTION__,__LINE__, outputfile);
	perror(NULL);
	return EXIT_FAILURE;
  }

  /* Reserve disk-space for dest_npart particles, written as a fortran binary */
  const size_t pos_vel_itemsize = sizeof(float)*3;
  const off_t header_disk_size = 4 + sizeof(struct io_header) + 4;  //header
  const off_t pos_disk_size = 4 + pos_vel_itemsize*dest_npart + 4; //all 3 positions for each particle
  const off_t vel_disk_size = 4 + pos_vel_itemsize*dest_npart + 4; //all 3 velocities for each particle
  const off_t id_disk_size = 	4 + id_bytes*dest_npart + 4;//all particle IDs

  /* These are all of the fields to be written */
  const off_t outputfile_size = header_disk_size + pos_disk_size + vel_disk_size + id_disk_size;

  status = posix_fallocate(out_fd, 0, outputfile_size);
  if(status < 0) {
  	fprintf(stderr,"Error: Could not reserve disk space for %d particles (output file: `%s', expected file-size on disk = %zu bytes)\n",
  			dest_npart, outputfile, (size_t) outputfile_size);
  	perror(NULL);
  	return status;
  }

  struct io_header hdr;
  int dummy1=0,dummy2=0;
  const off_t pos_start_offset = header_disk_size + 4;
  const off_t vel_start_offset = header_disk_size + pos_disk_size + 4;
  const off_t id_start_offset = header_disk_size + pos_disk_size + vel_disk_size + 4;

  read(in_fd, &dummy1, 4);
  read(in_fd, &hdr, sizeof(struct io_header));
  read(in_fd, &dummy2, 4);
  if(dummy1 != 256 || dummy2 != 256) {
	fprintf(stderr,"Error: Padding bytes for header = %d (front) and %d (end) should be exactly 256\n",dummy1, dummy2);
	return EXIT_FAILURE;
  }
#ifdef USE_MMAP
  if(fraction == 1.0) {
	XRETURN(dest_npart == hdr.npart[1],EXIT_FAILURE,
			"for fraction = 1.0, input npart = %d must equal subsampled npart = %d\n",hdr.npart[1], dest_npart);
	/* XRETURN(sb.st_size == outputfile_size, EXIT_FAILURE, */
	/* 		"for fraction = 1.0, input (=%zu bytes) and output file sizes (=%zu bytes) must be the same",sb.st_size, outputfile_size); */
  }
#endif

  for(int type=0;type<6;type++) {
	if(type==1) {
	  continue;
	}

	/* There should only be dark-matter particles */
	if(hdr.npart[type] > 0) {
	  fprintf(stderr,"Error: Input file `%s' contains %d particles of type=%d. This code only works for dark-matter only simulations\n",
			  inputfile, hdr.npart[type], type);
	  return EXIT_FAILURE;
	}
  }

  /* The number of subsampled particles wanted can at most be the numbers present in the file*/
  if(dest_npart > hdr.npart[1]) {
	fprintf(stderr,"Error: Number of subsampled particles = %d exceeds the number of particles in the file = %d\n",
			dest_npart,  hdr.npart[1]);
	return EXIT_FAILURE;
  }

  struct io_header out_hdr = hdr;
  out_hdr.npart[1] = dest_npart;
  out_hdr.npartTotal[1] = nparttotal;
  out_hdr.npartTotalHighWord[1] = (nparttotal >> 32);
  out_hdr.mass[1] /= fraction;

  write(out_fd, &dummy1, sizeof(dummy1));
  write(out_fd, &out_hdr, sizeof(struct io_header));
  write(out_fd, &dummy2, sizeof(dummy2));

  //create an array of random indices
  size_t *random_indices = calloc(dest_npart, sizeof(*random_indices));
  status = gsl_ran_arr_index(rng, random_indices, (size_t) dest_npart, (size_t) hdr.npart[1]);

  //write positions
  current_utc_time(&t0);
  int posvel_disk_size=pos_vel_itemsize*dest_npart;
  write(out_fd, &posvel_disk_size, sizeof(posvel_disk_size));
  status = write_random_subsample_of_field(in_fd, out_fd, pos_start_offset, -1,dest_npart, pos_vel_itemsize, random_indices
#ifdef USE_MMAP
											   ,in_memblock
#endif
											   );
  if(status != EXIT_SUCCESS) {
	return status;
  }
  write(out_fd, &posvel_disk_size, sizeof(posvel_disk_size));
  current_utc_time(&t1);
  double pos_time = REALTIME_ELAPSED_NS(t0,t1)*1e-9;

  //write velocities
  write(out_fd, &posvel_disk_size, sizeof(posvel_disk_size));
  status = write_random_subsample_of_field(in_fd, out_fd, vel_start_offset, -1,dest_npart, pos_vel_itemsize, random_indices
#ifdef USE_MMAP
											   ,in_memblock
#endif
											   );
  if(status != EXIT_SUCCESS) {
	return status;
  }
  write(out_fd, &posvel_disk_size, sizeof(posvel_disk_size));
  current_utc_time(&t0);
  double vel_time = REALTIME_ELAPSED_NS(t1, t0)*1e-9;

  //write ids
  int id_size = id_bytes * dest_npart;
  write(out_fd, &id_size, sizeof(id_size));
  status = write_random_subsample_of_field(in_fd, out_fd, vel_start_offset, -1,dest_npart, id_bytes, random_indices
#ifdef USE_MMAP
                                           ,in_memblock
#endif
                                           );
  if(status != EXIT_SUCCESS) {
	return status;
  }
  write(out_fd, &id_size, sizeof(id_size));
  current_utc_time(&t1);
  double id_time = REALTIME_ELAPSED_NS(t0,t1)*1e-9;

  free(random_indices);

  //close the input file -> we are only reading, unlikely to be error
  close(in_fd);

#ifdef USE_MMAP
  munmap(in_memblock, sb.st_size);
#endif

  //check for error code here since disk quota might be hit
  status = close(out_fd);
  if(status != EXIT_SUCCESS){
	fprintf(stderr,"Error while closing output file = `%s'\n",outputfile);
	perror(NULL);
	return status;
  }
  current_utc_time(&t1);
  /* fprintf(stderr,"Done with file. Total time taken = %8.4lf seconds (pos_time = %6.3e vel_time = %6.3e id_time = %6.3e seconds)\n",REALTIME_ELAPSED_NS(tstart,t1)*1e-9, */
  /*   	  pos_time,vel_time,id_time); */

  return EXIT_SUCCESS;
}

int main(int argc,char **argv)
{
  const char argnames[][100]={"fraction","input file","output file"};
  int nargs=sizeof(argnames)/(sizeof(char)*100);
  char input_filename[MAXLEN];
  char output_filename[MAXLEN];
  struct timespec tstart,t0,t1;
  int64_t nparticles_written=0,nparticles_withmass=0;
  const gsl_rng_type * rng_type = gsl_rng_ranlxd1;
  gsl_rng * all_procs_rng = gsl_rng_alloc(rng_type);
  unsigned long seed = 42;
  int64_t TotNumPart;
  current_utc_time(&tstart);

  if (argc != 4 )  {
	fprintf(stderr,"ERROR: %s usage - <fraction>  <gadget snapshot name>  <output filename>\n",argv[0]);
	fprintf(stderr,"Each file will be subsampled to get (roughly) that fraction for each particle-type\n");
    fprintf(stderr,"\nFound: %d parameters\n ",argc-1);
	int i;
    for(i=1;i<argc;i++) {
      if(i <= nargs)
		fprintf(stderr,"\t\t %s = `%s' \n",argnames[i-1],argv[i]);
      else
		fprintf(stderr,"\t\t <> = `%s' \n",argv[i]);
    }
    if(i <= nargs) {
      fprintf(stderr,"\nMissing required parameters: \n");
      for(i=argc;i<=nargs;i++)
		fprintf(stderr,"\t\t %20s = `?'\n",argnames[i-1]);
    }
	fprintf(stderr,"\n\n");
    exit(EXIT_FAILURE);
  }

  double fraction = atof(argv[1]);
  XRETURN(fraction > 0 && fraction <= 1.0, EXIT_FAILURE, "Subsample fraction = %lf needs to be in (0,1]", fraction);
  strncpy(input_filename,argv[2],MAXLEN);
  strncpy(output_filename,argv[3],MAXLEN);
  XRETURN(strncmp(input_filename,output_filename,MAXLEN) != 0, EXIT_FAILURE, "Input filename = `%s' and output filename = `%s' are the same",input_filename, output_filename);
  const int nfiles = get_gadget_nfiles(input_filename);
  struct io_header header = get_gadget_header(input_filename);
  TotNumPart = get_Numpart(&header);
  XRETURN(header.npartTotal[0] == 0 && header.npartTotalHighWord[0]  == 0, EXIT_FAILURE, "Subsampling will not work with gas particles");

  fprintf(stderr,"Running `%s' on %d files with the following parameters \n",argv[0],nfiles);
  fprintf(stderr,"\n\t\t ---------------------------------------------\n");
  for(int i=1;i<=nargs;i++) {
	fprintf(stderr,"\t\t %-25s = %s \n",argnames[i-1],argv[i]);
  }
#ifdef _OPENMP
#pragma omp parallel
  {
      int nthreads = omp_get_num_threads();
      int tid = omp_get_thread_num();
      if(tid == 0)
          fprintf(stderr,"\t\t %-25s = %d \n","nthreads", nthreads);
  }
#endif
  fprintf(stderr,"\t\t ---------------------------------------------\n\n");

  if(SIZE_MAX != 0xFFFFFFFFFFFFFFFF) {
      fprintf(stderr,"Error: SIZE_MAX=%zu should be (2^64-1)\n", SIZE_MAX);
      return EXIT_FAILURE;
  }

  size_t id_bytes;
  int interrupted=0;
  size_t seedtable[nfiles];
  int64_t nparttotal=0;
  fprintf(stderr,"Checking all input files ...\n");
  init_my_progressbar(nfiles, &interrupted);
  for(int ifile=0;ifile<nfiles;ifile++) {
      seedtable[ifile] = SIZE_MAX * gsl_rng_uniform(all_procs_rng);
      char inputfile[MAXLEN];
      my_snprintf(inputfile,MAXLEN,"%s.%d",input_filename, ifile);
      if(ifile == 0) {
          id_bytes = get_gadget_id_bytes(inputfile);
          fprintf(stderr,"Gadget ID bytes = %zu\n",id_bytes);
          interrupted=1;
      }
      struct io_header hdr = get_gadget_header(inputfile);
      const int dest_npart = fraction * hdr.npart[1];
      XRETURN((int) (dest_npart*3*sizeof(float))  < INT_MAX, EXIT_FAILURE,
              "Padding bytes will overflow, please reduce the value of fraction (currently, fraction = %lf)\n",fraction);
      my_progressbar(ifile,&interrupted);
      nparttotal += dest_npart;
  }
  finish_myprogressbar(&interrupted);
  fprintf(stderr,"Checking all input files .....done\n\n");

  int numdone=0, errorflag=0, savestatus=0;

  init_my_progressbar(nfiles, &interrupted);
#ifdef _OPENMP
#pragma omp parallel shared(numdone, savestatus)
  {
      int nthreads = omp_get_num_threads();
      int tid = omp_get_thread_num();
#pragma omp for schedule(dynamic)
#endif
      for(int ifile=0;ifile<nfiles;ifile++) {
          if(errorflag == 0) {
              //Setup the rng
              gsl_rng * rng = gsl_rng_alloc(rng_type);
              gsl_rng_set(rng, seedtable[ifile]);

#ifdef _OPENMP
#pragma omp atomic
              numdone++;

              if(tid == 0)
#endif
                  my_progressbar(numdone,&interrupted);

              char inputfile[MAXLEN],outputfile[MAXLEN];
              my_snprintf(inputfile,MAXLEN,"%s.%d",input_filename,ifile);
              my_snprintf(outputfile, MAXLEN,"%s.%d",output_filename,ifile);
              struct io_header hdr = get_gadget_header(inputfile);
              const int dest_npart = fraction * hdr.npart[1];
              int status = subsample_single_gadgetfile(dest_npart, inputfile, outputfile, id_bytes, rng, fraction, nparttotal);
              if(status != EXIT_SUCCESS) {
                  savestatus = status;
                  errorflag = 1;
              }
              gsl_rng_free(rng);

              /* #ifdef _OPENMP           */
              /* #pragma omp critical(info) */
              /*               { */
              /* #endif */
              /*                   fprintf(stderr,"Wrote %d subsampled particles out of %d particles (%d/%d files).\n", */
              /*                           dest_npart, hdr.npart[1], ifile+1,nfiles); */

              /*                   interrupted=1; */
              /* #ifdef _OPENMP */
              /*               }//critical region */
              /* #endif */

          }//errorflag == 0 condition
      }//for loop over nfiles

#ifdef _OPENMP
  }//omp parallel region
#endif

  gsl_rng_free(all_procs_rng);
  finish_myprogressbar(&interrupted);

  if(errorflag != 0) {
      return savestatus;
  }

  current_utc_time(&t1);
  fprintf(stderr,"subsample_Gadget> Done. Wrote %"PRId64" particles to file `%s'. Time taken = %6.2lf mins\n",
		  nparttotal,output_filename,REALTIME_ELAPSED_NS(tstart, t1)*1e-9/60.0);

  return EXIT_SUCCESS;
}

max_stars_repo_name: manodeep/subsampleGadget

max_issues_repo_name: manodeep/subsampleGadget

max_forks_repo_name: manodeep/subsampleGadget


text:
#pragma once

#include "bits_of_matcha/engine/memory/Buffer.h"
#include "bits_of_matcha/engine/iterations/MatrixwiseBinaryCtx.h"

#include <cblas.h>

namespace matcha::engine::cpu {

void mm(engine::Buffer* a, engine::Buffer* b, engine::Buffer* c, const MatrixwiseBinaryCtx& ctx) {
  size_t sizeA = ctx.rowsA * ctx.colsA;
  size_t sizeB = ctx.rowsB * ctx.colsB;
  size_t sizeC = ctx.rowsA * ctx.colsB;

  dynamic_cast<Buffer*>(a);
  dynamic_cast<Buffer*>(b);
  dynamic_cast<Buffer*>(c);

  auto valsA = a->as<float*>();
  auto valsB = b->as<float*>();
  auto valsC = c->as<float*>();

  int strides = (int) ctx.prefixStridesA.size();

  std::vector<unsigned> beginA(strides, 0);
  std::vector<unsigned> beginB(strides, 0);
  std::vector<unsigned> beginC(strides, 0);
  auto iterC = valsC;

  int counter = 0;
  int axis = 1;
  while (true) {
//    if (counter++ > 20) exit(69);
    if (axis < strides - 1) {
      print("==========");
      print("axis ", axis, " beginC ", beginC[axis], " stridesC ", ctx.prefixStridesC[axis - 1]);
//    print("axis: ", axis, " (strides ", strides, ")");
//      print("strides: ", ctx.prefixStridesA[axis], " ", ctx.prefixStridesB[axis], " -> ", ctx.prefixStridesC[axis]);
//      print("begins: ", beginA[axis], " ", beginB[axis], " -> ", beginC[axis]);
//      print("BC: ", beginC[axis], " ", beginC[axis - 1] + ctx.prefixStridesC[axis - 1]);
//      print("beginC ", beginC[axis]);
      if (beginC[axis] == beginC[axis - 1] + ctx.prefixStridesC[axis - 1]) {
//        print("-> dec axis");
        if (axis == 1) break;
        axis--;
      } else {
        beginA[axis + 1] = beginA[axis];
        beginB[axis + 1] = beginB[axis];
        beginC[axis + 1] = beginC[axis];
        beginA[axis] += ctx.prefixStridesA[axis];
        beginB[axis] += ctx.prefixStridesB[axis];
        beginC[axis] += ctx.prefixStridesC[axis];
//        print("-> inc axis");
        axis++;
      }
    } else {
      if (axis >= strides) {
        axis = strides - 1;
      }

//      print(ctx.prefixStridesC[axis - 1]);
//      print(beginA[axis]);
      auto matA = valsA + beginA[axis] * sizeA;
      auto matB = valsB + beginB[axis] * sizeB;
      auto loops = axis - 1 >= 0 ? ctx.prefixStridesC[axis - 1] : 1;
      for (auto itC = 0; itC != loops; itC++) {
        auto matC = valsC + (beginC[axis] + itC) * sizeC;

//        print(matA - valsA, " ", matB - valsB, " -> ", matC - valsC);
///*
        cblas_sgemm(
          CblasRowMajor,
          CblasNoTrans,
          CblasNoTrans,
          (int) ctx.rowsA,
          (int) ctx.colsB,
          (int) ctx.colsA,
          1,
          matA,
          (int) ctx.colsA,
          matB,
          (int) ctx.colsB,
          0,
          matC,
          (int) ctx.colsB
        );
//        */

        matA += sizeA * ctx.prefixStridesA[axis];
        matB += sizeB * ctx.prefixStridesB[axis];
      }
//      print("blas");

      if (axis <= 1) break;
      axis--;
    }

  }

}

}

ext: h

max_stars_repo_name: matcha-ai/matcha

max_issues_repo_name: matcha-ai/matcha

max_forks_repo_name: matcha-ai/matcha


text:
/* Pre-requisites for the package regression */
#include "config.h"

#ifdef HAVE_GSL_HEADER

#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_vector.h>

#include "R.h"

#if _CUDA_
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cublas.h>
#include <cula_blas.h>
#include <cula.h>
#endif

#if _CUDA_
	#define PREC float
	#define PREC_EPS FLT_EPSILON
	#define GSL_TYPE(type) gsl_ ## type ## _float
	#define	GSL_FUNCTION(type,name) gsl_ ## type ## _float_ ## name
	#define GSL_STATS_FUNCTION(name) gsl_stats ## _float_ ## name
	#define MATHS_FUNCTION(name) name ## f
	#define SVD_FUNCTION svdAnyMatCuda
	#define BLAS_FUNCTION(name) my_cula_s ## name
	#define MY_FUNCTION(name) my_cula_ ## name
	#define UGAUSSIAN_FUNCTION my_ugaussian_function
	#define GSL_BLAS_FUNCTION(name) gsl_blas_s ## name
	#define PREPARE_FUNCTION(name) prepare ## name ## ForCoordinateDescentCuda
	#define PREC_DIFF -4
#else
	#define PREC double
	#define PREC_EPS DBL_EPSILON
	#define GSL_TYPE(type) gsl_ ## type
	#define	GSL_FUNCTION(type,name) gsl_ ## type ## _ ## name
	#define GSL_STATS_FUNCTION(name) gsl_stats ## _ ## name
	#define MATHS_FUNCTION(name) name
	#define SVD_FUNCTION svdAnyMat
	#define BLAS_FUNCTION(name) gsl_blas_d ## name
	#define MY_FUNCTION(name) my_gsl_ ## name
	#define UGAUSSIAN_FUNCTION gsl_cdf_ugaussian_P
	#define GSL_BLAS_FUNCTION(name) gsl_blas_d ## name
	#define PREPARE_FUNCTION(name) prepare ## name ## ForCoordinateDescent
	#define PREC_DIFF -6
#endif

#endif

max_stars_repo_name: dfrankow/ridge-cran

max_issues_repo_name: dfrankow/ridge-cran

max_forks_repo_name: dfrankow/ridge-cran


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 * @brief scale codec
 * @file Scale.h
 */
#pragma once
#include "Common.h"
#include "ScaleDecoderStream.h"
#include "ScaleEncoderStream.h"
#include <boost/system/system_error.hpp>
#include <boost/throw_exception.hpp>
#include <gsl/span>
#include <vector>

namespace bcos
{
namespace codec
{
namespace scale
{
/**
 * @brief convenience function for encoding primitives data to stream
 * @tparam Args primitive types to be encoded
 * @param args data to encode
 * @return encoded data
 */
template <typename... Args>
void encode(std::shared_ptr<bytes> _encodeData, Args&&... _args)
{
    ScaleEncoderStream s;
    (s << ... << std::forward<Args>(_args));
    *_encodeData = s.data();
}

template <typename... Args>
bytes encode(Args&&... _args)
{
    ScaleEncoderStream s;
    (s << ... << std::forward<Args>(_args));
    return s.data();
}

/**
 * @brief convenience function for decoding primitives data from stream
 * @tparam T primitive type that is decoded from provided span
 * @param span of bytes with encoded data
 * @return decoded T
 */
template <class T>
void decode(T& _decodedObject, gsl::span<byte const> _span)
{
    ScaleDecoderStream s(_span);
    s >> _decodedObject;
}

template <class T>
T decode(gsl::span<byte const> _span)
{
    T t;
    decode(t, _span);
    return t;
}
}  // namespace scale
}  // namespace codec
}  // namespace bcos


text:
/**
 *
 * @file core_dsetvar.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 * CORE_dsetvar sets a single variable, x := alpha.
 *
 *******************************************************************************
 *
 * @param[in] alpha
 *         Scalar to set x to, passed by pointer so it can depend on runtime value.
 *
 * @param[out] x
 *         On exit, x = alpha.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dsetvar = PCORE_dsetvar
#define CORE_dsetvar PCORE_dsetvar
#endif
void CORE_dsetvar(const double *alpha, double *x)
{
    *x = *alpha;
}

ext: c


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

int main()
{

    //set up some data
    int n=300;
    float* x = malloc( n*sizeof(float) );
    float* y = malloc( n*sizeof(float) );
    for( int i=0; i<n; ++i )
    {
        x[i]=i+1;
        y[i]=(i+1)*(i+1);
    }

    //calculate y = a*x + y
    float a=7.0f;
    int incx=1;
    int incy=1;
    cblas_saxpy( n, a, x, incx, y, incy );

    //check it
    int status = 0;
    for( int i=0; i<n; ++i )
    {
        float expected = 7.0f*(i+1) + (i+1)*(i+1);
        if( fabs( y[i]/expected - 1 )>1e-4 ){
            printf("error matching element i=%d expected=%g received=%g\n",i+1,expected,y[i]);
            status=-1;
            break;
        }
    }

    //free data
    free(x);
    free(y);

    return status;

}

max_stars_repo_name: murraypurves/BootsOnTheGround

max_issues_repo_name: murraypurves/BootsOnTheGround

max_forks_repo_name: murraypurves/BootsOnTheGround


text:
#ifndef QDM_MCMC_H
#define QDM_MCMC_H 1

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <hdf5.h>

#include "ijk.h"
#include "tau.h"

typedef struct {
  const gsl_rng *rng;

  size_t burn;
  size_t iter;
  size_t thin;

  size_t acc_check;
  size_t spline_df;

  const gsl_vector *xi;
  const gsl_vector *ll;
  const gsl_vector *tau;
  const gsl_matrix *theta;

  double xi_prior_mean;
  double xi_prior_var;
  double xi_tune_sd;

  double theta_min;
  double theta_tune_sd;

  const gsl_vector *x;
  const gsl_vector *y;

  const qdm_tau *t;
  const gsl_vector *m_knots;

  bool truncate;
} qdm_mcmc_parameters;

typedef struct {
  gsl_vector *xi;
  gsl_vector *xi_acc;
  gsl_vector *xi_p;
  gsl_vector *xi_tune;

  gsl_vector *ll;
  gsl_vector *ll_p;

  gsl_vector *tau;
  gsl_vector *tau_p;

  gsl_matrix *theta;
  gsl_matrix *theta_acc;
  gsl_matrix *theta_p;
  gsl_matrix *theta_star;
  gsl_matrix *theta_star_p;
  gsl_matrix *theta_tune;
} qdm_mcmc_workspace;

typedef struct {
  size_t s;

  qdm_ijk *theta;
  qdm_ijk *theta_star;

  qdm_ijk *ll;
  qdm_ijk *tau;

  qdm_ijk *xi;
} qdm_mcmc_results;

typedef struct {
  qdm_mcmc_parameters p;
  qdm_mcmc_workspace w;
  qdm_mcmc_results r;
} qdm_mcmc;

qdm_mcmc *
qdm_mcmc_alloc(
    qdm_mcmc_parameters p
);

void
qdm_mcmc_free(
    qdm_mcmc *mcmc
);

int
qdm_mcmc_run(
    qdm_mcmc *mcmc
);

int
qdm_mcmc_next(
    qdm_mcmc *mcmc
);

int
qdm_mcmc_update_theta(
    qdm_mcmc *mcmc
);

int
qdm_mcmc_update_xi(
    qdm_mcmc *mcmc
);

void
qdm_mcmc_update_tune(
    qdm_mcmc *mcmc
);

int
qdm_mcmc_save(
    qdm_mcmc *mcmc,
    size_t k
);

int
qdm_mcmc_write(
    hid_t id,
    const qdm_mcmc *mcmc
);

int
qdm_mcmc_read(
    hid_t id,
    qdm_mcmc **mcmc
);

#endif /* QDM_MCMC_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

#include "su3_matrix.c"

#ifndef gauge_config
#include "config.h"
#define gauge_config 1

#define sq(x) ((x)*(x))
#define cb(x) ((x)*(x)*(x))
#endif

int acc, tot; // acceptance rate

#define FORWARD 1 //true
#define BACKWARD 0 //false

// Access pattern methods --------------------------------------------------

// index of link site=(t,x,y,z) in direction dir
int index(const int* site, int dir){
  return dir+site[0]*4+site[1]*N*4+site[2]*N*N*4+site[3]*N*N*N*4;
}

// shift site array by one to direction dir
void shift(int* site, int dir, int forward){
  if(forward == 1){
    site[dir] = (site[dir]+1)%N;
  }else{
    site[dir] = (site[dir]-1+N)%N;
  }
}

// spatial rotation by 90deg of shape array
void rotate(int n, int* shape){
  for(int i=0;i<n;i++){
    if(shape[i]==1) shape[i]=2;
    if(shape[i]==2) shape[i]=3;
    if(shape[i]==3) shape[i]=1;
  }
}

// end of access pattern methods -------------------------------------------

// Plaquette methods -------------------------------------------------------

// plaquette operator at site (t,x,y,z) in plane (mu,nu)
// normally mu>nu
double plaquette(const su3_matrix* links, int mu, int nu, int* site){
  su3_matrix A,B,C,D;
  A = links[index(site,mu)];
  shift(site,mu,FORWARD);
  B = links[index(site,nu)];
  shift(site,mu,BACKWARD);
  shift(site,nu,FORWARD);
  C = su3_inv(links[index(site,mu)]);
  shift(site,nu,BACKWARD); // restore site position
  D = su3_inv(links[index(site,nu)]);

  A = su3_mul(A,B);
  A = su3_mul(A,C);
  A = su3_mul(A,D);

  return su3_trace(A)/3.0;
}

// plaquette spatial mean
double plaq_mean_s(int t, const su3_matrix* links){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
          site[0]=t; site[1]=x; site[2]=y; site[3]=z;
          for(int mu=1;mu<4;mu++){
            for(int nu=mu+1;nu<4;nu++){
              mean += plaquette(links,mu,nu,site);
            }
          }
      }
    }
  }
  return mean/ (3*N*N*N); // 3 directions & N^3 sites
}

// end of plaquette methods ------------------------------------------------

// Metropolis update --------------------------------------------------------

void do_hb_update(su3_matrix* links, const gsl_rng* r, const su3_matrix* rands, int* site, int dir);
double action(su3_matrix* links, int* site, int dir);

// the big loop
void hb_update(su3_matrix* links, const gsl_rng* r, const su3_matrix* rands){
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
          for(int dir=0;dir<4;dir++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            do_hb_update(links,r,rands,site,dir);
          }
        }
      }
    }
  }
}

// update
void do_hb_update(su3_matrix* links, const gsl_rng* r, const su3_matrix* rands, int* site, int dir){
  if(verbose) {printf("Updating {%d, %d, %d, %d} direction %d",site[0],site[1],site[2],site[3],dir);getchar();}
  double dS,action_old;
  su3_matrix link_old;
  for(int times=0; times<Nk; times++){
    action_old = action(links,site,dir);
    link_old = links[index(site,dir)];
    links[index(site,dir)]=su3_mul(links[index(site,dir)],rands[gsl_rng_uniform_int(r, 50)]);

    dS = action(links,site,dir)-action_old;

    if(verbose) {printf("Old action = %f, ",action_old); su3_print(link_old); getchar();}
    if(verbose) {printf("New action = %f, ",action_old+dS); su3_print(links[index(site,dir)]); getchar();}
    acc++;tot++;
    if(dS>0){
      if(exp(-dS)<gsl_rng_uniform(r)){
        acc--;
        links[index(site,dir)]=link_old; // revert update
        if(verbose) {printf("Reverted"); getchar();}
      }
    }
  }

}

// calculating the action
// for each direction we need to find plaquette at 2 sites in each of the 3 planes
// mu > nu
double action(su3_matrix* links, int* site, int dir){
  int mu,nu;
  double result=0;
  for(int k=0; k<4; k++){
    if(k==dir) continue; // loop over planes
    if(k>dir){mu=k; nu=dir;} // mu=larger number
    if(k<dir){mu=dir; nu=k;}

    result += plaquette(links, mu,nu,site);
    shift(site,k,BACKWARD); // shift to the perpendicular direction backwards
    result += plaquette(links, mu,nu,site);
    shift(site,k,FORWARD);  // shift it back
  }

  return -beta*result;
}

// end of heat bath update -------------------------------------------------

// Loops -------------------------------------------------------------------
// Generic loop ------------------------------------------------------------

// TODO: not working
double gloop(int n, int* shape, const su3_matrix* links, int* site){
  su3_matrix A[2*n];

  // forward
  for(int i=0;i<n;i++){
    A[i]=links[index(site,shape[i])];
    shift(site,shape[i],FORWARD);
  }

  // backward
  for(int i=0;i<n;i++){
    shift(site,shape[i],BACKWARD);
    A[n+i]=su3_inv(links[index(site,shape[i])]);
  }

  // multiply
  for(int i=1;i<2*n;i++){
    A[0]=su3_mul(A[0],A[i]);
  }

  return su3_trace(A[0])/3.0;
}

double loop(int r, int t, int mu, int nu, const su3_matrix* links, int* site){
  int shape[r+t];
  for(int i=0;i<t;i++)shape[i]=nu;
  for(int i=0;i<r;i++)shape[t+i]=mu;

  return gloop(r+t,shape,links,site);
  /*
  su3_matrix A[(r+t)*2];

  // bottom side
  for(int i=0;i<t;i++){
    A[i]=links[index(site,nu)];
    shift(site,nu,FORWARD);
  }

  // right side
  for(int i=0;i<r;i++){
    A[t+i]=links[index(site,mu)];
    shift(site,mu,FORWARD);
  }

  // top side
  for(int i=0;i<t;i++){
    shift(site,nu,BACKWARD);
    A[t+r+i]=su3_inv(links[index(site,nu)]);
  }

  // left side
  for(int i=0;i<r;i++){
    shift(site,mu,BACKWARD);
    A[t+r+t+i]=su3_inv(links[index(site,mu)]);
  }

  // multiply
  for(int i=1;i<2*(r+t);i++){
    A[0]=su3_mul(A[0],A[i]);
  }

  return su3_trace(A[0])/3.0;*/
}

// end of Loop -------------------------------------------------------------
// Wilson loop -------------------------------------------------------------

double wloop(int r, int t, int mu,const su3_matrix* links, int* site){
  return loop(r, t, mu, 0, links, site);
}

double wloop_mean(int wR, int wT, const su3_matrix* links){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            for(int mu=1;mu<4;mu++){
              mean += wloop(wR,wT,mu,links,site);
            }
        }
      }
    }
  }
  return mean/ (3*N*N*N*N); // 3 directions & N^3 sites
}
// end of Wilson loop ------------------------------------------------------
// Polyakov loop -----------------------------------------------------------
double polyakov(const su3_matrix* links, int* site){
  su3_matrix A[N];

  for(int i=0;i<N;i++){
    A[i]=links[index(site,0)];
    shift(site,0,FORWARD);
  }

  // multiply
  for(int i=1;i<N;i++){
    A[0]=su3_mul(A[0],A[i]);
  }

  return su3_trace(A[0])/3.0;
}

double polyakov_mean(const su3_matrix* links){
  double mean = 0;
  int site[4];
  site[0]=0;
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        site[1]=x; site[2]=y; site[3]=z;
        mean += polyakov(links,site);
      }
    }
  }
  return mean/ (N*N*N);
}
// end of Polyakov loop ----------------------------------------------------

// Twisted Wilson loop -----------------------------------------------------
// Ref: Lepage p.33

// r=sqrt(2)
double twloop2(int t, int mu, int nu ,const su3_matrix* links, int* site){
  su3_matrix A[(2+t)*2];

  // bottom side
  for(int i=0;i<t;i++){
    A[i]=links[index(site,0)];
    shift(site,0,FORWARD);
  }

  // right side, twisted
  A[t+0]=links[index(site,mu)];
  shift(site,mu,FORWARD);
  A[t+1]=links[index(site,nu)];
  shift(site,nu,FORWARD);

  // top side
  for(int i=0;i<t;i++){
    shift(site,0,BACKWARD);
    A[t+2+i]=su3_inv(links[index(site,0)]);
  }

  // left side, twisted
  shift(site,mu,BACKWARD);
  A[t+2+t+0]=su3_inv(links[index(site,mu)]);
  shift(site,nu,BACKWARD);
  A[t+2+t+1]=su3_inv(links[index(site,nu)]);

  // multiply
  for(int i=1;i<2*(2+t);i++){
    A[0]=su3_mul(A[0],A[i]);
  }

  return su3_trace(A[0])/3.0;
}

double twloop2_mean(int wT, const su3_matrix* links){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            for(int mu=1;mu<4;mu++){
              for(int nu=1;nu<4;nu++){
                if(mu!=nu){
                  mean += twloop2(wT,mu,nu,links,site);
                }
              }
            }
        }
      }
    }
  }
  return mean/ (6*N*N*N*N); // 3 x 2 directions & N^3 sites
}

// r=sqrt(3)
double twloop3(int t, const su3_matrix* links, int* site){
  su3_matrix A[(3+t)*2];

  // bottom side
  for(int i=0;i<t;i++){
    A[i]=links[index(site,0)];
    shift(site,0,FORWARD);
  }

  // right side, twisted
  A[t+0]=links[index(site,1)];
  shift(site,1,FORWARD);
  A[t+1]=links[index(site,2)];
  shift(site,2,FORWARD);
  A[t+2]=links[index(site,3)];
  shift(site,3,FORWARD);

  // top side
  for(int i=0;i<t;i++){
    shift(site,0,BACKWARD);
    A[t+3+i]=su3_inv(links[index(site,0)]);
  }

  // left side, twisted
  shift(site,1,BACKWARD);
  A[t+3+t+0]=su3_inv(links[index(site,1)]);
  shift(site,2,BACKWARD);
  A[t+3+t+1]=su3_inv(links[index(site,2)]);
  shift(site,3,BACKWARD);
  A[t+3+t+2]=su3_inv(links[index(site,3)]);

  // multiply
  for(int i=1;i<2*(3+t);i++){
    A[0]=su3_mul(A[0],A[i]);
  }

  return su3_trace(A[0])/3.0;
}

double twloop3_mean(int wT, const su3_matrix* links){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            mean += twloop3(wT,links,site);
        }
      }
    }
  }
  return mean/ (N*N*N*N); // N^3 sites
}
// end of twisted Wilson loop --------------------------------------------

// Plaquette correlator
double plaq_corr(int T, const su3_matrix* links){
  double mean=0;
  for(int i=0;i<N;i++){
    mean+=plaq_mean_s(i,links)*plaq_mean_s(i+T,links);
  }
  return mean/N;
}
// end of plaquette correlator

// TODO: vvvvvvvv

// Smeared Wilson loop ---------------------------------------------------

// end of smeared Wilson loop --------------------------------------------

max_stars_repo_name: uchuutamashi/lqcd

max_issues_repo_name: uchuutamashi/lqcd

max_forks_repo_name: uchuutamashi/lqcd


text:
#include <pygsl/solver.h>
#include <pygsl/rng.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

const char  * filename = __FILE__;
PyObject *module = NULL;
static const char monte_f_type_name[] = "F-monte";
static const char module_doc[] = "XXX Missing";
static const char PyGSL_monte_type[] = "monte";

enum PyGSL_gsl_monte_type{
     PyGSL_MONTE_plain = 1,
     PyGSL_MONTE_miser = 2,
     PyGSL_MONTE_vegas = 3
};

struct _monte_csys{
     enum PyGSL_gsl_monte_type type;
};

typedef struct _monte_csys monte_csys;

static double
PyGSL_g_test(double *k, size_t dim, void *params)
{
     double A = 1.0 / (M_PI * M_PI * M_PI);
     return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
}

static double
PyGSL_monte_function_wrap(double *x,  size_t dim, void *params)
{
     double tmp, *dresult;
     int dimension;
     PyGSL_solver *s;
     PyArrayObject *a_array = NULL;
     PyObject *result = NULL, *arglist=NULL, *callback, *retval;
     PyGSL_error_info  info;
     /* the line number to appear in the traceback */
     int trb_lineno = -1, i;
     struct pygsl_array_cache * cache_ptr;
     gsl_vector_view view;

     FUNC_MESS_BEGIN();
     s = (PyGSL_solver  *) params;

     /*
     flag = PyGSL_function_wrap_On_O(&view.vector, s->cbs[0], s->args,
                                     &tmp, NULL, view.vector.size,
                                     "monte_wrap");
     */

     FUNC_MESS_BEGIN();
     callback = s->cbs[0];
     assert(s->args != NULL);
     assert(callback != NULL);

     /* Do I need to copy the array ??? */
     for(i = 2; i < PyGSL_SOLVER_N_ARRAYS; ++i){
	  cache_ptr = &(s->cache[i]);
	  if(cache_ptr->ref == NULL)
	       /* no array found */
	       break;
	  if(x == cache_ptr->data){
	       a_array = cache_ptr->ref;
	       break;
	  }
     }
     if(i >= PyGSL_SOLVER_N_ARRAYS){
	  trb_lineno = __LINE__ -1;
	  pygsl_error("Not enough space to cache arrays...", filename, trb_lineno, GSL_ESANITY);
	  goto fail;
     }
     if(a_array == NULL){
	  dimension = dim;
	  a_array = (PyArrayObject *) PyArray_FromDimsAndData(1, &dimension, PyArray_DOUBLE, (char *)x);
	  cache_ptr->ref = a_array;
	  cache_ptr->data = x;
	  /*
	   * Required for deallocation.
	   * Tuples steal the references. But this array will be dereferenced by the
	   * normal procedure!
	   */
	  Py_INCREF(a_array);
     }
     if (a_array == NULL){
	  trb_lineno = __LINE__ - 2;
	  goto fail;
     }
     /* arglist was already set up */
     result = (PyObject *) s->cache[1].ref;
     dresult = (double *) s->cache[1].data;
     arglist = (PyObject *) s->cache[0].ref;

     Py_INCREF(s->args);
     PyTuple_SET_ITEM(arglist, 0, (PyObject *) a_array);
     PyTuple_SET_ITEM(arglist, 1, s->args);

     /*
     arglist = Py_BuildValue("(OOO)", a_array, s->args, result);
     if(DEBUG > 2){
	  fprintf(stderr, "callback = %p, arglist = %p\n", callback, arglist);
     }
     */
     FUNC_MESS("    Call Python Object BEGIN");
     retval = PyEval_CallObject(callback, arglist);
     FUNC_MESS("    Call Python Object END");

     /*
     info.callback = callback;
     info.message  = __FUNCTION__;
     if(PyGSL_CHECK_PYTHON_RETURN(retval, 0, &info) != GSL_SUCCESS){
	  trb_lineno = __LINE__ - 1;
	  goto fail;
     }
     info.argnum = 1;
     DEBUG_MESS(3, "result was %e", *dresult);
     */
     FUNC_MESS_END();
     return *dresult;

 fail:
     PyGSL_add_traceback(NULL, __FILE__, __FUNCTION__, trb_lineno);
     FUNC_MESS("Failure");

     if(s->isset == 1){
	  longjmp(s->buffer, GSL_EFAILED);
	  FUNC_MESS("\t\t Using jump buffer");
     } else {
	  FUNC_MESS("\t\t Jump buffer was not defined!");
     }
     tmp = gsl_nan();
     return tmp;
}

static PyObject *
PyGSL_monte_init(PyGSL_solver *self, PyObject *args)
{
     int flag = GSL_EFAILED;
     monte_csys * csys;

     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(self));
     csys = (monte_csys *)self->c_sys;
     assert(csys);
     switch(csys->type){
     case PyGSL_MONTE_plain:
	  flag = gsl_monte_plain_init(self->solver);
	  break;
     case PyGSL_MONTE_miser:
	  flag = gsl_monte_miser_init(self->solver);
	  break;
     case PyGSL_MONTE_vegas:
	  flag = gsl_monte_vegas_init(self->solver);
	  break;
     default:
	  DEBUG_MESS(2, "Monte type %d unknown",flag);
	  PyGSL_ERROR_NULL("Unknown monte type!", GSL_ESANITY);
     }

     if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
	  return NULL;
     }
     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;
}

typedef int (*pygsl_monte_fptr_t)(gsl_monte_function * F, double XL[], double XU[], size_t DIM,
		     size_t CALLS, gsl_rng * R, gsl_monte_plain_state * S,
		     double * RESULT, double * ABSERR);

static PyObject *
PyGSL_monte_integrate(PyGSL_solver *self, PyObject *args)
{

     int flag = GSL_EFAILED, line, dim;
     unsigned int ncalls;
     monte_csys * csys;
     gsl_rng * rng;
     PyObject  *func=NULL, *xlo=NULL, *xuo=NULL, *rng_obj=NULL, *ncallso=NULL, *argso=NULL;
     PyArrayObject *xla=NULL;
     PyArrayObject *xua=NULL;
     double result, abserr;
     gsl_monte_function mfunc;
     pygsl_monte_fptr_t fptr;

     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(self));

     if(!PyArg_ParseTuple(args, "OOOOO|O", &func, &xlo, &xuo, &ncallso, &rng_obj, &argso)){
	  line = __LINE__ - 1;
	  goto fail;
     }

     if(!PyCallable_Check(func)){
	  line = __LINE__ - 1;
	  pygsl_error("The function argument must be callable!", filename, line, GSL_EINVAL);
	  goto fail;
     }

     if(PyGSL_PYLONG_TO_UINT(ncallso, &ncalls, NULL) != GSL_SUCCESS){
	  line = __LINE__ - 1;
	  goto fail;
     }

     if(!PyGSL_RNG_Check(rng_obj)){
	  line = __LINE__ - 1;
	  pygsl_error("The rng object must be a rng!", filename, line, GSL_EINVAL);
	  goto fail;
     }

     rng = ((PyGSL_rng *) rng_obj)->rng;

     xla = PyGSL_vector_check(xlo, -1, PyGSL_DARRAY_CINPUT(2), NULL, NULL);
     if(xla == NULL){
	  line = __LINE__ - 2;
	  goto fail;
     }

     dim = xla->dimensions[0];
     xua = PyGSL_vector_check(xuo, dim, PyGSL_DARRAY_CINPUT(3), NULL, NULL);
     if(xua == NULL){
	  line = __LINE__ - 2;
	  goto fail;
     }

     if(argso == NULL){
	  Py_INCREF(Py_None);
	  argso = Py_None;
     }
     Py_XDECREF(self->cbs[0]);
     Py_INCREF(func);
     self->cbs[0] = func;
     self->args = argso;
     csys = (monte_csys *)self->c_sys;

     switch(csys->type){
     case PyGSL_MONTE_plain:	  fptr = (pygsl_monte_fptr_t)&gsl_monte_plain_integrate; break;
     case PyGSL_MONTE_miser:	  fptr = (pygsl_monte_fptr_t)&gsl_monte_miser_integrate; break;
     case PyGSL_MONTE_vegas:	  fptr = (pygsl_monte_fptr_t)&gsl_monte_vegas_integrate; break;
     default:
	  line = __LINE__ - 5;
	  DEBUG_MESS(2, "Monte type %d unknown",flag);
	  pygsl_error("Unknown monte type!", filename, line, GSL_ESANITY);
	  goto fail;
     }

     assert(fptr);
     mfunc.f = &PyGSL_monte_function_wrap;
     mfunc.f = &PyGSL_g_test;
     mfunc.dim = dim;
     mfunc.params = self;

     if(setjmp(self->buffer) != 0){
	  self->isset = 0;
	  line = __LINE__ - 2;
	  goto fail;
     }else{
	  self->isset = 1;
     }
     flag = fptr(&mfunc, (double *) xla->data, (double *)xua->data,
		 dim, ncalls, rng, self->solver, &result, &abserr);
     self->isset = 0;
     if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	  line = __LINE__ - 2;
	  return NULL;
     }

     Py_DECREF(xla);     xla = NULL;
     Py_DECREF(xua);     xua = NULL;
     Py_DECREF(self->cbs[0]);
     Py_DECREF(self->args);
     self->cbs[0] = NULL;
     self->args = NULL;
     FUNC_MESS_END();

     return Py_BuildValue("dd", result, abserr);

 fail:
     FUNC_MESS("FAIL");
     Py_XDECREF(xla);
     Py_XDECREF(xua);
     Py_XDECREF(self->cbs[0]);
     Py_XDECREF(self->args);
     self->cbs[0] = NULL;
     self->args = NULL;
     PyGSL_add_traceback(module, filename, __FUNCTION__, line);
     return NULL;
}

#define GET_SET(montetype, name, mode) \
static PyObject * GetSet_ ## montetype ## _ ## name(PyGSL_solver *self, PyObject *args) \
{ \
     gsl_monte_## montetype ## _state *state = self->solver; \
     return PyGSL_solver_GetSet((PyObject *)self, args, &(state->name), mode); \
}

#define GET_SET_DOUBLE(mintype, name) GET_SET(mintype, name, PyGSL_MODE_DOUBLE)
#define GET_SET_SIZE_T(mintype, name) GET_SET(mintype, name, PyGSL_MODE_SIZE_T)
#define GET_SET_INT(mintype, name)    GET_SET(mintype, name, PyGSL_MODE_INT)

GET_SET_DOUBLE(miser, estimate_frac)
GET_SET_SIZE_T(miser, min_calls)
GET_SET_SIZE_T(miser, min_calls_per_bisection)
GET_SET_DOUBLE(miser, alpha)
GET_SET_DOUBLE(miser, dither)

GET_SET_DOUBLE(vegas, result)
GET_SET_DOUBLE(vegas, sigma)
GET_SET_DOUBLE(vegas, chisq)
GET_SET_DOUBLE(vegas, alpha)
GET_SET_SIZE_T(vegas, iterations)
GET_SET_INT(vegas, stage)
GET_SET_INT(vegas, mode)
GET_SET_INT(vegas, verbose)

#undef GET_SET_DOUBLE
#undef GET_SET_SIZE_T
#undef GET_SET_INT
#undef GET_SET
#define GETSET(montetype, name, mode) {#name, (PyCFunction) GetSet_ ## montetype ## _ ## name, METH_VARARGS, NULL},
#define GET_SET_DOUBLE(montetype, name) GETSET(montetype, name, 0)
#define GET_SET_INT(montetype, name) GETSET(montetype, name, 0)
#define GET_SET_SIZE_T(montetype, name) GETSET(montetype, name, 0)

#define MONTE_STANDARD_METHODS \
     {"init",      (PyCFunction)PyGSL_monte_init,      METH_VARARGS, NULL}, \
     {"integrate", (PyCFunction)PyGSL_monte_integrate, METH_VARARGS, NULL}, \

static PyMethodDef PyGSL_monte_plain_methods[] = {
     MONTE_STANDARD_METHODS
     {NULL, NULL, 0, NULL}
};

static PyMethodDef PyGSL_monte_miser_methods[] = {
     MONTE_STANDARD_METHODS
     GET_SET_DOUBLE(miser, estimate_frac)
     GET_SET_SIZE_T(miser, min_calls)
     GET_SET_SIZE_T(miser, min_calls_per_bisection)
     GET_SET_DOUBLE(miser, alpha)
     GET_SET_DOUBLE(miser, dither)
     {NULL, NULL, 0, NULL}
};

static PyMethodDef PyGSL_monte_vegas_methods[] = {
     MONTE_STANDARD_METHODS
     GET_SET_DOUBLE(vegas, result)
     GET_SET_DOUBLE(vegas, sigma)
     GET_SET_DOUBLE(vegas, chisq)
     GET_SET_DOUBLE(vegas, alpha)
     GET_SET_SIZE_T(vegas, iterations)
     GET_SET_INT(vegas, stage)
     GET_SET_INT(vegas, mode)
     GET_SET_INT(vegas, verbose)
     {NULL, NULL, 0, NULL}
};

/* make the init look similar to the one */
static void *
PyGSL_gsl_monte_alloc(void * type, int n)
{
     enum PyGSL_gsl_monte_type flag = (enum PyGSL_gsl_monte_type) type;
     void * result = NULL;

     FUNC_MESS_BEGIN();
     switch(flag){
     case PyGSL_MONTE_plain:
	  result = gsl_monte_plain_alloc(n);
	  break;
     case PyGSL_MONTE_miser:
	  result = gsl_monte_miser_alloc(n);
	  break;
     case PyGSL_MONTE_vegas:
	  result = gsl_monte_vegas_alloc(n);
	  break;
     default:
	  DEBUG_MESS(2, "Monte type %d unknown",flag);
	  PyGSL_ERROR_NULL("Unknown monte type!", GSL_ESANITY);
     }
     FUNC_MESS_END();
     return result;
}

const struct _SolverStatic
monte_plain_solver_f   = {{(void_m_t) gsl_monte_plain_free,
			  (void_m_t) NULL,
			  NULL,
			  NULL},
			  1, PyGSL_monte_plain_methods,  PyGSL_monte_type},
monte_miser_solver_f   = {{(void_m_t) gsl_monte_miser_free,
			  (void_m_t) NULL,
			  NULL,
			   NULL},
			  1, PyGSL_monte_miser_methods,  PyGSL_monte_type},
monte_vegas_solver_f   = {{(void_m_t) gsl_monte_vegas_free,
			  (void_m_t) NULL,
			  NULL,
			  NULL},
			  1, PyGSL_monte_vegas_methods,  PyGSL_monte_type};

static PyObject*
PyGSL_init_monte(PyObject *self, PyObject *args, int type, const struct _SolverStatic *sost)
{
     PyGSL_solver *result = NULL;
     struct pygsl_array_cache *tmp;
     int line = -1, dims;
     monte_csys * csys;
     PyObject *mytuple;

     FUNC_MESS_BEGIN();
     const solver_alloc_struct alloc = {(const void *) type, (void *) PyGSL_gsl_monte_alloc, sost};
     result = (PyGSL_solver *) PyGSL_solver_dn_init(self, args, &alloc, 1);
     if(result == NULL){
	  line = __LINE__ - 2;
	  goto fail;
     }

     csys = (monte_csys *) calloc(1, sizeof(monte_csys));
     if(csys == NULL){
	  PyErr_NoMemory();
	  line = __LINE__ - 1;
	  goto fail;
     }

     dims = 1;
     tmp = result->cache;

     tmp[1].ref = PyGSL_New_Array(1, &dims, PyArray_DOUBLE);
     tmp[1].data = (double *) tmp[1].ref->data;
     mytuple =  PyTuple_New(3);
     PyTuple_SetItem(mytuple, 2, (PyObject *) tmp[1].ref);
     /*
      * Required for deallocation.
      * Tuples steal the references. But this array will be dereferenced by the
      * normal procedure!
      */
     Py_INCREF(tmp[1].ref);
     Py_INCREF(tmp[1].ref);

     /*
      * Initalise the others
      */
     Py_INCREF(Py_None);
     Py_INCREF(Py_None);
     Py_INCREF(Py_None);
     Py_INCREF(Py_None);
     PyTuple_SetItem(mytuple, 0, Py_None);
     PyTuple_SetItem(mytuple, 1, Py_None);

     tmp[0].ref = (PyArrayObject*) mytuple;

     result->cache = tmp;
     csys->type = type;
     result->c_sys = csys;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     FUNC_MESS("Fail");
     Py_XDECREF(result);
     PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
     return NULL;
}

#define MONTE_INIT(name) \
static PyObject* PyGSL_init_ ## name(PyObject * self, PyObject *args) \
{ return PyGSL_init_monte(self, args, PyGSL_MONTE_ ##name, &monte_ ## name ## _solver_f); }

MONTE_INIT(plain)
MONTE_INIT(miser)
MONTE_INIT(vegas)

static PyMethodDef mMethods[] = {
     {"plain", PyGSL_init_plain, METH_VARARGS, NULL},
     {"miser", PyGSL_init_miser, METH_VARARGS, NULL},
     {"vegas", PyGSL_init_vegas, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

void
initmonte(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("monte", mMethods);
     module = m;
     assert(m);
     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     init_pygsl()
     import_pygsl_solver();
     assert(PyGSL_API);

     if (!(item = PyString_FromString((char*)module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }
     PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_IMPORTANCE",      GSL_VEGAS_MODE_IMPORTANCE);
     PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_IMPORTANCE_ONLY", GSL_VEGAS_MODE_IMPORTANCE_ONLY);
     PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_STRATIFIED",      GSL_VEGAS_MODE_STRATIFIED);

     FUNC_MESS_END();
     return;

 fail:
     FUNC_MESS("FAIL");
     return;

}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#if !defined(XFORMCVSSTATUS_H_INCLUDED)
#define XFORMCVSSTATUS_H_INCLUDED

#include <filesystem>
#include <gsl/gsl>
#include <iosfwd>
#include <map>
#include <string>

class XformCvsStatus
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName,
		const char* pMsg);

	XformCvsStatus(::gsl::span<const char*const> args);
	int run();

	XformCvsStatus(const XformCvsStatus&) = delete;
	XformCvsStatus& operator=(const XformCvsStatus&) = delete;
	XformCvsStatus(XformCvsStatus&&) = delete;
	XformCvsStatus& operator=(XformCvsStatus&&) = delete;

private:
	using Path = ::std::filesystem::path;
	using FileToStatusMap = ::std::map<Path, ::std::string>;

	void insertInMap(const ::std::string& file, const ::std::string& status);
	Path buildPath(const ::std::string& file) const;
	static ::std::string xlateStatus(const ::std::string& status);
	static ::std::string pathStringToConsoleString(const Path::string_type& str);

	bool					m_suppressNew;
	bool					m_suppressUpToDate;
	bool					m_suppressLocal;
	Path					m_lastWorkingDir;
	FileToStatusMap		m_map;
};

#endif // XFORMCVSSTATUS_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
// A rectangular block of data of type float.  A transparent block
// cache is used for each instance of this class if necessary, so the
// image can be much bigger than will fit in memory, and pixels can
// still be accessed and written efficiently, provided subsequent
// accesses are spatially correlated.  A variety of useful methods are
// implemented (filtering, subsetting, interpolating, etc.)
//
// Don't try to access the same instance concurrently.  Split your
// images up into separate instances if you must parallelize things.
//
// For many methods, arguments of type ssize_t are used, but are not
// allowed to be negative.  This is to help prevent people from
// shooting themselves in the foor by accidently passing negative
// values which don't yield warnings and can't be caught with
// assertions.

#ifndef FLOAT_IMAGE_H
#define FLOAT_IMAGE_H

#ifndef solaris
#  include <stdint.h>
#endif
#include "asf_meta.h"
#include <stdio.h>
#include <sys/types.h>

#include <glib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_histogram.h>

// sometimes we don't have this - choose a conservative value
#ifndef SSIZE_MAX
#define SSIZE_MAX 32767
#endif

// Instance structure.  Everything here is private and need not be
// used or understood by client code, except for the size_x and size_y
// fields.
typedef struct {
  size_t size_x, size_y;    // Image dimensions.
  size_t cache_space;       // Memory cache space in bytes.
  size_t cache_area;        // Memory cache area in pixels.
  size_t tile_size;         // Tile size in pixels on a side.
  size_t cache_size_in_tiles;   // Number of tiles in cache.
  size_t tile_count_x;      // Number of tiles in image in x direction.
  size_t tile_count_y;      // Number of tiles in image in y direction.
  size_t tile_count;        // Total number of tiles in image.
  size_t tile_area;         // Area of a tile, in pixels.
  float *cache;             // Memory cache.
  float **tile_addresses;   // Addresss of individual tiles in the cache.
  GQueue *tile_queue;       // Queue of tile offsets kept in load order.
  FILE *tile_file;          // File with tiles stored contiguously.
  GString *tile_file_name;  // Name of the tile file
  int reference_count;      // For optional reference counting.
} FloatImage;

///////////////////////////////////////////////////////////////////////////////
//
// Creating New Instances
//
// Includeing methods which create new instances by copying existing
// ones.
//
///////////////////////////////////////////////////////////////////////////////

// Thaw out a previously frozen instance (produced with
// float_image_freeze) using data pointed to by file_pointer.  Frozen
// instances aren't portable between platforms.  After thawing
// file_pointer points to the data immediately following the data from
// the thawed instance (or to the end of the file).
FloatImage *
float_image_thaw (FILE *file_pointer);

// Create a new image filled with zero pixel values.
FloatImage *
float_image_new (ssize_t size_x, ssize_t size_y);

// Create a new image with pixels initialized to value.
FloatImage *
float_image_new_with_value (ssize_t size_x, ssize_t size_y, float value);

// Create a new image from memory.  This pixels are assumed to be
// layed out in memory in the usual way, i.e. contiguous rows of
// pixels in the x direction are contiguous in memory.
FloatImage *
float_image_new_from_memory (ssize_t size_x, ssize_t size_y, float *buffer);

// Create a new independent copy of model.
FloatImage *
float_image_copy (FloatImage *model);

// Form reduced resolution version of the model.  The scale_factor
// must be positive and odd.  The new image will be round ((double)
// model->size_x / scale_factor) pixels by round ((double)
// model->size_y / scale_factor) pixels.  Scaling is performed by
// averaging blocks of pixels together, using odd pixel reflection
// around the image edges (see the description of the apply_kernel
// method).  The upper and leftmost blocks of pixels averaged together
// are always centered at the 0 index in the direction in question, so
// reflection is always used for these edges.  Whether reflection is
// used for the right and lower edges depends on the relationship
// between the model dimensions and the scale factor.
FloatImage *
float_image_new_from_model_scaled (FloatImage *model, ssize_t scale_factor);

// Create a new image by copying the portion of model with upper left
// corner at model coordinates (x, y), width size_x, and height
// size_y.
FloatImage *
float_image_new_subimage (FloatImage *model, ssize_t x, ssize_t y,
              ssize_t size_x, ssize_t size_y);

// Type used to specify whether disk files should be in big or little
// endian byte order.
typedef enum {
  FLOAT_IMAGE_BYTE_ORDER_LITTLE_ENDIAN=1,
  FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN
} float_image_byte_order_t;

// Create a new image from data at byte offset in file.  The pixel
// layout in the file is assumed to be the same as for the
// float_image_new_from_memory method.  The byte order of individual
// pixels in the file should be byte_order.
FloatImage *
float_image_new_from_file (ssize_t size_x, ssize_t size_y, const char *file,
               off_t offset, float_image_byte_order_t byte_order);

// The method is like new_from_file, but takes a file pointer instead
// of a file name, and the offset argument is with respect to the
// current position in the file_pointer stream.
FloatImage *
float_image_new_from_file_pointer (ssize_t size_x, ssize_t size_y,
                   FILE *file_pointer, off_t offset,
                   float_image_byte_order_t byte_order);

// Form a low quality reduced resolution version of the
// original_size_x by original_size_y image in file.  The new image
// will be size_x by size_y pixels.  This method is like new_from_file
// method, but gets its data by sampling in each dimension using
// bilinear interpolation.  This is a decent way of forming quick
// thumbnails of images, or of scaling images down just slightly (by a
// factor of say 1.5 or less), but not much else.
FloatImage *
float_image_new_from_file_scaled (ssize_t size_x, ssize_t size_y,
                  ssize_t original_size_x,
                  ssize_t original_size_y,
                  const char *file, off_t offset,
                  float_image_byte_order_t byte_order);

// The function that does it all, generating an instance of FloatImage
// from a file and the metadata
FloatImage *
float_image_new_from_metadata(meta_parameters *meta, const char *file);

// For multi-band imagery the previous function needs to be more specific.
FloatImage *
float_image_band_new_from_metadata(meta_parameters *meta,
                   int band, const char *file);

// Sample type of an image that is to be used to create a float_image
// instance.  For example, floating point image can be created from
// signed sixteen bit integer data.
typedef enum {
  FLOAT_IMAGE_SAMPLE_TYPE_SIGNED_TWO_BYTE_INTEGER=1,
  FLOAT_IMAGE_SAMPLE_TYPE_UNSIGNED_BYTE
} float_image_sample_type;

// Form a new image by reading a data file full of sample_type
// samples.  The integer sample type are converted to floating point
// values by simple assignment (i.e. C assignment semantics apply,
// which for the smaller integer types should mean than an exact
// floating point representation is possible).  The other arguments
// are like those of the new_from_file method.
FloatImage *
float_image_new_from_file_with_sample_type
  (ssize_t size_x, ssize_t size_y, const char *file, off_t offset,
   float_image_byte_order_t byte_order, float_image_sample_type sample_type);

// This method is to new_from_file_with_sample_type as
// new_from_file_pointer is to new_from_file.
FloatImage *
float_image_new_from_file_pointer_with_sample_type
  (ssize_t size_x, ssize_t size_y, FILE *file_pointer, off_t offset,
   float_image_byte_order_t byte_order, float_image_sample_type sample_type);

///////////////////////////////////////////////////////////////////////////////
//
// Getting and Setting Image Pixels and Regions
//
///////////////////////////////////////////////////////////////////////////////

// Get pixel at 0-indexed position x, y.  A cache is used so pixel
// access does not usually involve the hard disk, provided subsequent
// pixel lookups are spatially close together.  You can probably get
// away with trating this function just as if everything was in
// memory.  For details, see the cache control methods below.
//
// The x and y arguments should always be positive, ssize_t is used
// only so people can't as easily shoot themselves in the foot by
// accidently supplying negative arguments which get promoted in some
// strange way.
float
float_image_get_pixel (FloatImage *self, ssize_t x, ssize_t y);

// Set pixel at 0-indexed position x, y to value.  A cache is used to
// make this fast, as for the float_image_get_pixel method.
void
float_image_set_pixel (FloatImage *self, ssize_t x, ssize_t y, float value);

// Get rectangular image region of size_x, size_y having upper left
// corner at x, y and copy it into already allocated buffer.  There is
// not necessarily any caching help for this method, i.e. it may
// always involve disk access and always be slow.
void
float_image_get_region (FloatImage *self, ssize_t x, ssize_t y,
            ssize_t size_x, ssize_t size_y, float *buffer);

// This method is analogous to float_image_get_region.
void
float_image_set_region (FloatImage *self, size_t x, size_t y, size_t size_x,
            size_t size_y, float *buffer);

// Get a full row of pixels, copying the data into already allocated
// buffer.  This method will be fast on the average for calls with
// sequential row numbers.
void
float_image_get_row (FloatImage *self, size_t row, float *buffer);

// Get a pixel, performing odd reflection at image edges if the pixel
// indicies fall outside the image.  See the description of the
// apply_kernel method for an explanation of reflection.
float
float_image_get_pixel_with_reflection (FloatImage *self, ssize_t x, ssize_t y);

///////////////////////////////////////////////////////////////////////////////
//
// Image Analysis and Statistics
//
///////////////////////////////////////////////////////////////////////////////

// Default mask value when figuring image stats
#define FLOAT_IMAGE_DEFAULT_MASK (0.0)

// Finds the minimum and maximum pixel values in the image, and the
// mean and standard deviation of all pixels.  This function considers
// every pixel in the image when mask is NAN, otherwise it discounts
// all values within .00000000001 of the mask
void
float_image_statistics (FloatImage *self, float *min, float *max, float *mean,
                        float *standard_deviation, float mask);

//  Does the same thing as float_image_statistics() but for only one selected
// band within a multi-band float image.  Returns 1 on error, otherwise 0
int
float_image_band_statistics(FloatImage *self, meta_stats *stats,
                            int line_count, int band_no,
                            float mask);

// This method works like the statistics method, except values in the
// interval [interval_start, interval_end] are not considered at all
// for the purposes of determining any of the outputs.
void
float_image_statistics_with_mask_interval (FloatImage *self, float *min,
                       float *max, float *mean,
                       float *standard_deviation,
                       double interval_start,
                       double interval_end);

// Compute an efficient estimate of the mean and standard deviation of
// the pixels in the image, by sampling every stride th pixel in each
// dimension, beginning with pixel (0, 0). If the mask is a non-NAN
// value, this function will discount all values within .00000000001
// of the mask
void
float_image_approximate_statistics (FloatImage *self, size_t stride,
                                    float *mean, float *standard_deviation,
                                    float mask);

// This method is a logical combination of the
// statistics_with_mask_interval and approximate_statistics methods.
void
float_image_approximate_statistics_with_mask_interval
  (FloatImage *self, size_t stride, float *mean, float *standard_deviation,
   double interval_start, double interval_end);

// Creates a gsl_histogram with 'num_bins' bins evenly spaced between
// 'min' and 'max'.  This function considers every pixel in the image.
gsl_histogram *
float_image_gsl_histogram (FloatImage *self, float min, float max,
                           size_t num_bins);

///////////////////////////////////////////////////////////////////////////////
//
// Kernels, Interpolation, and Sampling
//
///////////////////////////////////////////////////////////////////////////////

// Apply kernel centerd at pixel x, y and return the value.  The
// kernel matrix must be have equal odd dimensions.  The values in the
// kernel are multiplied by the pixels, and the sum of the products
// returned.  When part of the kernel would fall outside the image
// extents, the values used for the out-of-image pixels are the mirror
// images of the corresponding in-image pixels, with the edge pixels
// not duplicated, i.e. reflection about the middle of the edge pixels
// is used.
float
float_image_apply_kernel (FloatImage *self, ssize_t x, ssize_t y,
              gsl_matrix_float *kernel);

// Type used to specify whether disk files should be in big or little
// endian byte order.

// Sample method types.  These dictate how nearby pixels are
// considered when we want to find the approximate value for a point
// which falls between pixel indicies.
typedef enum {
  // Nearest pixel.
  FLOAT_IMAGE_SAMPLE_METHOD_NEAREST_NEIGHBOR,
  // Linearly weited average of four nearest pixels
  FLOAT_IMAGE_SAMPLE_METHOD_BILINEAR,
  // Bicubic spline interpolation (which consideres the nearest 16 pixels).
  FLOAT_IMAGE_SAMPLE_METHOD_BICUBIC
} float_image_sample_method_t;

float
float_image_sample (FloatImage *self, float x, float y,
            float_image_sample_method_t sample_method);

///////////////////////////////////////////////////////////////////////////////
//
// Comparing Images
//
///////////////////////////////////////////////////////////////////////////////

// Return true iff self and other have identical sizes and pixels that
// are all approximately equal to relative accuracy epsilon, as
// understood by the GNU Scientific Library function gsl_fcmp.
gboolean
float_image_equals (FloatImage *self, FloatImage *other, float epsilon);

///////////////////////////////////////////////////////////////////////////////
//
// Manipulating Images
//
///////////////////////////////////////////////////////////////////////////////

// Flip an image about a horizontal line through the center of the image
void
float_image_flip_y(FloatImage *self);

// Flip an image about a vertical line through the center of the image
void
float_image_flip_x(FloatImage *self);

///////////////////////////////////////////////////////////////////////////////
//
// Storing Images in Files
//
///////////////////////////////////////////////////////////////////////////////

// Store instance self at position pointed to by file_pointer, for
// later retrieval using float_image_thaw.  The serialized version of
// self is not portable between platforms.
void
float_image_freeze (FloatImage *self, FILE *file_pointer);

// Store image pixels in file.  The image is stored in the usual
// order, i.e. contiguous rows of pixels in the x direction are stored
// contiguously in memory.  Individual pixels are stored in byte order
// byte_order.  Returns 0 on success, nonzero on error.
int
float_image_store (FloatImage *self, const char *file,
           float_image_byte_order_t byte_order);
int
float_image_store_ext(FloatImage *self, const char *file,
                      float_image_byte_order_t byte_order, int append_flag);

int
float_image_band_store(FloatImage *self, const char *file,
               meta_parameters *meta, int append_flag);

///////////////////////////////////////////////////////////////////////////////
//
// Exporting Images in Various Image File Formats
//
///////////////////////////////////////////////////////////////////////////////

// Export image to file as a gray scaled jpeg image, with largest
// dimension no larger than max_dimension.  The max_dimension argument
// must be less than or equal to the largest dimension of the image.
// The image may be scaled st its largest dimension is considerably
// less than max_dimension.  Scaling is performed by averaging blocks
// of pixels together, using odd pixel reflection around the image
// edges (see the description of the apply_kernel method).  The JPEG
// format uses byte-valued gray scale samples, so the dynamic range of
// gray scale output pixels is limited to [0, 255].  Image pixel
// values inside two standard deviations of the mean pixel value are
// mapped linearly into this range; image pixel values outside two
// standard are clamped at the appropriate limit.  In determining the
// image statistics (mean and standard deviation), values equal to mask
// are not considered, unless mask is NAN, in which case mask has no
// effect.  If all image pixels have the same value, the output is
// made black if the pixels have value 0.0, and white otherwise.  This
// routine slurps the whole image into memory, so beware.  Returns 0
// on success, nonzero on error.
int
float_image_export_as_jpeg (FloatImage *self, const char *file,
                size_t max_dimension, double mask);

// This method works like the export_as_jpeg method, but all values in
// the interval [interval_start, interval_end] are considered to be
// uninteresting low values which should be mapped to zero in the
// output image, and should not be included in the calculations which
// determine the image statistics used to map the other floating point
// values into bytes.  For example, if one has an a bright island
// surrounded by very dark water in a radar image, setting the dark
// water covered areas to zero and not including them in the mean and
// standard deviation calculations will prevent the more interesting
// land areas from being driven into saturation in the generated
// image.
int
float_image_export_as_jpeg_with_mask_interval (FloatImage *self,
                           const char *file,
                           ssize_t max_dimension,
                           double interval_start,
                           double interval_end);

// Export image to file as a gray scaled tiff image, with largest
// dimension no larger than max_dimension.  The max_dimension argument
// must be less than or equal to the largest dimension of the image.
// The image may be scaled st its largest dimension is considerably
// less than max_dimension.  Scaling is performed by averaging blocks
// of pixels together, using odd pixel reflection around the image
// edges (see the description of the apply_kernel method).  The TIFF
// format (here) uses byte-valued gray scale samples, so the dynamic range
// of gray scale output pixels is limited to [0, 255].  Image pixel
// values inside two standard deviations of the mean pixel value are
// mapped linearly into this range; image pixel values outside two
// standard are clamped at the appropriate limit.  In determining the
// image statistics (mean and standard deviation), values equal to mask
// are not considered, unless mask is NAN, in which case mask has no
// effect.  If all image pixels have the same value, the output is
// made black if the pixels have value 0.0, and white otherwise.  This
// routine slurps the whole image into memory, so beware.  Returns 0
// on success, nonzero on error.
int
float_image_export_as_tiff (FloatImage *self, const char *file,
                            size_t max_dimension, double mask);

// This method exports the float image data as an ascii CSV file.
// Don't use this method on large sections of data - it will give
// an assertion failure if either dimension is larger than 255.
int
float_image_export_as_csv (FloatImage *self, const char *file);

///////////////////////////////////////////////////////////////////////////////
//
// Controlling Image Data Caching
//
// It probably isn't necessary to use these methods.  They are
// provided largely to make it clear how the cache works.  The major
// tunable parameter is the size of the in-memory cache to use.
//
// When a new image is created, the following steps are performed:
//
//      1. The image is divided up into square tiles st two full rows
//         or columns of tiles will fit in the memory cache.
//
//      2. A copy of the image is created on disk with the memory
//         layout rearranged st individual tiles are contiguous in
//         memory.  This allows tiles to be quickly retrieved later.
//
// When a pixel is accessed (read or set), the following happens:
//
//      1. If the pixel is in a tile already loaded into the cache,
//         it is simply fetched or set.
//
//      2. Otherwise, the tile containing the pixel is loaded,
//         possibly displacing an already loaded tile, and then the
//         pixel is fetched or set.  The tile displaced is the one
//         loaded longest ago (there is no most-recently-accessed
//         heuristic, as this would make pixel access too slow).
//
// Thus, using a larger memory cache will result in larger tiles being
// used, and fewer tile loads being needed.  In general, the default
// behavior is pretty good, but if you know will be performing lots of
// widely (but not too widely) scattered accesses, you might want to
// make it bigger.
//
///////////////////////////////////////////////////////////////////////////////

// Get the image memory cache size setting, in bytes.  Note that this
// is the memory cache used per image, not the class-wide cache usage.
// If you will have a lot of objects instantiated simultaneously, you
// may find it necessary to use a smaller cache for each image.
size_t
float_image_get_cache_size (FloatImage *self);

// Set the image memory cache to size bytes.  Changing the cache size
// requires the tiling to be recomputed, the on-disk tile cache to be
// regenerated, and the in memory cache to be flushed, so its slow.
void
float_image_set_cache_size (FloatImage *self, size_t size);

///////////////////////////////////////////////////////////////////////////////
//
// Reference Counting or Freeing Instances
//
///////////////////////////////////////////////////////////////////////////////

// Increment reference count.  Return pointer to self for convenience.
FloatImage *
float_image_ref (FloatImage *self);

// Decrement reference count, freeing instance if count falls to 0.
void
float_image_unref (FloatImage *self);

// Destroy self, regardless of reference count.
void
float_image_free (FloatImage *self);

#endif // #ifndef FLOAT_IMAGE_H

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/**
 *
 * @file core_sgemv.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:49 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 * CORE_sgemv performs one of the matrix-vector operations
 *
 *    y := alpha*A*x    + beta*y,   or
 *    y := alpha*A**T*x + beta*y,   or
 *    y := alpha*A**T*x + beta*y,
 *
 * where alpha and beta are scalars, x and y are vectors,
 * and A is an m by n matrix.
 *
 *******************************************************************************
 *
 * @param[in] trans
 *         @arg PlasmaNoTrans:   y := alpha*A*x    + beta*y
 *         @arg PlasmaTrans:     y := alpha*A**T*x + beta*y
 *         @arg PlasmaTrans: y := alpha*A**T*x + beta*y
 *
 * @param[in] m
 *         Number of rows of matrix A.
 *
 * @param[in] n
 *         Number of columns of matrix A.
 *
 * @param[in] alpha
 *         Scalar alpha.
 *
 * @param[in] A
 *         On entry, m by n matrix A. Dimension (lda,n).
 *
 * @param[in] lda
 *         Leading dimension of array A. lda >= max(1,m).
 *
 * @param[in] x
 *         On entry, vector x.
 *         If trans == PlasmaNoTrans, the n vector x has dimension 1 + (n-1)*abs(incx).
 *         Else, the m vector x has dimension 1 + (m-1)*abs(incx).
 *
 * @param[in] incx
 *         Increment between elements of x. incx must not be zero.
 *
 * @param[in] beta
 *         Scalar beta.
 *
 * @param[in,out] y
 *         On entry, vector y. On exit, y  is overwritten by updated vector y.
 *         If trans == PlasmaNoTrans, the m vector y has dimension 1 + (m-1)*abs(incy).
 *         Else, the n vector y has dimension 1 + (n-1)*abs(incy).
 *
 * @param[in] incy
 *         Increment between elements of y. incy must not be zero.
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgemv = PCORE_sgemv
#define CORE_sgemv PCORE_sgemv
#endif
void CORE_sgemv(PLASMA_enum trans,
                int m, int n,
                float alpha, const float *A, int lda,
                                          const float *x, int incx,
                float beta,        float *y, int incy)
{
    cblas_sgemv(
        CblasColMajor,
        (CBLAS_TRANSPOSE)trans,
        m, n,
        (alpha), A, lda,
        x, incx,
        (beta), y, incy );
}


text:
/**
 * Copyright (C) 2019 Samsung Electronics Co., Ltd. All Rights Reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *   http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *
 * @file	tensor.h
 * @date	04 December 2019
 * @brief	This is Tensor class for calculation
 * @see		https://github.com/nnstreamer/nntrainer
 * @author	Jijoong Moon <jijoong.moon@samsung.com>
 * @bug		No known bugs except for NYI items
 *
 */

#ifndef __TENSOR_H__
#define __TENSOR_H__
#ifdef __cplusplus

#ifdef USE_BLAS
extern "C" {
#include <cblas.h>
}
#endif

#include <cmath>
#include <fstream>
#include <iostream>
#include <memory>
#include <regex>
#include <vector>

namespace nntrainer {

#define MAXDIM 4

class TensorDim {
public:
  TensorDim() {
    for (int i = 0; i < MAXDIM; ++i) {
      dim[i] = 1;
    }
  }
  ~TensorDim(){};
  unsigned int batch() { return dim[0]; };
  unsigned int channel() { return dim[1]; };
  unsigned int height() { return dim[2]; };
  unsigned int width() { return dim[3]; };

  void batch(unsigned int b) { dim[0] = b; };
  void channel(unsigned int c) { dim[1] = c; };
  void height(unsigned int h) { dim[2] = h; };
  void width(unsigned int w) { dim[3] = w; };

  unsigned int *getDim() { return dim; }

  int setTensorDim(std::string input_shape);

private:
  unsigned int dim[4];
};

/**
 * @class   Tensor Class for Calculation
 * @brief   Tensor Class for Calculation
 */
class Tensor {
public:
  /**
   * @brief     Constructor of Tensor
   */
  Tensor() : height(0), width(0), batch(0), ndim(0), len(0){};

  /**
   * @brief     Constructor of Tensor with batch size one
   * @param[in] heihgt Height of Tensor
   * @param[in] width Width of Tensor
   */
  Tensor(int height, int width);

  /**
   * @brief     Constructor of Tensor
   * @param[in] batch Batch of Tensor
   * @param[in] heihgt Height of Tensor
   * @param[in] width Width of Tensor
   */
  Tensor(int batch, int height, int width);

  /**
   * @brief   Constructor of Tensor
   * @param[in] d data for the Tensor with batch size one
   */
  Tensor(std::vector<std::vector<float>> const &d);

  /**
   * @brief     Constructor of Tensor
   * @param[in] d data for the Tensor
   */
  Tensor(std::vector<std::vector<std::vector<float>>> const &d);

  /**
   * @brief     return value at specific location
   * @param[in] batch batch location
   * @param[in] h height location
   * @param[in] w width location
   */
  float getValue(int batch, int h, int w);

  /**
   * @brief     Multiply value element by element
   * @param[in] value multiplier
   * @retval    Calculated Tensor
   */
  Tensor multiply(float const &value);

  /**
   * @brief     Divide value element by element
   * @param[in] value Divisor
   * @retval    Calculated Tensor
   */
  Tensor divide(float const &value);

  /**
   * @brief     Add Tensor Element by Element
   * @param[in] m Tensor to be added
   * @retval    Calculated Tensor
   */
  Tensor add(Tensor const &m) const;

  /**
   * @brief     Add value Element by Element
   * @param[in] value value to be added
   * @retval    Calculated Tensor
   */
  Tensor add(float const &value);

  /**
   * @brief     Substract Tensor Element by Element
   * @param[in] m Tensor to be added
   * @retval    Calculated Tensor
   */
  Tensor subtract(Tensor const &m) const;

  /**
   * @brief     subtract value Element by Element
   * @param[in] value value to be added
   * @retval    Calculated Tensor
   */
  Tensor subtract(float const &value);

  /**
   * @brief     Multiply Tensor Element by Element ( Not the MxM )
   * @param[in] m Tensor to be multiplied
   * @retval    Calculated Tensor
   */
  Tensor multiply(Tensor const &m) const;

  /**
   * @brief     Divide Tensor Element by Element
   * @param[in] m Divisor Tensor
   * @retval    Calculated Tensor
   */
  Tensor divide(Tensor const &m) const;

  /**
   * @brief     Dot Product of Tensor ( equal MxM )
   * @param[in] m Tensor
   * @retval    Calculated Tensor
   */
  Tensor dot(Tensor const &m) const;

  /**
   * @brief     Transpose Tensor
   * @retval    Calculated Tensor
   */
  Tensor transpose() const;

  /**
   * @brief     sum all the Tensor elements according to the batch
   * @retval    Calculated Tensor(batch, 1, 1)
   */
  Tensor sum() const;

  /**
   * @brief     sum all the Tensor elements according to the axis
   * @retval    Calculated Tensor
   */
  Tensor sum(int axis) const;

  /**
   * @brief     Averaging the Tensor elements according to the batch
   * @retval    Calculated Tensor(1, height, width)
   */
  Tensor average() const;

  /**
   * @brief     Softmax the Tensor elements
   * @retval    Calculated Tensor
   */
  Tensor softmax() const;

  /**
   * @brief     l2norm the Tensor elements
   * @retval    Calculated l2norm
   */
  float l2norm() const;

  /**
   * @brief     Normalize the Tensor elements
   * @retval    Calculated Tensor
   */
  Tensor normalization() const;

  /**
   * @brief     Standardize the Tensor elements
   * @retval    Calculated Tensor
   */
  Tensor standardization() const;

  /**
   * @brief     Fill the Tensor elements with zero
   */
  void setZero();

  /**
   * @brief     Reduce Rank ( Tensor to Vector )
   * @retval    Saved vector
   */
  std::vector<float> mat2vec();

  /**
   * @brief     Apply function element by element
   * @param[in] *function function pointer applied
   * @retval    Tensor
   */
  Tensor apply(float (*function)(float)) const;

  /**
   * @brief     Apply function to Tensor
   * @param[in] *function function pointer applied
   * @retval    Tensor
   */
  Tensor apply(Tensor (*function)(Tensor)) const;

  /**
   * @brief     Print element
   * @param[in] out out stream
   * @retval    Tensor
   */
  void print(std::ostream &out) const;

  /**
   * @brief     Get Width of Tensor
   * @retval    int Width
   */
  int getWidth() { return width; };

  /**
   * @brief     Get Height of Tensor
   * @retval    int Height
   */
  int getHeight() { return height; };

  /**
   * @brief     Get Batch of Tensor
   * @retval    int Batch
   */
  int getBatch() { return batch; };

  /**
   * @brief     Set the elelemnt value
   * @param[in] batch batch location
   * @param[in] i height location
   * @param[in] j width location
   * @param[in] value value to be stored
   */
  void setValue(int batch, int i, int j, float value);

  /**
   * @brief     Copy the Tensor
   * @param[in] from Tensor to be Copyed
   * @retval    Matix
   */
  Tensor &copy(Tensor const &from);

  /**
   * @brief     Save the Tensor into file
   * @param[in] file output file stream
   */
  void save(std::ofstream &file);

  /**
   * @brief     Read the Tensor from file
   * @param[in] file input file stream
   */
  void read(std::ifstream &file);

  /**
   * @brief     return argument index which value is max
   * @retval    int argument index
   */
  int argmax();

  /**
   * @brief     return Data pointer of Tensor
   * @retval    float pointer
   */
  float *getData() { return data.data(); }

private:
  /**< handle the data as a std::vector type */
  std::vector<float> data;
  int height;
  int width;
  int batch;
  int ndim;
  int len;
};

/**
 * @brief   Overriding output stream
 */
std::ostream &operator<<(std::ostream &out, Tensor const &m);

} /* namespace nntrainer */

#endif /* __cplusplus */
#endif /* __TENSOR_H__ */

max_stars_repo_name: zhoonit/nntrainer

max_issues_repo_name: zhoonit/nntrainer

max_forks_repo_name: zhoonit/nntrainer


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#include <string.h>
#include "cconfigspace_internal.h"
#include "distribution_internal.h"

struct _ccs_distribution_multivariate_data_s {
	_ccs_distribution_common_data_t  common_data;
	size_t                           num_distributions;
	ccs_distribution_t              *distributions;
	size_t                          *dimensions;
	ccs_interval_t                  *bounds;
};
typedef struct _ccs_distribution_multivariate_data_s _ccs_distribution_multivariate_data_t;

static ccs_result_t
_ccs_distribution_multivariate_del(ccs_object_t o) {
	struct _ccs_distribution_multivariate_data_s *data =
		(struct _ccs_distribution_multivariate_data_s *)(((ccs_distribution_t)o)->data);
	for (size_t i = 0; i < data->num_distributions; i++) {
		ccs_release_object(data->distributions[i]);
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_multivariate_get_bounds(_ccs_distribution_data_t *data,
                                          ccs_interval_t           *interval_ret);

static ccs_result_t
_ccs_distribution_multivariate_samples(_ccs_distribution_data_t *data,
                                       ccs_rng_t                 rng,
                                       size_t                    num_values,
                                       ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_multivariate_strided_samples(_ccs_distribution_data_t *data,
                                               ccs_rng_t                 rng,
                                               size_t                    num_values,
                                               size_t                    stride,
                                               ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_multivariate_soa_samples(_ccs_distribution_data_t  *data,
                                           ccs_rng_t                  rng,
                                           size_t                     num_values,
                                           ccs_numeric_t            **values);

static _ccs_distribution_ops_t _ccs_distribution_multivariate_ops = {
	{ &_ccs_distribution_multivariate_del },
	&_ccs_distribution_multivariate_samples,
	&_ccs_distribution_multivariate_get_bounds,
	&_ccs_distribution_multivariate_strided_samples,
	&_ccs_distribution_multivariate_soa_samples
};

ccs_result_t
ccs_create_multivariate_distribution(size_t              num_distributions,
                                     ccs_distribution_t *distributions,
                                     ccs_distribution_t *distribution_ret) {
	CCS_CHECK_ARY(num_distributions, distributions);
	CCS_CHECK_PTR(distribution_ret);
	if (!num_distributions || num_distributions > INT64_MAX)
		return -CCS_INVALID_VALUE;

	ccs_result_t err;
	size_t       i = 0;
	size_t       dimension = 0;
	ccs_distribution_t distrib;
	_ccs_distribution_multivariate_data_t *distrib_data;

        for (i = 0; i < num_distributions; i++) {
		size_t dim;
		CCS_VALIDATE(ccs_distribution_get_dimension(distributions[i], &dim));
		dimension += dim;
	}

	uintptr_t mem, cur_mem;
	mem = (uintptr_t)calloc(1,
		sizeof(struct _ccs_distribution_s) +
		sizeof(_ccs_distribution_multivariate_data_t) +
		sizeof(ccs_distribution_t)*num_distributions +
		sizeof(size_t)*num_distributions +
		sizeof(ccs_interval_t)*dimension +
		sizeof(ccs_numeric_type_t)*dimension);

	if (!mem)
		return -CCS_OUT_OF_MEMORY;
        cur_mem = mem;

	distrib = (ccs_distribution_t)cur_mem;
	cur_mem += sizeof(struct _ccs_distribution_s);
	_ccs_object_init(&(distrib->obj), CCS_DISTRIBUTION, (_ccs_object_ops_t *)&_ccs_distribution_multivariate_ops);
	distrib_data = (_ccs_distribution_multivariate_data_t *)(cur_mem);
	cur_mem += sizeof(_ccs_distribution_multivariate_data_t);
	distrib_data->common_data.type        = CCS_MULTIVARIATE;
	distrib_data->common_data.dimension   = dimension;
	distrib_data->num_distributions       = num_distributions;
	distrib_data->distributions           = (ccs_distribution_t *)(cur_mem);
	cur_mem += sizeof(ccs_distribution_t)*num_distributions;
	distrib_data->dimensions              = (size_t *)(cur_mem);
	cur_mem += sizeof(size_t)*num_distributions;
	distrib_data->bounds                  = (ccs_interval_t *)(cur_mem);
	cur_mem += sizeof(ccs_interval_t)*dimension;
	distrib_data->common_data.data_types  = (ccs_numeric_type_t *)(cur_mem);
	cur_mem += sizeof(ccs_numeric_type_t)*dimension;

	dimension = 0;
        for (i = 0; i < num_distributions; i++) {
		size_t dim;
		CCS_VALIDATE_ERR_GOTO(err,
		  ccs_distribution_get_data_types(distributions[i], distrib_data->common_data.data_types + dimension),
		  errmemory);
		CCS_VALIDATE_ERR_GOTO(err,
		  ccs_distribution_get_bounds(distributions[i], distrib_data->bounds + dimension),
		  errmemory);
		CCS_VALIDATE_ERR_GOTO(err, ccs_distribution_get_dimension(distributions[i], &dim), errmemory);
		distrib_data->dimensions[i] = dim;
		dimension += dim;
	}

	for (i = 0; i < num_distributions; i++) {
		CCS_VALIDATE_ERR_GOTO(err, ccs_retain_object(distributions[i]), distrib);
		distrib_data->distributions[i] = distributions[i];
	}
	distrib->data = (_ccs_distribution_data_t *)distrib_data;
	*distribution_ret = distrib;

	return CCS_SUCCESS;
distrib:
	for (i = 0; i < num_distributions; i++) {
		if (distrib_data->distributions[i])
			ccs_release_object(distributions[i]);
	}
errmemory:
	free((void *)mem);
	return err;
}

static ccs_result_t
_ccs_distribution_multivariate_get_bounds(_ccs_distribution_data_t *data,
                                     ccs_interval_t           *interval_ret) {
	_ccs_distribution_multivariate_data_t *d =
		(_ccs_distribution_multivariate_data_t *)data;
	memcpy(interval_ret, d->bounds, d->common_data.dimension*sizeof(ccs_interval_t));
	return CCS_SUCCESS;
}

static inline _ccs_distribution_ops_t *
ccs_distribution_get_ops(ccs_distribution_t distribution) {
	return (_ccs_distribution_ops_t *)distribution->obj.ops;
}

static ccs_result_t
_ccs_distribution_multivariate_samples(_ccs_distribution_data_t *data,
                                       ccs_rng_t                 rng,
                                       size_t                    num_values,
                                       ccs_numeric_t            *values) {
	_ccs_distribution_multivariate_data_t *d =
		(_ccs_distribution_multivariate_data_t *)data;

	for (size_t i = 0; i < d->num_distributions; i++) {
		CCS_VALIDATE(ccs_distribution_get_ops(d->distributions[i])->strided_samples(
			d->distributions[i]->data, rng, num_values,
			d->common_data.dimension, values));
		values += d->dimensions[i];
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_multivariate_strided_samples(_ccs_distribution_data_t *data,
                                               ccs_rng_t                 rng,
                                               size_t                    num_values,
                                               size_t                    stride,
                                               ccs_numeric_t            *values) {
	_ccs_distribution_multivariate_data_t *d =
		(_ccs_distribution_multivariate_data_t *)data;

	for (size_t i = 0; i < d->num_distributions; i++) {
		CCS_VALIDATE(ccs_distribution_get_ops(d->distributions[i])->strided_samples(
			d->distributions[i]->data, rng, num_values,
			stride, values));
		values += d->dimensions[i];
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_multivariate_soa_samples(_ccs_distribution_data_t  *data,
                                           ccs_rng_t                  rng,
                                           size_t                     num_values,
                                           ccs_numeric_t            **values) {
	_ccs_distribution_multivariate_data_t *d =
		(_ccs_distribution_multivariate_data_t *)data;

	for (size_t i = 0; i < d->num_distributions; i++) {
		CCS_VALIDATE(ccs_distribution_get_ops(d->distributions[i])->soa_samples(
			d->distributions[i]->data, rng, num_values, values));
		values += d->dimensions[i];
	}
	return CCS_SUCCESS;
}

ccs_result_t
ccs_multivariate_distribution_get_num_distributions(ccs_distribution_t  distribution,
                                               size_t             *num_distributions_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_MULTIVARIATE);
	CCS_CHECK_PTR(num_distributions_ret);
	_ccs_distribution_multivariate_data_t * data = (_ccs_distribution_multivariate_data_t *)distribution->data;
	*num_distributions_ret = data->num_distributions;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_multivariate_distribution_get_distributions(ccs_distribution_t  distribution,
                                           size_t              num_distributions,
                                           ccs_distribution_t *distributions,
                                           size_t             *num_distributions_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_MULTIVARIATE);
	CCS_CHECK_ARY(num_distributions, distributions);
	if (!distributions && !num_distributions_ret)
		return -CCS_INVALID_VALUE;
	_ccs_distribution_multivariate_data_t * data = (_ccs_distribution_multivariate_data_t *)distribution->data;
	if (distributions) {
		if (num_distributions < data->num_distributions)
			return -CCS_INVALID_VALUE;
		for (size_t i = 0; i < data->num_distributions; i++)
			distributions[i] = data->distributions[i];
		for (size_t i = data->num_distributions; i < num_distributions; i++)
			distributions[i] = NULL;
	}
	if (num_distributions_ret)
		*num_distributions_ret = data->num_distributions;
	return CCS_SUCCESS;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS