apop_asst.c 19.8 KB
 Jerome Benoit committed Aug 15, 2014 1 2  /** \file apop_asst.c The odds and ends bin.  Jerome Benoit committed Jul 10, 2015 3 Copyright (c) 2005--2007, 2010 by Ben Klemens. Licensed under the GPLv2; see COPYING. */  Jerome Benoit committed Aug 15, 2014 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29  #include "apop_internal.h" #include #include #include #ifdef _OPENMP #include #define omp_threadnum omp_get_thread_num() #else #define omp_threadnum 0 #endif extern char *apop_nul_string; //more efficient than xprintf, but a little less versatile. static void apop_tack_on(char **in, char *addme){ if (!addme) return; size_t inlen = *in? strlen(*in): 0; size_t total_len= inlen + strlen(addme); *in = realloc(*in, total_len+1); strcpy(*in+inlen, addme); } typedef int (*apop_fn_riip)(apop_data*, int, int, void*);  Jerome Benoit committed Jul 10, 2015 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 /** Join together the \c text grid of an \ref apop_data set into a single string. For example, say that we have a data set with some text: row 0 has \c "a0", \c "b0", \c "c0"; row 2 has \c "a1", \c "b1", \c "c1"; and so on. We would like to produce \code insert into tab values ('a0', 'b0', 'c0'); insert into tab values ('a1', 'b1', 'c1'); ... \endcode This could be sent to an SQL engine to copy the data to a database (but this is just an example for demonstration---use \ref apop_data_print to write to a database table). To construct this single string from the text grid, we would need to add: \li before the text, Insert into tab values ('. \li between each element on a row: ', ' \li between rows: '); \\ninsert into tab values(' \li at the tail end: ');' Thus, do the conversion via: \code char *insert_string = apop_text_paste(indata, .before="Insert into tab values ('", .between="', '", .between_cols="'); \\ninsert into tab values(', .after="');'" ); \endcode  Jerome Benoit committed Aug 15, 2014 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80  \param strings An \ref apop_data set with a grid of text to be combined into a single string \param between The text to put in between the rows of the table, such as ", ". (Default is a single space: " ") \param before The text to put at the head of the string. For the query example, this would be .before="select ". (Default: NULL) \param after The text to put at the tail of the string. For the query example, .after=" from data_table". (Default: NULL) \param between_cols The text to insert between columns of text. See below for an example (Default is set to equal .between) \param prune If you don't want to use the entire text set, you can provide a function to indicate which elements should be pruned out. Some examples: \code //Just use column 3 int is_not_col_3(apop_data *indata, int row, int col, void *ignore){ return col!=3; } //Jump over blanks as if they don't exist. int is_blank(apop_data *indata, int row, int col, void *ignore){ return strlen(indata->text[row][col])==0; } \endcode \param prune_parameter A void pointer to pass to your \c prune function.  Jerome Benoit committed Jul 10, 2015 81 82 \return A single string with the elements of the \c strings table joined as per your specification. Allocated by the function, to be freed by you if desired.  Jerome Benoit committed Aug 15, 2014 83   Jerome Benoit committed Jul 10, 2015 84 85 86 87 88  \li If the table of strings is \c NULL or has no text, the output string will have only the .before and .after parts with nothing in between. \li if apop_opts.verbose >=3, then print the pasted text to stderr. \li It is sometimes useful to use \c Apop_r and \c Apop_rs to get a view of only one or a few rows in conjunction with this function.  Jerome Benoit committed Aug 15, 2014 89   Jerome Benoit committed Jul 10, 2015 90 91 92  \li This function uses the \ref designated syntax for inputs. This sample snippet generates the SQL for a query using a list of column names (where  Jerome Benoit committed Aug 15, 2014 93 94 the query begins with select , ends with from datatab, and has commas in between each element), re-processes the same list to produce the head of an HTML  Jerome Benoit committed Jul 10, 2015 95 table, then produces the body of the table with the query result.  Jerome Benoit committed Aug 15, 2014 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115  \include sql_to_html.c */ #ifdef APOP_NO_VARIADIC char * apop_text_paste(apop_data const *strings, char *between, char *before, char *after, char *between_cols, apop_fn_riip prune, void *prune_parameter){ #else apop_varad_head(char *, apop_text_paste){ apop_data const *apop_varad_var(strings, NULL); char *apop_varad_var(between, " "); char *apop_varad_var(before, NULL); char *apop_varad_var(after, NULL); char *apop_varad_var(between_cols, between); apop_fn_riip apop_varad_var(prune, NULL); void *apop_varad_var(prune_parameter, NULL); return apop_text_paste_base(strings, between, before, after, between_cols, prune, prune_parameter); } char * apop_text_paste_base(apop_data const *strings, char *between, char *before, char *after, char *between_cols, apop_fn_riip prune, void *prune_parameter){ #endif char *prior_line=NULL, *oneline=NULL, *out = before ? strdup(before) : NULL;  Jerome Benoit committed Jul 10, 2015 116  for (int i=0; i< ((!strings || !*strings->textsize)? 0 : *strings->textsize); i++){  Jerome Benoit committed Aug 15, 2014 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139  free(oneline); oneline = NULL; for (int j=0; j< strings->textsize[1]; j++){ if (prune && !prune((apop_data*)strings, i, j, prune_parameter)) continue; apop_tack_on(&oneline, strings->text[i][j]); if (j textsize[1]-1) apop_tack_on(&oneline, between_cols); } apop_tack_on(&out, prior_line); if (prior_line && oneline) apop_tack_on(&out, between); free(prior_line); prior_line=oneline ? strdup(oneline): NULL; //if (i textsize[0]-1) apop_tack_on(&out, between); //if (oneline) apop_tack_on(&out, oneline); } apop_tack_on(&out, oneline); //the final one never got a chance to be prior_line apop_tack_on(&out, after); Apop_notify(3, "%s", out); return out; } /** Calculate \f$\sum_{n=1}^N {1\over n^s}\f$ \li There are no doubt efficient shortcuts do doing this, but I use brute force. [Though Knuth's Art of Programming v1 doesn't offer anything, which is strong indication of nonexistence.] To speed things along, I save the results so that they can just be looked up should you request the same calculation.  Jerome Benoit committed Jul 10, 2015 140 \li If \c N is zero or negative, return NaN. Notify the user if apop_opts.verbosity >=0  Jerome Benoit committed Aug 15, 2014 141 142 143 144 145  For example: \include test_harmonic.c */  Jerome Benoit committed Jul 10, 2015 146 long double apop_generalized_harmonic(int N, double s){  Jerome Benoit committed Aug 15, 2014 147 148 149 150 151 152 153 154 155 /* Each row in the saved-results structure is an \f$s\f$, and each column is \f$1\dots n\f$, up to the largest \f$n\f$ calculated to date. When reading the code, remember that the zeroth element holds the value for N=1, and so on. */ Apop_stopif(N<=0, return GSL_NAN, 0, "N is %i, but must be greater than 0.", N); static double * eses = NULL; static int * lengths= NULL; static int count = 0;  Jerome Benoit committed Jul 10, 2015 156  static long double ** precalced=NULL;  Jerome Benoit committed Aug 15, 2014 157 158 159 160 161 162 163 164 165  int old_len, i; OMP_critical(generalized_harmonic) { //Due to memoization, this can't parallelize. for (i=0; i< count; i++) if (eses == NULL || eses[i] == s) break; if (i == count){ //you need to build the vector from scratch. count ++; i = count - 1;  Jerome Benoit committed Jul 10, 2015 166  precalced = realloc(precalced, sizeof (long double*) * count);  Jerome Benoit committed Aug 15, 2014 167 168  lengths = realloc(lengths, sizeof (int*) * count); eses = realloc(eses, sizeof (double) * count);  Jerome Benoit committed Jul 10, 2015 169  precalced[i] = malloc(sizeof(long double) * N);  Jerome Benoit committed Aug 15, 2014 170 171 172 173 174 175 176 177 178  lengths[i] = N; eses[i] = s; precalced[i][0] = 1; old_len = 1; } else { //then you found it. old_len = lengths[i]; } if (N-1 >= old_len){ //It's there, but you need to extend what you have.  Jerome Benoit committed Jul 10, 2015 179  precalced[i] = realloc(precalced[i], sizeof(long double) * N);  Jerome Benoit committed Aug 15, 2014 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221  for (int j = old_len; jprintf-style arguments. E.g., \code char filenames[] = "apop_asst.c apop_asst.o" apop_system("ls -l %s", filenames); \endcode \return The return value of the \c system() call. */ int apop_system(const char *fmt, ...){ char *q; va_list argp; va_start(argp, fmt); Apop_stopif(vasprintf(&q, fmt, argp)==-1, return -1, 0, "Trouble writing to a string."); va_end(argp); int out = system(q); free(q); return out; } static int count_parens(const char *string){ int out = 0; int last_was_backslash = 0; for(const char *step =string; *step !='\0'; step++){ if (*step == '\\' && !last_was_backslash){ last_was_backslash = 1; continue; } if (*step == ')' && !last_was_backslash) out++; last_was_backslash = 0; } return out; }  Jerome Benoit committed Jul 10, 2015 222 223 224 /** Extract subsets from a string via regular expressions. This function takes a regular expression and repeatedly applies it to an input string. It returns the count of matches, and optionally returns the matches themselves organized into the \c text grid of an \ref apop_data set.  Jerome Benoit committed Aug 15, 2014 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246  \li There are three common flavors of regular expression: Basic, Extended, and Perl-compatible (BRE, ERE, PCRE). I use EREs, as per the specs of your C library, which should match POSIX's ERE specification. For example, "p.val" will match "P value", "p.value", "p values" (and even "tempeval", so be careful). If you give a non-\c NULL address in which to place a table of paren-delimited substrings, I'll return them as a row in the text element of the returned \ref apop_data set. I'll return all the matches, filling the first row with substrings from the first application of your regex, then filling the next row with another set of matches (if any), and so on to the end of the string. Useful when parsing a list of items, for example. \param string The string to search (no default) \param regex The regular expression (no default) \param substrings Parens in the regex indicate that I should return matching substrings. Give me the _address_ of an \ref apop_data* set, and I will allocate and fill the text portion with matches. Default= \c NULL, meaning do not return substrings (even if parens exist in the regex). If no match, return an empty \ref apop_data set, so output->textsize[0]==0. \param use_case Should I be case sensitive, \c 'y' or \c 'n'? (default = \c 'n', which is not the POSIX default.) \return Count of matches found. 0 == no match. \c substrings may be allocated and filled if needed. \li If apop_opts.stop_on_warning='n' returns -1 on error (e.g., regex \c NULL or didn't compile). \li If strings==NULL, I return 0---no match---and if \c substrings is provided, set it to \c NULL. \li Here is the test function. Notice that the substring-pulling  Jerome Benoit committed Jul 10, 2015 247 248 function call passes \c &subs, not plain \c subs.  Jerome Benoit committed Aug 15, 2014 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310  \include test_regex.c \li Each set of matches will be one row of the output data. E.g., given the regex ([A-Za-z])([0-9]), the column zero of outdata will hold letters, and column one will hold numbers. Use \ref apop_data_transpose to reverse this so that the letters are in outdata->text[0] and numbers in outdata->text[1]. */ #ifdef APOP_NO_VARIADIC int apop_regex(const char *string, const char* regex, apop_data **substrings, const char use_case){ #else apop_varad_head(int, apop_regex){ const char * apop_varad_var(string, NULL); apop_data **apop_varad_var(substrings, NULL); if (!string) { if (substrings) *substrings=NULL; return 0; } const char * apop_varad_var(regex, NULL); Apop_stopif(!regex, return -1, 0, "You gave me a NULL regex."); const char apop_varad_var(use_case, 'n'); return apop_regex_base(string, regex, substrings, use_case); } int apop_regex_base(const char *string, const char* regex, apop_data **substrings, const char use_case){ #endif regex_t re; int matchcount=count_parens(regex); int found, found_ct=0; regmatch_t result[matchcount+1]; int compiled_ok = !regcomp(&re, regex, REG_EXTENDED + (use_case=='y' ? 0 : REG_ICASE) + (substrings ? 0 : REG_NOSUB) ); Apop_stopif(!compiled_ok, return -1, 0, "This regular expression didn't compile: \"%s\"", regex); int matchrow = 0; if (substrings) *substrings = apop_data_alloc(); do { found_ct+= found = !regexec(&re, string, matchcount+1, result, matchrow ? REG_NOTBOL : 0); if (substrings && found){ *substrings = apop_text_alloc(*substrings, matchrow+1, matchcount); //match zero is the whole string; ignore. for (int i=0; i< matchcount; i++){ if (result[i+1].rm_eo > 0){//GNU peculiarity: match-to-empty marked with -1. int length_of_match = result[i+1].rm_eo - result[i+1].rm_so; if ((*substrings)->text[matchrow][i] != apop_nul_string) free((*substrings)->text[matchrow][i]); (*substrings)->text[matchrow][i] = malloc(strlen(string)+1); memcpy((*substrings)->text[matchrow][i], string + result[i+1].rm_so, length_of_match); (*substrings)->text[matchrow][i][length_of_match] = '\0'; } //else matches nothing; apop_text_alloc already made this cell this NULL. } string += result[0].rm_eo; //end of whole match; matchrow++; } } while (substrings && found && string[0]!='\0'); regfree(&re); return found_ct; } /** RNG from a Generalized Hypergeometric type B3. Devroye uses this as the base for many of his distribution-generators, including the Waring.  Jerome Benoit committed Jul 10, 2015 311 \li If one of the inputs is <=0, error; return NaN and print a warning.  Jerome Benoit committed Aug 15, 2014 312 313 */ //Header in stats.h double apop_rng_GHgB3(gsl_rng * r, double* a){  Jerome Benoit committed Jul 10, 2015 314  Apop_stopif(!((a[0]>0) && (a[1] > 0) && (a[2] > 0)), return NAN, 0, "all inputs must be positive.");  Jerome Benoit committed Aug 15, 2014 315 316 317 318 319 320 321 322 323  double aa = gsl_ran_gamma(r, a[0], 1), b = gsl_ran_gamma(r, a[1], 1), c = gsl_ran_gamma(r, a[2], 1); int p = gsl_ran_poisson(r, aa*b/c); return p; } /** The Beta distribution is useful for modeling because it is bounded between zero and one, and can be either unimodal (if the variance is low) or bimodal (if the variance is high), and can have either a slant toward the bottom or top of the range (depending on the mean).  Jerome Benoit committed Jul 10, 2015 324 The distribution has two parameters, typically named \f$\alpha\f$ and \f$\beta\f$, which can be difficult to interpret. However, there is a one-to-one mapping between (alpha, beta) pairs and (mean, variance) pairs. Since we have good intuition about the meaning of means and variances, this function takes in a mean and variance, calculates alpha and beta behind the scenes, and returns the appropriate Beta distribution.  Jerome Benoit committed Aug 15, 2014 325 326 327 328 329 330 331 332  \param m The mean the Beta distribution should have. Notice that m is in [0,1]. \param v The variance which the Beta distribution should have. It is in (0, 1/12), where (1/12) is the variance of a Uniform(0,1) distribution. Funny things happen with variance near 1/12 and mean far from 1/2.  Jerome Benoit committed Jul 10, 2015 333 334 335 \return Returns an \ref apop_model produced by copying the \c apop_beta model and setting its parameters appropriately.  Jerome Benoit committed Aug 15, 2014 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 \exception out->error=='r' Range error: mean is not within [0, 1]. */ apop_model *apop_beta_from_mean_var(double m, double v){ Apop_stopif(m>=1|| m<=0, apop_model *out = apop_model_copy(apop_beta); out->error='r'; return out, 0, "You asked for a beta distribution " "with mean %g, but the mean of the beta will always " "be strictly between zero and one.", m); double k = (m * (1- m)/ v) -1; double alpha = m*k; double beta = k * (1-m); return apop_model_set_parameters(apop_beta, alpha, beta); } /** \def apop_rng_get_thread The \c gsl_rng is not itself thread-safe, in the sense that it can not be used  Jerome Benoit committed Jul 10, 2015 352 353 simultaneously by multiple threads. However, if each thread has its own \c gsl_rng, then each will safely operate independently.  Jerome Benoit committed Aug 15, 2014 354 355 356 357 358 359  Thus, Apophenia keeps an internal store of RNGs for use by threaded functions. If the input to this function, \c thread, is greater than any previous input, then the array of gsl_rngs is extended to length \c thread, and each element extended using ++apop_opts.rng_seed (i.e., the seed is incremented before use).  Jerome Benoit committed Jul 10, 2015 360 361 362 363 364 365 366 This function can be used anywhere a \c gsl_rng would be used. \param thread_in The number of the RNG to retrieve, starting at zero (which is how OpenMP numbers its threads). If -1, I'll look up the current thread (via \c omp_get_thread_num) for you. See \ref threading for additional notes. In most cases, you want to use apop_rng_get_thread(-1).  Jerome Benoit committed Aug 15, 2014 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416  \return The appropriate RNG, initialized if necessary. \hideinitializer */ gsl_rng *apop_rng_get_thread_base(int thread){ static gsl_rng **rngs; static int rng_ct = -1; if (thread==-1){ #ifdef OpenMP thread = omp_get_thread_num(); #else thread = 0; #endif } OMP_critical(rng_get_thread) if (thread > rng_ct) { rngs = realloc(rngs, sizeof(gsl_rng*)*(thread+1)); for (int i=rng_ct+1; i<= thread; i++) rngs[i] = apop_rng_alloc(++apop_opts.rng_seed); rng_ct = thread; } return rngs[thread]; } /** Make a set of random draws from a model and write them to an \ref apop_data set. \param model The model from which draws will be made. Must already be prepared and/or estimated. \param count The number of draws to make. If \c draw_matrix is not \c NULL, then this is ignored and count=draw_matrix->matrix->size1. default=1000. \param draws If not \c NULL, a pre-allocated data set whose \c matrix element will be filled with draws. \return An \ref apop_data set with the matrix filled with \c size draws. If draw_matrix!=NULL, then return a pointer to it. \exception out->error=='m' Input model isn't good for making draws: it is \c NULL, or m->dsize=0. \exception out->error=='s' You gave me a \c draws matrix, but its size is less than the size of a single draw from the data, model->dsize. \exception out->error=='d' Trouble drawing from the distribution for at least one row. That row is set to all \c NAN. \li Prints a warning if you send in a non-NULL apop_data set, but its \c matrix element is \c NULL, when apop_opts.verbose>=1. \li See also \ref apop_draw, which makes a single draw. \li Random numbers are generated using RNGs from \ref apop_rng_get_thread, qv. Here is a two-line program to draw a different set of ten Standard Normals on every run (provided runs are more than a second apart): \include draw_some_normals.c  Jerome Benoit committed Jul 10, 2015 417 418 419  \li This function uses the \ref designated syntax for inputs. */  Jerome Benoit committed Aug 15, 2014 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 #ifdef APOP_NO_VARIADIC apop_data * apop_model_draws(apop_model *model, int count, apop_data *draws){ #else apop_varad_head(apop_data *, apop_model_draws){ apop_model * apop_varad_var(model, NULL); Apop_stopif(!model, apop_return_data_error(n), 0, "Input model is NULL."); Apop_stopif(!model->dsize, apop_return_data_error(n), 0, "Input model has dsize==0."); apop_data * apop_varad_var(draws, NULL); int apop_varad_var(count, 1000); if (draws) { Apop_stopif(!draws->matrix, draws->error='m'; return draws, 1, "Input data set's matrix is NULL."); Apop_stopif((int)draws->matrix->size2 < model->dsize, draws->error='s'; draws->error='m'; return draws, 1, "Input data set's matrix column count is less than model->dsize."); count = draws->matrix->size1; } else Apop_stopif(model->dsize<=0, apop_return_data_error(n), 0, "model->dsize<=0, so I don't know the size of matrix to allocate."); return apop_model_draws_base(model, count, draws); } apop_data * apop_model_draws_base(apop_model *model, int count, apop_data *draws){ #endif apop_data *out = draws ? draws : apop_data_alloc(count, model->dsize); OMP_for (int i=0; i< count; i++){ apop_data *onerow = Apop_r(out, i); Apop_stopif(apop_draw(onerow->matrix->data, apop_rng_get_thread(omp_threadnum), model), gsl_matrix_set_all(onerow->matrix, GSL_NAN); out->error='d', 0, "Trouble drawing for row %i. " "I set it to all NANs and set out->error='d'.", i); } return out; }