Adopt use of gnulib for portability.
[pspp] / src / sample.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18    02110-1301, USA. */
19
20 #include <config.h>
21 #include <gsl/gsl_rng.h>
22 #include <limits.h>
23 #include <stdio.h>
24 #include <math.h>
25 #include "alloc.h"
26 #include "command.h"
27 #include "error.h"
28 #include "lexer.h"
29 #include "settings.h"
30 #include "str.h"
31 #include "var.h"
32
33 #include "gettext.h"
34 #define _(msgid) gettext (msgid)
35
36 #include "debug-print.h"
37
38 /* The two different types of samples. */
39 enum
40   {
41     TYPE_A_FROM_B,              /* 5 FROM 10 */
42     TYPE_FRACTION               /* 0.5 */
43   };
44
45 /* SAMPLE transformation. */
46 struct sample_trns
47   {
48     struct trns_header h;
49     int type;                   /* One of TYPE_*. */
50     int n, N;                   /* TYPE_A_FROM_B: n from N. */
51     int m, t;                   /* TYPE_A_FROM_B: # picked so far; # so far. */
52     unsigned frac;              /* TYPE_FRACTION: a fraction of UINT_MAX. */
53   };
54
55 static trns_proc_func sample_trns_proc;
56
57 int
58 cmd_sample (void)
59 {
60   struct sample_trns *trns;
61
62   int type;
63   int a, b;
64   unsigned frac;
65
66   if (!lex_force_num ())
67     return CMD_FAILURE;
68   if (!lex_is_integer ())
69     {
70       unsigned long min = gsl_rng_min (get_rng ());
71       unsigned long max = gsl_rng_max (get_rng ());
72
73       type = TYPE_FRACTION;
74       if (tokval <= 0 || tokval >= 1)
75         {
76           msg (SE, _("The sampling factor must be between 0 and 1 "
77                      "exclusive."));
78           return CMD_FAILURE;
79         }
80           
81       frac = tokval * (max - min) + min;
82       a = b = 0;
83     }
84   else
85     {
86       type = TYPE_A_FROM_B;
87       a = lex_integer ();
88       lex_get ();
89       if (!lex_force_match_id ("FROM"))
90         return CMD_FAILURE;
91       if (!lex_force_int ())
92         return CMD_FAILURE;
93       b = lex_integer ();
94       if (a >= b)
95         {
96           msg (SE, _("Cannot sample %d observations from a population of "
97                      "%d."),
98                a, b);
99           return CMD_FAILURE;
100         }
101       
102       frac = 0;
103     }
104   lex_get ();
105
106   trns = xmalloc (sizeof *trns);
107   trns->h.proc = sample_trns_proc;
108   trns->h.free = NULL;
109   trns->type = type;
110   trns->n = a;
111   trns->N = b;
112   trns->m = trns->t = 0;
113   trns->frac = frac;
114   add_transformation ((struct trns_header *) trns);
115
116   return lex_end_of_command ();
117 }
118
119 /* Executes a SAMPLE transformation. */
120 static int
121 sample_trns_proc (struct trns_header * trns, struct ccase *c UNUSED,
122                   int case_num UNUSED)
123 {
124   struct sample_trns *t = (struct sample_trns *) trns;
125   double U;
126
127   if (t->type == TYPE_FRACTION) 
128     {
129       if (gsl_rng_get (get_rng ()) <= t->frac)
130         return -1;
131       else
132         return -2;
133     }
134
135   if (t->m >= t->n)
136     return -2;
137
138   U = gsl_rng_uniform (get_rng ());
139   if ((t->N - t->t) * U >= t->n - t->m)
140     {
141       t->t++;
142       return -2;
143     }
144   else
145     {
146       t->m++;
147       t->t++;
148       return -1;
149     }
150 }