Add new dependency to gui-tools
[pspp-builds.git] / src / math / order-stats.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2008, 2009 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18 #include "order-stats.h"
19 #include <libpspp/assertion.h>
20 #include <data/val-type.h>
21 #include <gl/xalloc.h>
22 #include <data/variable.h>
23 #include <data/casereader.h>
24 #include <string.h>
25
26 #if 0
27
28 #include <stdio.h>
29
30 static void
31 order_stats_dump_k1 (const struct order_stats *os)
32 {
33   struct k *k = &os->k[0];
34   printf ("K1: tc %g; c %g cc %g ccp %g\n",
35           k->tc, k->c, k->cc, k->cc_p1);
36
37 }
38
39 static void
40 order_stats_dump_k2 (const struct order_stats *os)
41 {
42   struct k *k = &os->k[1];
43   printf ("K2: tc %g; c %g cc %g ccp %g\n",
44           k->tc, k->c, k->cc, k->cc_p1);
45 }
46
47
48 void
49 order_stats_dump (const struct order_stats *os)
50 {
51   order_stats_dump_k1 (os);
52   order_stats_dump_k2 (os);
53 }
54
55 #endif
56
57 static void
58 update_k_lower (struct k *kk,
59                 double y_i, double c_i, double cc_i)
60 {
61   if ( cc_i <= kk->tc )
62     {
63       kk->cc = cc_i;
64       kk->c = c_i;
65       kk->y = y_i;
66     }
67 }
68
69
70 static void
71 update_k_upper (struct k *kk,
72                 double y_i, double c_i, double cc_i)
73 {
74   if ( cc_i > kk->tc && kk->c_p1 == 0)
75     {
76       kk->cc_p1 = cc_i;
77       kk->c_p1 = c_i;
78       kk->y_p1 = y_i;
79     }
80 }
81
82
83 static void
84 update_k_values (const struct ccase *cx, double y_i, double c_i, double cc_i,
85                  struct order_stats **os, size_t n_os)
86 {
87   int j;
88
89   for (j = 0 ; j < n_os ; ++j)
90     {
91       int k;
92       struct order_stats *tos = os[j];
93       struct statistic  *stat = (struct statistic *) tos;
94       for (k = 0 ; k < tos->n_k; ++k)
95         {
96           struct k *myk = &tos->k[k];
97           update_k_lower (myk, y_i, c_i, cc_i);
98           update_k_upper (myk, y_i, c_i, cc_i);
99         }
100
101       if ( stat->accumulate )
102         stat->accumulate (stat, cx, c_i, cc_i, y_i);
103
104       tos->cc = cc_i;
105     }
106 }
107
108
109 void
110 order_stats_accumulate (struct order_stats **os, size_t nos,
111                         struct casereader *reader,
112                         const struct variable *wv,
113                         const struct variable *var,
114                         enum mv_class exclude)
115 {
116   struct ccase *cx;
117   struct ccase *prev_cx = NULL;
118   double prev_value = -DBL_MAX;
119
120   double cc_i = 0;
121   double c_i = 0;
122
123   for (; (cx = casereader_read (reader)) != NULL; case_unref (cx))
124     {
125       const double weight = wv ? case_data (cx, wv)->f : 1.0;
126       const double this_value = case_data (cx, var)->f;
127
128       /* The casereader MUST be sorted */
129       assert (this_value >= prev_value);
130
131       if ( var_is_value_missing (var, case_data (cx, var), exclude))
132         continue;
133
134       case_unref (prev_cx);
135
136       if ( prev_value == -DBL_MAX || prev_value == this_value)
137         c_i += weight;
138
139       if ( prev_value > -DBL_MAX && this_value > prev_value)
140         {
141           update_k_values (prev_cx, prev_value, c_i, cc_i, os, nos);
142           c_i = weight;
143         }
144
145       cc_i += weight;
146       prev_value = this_value;
147       prev_cx = case_ref (cx);
148     }
149
150   update_k_values (prev_cx, prev_value, c_i, cc_i, os, nos);
151   case_unref (prev_cx);
152
153   casereader_destroy (reader);
154 }
155
156
157