# HG changeset patch # User Andrew L Janke # Date 1331674449 -36000 # Node ID 72d418b0012633ce26a01339b7237f965d58265e # Parent aa328918dc06c8eb4d74933be8284453d6f4728a Added mincsample and bumped version number diff --git a/progs/mincsample/README b/progs/mincsample/README new file mode 100644 --- /dev/null +++ b/progs/mincsample/README @@ -0,0 +1,2 @@ +Andrew Janke - a.janke@gmail.com +http://a.janke.googlepages.com/ diff --git a/progs/mincsample/mincsample.c b/progs/mincsample/mincsample.c new file mode 100644 --- /dev/null +++ b/progs/mincsample/mincsample.c @@ -0,0 +1,606 @@ +/* mincsample.c + * + * Generate samplings from MINC files + * + * Andrew Janke - a.janke@gmail.com + * Mark Griffin + * + * Copyright Andrew Janke and Mark Griffin, McConnell Brain Imaging Centre + * Permission to use, copy, modify, and distribute this software and its + * documentation for any purpose and without fee is hereby granted, + * provided that the above copyright notice appear in all copies. The + * author and the University make no representations about the + * suitability of this software for any purpose. It is provided "as is" + * without express or implied warranty. + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "mt19937ar.h" + +#ifndef FALSE +#define FALSE 0 +#endif + +#ifndef TRUE +#define TRUE 1 +#endif + +#define WORLD_NDIMS 3 +#define DEFAULT_INT -1 + +/* typedefs */ +typedef enum { SAMPLE_ALL, SAMPLE_RND } Sample_enum; +typedef enum { OUTPUT_ASCII, OUTPUT_DOUBLE } Output_enum; + +typedef struct { + Sample_enum sample_type; + + /* input parameters */ + int masking; + double mask_val; + int mask_idx; + + /* sampling */ + int rand_samples; + int max_samples; + + /* output parameters */ + int sample_mask; + int sample_mask_idx; + + Output_enum output_type; + int output_coords; + FILE *outFP; + } Loop_Data; + +/* function prototypes */ +void count_points(void *caller_data, long num_voxels, int input_num_buffers, + int input_vector_length, double *input_data[], + int output_num_buffers, int output_vector_length, + double *output_data[], Loop_Info * loop_info); +void get_points(void *caller_data, long num_voxels, int input_num_buffers, + int input_vector_length, double *input_data[], int output_num_buffers, + int output_vector_length, double *output_data[], + Loop_Info * loop_info); +void write_data(FILE * fp, double value, Output_enum ot); +void get_minc_attribute(int mincid, char *varname, char *attname, + int maxvals, double vals[]); +int get_minc_ndims(int mincid); +void find_minc_spatial_dims(int mincid, int space_to_dim[], int dim_to_space[]); +void get_minc_voxel_to_world(int mincid, + double voxel_to_world[WORLD_NDIMS][WORLD_NDIMS + 1]); +void normalize_vector(double vector[]); +void transform_coord(double out_coord[], + double transform[WORLD_NDIMS][WORLD_NDIMS + 1], + double in_coord[]); + +/* global vars for printing coordinates */ +int space_to_dim[WORLD_NDIMS] = { -1, -1, -1 }; +int dim_to_space[MAX_VAR_DIMS]; +int file_ndims = 0; +double voxel_to_world[WORLD_NDIMS][WORLD_NDIMS + 1]; + +/* Argument variables and table */ +static int verbose = FALSE; +static int quiet = FALSE; +static int clobber = FALSE; +static int max_buffer = 4 * 1024; +static char *mask_fname = NULL; +static char *sample_fname = NULL; +static char *out_fname = NULL; +static int append_output = FALSE; +static int rand_seed = DEFAULT_INT; +static Loop_Data md = { + SAMPLE_ALL, + FALSE, 1.0, 0, + 0, 0, + FALSE, 0, + OUTPUT_ASCII, FALSE, + NULL + }; + +static ArgvInfo argTable[] = { + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, + "General options:"}, + {"-verbose", ARGV_CONSTANT, (char *)TRUE, (char *)&verbose, + "print out extra information."}, + {"-quiet", ARGV_CONSTANT, (char *)TRUE, (char *)&quiet, + "be very quiet."}, + {"-clobber", ARGV_CONSTANT, (char *)TRUE, (char *)&clobber, + "clobber existing files."}, + {"-max_buffer", ARGV_INT, (char *)1, (char *)&max_buffer, + "maximum size of buffers (in kbytes)"}, + {"-mask", ARGV_STRING, (char *)1, (char *)&mask_fname, + "select voxels within the specified mask"}, + {"-mask_val", ARGV_FLOAT, (char *)1, (char *)&(md.mask_val), + "mask value to use"}, + + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, + "\nSampling Types:"}, + {"-all", ARGV_CONSTANT, (char *)SAMPLE_ALL, (char *)&md.sample_type, + "sample all the data (Default)"}, + {"-random_seed", ARGV_INT, (char *)1, (char *)&rand_seed, + "Random seed to use (use to get reproducible runs) Default: use tv_usec"}, + {"-random_samples", ARGV_INT, (char *)1, (char *)&md.rand_samples, + "take # random samples from the input data"}, + + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, + "\nOutput Options:"}, + {"-sample", ARGV_STRING, (char *)1, (char *)&sample_fname, + "Output a file of chosen points"}, + {"-outfile", ARGV_STRING, (char *)1, (char *)&out_fname, + " for output data (Default: stdout)"}, + {"-append", ARGV_CONSTANT, (char *)TRUE, (char *)&append_output, + "append output data to existing file"}, + {"-ascii", ARGV_CONSTANT, (char *)OUTPUT_ASCII, (char *)&md.output_type, + "Write out data as ascii strings (default)"}, + {"-double", ARGV_CONSTANT, (char *)OUTPUT_DOUBLE, (char *)&md.output_type, + "Write out data as double precision floating-point values"}, + {"-coords", ARGV_CONSTANT, (char *)TRUE, (char *)&md.output_coords, + "Write out world co-ordinates as well as values"}, + + {NULL, ARGV_HELP, NULL, NULL, ""}, + {NULL, ARGV_END, NULL, NULL, NULL} + }; + +int main(int argc, char *argv[]) +{ + char **infiles; + char *outfiles[1]; + int n_infiles, n_outfiles; + char *arg_string; + Loop_Options *loop_opts; + int mincid; + struct timeval timer; + int i; + + /* Save time stamp and args */ + arg_string = time_stamp(argc, argv); + + /* get arguments */ + if(ParseArgv(&argc, argv, argTable, 0) || (argc < 2)){ + fprintf(stderr, "\nUsage: %s [options] ... \n", argv[0]); + fprintf(stderr, " %s -help\n\n", argv[0]); + exit(EXIT_FAILURE); + } + + /* check arguments */ + if(md.rand_samples != 0){ + md.sample_type = SAMPLE_RND; + if(md.rand_samples < 0){ + fprintf(stderr, "%s: -rand_samples (%d) must be greater than 0\n\n", + argv[0], md.rand_samples); + exit(EXIT_FAILURE); + } + } + + /* check arguments */ + if(rand_seed != DEFAULT_INT && rand_seed < 0){ + fprintf(stderr, "%s: -rand_seed (%d) must be 0 or greater\n\n", argv[0], rand_seed); + exit(EXIT_FAILURE); + } + + /* get infile names */ + n_infiles = argc - 1; + infiles = (char **)malloc(sizeof(char *) * (n_infiles + 1)); /* + 1 for mask */ + for(i = 0; i < n_infiles; i++){ + infiles[i] = argv[i + 1]; + } + if(mask_fname != NULL){ + infiles[n_infiles] = mask_fname; + md.masking = TRUE; + md.mask_idx = n_infiles; + n_infiles++; + } + + /* check for the infile(s) */ + for(i = 0; i < n_infiles; i++){ + if(access(infiles[i], F_OK) != 0){ + fprintf(stderr, "%s: Couldn't find input file: %s\n\n", argv[0], infiles[i]); + exit(EXIT_FAILURE); + } + } + + /* set up and check for voxel_loop outfiles */ + if(sample_fname != NULL){ + if(access(sample_fname, F_OK) == 0 && !clobber){ + fprintf(stderr, "%s: %s exists, use -clobber to overwrite\n\n", argv[0], + sample_fname); + exit(EXIT_FAILURE); + } + md.sample_mask = TRUE; + md.sample_mask_idx = 0; + n_outfiles = 1; + outfiles[0] = sample_fname; + } + else { + n_outfiles = 0; + } + + /* set up data outfile */ + if(out_fname == NULL || strcmp(out_fname, "-") == 0){ + md.outFP = stdout; + } + else { + if(!append_output && access(out_fname, F_OK) == 0 && !clobber){ + fprintf(stderr, "%s: %s exists, use -clobber to overwrite\n\n", argv[0], + out_fname); + exit(EXIT_FAILURE); + } + + if((md.outFP = fopen(out_fname, (append_output) ? "a" : "w")) == NULL){ + fprintf(stderr, "%s: problems opening %s\n", argv[0], out_fname); + exit(EXIT_FAILURE); + } + } + + /* Get some information from the first file for printing co-ordinates */ + mincid = miopen(infiles[0], NC_NOWRITE | 0x8000); + file_ndims = get_minc_ndims(mincid); + find_minc_spatial_dims(mincid, space_to_dim, dim_to_space); + get_minc_voxel_to_world(mincid, voxel_to_world); + + /* set up voxel loop options */ + loop_opts = create_loop_options(); + set_loop_verbose(loop_opts, FALSE); + set_loop_clobber(loop_opts, clobber); + set_loop_buffer_size(loop_opts, (long)1024 * max_buffer); + + /* set up random sampling if required */ + if(md.sample_type == SAMPLE_RND){ + void *tmp = NULL; /* for gettimeofday */ + + /* get max number of samples */ + voxel_loop(1, (md.masking) ? &mask_fname : infiles, 0, NULL, NULL, + loop_opts, count_points, (void *)&md); + if(verbose){ + fprintf(stderr, " | Got max # of points: %d\n", md.max_samples); + } + + if(md.rand_samples > md.max_samples){ + fprintf(stderr, "%s: -rand_samples (%d) must be less than max samples (%d)\n\n", + argv[0], md.rand_samples, md.max_samples); + exit(EXIT_FAILURE); + } + + /* initialise random number generator */ + if(rand_seed == DEFAULT_INT){ + gettimeofday(&timer, tmp); + rand_seed = timer.tv_usec; + } + + if(verbose){ + fprintf(stderr, " | Using random seed: %d\n", rand_seed); + } + + init_genrand((unsigned long)rand_seed); + } + + /* do the sampling */ + voxel_loop(n_infiles, infiles, n_outfiles, outfiles, arg_string, + loop_opts, get_points, (void *)&md); + + /* tidy up */ + fclose(md.outFP); + free_loop_options(loop_opts); + + return (EXIT_SUCCESS); + } + +/* get points from file(s), write out to an input FP */ +void get_points(void *caller_data, long num_voxels, int input_num_buffers, + int input_vector_length, double *input_data[], int output_num_buffers, + int output_vector_length, double *output_data[], Loop_Info * loop_info) +{ + Loop_Data *md = (Loop_Data *) caller_data; + int i, idim, ivox; + int n_infiles; + int do_sample; + double mask_value; + int dim_index; + long index[MAX_VAR_DIMS]; + double voxel_coord[WORLD_NDIMS]; + double world_coord[WORLD_NDIMS]; + + /* shut the compiler up */ + (void)output_num_buffers; + (void)output_vector_length; + + n_infiles = (md->masking) ? input_num_buffers - 1 : input_num_buffers; + + /* for each voxel */ + for(ivox = 0; ivox < num_voxels * input_vector_length; ivox++){ + + /* nasty way that works for masking or not */ + mask_value = 0; + if(!md->masking || + (md->masking && fabs(input_data[md->mask_idx][ivox] - md->mask_val) < 0.5)){ + + /* flag to sample */ + do_sample = FALSE; + + switch (md->sample_type){ + case SAMPLE_ALL: + do_sample = TRUE; + mask_value = 1.0; + break; + + case SAMPLE_RND: + /* if this voxel 'qualifies', write out data */ + if(genrand_res53() < ((double)md->rand_samples / md->max_samples)){ + md->rand_samples--; + do_sample = TRUE; + mask_value = 1.0; + } + + md->max_samples--; + break; + + default: + fprintf(stderr, "ERROR - Sample type is undefined (%d)\n", md->sample_type); + exit(EXIT_FAILURE); + } + + /* now write out the data */ + if(do_sample){ + + /* get and convert voxel to world coordinates */ + if(md->output_coords){ + get_info_voxel_index(loop_info, ivox, file_ndims, index); + for(idim = 0; idim < WORLD_NDIMS; idim++){ + dim_index = space_to_dim[idim]; + if(dim_index >= 0){ + voxel_coord[idim] = index[dim_index]; + } + } + transform_coord(world_coord, voxel_to_world, voxel_coord); + } + + switch (md->output_type){ + case OUTPUT_ASCII: + if(md->output_coords){ + fprintf(md->outFP, "%.20g\t%.20g\t%.20g\t", world_coord[0], + world_coord[1], world_coord[2]); + } + + for(i = 0; i < n_infiles; i++){ + fprintf(md->outFP, "%.20g\t", input_data[i][ivox]); + } + fprintf(md->outFP, "\n"); + break; + + case OUTPUT_DOUBLE: + if(md->output_coords){ + fwrite(&(world_coord[0]), sizeof(double), 1, md->outFP); + fwrite(&(world_coord[1]), sizeof(double), 1, md->outFP); + fwrite(&(world_coord[2]), sizeof(double), 1, md->outFP); + } + + for(i = 0; i < n_infiles; i++){ + fwrite(&(input_data[i][ivox]), sizeof(double), 1, md->outFP); + } + + break; + + default: + fprintf(stderr, "ERROR - Output type is undefined (%d)\n", + md->output_type); + exit(EXIT_FAILURE); + } + } + + } + + /* output sampling mask */ + if(md->sample_mask){ + output_data[md->sample_mask_idx][ivox] = mask_value; + } + } + + } + +/* count points that are within a mask */ +void count_points(void *caller_data, long num_voxels, int input_num_buffers, + int input_vector_length, double *input_data[], + int output_num_buffers, int output_vector_length, + double *output_data[], Loop_Info * loop_info) +{ + Loop_Data *md = (Loop_Data *) caller_data; + int ivox; + + /* shut the compiler up */ + (void)input_num_buffers; + (void)output_num_buffers; + (void)output_vector_length; + (void)output_data; + (void)loop_info; + + for(ivox = 0; ivox < num_voxels * input_vector_length; ivox++){ + if(md->masking){ + if(fabs(input_data[0][ivox] - md->mask_val) < 0.5){ + md->max_samples++; + } + } + else { + md->max_samples++; + } + } + } + +inline void write_data(FILE * fp, double value, Output_enum ot) +{ + switch (ot){ + case OUTPUT_ASCII: + fprintf(fp, "%.20g\n", value); + break; + + case OUTPUT_DOUBLE: + fwrite(&(value), sizeof(double), 1, fp); + break; + + default: + fprintf(stderr, "ERROR - Output type is undefined (%d)\n", ot); + exit(EXIT_FAILURE); + } + } + +void normalize_vector(double vector[]) +{ + int idim; + double magnitude; + + magnitude = 0.0; + for(idim = 0; idim < WORLD_NDIMS; idim++){ + magnitude += (vector[idim] * vector[idim]); + } + magnitude = sqrt(magnitude); + if(magnitude > 0.0){ + for(idim = 0; idim < WORLD_NDIMS; idim++){ + vector[idim] /= magnitude; + } + } + } + +/* Transforms a coordinate through a linear transform -- from mincstats */ +void transform_coord(double out_coord[], + double transform[WORLD_NDIMS][WORLD_NDIMS + 1], double in_coord[]) +{ + int idim, jdim; + double homogeneous_coord[WORLD_NDIMS + 1]; + + for(idim = 0; idim < WORLD_NDIMS; idim++){ + homogeneous_coord[idim] = in_coord[idim]; + } + homogeneous_coord[WORLD_NDIMS] = 1.0; + + for(idim = 0; idim < WORLD_NDIMS; idim++){ + out_coord[idim] = 0.0; + for(jdim = 0; jdim < WORLD_NDIMS + 1; jdim++){ + out_coord[idim] += transform[idim][jdim] * homogeneous_coord[jdim]; + } + } + } + +/* Get the voxel to world transform (for column vectors) -- from mincstats */ +void get_minc_voxel_to_world(int mincid, + double voxel_to_world[WORLD_NDIMS][WORLD_NDIMS + 1]) +{ + int idim, jdim; + double dircos[WORLD_NDIMS]; + double step, start; + char *dimensions[] = { MIxspace, MIyspace, MIzspace }; + + /* Zero the matrix */ + for(idim = 0; idim < WORLD_NDIMS; idim++){ + for(jdim = 0; jdim < WORLD_NDIMS + 1; jdim++) + voxel_to_world[idim][jdim] = 0.0; + } + + for(jdim = 0; jdim < WORLD_NDIMS; jdim++){ + + /* Set default values */ + step = 1.0; + start = 0.0; + for(idim = 0; idim < WORLD_NDIMS; idim++) + dircos[idim] = 0.0; + dircos[jdim] = 1.0; + + /* Get the attributes */ + get_minc_attribute(mincid, dimensions[jdim], MIstart, 1, &start); + get_minc_attribute(mincid, dimensions[jdim], MIstep, 1, &step); + get_minc_attribute(mincid, dimensions[jdim], MIdirection_cosines, + WORLD_NDIMS, dircos); + normalize_vector(dircos); + + /* Put them in the matrix */ + for(idim = 0; idim < WORLD_NDIMS; idim++){ + voxel_to_world[idim][jdim] = step * dircos[idim]; + voxel_to_world[idim][WORLD_NDIMS] += start * dircos[idim]; + } + } + } + +/* Get the mapping from spatial dimension - x, y, z - to file dimensions + and vice-versa. -- from mincstats */ +void find_minc_spatial_dims(int mincid, int space_to_dim[], int dim_to_space[]) +{ + int imgid, dim[MAX_VAR_DIMS]; + int idim, ndims, world_index; + char dimname[MAX_NC_NAME]; + + /* Set default values */ + for(idim = 0; idim < 3; idim++) + space_to_dim[idim] = -1; + for(idim = 0; idim < MAX_VAR_DIMS; idim++) + dim_to_space[idim] = -1; + + /* Get the dimension ids for the image variable */ + imgid = ncvarid(mincid, MIimage); + (void)ncvarinq(mincid, imgid, NULL, NULL, &ndims, dim, NULL); + + /* Loop over them to find the spatial ones */ + for(idim = 0; idim < ndims; idim++){ + + /* Get the name and check that this is a spatial dimension */ + (void)ncdiminq(mincid, dim[idim], dimname, NULL); + if((dimname[0] == '\0') || (strcmp(&dimname[1], "space") != 0)){ + continue; + } + + /* Look for the spatial dimensions */ + switch (dimname[0]){ + case 'x': + world_index = 0; + break; + case 'y': + world_index = 1; + break; + case 'z': + world_index = 2; + break; + default: + world_index = 0; + break; + } + space_to_dim[world_index] = idim; + dim_to_space[idim] = world_index; + } + } + +/* Get a double attribute from a minc file -- from mincstats */ +void get_minc_attribute(int mincid, char *varname, char *attname, + int maxvals, double vals[]) +{ + int varid; + int old_ncopts; + int att_length; + + if(!mivar_exists(mincid, varname)) + return; + varid = ncvarid(mincid, varname); + old_ncopts = ncopts; + ncopts = 0; + (void)miattget(mincid, varid, attname, NC_DOUBLE, maxvals, vals, &att_length); + ncopts = old_ncopts; + } + +/* Get the total number of image dimensions in a minc file -- from mincstats */ +int get_minc_ndims(int mincid) +{ + int imgid; + int ndims; + + imgid = ncvarid(mincid, MIimage); + (void)ncvarinq(mincid, imgid, NULL, NULL, &ndims, NULL, NULL); + + return ndims; + } diff --git a/progs/mincsample/mincsample.man1 b/progs/mincsample/mincsample.man1 new file mode 100644 --- /dev/null +++ b/progs/mincsample/mincsample.man1 @@ -0,0 +1,99 @@ +.\" Hey, EMACS: -*- nroff -*- +.\" Copyright 2004 Andrew Janke, McConnell Brain Imaging Centre, +.\" Montreal Neurological Institute, McGill University. +.\" Permission to use, copy, modify, and distribute this +.\" software and its documentation for any purpose and without +.\" fee is hereby granted, provided that the above copyright +.\" notice appear in all copies. The author and McGill University +.\" make no representations about the suitability of this +.\" software for any purpose. It is provided "as is" without +.\" express or implied warranty. +.\" +.\" +.TH MINCSAMPLE 1 "MINC User's Guide" + +.SH NAME +mincsample - generate samplings from minc files. + +.SH SYNOPSIS +.B mincsample +[] [ [<..>]] + +.SH DESCRIPTION +\fIMincsample\fR produces a data sampling on STDOUT from an input +series of minc files. The output can be either ascii (-ascii) or as +a raw binary stream of doubles (-double). The output data is ordered +first by file then voxel. When -ascii is used the data values from +each file are separated by a tab and the sampling points with a newline. +When using -double, no separators are used. + +If -coords is also specified, the world co-ordinate at each sampling point +will precede the data from each of the files. An optional -outfile +argument can also be used to direct the output to a file and -append +used to append the data to the file as opposed to overwriting the file. + +By default all data points are written out (-all) the output of points can +also be constrained to be points within a mask (-mask and -mask_val) and further +by a random sampling of a sub-set of points via the -random_samples and +-random_seed arguments + +.SH OPTIONS +.TP +\fB\-verbose\fR +Print extra information during processing. +.TP +\fB\-quiet\fR +Print only the ouput sampling data. +.TP +\fB\-clobber\fR +Overwrite existing output files (when used with -outfile) +.TP +\fB\-max_buffer\fR \fIsize\fR +Specify the maximum size of the internal buffers (in kbytes). Default is 4096 (4MB). +.TP +\fB\-mask\fR \fImask.mnc\fR +Specify and input mask, only sampling points within this mask will be used. +.TP +\fB\-mask_val\fR \fIvalue\fR +Specify the value to use from the mask (Default: 1). +.TP +\fB\-all\fR +Sample all the data points. +.TP +\fB\-random_seed\fR \fIvalue\fR +Specify the random seed value to use during random sampling, this is to enable +reproducible runs. If no seed is given a semi-random seed will be chosen (from time). +.TP +\fB\-random_samples\fR \fIvalue\fR +Specify the number of random samples to take from the input files. This value must be smaller +than the maximum possible number of samples. +.TP +\fB\-sample\fR \fIsample.mnc\fR +Output a mask file that corresponds to where samples were taken from. +.TP +\fB\-outfile\fR \fIfile\fR +Output sampling data to a file. (Default: STDOUT). +.TP +\fB\-append\fR +Append output data to an existing file. +.TP +\fB\-ascii\fR +Write out data as ascii strings (Default). +.TP +\fB\-ascii\fR +Write out data as double precision floating-point values. +.TP +\fB\-coords\fR +Write out world co-ordinates as well as sampling values. +.TP +\fB-help\fR +Print summary of command-line options and exit. +.TP +\fB\-version\fR +Print the program's version number and exit. + +.SH AUTHOR +Andrew Janke and Mark Griffin + +.SH COPYRIGHTS +Copyright \(co 2004 by Andrew Janke and Mark Griffin diff --git a/progs/mincsample/mt19937ar.c b/progs/mincsample/mt19937ar.c new file mode 100644 --- /dev/null +++ b/progs/mincsample/mt19937ar.c @@ -0,0 +1,190 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length) +{ + int i, j, k; + init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,0x7fffffff]-interval */ +long genrand_int31(void) +{ + return (long)(genrand_int32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +double genrand_real1(void) +{ + return genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +double genrand_real2(void) +{ + return genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +double genrand_real3(void) +{ + return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void) +{ + unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} + + +/* These real versions are due to Isaku Wada, 2002/01/09 added */ + +// int main(void) +// { +// int i; +// unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4; +// init_by_array(init, length); +// printf("1000 outputs of genrand_int32()\n"); +// for (i=0; i<1000; i++) { +// printf("%10lu ", genrand_int32()); +// if (i%5==4) printf("\n"); +// } +// printf("\n1000 outputs of genrand_real2()\n"); +// for (i=0; i<1000; i++) { +// printf("%10.8f ", genrand_real2()); +// if (i%5==4) printf("\n"); +// } +// return 0; +// } diff --git a/progs/mincsample/mt19937ar.h b/progs/mincsample/mt19937ar.h new file mode 100644 --- /dev/null +++ b/progs/mincsample/mt19937ar.h @@ -0,0 +1,70 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s); + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length); + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void); + +/* generates a random number on [0,0x7fffffff]-interval */ +long genrand_int31(void); + +/* generates a random number on [0,1]-real-interval */ +double genrand_real1(void); + +/* generates a random number on [0,1)-real-interval */ +double genrand_real2(void); + +/* generates a random number on (0,1)-real-interval */ +double genrand_real3(void); + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void); +