1 rick  1.1 /*

  2 rick  1.2  *  smpl_07.c						$DRMS/proj/cookbook/

  3 rick  1.1  *
  4            *  Populates a data series consisting of variable size images with random

  5 rick  1.2  *    data from a selected distribution function

  6 rick  1.1  *  Illustrates features of the DRMS_Array and DRMS_Segment structs, and
  7            *    writing to DRMS record segments
  8            *
  9            *  Usage:

 10 rick  1.2  *    smpl_07 [ds= cols= rows= ]

 11 rick  1.1  *
 12            *  Revision history is at end of file.
 13            */
 14           

 15 rick  1.2 char *module_name = "CookbookRecipe:07";

 16 rick  1.1 char *version_id = "1.0";
 17           
 18           #include <jsoc_main.h>
 19           
 20           ModuleArgs_t module_args[] = {
 21             {ARG_STRING,	"ds",   "drms.images", "name of data series"},
 22             {ARG_INT,	"cols", "0",  "number of columns in image (default (rows)"},
 23             {ARG_INT,	"rows", "0",  "number of rows in image (default: cols)"},
 24             {ARG_NUME,	"dist", "uniform", "distribution for pseudo-randoms",
 25                 "uniform, dblexp, normal"},
 26             {ARG_INT,	"seed", "0",  "initial seed for pseudo-random number generator"},
 27             {ARG_FLAG, "v", "", "verbose flag"}, 
 28             {ARG_END}
 29           };
 30           
 31           enum dist_type {UNIFORM, DBL_EXP, NORMAL};
 32           
 33           float next_rand (int distrib) {
 34             double val0, val1, rmod, gf;
 35             static double sav;
 36             static int avail = 0;
 37 rick  1.1 
 38             switch (distrib) {
 39               case (UNIFORM):
 40                 return 2 * drand48 () - 1.0;
 41               case (DBL_EXP):
 42                 val0 = -log (drand48 ());
 43                 val1 = (drand48 () < 0.5) ? -0.5 : +0.5;
 44                 return val0 * val1;
 45               case (NORMAL):
 46                 if (avail) {
 47           	avail = 0;
 48           	return sav;
 49                 }
 50                 rmod = 2.0;
 51                 while (rmod >= 1.0) {
 52           	val0 = 2 * drand48 () - 1.0;
 53           	val1 = 2 * drand48 () - 1.0;
 54           	rmod = val0 * val0 + val1 * val1;
 55                 }
 56                 gf = sqrt (-2.0 * log (rmod) / rmod);
 57                 sav = gf * val0;
 58 rick  1.1       avail = 1;
 59                 return gf * val1;
 60               default:
 61                 return 0.0;
 62             }
 63           }
 64           
 65           void rand_fill_array (float *data, int cols, int rows, int df) {
 66             double xr;
 67             int n, col, row, pr;
 68             int colsM1 = cols - 1;
 69             double kf1 = 1.0, kf2 = 0.5, kf3 = 1.0 / 3.0, kf4 = 0.25, kf5 = 0.2;
 70           
 71             n = row = 0;
 72             data[n++] = next_rand (df);
 73             for (col = 1; col < colsM1; col++, n++) {
 74               data[n] = kf1 * data[n-1] + next_rand (df);
 75             }
 76             data[n++] = kf2 * (data[n-1] + data[n+1-cols]) + next_rand (df);
 77           
 78             row++;
 79 rick  1.1   for (; row < rows; row++) {
 80               pr = n - cols;
 81               data[n++] = kf3 * (data[n-1] + data[pr] + data[pr+1]) + next_rand (df);
 82               pr++;
 83               for (col = 1; col < colsM1; col++, n++, pr++) {
 84                 data[n] = kf4 * (data[pr-1] + data[pr] + data[pr+1] + data[n-1]) +
 85           	  next_rand (df);
 86               }
 87               data[n++] = kf5 * (data[pr-1] + data[pr] + data[pr+1] + data[n-1] + data[n+1-cols]) +
 88           	next_rand (df);
 89             }
 90           }
 91           
 92           int DoIt (void) {
 93             CmdParams_t *params = &cmdparams;
 94             DRMS_Record_t *record;
 95             DRMS_Segment_t *record_segment;
 96             DRMS_Array_t *data_array;
 97             float *data;
 98             int axes[2];
 99             int status = 0;
100 rick  1.1 

101 rick  1.4   char *series = strdup (params_get_str (params, "ds"));

102 rick  1.1   int cols = params_get_int (params, "cols");
103             int rows = params_get_int (params, "rows");
104             int distrib = params_get_int (params, "dist");
105             long int seed = params_get_int (params, "seed");
106             int verbose = params_isflagset (params, "v");
107           	/*  make sure that at least one of cols and rows has been specified  */
108             if (cols <= 0) cols = rows;
109             if (rows <= 0) rows = cols;
110             if (rows <= 0) {
111               fprintf (stderr, "Error: at least one of cols and rows must be specified\n");
112               return 1;
113             }
114           					   /*  create the output data array  */
115             axes[0] = cols;
116             axes[1] = rows;
117             data = (float *)malloc (cols * rows * sizeof (float));
118             data_array = drms_array_create (DRMS_TYPE_FLOAT, 2, axes, (void *)data,
119                 &status);
120           			       /*  fill the output array (with random data)  */
121             srand48 (seed);
122             rand_fill_array (data, cols, rows, distrib);
123 rick  1.1 
124             if (verbose) {
125               int col, row, n;
126               int mincol, maxcol, minrow, maxrow;
127               float minval, maxval;
128           
129               minval = maxval = data[0];
130               mincol = minrow = maxcol = maxrow = 0;
131               for (n = 0, row = 0; row < rows; row++) {
132                 for (col = 0; col < cols; col++, n++) {
133           	if (data[n] < minval) {
134           	  minval = data[n];
135           	  mincol = col;
136           	  minrow = row;
137           	}
138           	if (data[n] > maxval) {
139           	  maxval = data[n];
140           	  maxcol = col;
141           	  maxrow = row;
142           	}
143                 }
144 rick  1.1     }
145               printf ("min @ [%d, %d] : %f\n", mincol, minrow, minval);
146               printf ("max @ [%d, %d] : %f\n", maxcol, maxrow, maxval);
147             }
148           
149             record = drms_create_record (drms_env, series, DRMS_PERMANENT, &status);
150             if (!record) {
151               fprintf (stderr, "Error: Unable to create record in series %s\n", series);
152               fprintf (stderr, "       Does series exist and is it writeable by you?\n");
153               return 1;
154             }
155             record_segment = drms_segment_lookupnum (record, 0);
156           			/*  set the output scaling parameters for the array  */
157             data_array->bscale = 1.0;
158             data_array->bzero = 0.0;
159             status = drms_segment_write (record_segment, data_array, 0);
160             if (status) {
161               fprintf (stderr, "Error writing data to segment 0 of series %s\n", series);
162               fprintf (stderr, "      series may not have appropriate structure\n");
163               drms_free_array (data_array);
164               drms_close_record (record, DRMS_FREE_RECORD);
165 rick  1.1     return 1;
166             }
167             drms_close_record (record, DRMS_INSERT_RECORD);
168             return 0;
169           }
170           /*
171            *  Revision History
172            *
173            *  09.08.02	file created by R Bogart
174            */