tesseract  5.0.0-alpha-619-ge9db
cluster.cpp
Go to the documentation of this file.
1 /******************************************************************************
2  ** Filename: cluster.cpp
3  ** Purpose: Routines for clustering points in N-D space
4  ** Author: Dan Johnson
5  **
6  ** (c) Copyright Hewlett-Packard Company, 1988.
7  ** Licensed under the Apache License, Version 2.0 (the "License");
8  ** you may not use this file except in compliance with the License.
9  ** You may obtain a copy of the License at
10  ** http://www.apache.org/licenses/LICENSE-2.0
11  ** Unless required by applicable law or agreed to in writing, software
12  ** distributed under the License is distributed on an "AS IS" BASIS,
13  ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  ** See the License for the specific language governing permissions and
15  ** limitations under the License.
16  *****************************************************************************/
17 
18 #define _USE_MATH_DEFINES // for M_PI
19 #include <cfloat> // for FLT_MAX
20 #include <cmath> // for M_PI
21 #include <vector> // for std::vector
22 
23 #include "cluster.h"
24 #include "emalloc.h"
25 #include "genericheap.h"
26 #include <tesseract/helpers.h>
27 #include "kdpair.h"
28 #include "matrix.h"
29 #include "tprintf.h"
30 
31 #define HOTELLING 1 // If true use Hotelling's test to decide where to split.
32 #define FTABLE_X 10 // Size of FTable.
33 #define FTABLE_Y 100 // Size of FTable.
34 
35 // Table of values approximating the cumulative F-distribution for a confidence of 1%.
36 const double FTable[FTABLE_Y][FTABLE_X] = {
37  {4052.19, 4999.52, 5403.34, 5624.62, 5763.65, 5858.97, 5928.33, 5981.10, 6022.50, 6055.85,},
38  {98.502, 99.000, 99.166, 99.249, 99.300, 99.333, 99.356, 99.374, 99.388, 99.399,},
39  {34.116, 30.816, 29.457, 28.710, 28.237, 27.911, 27.672, 27.489, 27.345, 27.229,},
40  {21.198, 18.000, 16.694, 15.977, 15.522, 15.207, 14.976, 14.799, 14.659, 14.546,},
41  {16.258, 13.274, 12.060, 11.392, 10.967, 10.672, 10.456, 10.289, 10.158, 10.051,},
42  {13.745, 10.925, 9.780, 9.148, 8.746, 8.466, 8.260, 8.102, 7.976, 7.874,},
43  {12.246, 9.547, 8.451, 7.847, 7.460, 7.191, 6.993, 6.840, 6.719, 6.620,},
44  {11.259, 8.649, 7.591, 7.006, 6.632, 6.371, 6.178, 6.029, 5.911, 5.814,},
45  {10.561, 8.022, 6.992, 6.422, 6.057, 5.802, 5.613, 5.467, 5.351, 5.257,},
46  {10.044, 7.559, 6.552, 5.994, 5.636, 5.386, 5.200, 5.057, 4.942, 4.849,},
47  { 9.646, 7.206, 6.217, 5.668, 5.316, 5.069, 4.886, 4.744, 4.632, 4.539,},
48  { 9.330, 6.927, 5.953, 5.412, 5.064, 4.821, 4.640, 4.499, 4.388, 4.296,},
49  { 9.074, 6.701, 5.739, 5.205, 4.862, 4.620, 4.441, 4.302, 4.191, 4.100,},
50  { 8.862, 6.515, 5.564, 5.035, 4.695, 4.456, 4.278, 4.140, 4.030, 3.939,},
51  { 8.683, 6.359, 5.417, 4.893, 4.556, 4.318, 4.142, 4.004, 3.895, 3.805,},
52  { 8.531, 6.226, 5.292, 4.773, 4.437, 4.202, 4.026, 3.890, 3.780, 3.691,},
53  { 8.400, 6.112, 5.185, 4.669, 4.336, 4.102, 3.927, 3.791, 3.682, 3.593,},
54  { 8.285, 6.013, 5.092, 4.579, 4.248, 4.015, 3.841, 3.705, 3.597, 3.508,},
55  { 8.185, 5.926, 5.010, 4.500, 4.171, 3.939, 3.765, 3.631, 3.523, 3.434,},
56  { 8.096, 5.849, 4.938, 4.431, 4.103, 3.871, 3.699, 3.564, 3.457, 3.368,},
57  { 8.017, 5.780, 4.874, 4.369, 4.042, 3.812, 3.640, 3.506, 3.398, 3.310,},
58  { 7.945, 5.719, 4.817, 4.313, 3.988, 3.758, 3.587, 3.453, 3.346, 3.258,},
59  { 7.881, 5.664, 4.765, 4.264, 3.939, 3.710, 3.539, 3.406, 3.299, 3.211,},
60  { 7.823, 5.614, 4.718, 4.218, 3.895, 3.667, 3.496, 3.363, 3.256, 3.168,},
61  { 7.770, 5.568, 4.675, 4.177, 3.855, 3.627, 3.457, 3.324, 3.217, 3.129,},
62  { 7.721, 5.526, 4.637, 4.140, 3.818, 3.591, 3.421, 3.288, 3.182, 3.094,},
63  { 7.677, 5.488, 4.601, 4.106, 3.785, 3.558, 3.388, 3.256, 3.149, 3.062,},
64  { 7.636, 5.453, 4.568, 4.074, 3.754, 3.528, 3.358, 3.226, 3.120, 3.032,},
65  { 7.598, 5.420, 4.538, 4.045, 3.725, 3.499, 3.330, 3.198, 3.092, 3.005,},
66  { 7.562, 5.390, 4.510, 4.018, 3.699, 3.473, 3.305, 3.173, 3.067, 2.979,},
67  { 7.530, 5.362, 4.484, 3.993, 3.675, 3.449, 3.281, 3.149, 3.043, 2.955,},
68  { 7.499, 5.336, 4.459, 3.969, 3.652, 3.427, 3.258, 3.127, 3.021, 2.934,},
69  { 7.471, 5.312, 4.437, 3.948, 3.630, 3.406, 3.238, 3.106, 3.000, 2.913,},
70  { 7.444, 5.289, 4.416, 3.927, 3.611, 3.386, 3.218, 3.087, 2.981, 2.894,},
71  { 7.419, 5.268, 4.396, 3.908, 3.592, 3.368, 3.200, 3.069, 2.963, 2.876,},
72  { 7.396, 5.248, 4.377, 3.890, 3.574, 3.351, 3.183, 3.052, 2.946, 2.859,},
73  { 7.373, 5.229, 4.360, 3.873, 3.558, 3.334, 3.167, 3.036, 2.930, 2.843,},
74  { 7.353, 5.211, 4.343, 3.858, 3.542, 3.319, 3.152, 3.021, 2.915, 2.828,},
75  { 7.333, 5.194, 4.327, 3.843, 3.528, 3.305, 3.137, 3.006, 2.901, 2.814,},
76  { 7.314, 5.179, 4.313, 3.828, 3.514, 3.291, 3.124, 2.993, 2.888, 2.801,},
77  { 7.296, 5.163, 4.299, 3.815, 3.501, 3.278, 3.111, 2.980, 2.875, 2.788,},
78  { 7.280, 5.149, 4.285, 3.802, 3.488, 3.266, 3.099, 2.968, 2.863, 2.776,},
79  { 7.264, 5.136, 4.273, 3.790, 3.476, 3.254, 3.087, 2.957, 2.851, 2.764,},
80  { 7.248, 5.123, 4.261, 3.778, 3.465, 3.243, 3.076, 2.946, 2.840, 2.754,},
81  { 7.234, 5.110, 4.249, 3.767, 3.454, 3.232, 3.066, 2.935, 2.830, 2.743,},
82  { 7.220, 5.099, 4.238, 3.757, 3.444, 3.222, 3.056, 2.925, 2.820, 2.733,},
83  { 7.207, 5.087, 4.228, 3.747, 3.434, 3.213, 3.046, 2.916, 2.811, 2.724,},
84  { 7.194, 5.077, 4.218, 3.737, 3.425, 3.204, 3.037, 2.907, 2.802, 2.715,},
85  { 7.182, 5.066, 4.208, 3.728, 3.416, 3.195, 3.028, 2.898, 2.793, 2.706,},
86  { 7.171, 5.057, 4.199, 3.720, 3.408, 3.186, 3.020, 2.890, 2.785, 2.698,},
87  { 7.159, 5.047, 4.191, 3.711, 3.400, 3.178, 3.012, 2.882, 2.777, 2.690,},
88  { 7.149, 5.038, 4.182, 3.703, 3.392, 3.171, 3.005, 2.874, 2.769, 2.683,},
89  { 7.139, 5.030, 4.174, 3.695, 3.384, 3.163, 2.997, 2.867, 2.762, 2.675,},
90  { 7.129, 5.021, 4.167, 3.688, 3.377, 3.156, 2.990, 2.860, 2.755, 2.668,},
91  { 7.119, 5.013, 4.159, 3.681, 3.370, 3.149, 2.983, 2.853, 2.748, 2.662,},
92  { 7.110, 5.006, 4.152, 3.674, 3.363, 3.143, 2.977, 2.847, 2.742, 2.655,},
93  { 7.102, 4.998, 4.145, 3.667, 3.357, 3.136, 2.971, 2.841, 2.736, 2.649,},
94  { 7.093, 4.991, 4.138, 3.661, 3.351, 3.130, 2.965, 2.835, 2.730, 2.643,},
95  { 7.085, 4.984, 4.132, 3.655, 3.345, 3.124, 2.959, 2.829, 2.724, 2.637,},
96  { 7.077, 4.977, 4.126, 3.649, 3.339, 3.119, 2.953, 2.823, 2.718, 2.632,},
97  { 7.070, 4.971, 4.120, 3.643, 3.333, 3.113, 2.948, 2.818, 2.713, 2.626,},
98  { 7.062, 4.965, 4.114, 3.638, 3.328, 3.108, 2.942, 2.813, 2.708, 2.621,},
99  { 7.055, 4.959, 4.109, 3.632, 3.323, 3.103, 2.937, 2.808, 2.703, 2.616,},
100  { 7.048, 4.953, 4.103, 3.627, 3.318, 3.098, 2.932, 2.803, 2.698, 2.611,},
101  { 7.042, 4.947, 4.098, 3.622, 3.313, 3.093, 2.928, 2.798, 2.693, 2.607,},
102  { 7.035, 4.942, 4.093, 3.618, 3.308, 3.088, 2.923, 2.793, 2.689, 2.602,},
103  { 7.029, 4.937, 4.088, 3.613, 3.304, 3.084, 2.919, 2.789, 2.684, 2.598,},
104  { 7.023, 4.932, 4.083, 3.608, 3.299, 3.080, 2.914, 2.785, 2.680, 2.593,},
105  { 7.017, 4.927, 4.079, 3.604, 3.295, 3.075, 2.910, 2.781, 2.676, 2.589,},
106  { 7.011, 4.922, 4.074, 3.600, 3.291, 3.071, 2.906, 2.777, 2.672, 2.585,},
107  { 7.006, 4.917, 4.070, 3.596, 3.287, 3.067, 2.902, 2.773, 2.668, 2.581,},
108  { 7.001, 4.913, 4.066, 3.591, 3.283, 3.063, 2.898, 2.769, 2.664, 2.578,},
109  { 6.995, 4.908, 4.062, 3.588, 3.279, 3.060, 2.895, 2.765, 2.660, 2.574,},
110  { 6.990, 4.904, 4.058, 3.584, 3.275, 3.056, 2.891, 2.762, 2.657, 2.570,},
111  { 6.985, 4.900, 4.054, 3.580, 3.272, 3.052, 2.887, 2.758, 2.653, 2.567,},
112  { 6.981, 4.896, 4.050, 3.577, 3.268, 3.049, 2.884, 2.755, 2.650, 2.563,},
113  { 6.976, 4.892, 4.047, 3.573, 3.265, 3.046, 2.881, 2.751, 2.647, 2.560,},
114  { 6.971, 4.888, 4.043, 3.570, 3.261, 3.042, 2.877, 2.748, 2.644, 2.557,},
115  { 6.967, 4.884, 4.040, 3.566, 3.258, 3.039, 2.874, 2.745, 2.640, 2.554,},
116  { 6.963, 4.881, 4.036, 3.563, 3.255, 3.036, 2.871, 2.742, 2.637, 2.551,},
117  { 6.958, 4.877, 4.033, 3.560, 3.252, 3.033, 2.868, 2.739, 2.634, 2.548,},
118  { 6.954, 4.874, 4.030, 3.557, 3.249, 3.030, 2.865, 2.736, 2.632, 2.545,},
119  { 6.950, 4.870, 4.027, 3.554, 3.246, 3.027, 2.863, 2.733, 2.629, 2.542,},
120  { 6.947, 4.867, 4.024, 3.551, 3.243, 3.025, 2.860, 2.731, 2.626, 2.539,},
121  { 6.943, 4.864, 4.021, 3.548, 3.240, 3.022, 2.857, 2.728, 2.623, 2.537,},
122  { 6.939, 4.861, 4.018, 3.545, 3.238, 3.019, 2.854, 2.725, 2.621, 2.534,},
123  { 6.935, 4.858, 4.015, 3.543, 3.235, 3.017, 2.852, 2.723, 2.618, 2.532,},
124  { 6.932, 4.855, 4.012, 3.540, 3.233, 3.014, 2.849, 2.720, 2.616, 2.529,},
125  { 6.928, 4.852, 4.010, 3.538, 3.230, 3.012, 2.847, 2.718, 2.613, 2.527,},
126  { 6.925, 4.849, 4.007, 3.535, 3.228, 3.009, 2.845, 2.715, 2.611, 2.524,},
127  { 6.922, 4.846, 4.004, 3.533, 3.225, 3.007, 2.842, 2.713, 2.609, 2.522,},
128  { 6.919, 4.844, 4.002, 3.530, 3.223, 3.004, 2.840, 2.711, 2.606, 2.520,},
129  { 6.915, 4.841, 3.999, 3.528, 3.221, 3.002, 2.838, 2.709, 2.604, 2.518,},
130  { 6.912, 4.838, 3.997, 3.525, 3.218, 3.000, 2.835, 2.706, 2.602, 2.515,},
131  { 6.909, 4.836, 3.995, 3.523, 3.216, 2.998, 2.833, 2.704, 2.600, 2.513,},
132  { 6.906, 4.833, 3.992, 3.521, 3.214, 2.996, 2.831, 2.702, 2.598, 2.511,},
133  { 6.904, 4.831, 3.990, 3.519, 3.212, 2.994, 2.829, 2.700, 2.596, 2.509,},
134  { 6.901, 4.829, 3.988, 3.517, 3.210, 2.992, 2.827, 2.698, 2.594, 2.507,},
135  { 6.898, 4.826, 3.986, 3.515, 3.208, 2.990, 2.825, 2.696, 2.592, 2.505,},
136  { 6.895, 4.824, 3.984, 3.513, 3.206, 2.988, 2.823, 2.694, 2.590, 2.503}
137 };
138 
143 #define MINVARIANCE 0.0004
144 
151 #define MINSAMPLESPERBUCKET 5
152 #define MINSAMPLES (MINBUCKETS * MINSAMPLESPERBUCKET)
153 #define MINSAMPLESNEEDED 1
154 
161 #define BUCKETTABLESIZE 1024
162 #define NORMALEXTENT 3.0
163 
164 struct TEMPCLUSTER {
167 };
168 
171 
172 struct STATISTICS {
173  float AvgVariance;
174  float *CoVariance;
175  float *Min; // largest negative distance from the mean
176  float *Max; // largest positive distance from the mean
177 };
178 
179 struct BUCKETS {
180  DISTRIBUTION Distribution; // distribution being tested for
181  uint32_t SampleCount; // # of samples in histogram
182  double Confidence; // confidence level of test
183  double ChiSquared; // test threshold
184  uint16_t NumberOfBuckets; // number of cells in histogram
185  uint16_t Bucket[BUCKETTABLESIZE]; // mapping to histogram buckets
186  uint32_t *Count; // frequency of occurrence histogram
187  float *ExpectedCount; // expected histogram
188 };
189 
190 struct CHISTRUCT{
192  double Alpha;
193  double ChiSquared;
194 };
195 
196 // For use with KDWalk / MakePotentialClusters
198  ClusterHeap *heap; // heap used to hold temp clusters, "best" on top
199  TEMPCLUSTER *candidates; // array of potential clusters
200  KDTREE *tree; // kd-tree to be searched for neighbors
201  int32_t next; // next candidate to be used
202 };
203 
204 using DENSITYFUNC = double (*)(int32_t);
205 using SOLVEFUNC = double (*)(CHISTRUCT*, double);
206 
207 #define Odd(N) ((N)%2)
208 #define Mirror(N,R) ((R) - (N) - 1)
209 #define Abs(N) (((N) < 0) ? (-(N)) : (N))
210 
211 //--------------Global Data Definitions and Declarations----------------------
219 #define SqrtOf2Pi 2.506628275
220 static const double kNormalStdDev = BUCKETTABLESIZE / (2.0 * NORMALEXTENT);
221 static const double kNormalVariance =
223 static const double kNormalMagnitude =
224  (2.0 * NORMALEXTENT) / (SqrtOf2Pi * BUCKETTABLESIZE);
225 static const double kNormalMean = BUCKETTABLESIZE / 2;
226 
229 #define LOOKUPTABLESIZE 8
230 #define MAXDEGREESOFFREEDOM MAXBUCKETS
231 
232 static const uint32_t kCountTable[LOOKUPTABLESIZE] = {
233  MINSAMPLES, 200, 400, 600, 800, 1000, 1500, 2000
234 }; // number of samples
235 
236 static const uint16_t kBucketsTable[LOOKUPTABLESIZE] = {
237  MINBUCKETS, 16, 20, 24, 27, 30, 35, MAXBUCKETS
238 }; // number of buckets
239 
240 /*-------------------------------------------------------------------------
241  Private Function Prototypes
242 --------------------------------------------------------------------------*/
243 static void CreateClusterTree(CLUSTERER* Clusterer);
244 
245 static void MakePotentialClusters(ClusteringContext* context, CLUSTER* Cluster,
246  int32_t Level);
247 
248 static CLUSTER* FindNearestNeighbor(KDTREE*Tree, CLUSTER* Cluster,
249  float* Distance);
250 
251 static CLUSTER* MakeNewCluster(CLUSTERER* Clusterer, TEMPCLUSTER* TempCluster);
252 
253 static void ComputePrototypes(CLUSTERER* Clusterer, CLUSTERCONFIG* Config);
254 
255 static PROTOTYPE* MakePrototype(CLUSTERER* Clusterer, CLUSTERCONFIG* Config,
256  CLUSTER* Cluster);
257 
258 static PROTOTYPE* MakeDegenerateProto(uint16_t N,
259  CLUSTER* Cluster, STATISTICS* Statistics,
260  PROTOSTYLE Style, int32_t MinSamples);
261 
262 static PROTOTYPE* TestEllipticalProto(CLUSTERER* Clusterer,
263  CLUSTERCONFIG* Config, CLUSTER* Cluster,
264  STATISTICS* Statistics);
265 
266 static PROTOTYPE* MakeSphericalProto(CLUSTERER* Clusterer,
267  CLUSTER* Cluster, STATISTICS* Statistics,
268  BUCKETS* Buckets);
269 
270 static PROTOTYPE* MakeEllipticalProto(CLUSTERER* Clusterer,
271  CLUSTER* Cluster, STATISTICS* Statistics,
272  BUCKETS* Buckets);
273 
274 static PROTOTYPE* MakeMixedProto(CLUSTERER* Clusterer,
275  CLUSTER* Cluster, STATISTICS* Statistics,
276  BUCKETS* NormalBuckets, double Confidence);
277 
278 static void MakeDimRandom(uint16_t i, PROTOTYPE* Proto, PARAM_DESC* ParamDesc);
279 
280 static void MakeDimUniform(uint16_t i, PROTOTYPE* Proto, STATISTICS* Statistics);
281 
282 static STATISTICS* ComputeStatistics(int16_t N, PARAM_DESC ParamDesc[],
283  CLUSTER* Cluster);
284 
285 static PROTOTYPE* NewSphericalProto(uint16_t N, CLUSTER* Cluster,
286  STATISTICS* Statistics);
287 
288 static PROTOTYPE* NewEllipticalProto(int16_t N, CLUSTER* Cluster,
289  STATISTICS* Statistics);
290 
291 static PROTOTYPE* NewMixedProto(int16_t N, CLUSTER *Cluster, STATISTICS *Statistics);
292 
293 static PROTOTYPE* NewSimpleProto(int16_t N, CLUSTER *Cluster);
294 
295 static bool Independent(PARAM_DESC* ParamDesc,
296  int16_t N, float* CoVariance, float Independence);
297 
298 static BUCKETS *GetBuckets(CLUSTERER* clusterer,
299  DISTRIBUTION Distribution,
300  uint32_t SampleCount,
301  double Confidence);
302 
303 static BUCKETS *MakeBuckets(DISTRIBUTION Distribution,
304  uint32_t SampleCount,
305  double Confidence);
306 
307 static uint16_t OptimumNumberOfBuckets(uint32_t SampleCount);
308 
309 static double ComputeChiSquared(uint16_t DegreesOfFreedom, double Alpha);
310 
311 static double NormalDensity(int32_t x);
312 
313 static double UniformDensity(int32_t x);
314 
315 static double Integral(double f1, double f2, double Dx);
316 
317 static void FillBuckets(BUCKETS *Buckets,
318  CLUSTER *Cluster,
319  uint16_t Dim,
320  PARAM_DESC *ParamDesc,
321  float Mean,
322  float StdDev);
323 
324 static uint16_t NormalBucket(PARAM_DESC *ParamDesc,
325  float x,
326  float Mean,
327  float StdDev);
328 
329 static uint16_t UniformBucket(PARAM_DESC *ParamDesc,
330  float x,
331  float Mean,
332  float StdDev);
333 
334 static bool DistributionOK(BUCKETS* Buckets);
335 
336 static void FreeStatistics(STATISTICS *Statistics);
337 
338 static void FreeBuckets(BUCKETS *Buckets);
339 
340 static void FreeCluster(CLUSTER *Cluster);
341 
342 static uint16_t DegreesOfFreedom(DISTRIBUTION Distribution, uint16_t HistogramBuckets);
343 
344 static void AdjustBuckets(BUCKETS *Buckets, uint32_t NewSampleCount);
345 
346 static void InitBuckets(BUCKETS *Buckets);
347 
348 static int AlphaMatch(void *arg1, // CHISTRUCT *ChiStruct,
349  void *arg2); // CHISTRUCT *SearchKey);
350 
351 static CHISTRUCT *NewChiStruct(uint16_t DegreesOfFreedom, double Alpha);
352 
353 static double Solve(SOLVEFUNC Function,
354  void *FunctionParams,
355  double InitialGuess,
356  double Accuracy);
357 
358 static double ChiArea(CHISTRUCT *ChiParams, double x);
359 
360 static bool MultipleCharSamples(CLUSTERER* Clusterer,
361  CLUSTER* Cluster,
362  float MaxIllegal);
363 
364 static double InvertMatrix(const float* input, int size, float* inv);
365 
366 //--------------------------Public Code--------------------------------------
375 CLUSTERER *
376 MakeClusterer (int16_t SampleSize, const PARAM_DESC ParamDesc[]) {
377  CLUSTERER *Clusterer;
378  int i;
379 
380  // allocate main clusterer data structure and init simple fields
381  Clusterer = static_cast<CLUSTERER *>(Emalloc (sizeof (CLUSTERER)));
382  Clusterer->SampleSize = SampleSize;
383  Clusterer->NumberOfSamples = 0;
384  Clusterer->NumChar = 0;
385 
386  // init fields which will not be used initially
387  Clusterer->Root = nullptr;
388  Clusterer->ProtoList = NIL_LIST;
389 
390  // maintain a copy of param descriptors in the clusterer data structure
391  Clusterer->ParamDesc =
392  static_cast<PARAM_DESC *>(Emalloc (SampleSize * sizeof (PARAM_DESC)));
393  for (i = 0; i < SampleSize; i++) {
394  Clusterer->ParamDesc[i].Circular = ParamDesc[i].Circular;
395  Clusterer->ParamDesc[i].NonEssential = ParamDesc[i].NonEssential;
396  Clusterer->ParamDesc[i].Min = ParamDesc[i].Min;
397  Clusterer->ParamDesc[i].Max = ParamDesc[i].Max;
398  Clusterer->ParamDesc[i].Range = ParamDesc[i].Max - ParamDesc[i].Min;
399  Clusterer->ParamDesc[i].HalfRange = Clusterer->ParamDesc[i].Range / 2;
400  Clusterer->ParamDesc[i].MidRange =
401  (ParamDesc[i].Max + ParamDesc[i].Min) / 2;
402  }
403 
404  // allocate a kd tree to hold the samples
405  Clusterer->KDTree = MakeKDTree (SampleSize, ParamDesc);
406 
407  // Initialize cache of histogram buckets to minimize recomputing them.
408  for (auto & d : Clusterer->bucket_cache) {
409  for (auto & c : d)
410  c = nullptr;
411  }
412 
413  return Clusterer;
414 } // MakeClusterer
415 
429 SAMPLE* MakeSample(CLUSTERER * Clusterer, const float* Feature,
430  int32_t CharID) {
431  SAMPLE *Sample;
432  int i;
433 
434  // see if the samples have already been clustered - if so trap an error
435  // Can't add samples after they have been clustered.
436  ASSERT_HOST(Clusterer->Root == nullptr);
437 
438  // allocate the new sample and initialize it
439  Sample = static_cast<SAMPLE *>(Emalloc (sizeof (SAMPLE) +
440  (Clusterer->SampleSize -
441  1) * sizeof (float)));
442  Sample->Clustered = false;
443  Sample->Prototype = false;
444  Sample->SampleCount = 1;
445  Sample->Left = nullptr;
446  Sample->Right = nullptr;
447  Sample->CharID = CharID;
448 
449  for (i = 0; i < Clusterer->SampleSize; i++)
450  Sample->Mean[i] = Feature[i];
451 
452  // add the sample to the KD tree - keep track of the total # of samples
453  Clusterer->NumberOfSamples++;
454  KDStore(Clusterer->KDTree, Sample->Mean, Sample);
455  if (CharID >= Clusterer->NumChar)
456  Clusterer->NumChar = CharID + 1;
457 
458  // execute hook for monitoring clustering operation
459  // (*SampleCreationHook)(Sample);
460 
461  return (Sample);
462 } // MakeSample
463 
484  //only create cluster tree if samples have never been clustered before
485  if (Clusterer->Root == nullptr)
486  CreateClusterTree(Clusterer);
487 
488  //deallocate the old prototype list if one exists
489  FreeProtoList (&Clusterer->ProtoList);
490  Clusterer->ProtoList = NIL_LIST;
491 
492  //compute prototypes starting at the root node in the tree
493  ComputePrototypes(Clusterer, Config);
494  // We don't need the cluster pointers in the protos any more, so null them
495  // out, which makes it safe to delete the clusterer.
496  LIST proto_list = Clusterer->ProtoList;
497  iterate(proto_list) {
498  auto *proto = reinterpret_cast<PROTOTYPE *>(first_node(proto_list));
499  proto->Cluster = nullptr;
500  }
501  return Clusterer->ProtoList;
502 } // ClusterSamples
503 
514 void FreeClusterer(CLUSTERER *Clusterer) {
515  if (Clusterer != nullptr) {
516  free(Clusterer->ParamDesc);
517  if (Clusterer->KDTree != nullptr)
518  FreeKDTree (Clusterer->KDTree);
519  if (Clusterer->Root != nullptr)
520  FreeCluster (Clusterer->Root);
521  // Free up all used buckets structures.
522  for (auto & d : Clusterer->bucket_cache) {
523  for (auto & c : d)
524  if (c != nullptr)
525  FreeBuckets(c);
526  }
527 
528  free(Clusterer);
529  }
530 } // FreeClusterer
531 
538 void FreeProtoList(LIST *ProtoList) {
539  destroy_nodes(*ProtoList, FreePrototype);
540 } // FreeProtoList
541 
549 void FreePrototype(void *arg) { //PROTOTYPE *Prototype)
550  auto *Prototype = static_cast<PROTOTYPE *>(arg);
551 
552  // unmark the corresponding cluster (if there is one
553  if (Prototype->Cluster != nullptr)
554  Prototype->Cluster->Prototype = false;
555 
556  // deallocate the prototype statistics and then the prototype itself
557  free(Prototype->Distrib);
558  free(Prototype->Mean);
559  if (Prototype->Style != spherical) {
560  free(Prototype->Variance.Elliptical);
561  free(Prototype->Magnitude.Elliptical);
562  free(Prototype->Weight.Elliptical);
563  }
564  free(Prototype);
565 } // FreePrototype
566 
580 CLUSTER *NextSample(LIST *SearchState) {
581  CLUSTER *Cluster;
582 
583  if (*SearchState == NIL_LIST)
584  return (nullptr);
585  Cluster = reinterpret_cast<CLUSTER *>first_node (*SearchState);
586  *SearchState = pop (*SearchState);
587  for (;;) {
588  if (Cluster->Left == nullptr)
589  return (Cluster);
590  *SearchState = push (*SearchState, Cluster->Right);
591  Cluster = Cluster->Left;
592  }
593 } // NextSample
594 
602 float Mean(PROTOTYPE *Proto, uint16_t Dimension) {
603  return (Proto->Mean[Dimension]);
604 } // Mean
605 
613 float StandardDeviation(PROTOTYPE *Proto, uint16_t Dimension) {
614  switch (Proto->Style) {
615  case spherical:
616  return (static_cast<float>(sqrt (static_cast<double>(Proto->Variance.Spherical))));
617  case elliptical:
618  return (static_cast<float>(sqrt (static_cast<double>(Proto->Variance.Elliptical[Dimension]))));
619  case mixed:
620  switch (Proto->Distrib[Dimension]) {
621  case normal:
622  return (static_cast<float>(sqrt (static_cast<double>(Proto->Variance.Elliptical[Dimension]))));
623  case uniform:
624  case D_random:
625  return (Proto->Variance.Elliptical[Dimension]);
626  case DISTRIBUTION_COUNT:
627  ASSERT_HOST(!"Distribution count not allowed!");
628  }
629  }
630  return 0.0f;
631 } // StandardDeviation
632 
633 
634 /*---------------------------------------------------------------------------
635  Private Code
636 ----------------------------------------------------------------------------*/
650 static void CreateClusterTree(CLUSTERER *Clusterer) {
651  ClusteringContext context;
652  ClusterPair HeapEntry;
653  TEMPCLUSTER *PotentialCluster;
654 
655  // each sample and its nearest neighbor form a "potential" cluster
656  // save these in a heap with the "best" potential clusters on top
657  context.tree = Clusterer->KDTree;
658  context.candidates = static_cast<TEMPCLUSTER *>(Emalloc(Clusterer->NumberOfSamples * sizeof(TEMPCLUSTER)));
659  context.next = 0;
660  context.heap = new ClusterHeap(Clusterer->NumberOfSamples);
661  KDWalk(context.tree, reinterpret_cast<void_proc>(MakePotentialClusters), &context);
662 
663  // form potential clusters into actual clusters - always do "best" first
664  while (context.heap->Pop(&HeapEntry)) {
665  PotentialCluster = HeapEntry.data;
666 
667  // if main cluster of potential cluster is already in another cluster
668  // then we don't need to worry about it
669  if (PotentialCluster->Cluster->Clustered) {
670  continue;
671  }
672 
673  // if main cluster is not yet clustered, but its nearest neighbor is
674  // then we must find a new nearest neighbor
675  else if (PotentialCluster->Neighbor->Clustered) {
676  PotentialCluster->Neighbor =
677  FindNearestNeighbor(context.tree, PotentialCluster->Cluster,
678  &HeapEntry.key);
679  if (PotentialCluster->Neighbor != nullptr) {
680  context.heap->Push(&HeapEntry);
681  }
682  }
683 
684  // if neither cluster is already clustered, form permanent cluster
685  else {
686  PotentialCluster->Cluster =
687  MakeNewCluster(Clusterer, PotentialCluster);
688  PotentialCluster->Neighbor =
689  FindNearestNeighbor(context.tree, PotentialCluster->Cluster,
690  &HeapEntry.key);
691  if (PotentialCluster->Neighbor != nullptr) {
692  context.heap->Push(&HeapEntry);
693  }
694  }
695  }
696 
697  // the root node in the cluster tree is now the only node in the kd-tree
698  Clusterer->Root = static_cast<CLUSTER *>RootOf(Clusterer->KDTree);
699 
700  // free up the memory used by the K-D tree, heap, and temp clusters
701  FreeKDTree(context.tree);
702  Clusterer->KDTree = nullptr;
703  delete context.heap;
704  free(context.candidates);
705 } // CreateClusterTree
706 
716 static void MakePotentialClusters(ClusteringContext* context,
717  CLUSTER* Cluster, int32_t /*Level*/) {
718  ClusterPair HeapEntry;
719  int next = context->next;
720  context->candidates[next].Cluster = Cluster;
721  HeapEntry.data = &(context->candidates[next]);
722  context->candidates[next].Neighbor =
723  FindNearestNeighbor(context->tree,
724  context->candidates[next].Cluster,
725  &HeapEntry.key);
726  if (context->candidates[next].Neighbor != nullptr) {
727  context->heap->Push(&HeapEntry);
728  context->next++;
729  }
730 } // MakePotentialClusters
731 
745 static CLUSTER*
746 FindNearestNeighbor(KDTREE* Tree, CLUSTER* Cluster, float* Distance)
747 #define MAXNEIGHBORS 2
748 #define MAXDISTANCE FLT_MAX
749 {
750  CLUSTER *Neighbor[MAXNEIGHBORS];
751  float Dist[MAXNEIGHBORS];
752  int NumberOfNeighbors;
753  int32_t i;
754  CLUSTER *BestNeighbor;
755 
756  // find the 2 nearest neighbors of the cluster
758  &NumberOfNeighbors, reinterpret_cast<void **>(Neighbor), Dist);
759 
760  // search for the nearest neighbor that is not the cluster itself
761  *Distance = MAXDISTANCE;
762  BestNeighbor = nullptr;
763  for (i = 0; i < NumberOfNeighbors; i++) {
764  if ((Dist[i] < *Distance) && (Neighbor[i] != Cluster)) {
765  *Distance = Dist[i];
766  BestNeighbor = Neighbor[i];
767  }
768  }
769  return BestNeighbor;
770 } // FindNearestNeighbor
771 
781 static CLUSTER* MakeNewCluster(CLUSTERER* Clusterer,
782  TEMPCLUSTER* TempCluster) {
783  CLUSTER *Cluster;
784 
785  // allocate the new cluster and initialize it
786  Cluster = static_cast<CLUSTER *>(Emalloc(
787  sizeof(CLUSTER) + (Clusterer->SampleSize - 1) * sizeof(float)));
788  Cluster->Clustered = false;
789  Cluster->Prototype = false;
790  Cluster->Left = TempCluster->Cluster;
791  Cluster->Right = TempCluster->Neighbor;
792  Cluster->CharID = -1;
793 
794  // mark the old clusters as "clustered" and delete them from the kd-tree
795  Cluster->Left->Clustered = true;
796  Cluster->Right->Clustered = true;
797  KDDelete(Clusterer->KDTree, Cluster->Left->Mean, Cluster->Left);
798  KDDelete(Clusterer->KDTree, Cluster->Right->Mean, Cluster->Right);
799 
800  // compute the mean and sample count for the new cluster
801  Cluster->SampleCount =
802  MergeClusters(Clusterer->SampleSize, Clusterer->ParamDesc,
803  Cluster->Left->SampleCount, Cluster->Right->SampleCount,
804  Cluster->Mean, Cluster->Left->Mean, Cluster->Right->Mean);
805 
806  // add the new cluster to the KD tree
807  KDStore(Clusterer->KDTree, Cluster->Mean, Cluster);
808  return Cluster;
809 } // MakeNewCluster
810 
824 int32_t MergeClusters(int16_t N,
825  PARAM_DESC ParamDesc[],
826  int32_t n1,
827  int32_t n2,
828  float m[],
829  float m1[], float m2[]) {
830  int32_t i, n;
831 
832  n = n1 + n2;
833  for (i = N; i > 0; i--, ParamDesc++, m++, m1++, m2++) {
834  if (ParamDesc->Circular) {
835  // if distance between means is greater than allowed
836  // reduce upper point by one "rotation" to compute mean
837  // then normalize the mean back into the accepted range
838  if ((*m2 - *m1) > ParamDesc->HalfRange) {
839  *m = (n1 * *m1 + n2 * (*m2 - ParamDesc->Range)) / n;
840  if (*m < ParamDesc->Min)
841  *m += ParamDesc->Range;
842  }
843  else if ((*m1 - *m2) > ParamDesc->HalfRange) {
844  *m = (n1 * (*m1 - ParamDesc->Range) + n2 * *m2) / n;
845  if (*m < ParamDesc->Min)
846  *m += ParamDesc->Range;
847  }
848  else
849  *m = (n1 * *m1 + n2 * *m2) / n;
850  }
851  else
852  *m = (n1 * *m1 + n2 * *m2) / n;
853  }
854  return n;
855 } // MergeClusters
856 
865 static void ComputePrototypes(CLUSTERER* Clusterer, CLUSTERCONFIG* Config) {
866  LIST ClusterStack = NIL_LIST;
867  CLUSTER *Cluster;
868  PROTOTYPE *Prototype;
869 
870  // use a stack to keep track of clusters waiting to be processed
871  // initially the only cluster on the stack is the root cluster
872  if (Clusterer->Root != nullptr)
873  ClusterStack = push (NIL_LIST, Clusterer->Root);
874 
875  // loop until we have analyzed all clusters which are potential prototypes
876  while (ClusterStack != NIL_LIST) {
877  // remove the next cluster to be analyzed from the stack
878  // try to make a prototype from the cluster
879  // if successful, put it on the proto list, else split the cluster
880  Cluster = reinterpret_cast<CLUSTER *>first_node (ClusterStack);
881  ClusterStack = pop (ClusterStack);
882  Prototype = MakePrototype(Clusterer, Config, Cluster);
883  if (Prototype != nullptr) {
884  Clusterer->ProtoList = push (Clusterer->ProtoList, Prototype);
885  }
886  else {
887  ClusterStack = push (ClusterStack, Cluster->Right);
888  ClusterStack = push (ClusterStack, Cluster->Left);
889  }
890  }
891 } // ComputePrototypes
892 
908 static PROTOTYPE* MakePrototype(CLUSTERER* Clusterer, CLUSTERCONFIG* Config,
909  CLUSTER* Cluster) {
910  STATISTICS *Statistics;
911  PROTOTYPE *Proto;
912  BUCKETS *Buckets;
913 
914  // filter out clusters which contain samples from the same character
915  if (MultipleCharSamples (Clusterer, Cluster, Config->MaxIllegal))
916  return nullptr;
917 
918  // compute the covariance matrix and ranges for the cluster
919  Statistics =
920  ComputeStatistics(Clusterer->SampleSize, Clusterer->ParamDesc, Cluster);
921 
922  // check for degenerate clusters which need not be analyzed further
923  // note that the MinSamples test assumes that all clusters with multiple
924  // character samples have been removed (as above)
925  Proto = MakeDegenerateProto(
926  Clusterer->SampleSize, Cluster, Statistics, Config->ProtoStyle,
927  static_cast<int32_t>(Config->MinSamples * Clusterer->NumChar));
928  if (Proto != nullptr) {
929  FreeStatistics(Statistics);
930  return Proto;
931  }
932  // check to ensure that all dimensions are independent
933  if (!Independent(Clusterer->ParamDesc, Clusterer->SampleSize,
934  Statistics->CoVariance, Config->Independence)) {
935  FreeStatistics(Statistics);
936  return nullptr;
937  }
938 
939  if (HOTELLING && Config->ProtoStyle == elliptical) {
940  Proto = TestEllipticalProto(Clusterer, Config, Cluster, Statistics);
941  if (Proto != nullptr) {
942  FreeStatistics(Statistics);
943  return Proto;
944  }
945  }
946 
947  // create a histogram data structure used to evaluate distributions
948  Buckets = GetBuckets(Clusterer, normal, Cluster->SampleCount,
949  Config->Confidence);
950 
951  // create a prototype based on the statistics and test it
952  switch (Config->ProtoStyle) {
953  case spherical:
954  Proto = MakeSphericalProto(Clusterer, Cluster, Statistics, Buckets);
955  break;
956  case elliptical:
957  Proto = MakeEllipticalProto(Clusterer, Cluster, Statistics, Buckets);
958  break;
959  case mixed:
960  Proto = MakeMixedProto(Clusterer, Cluster, Statistics, Buckets,
961  Config->Confidence);
962  break;
963  case automatic:
964  Proto = MakeSphericalProto(Clusterer, Cluster, Statistics, Buckets);
965  if (Proto != nullptr)
966  break;
967  Proto = MakeEllipticalProto(Clusterer, Cluster, Statistics, Buckets);
968  if (Proto != nullptr)
969  break;
970  Proto = MakeMixedProto(Clusterer, Cluster, Statistics, Buckets,
971  Config->Confidence);
972  break;
973  }
974  FreeStatistics(Statistics);
975  return Proto;
976 } // MakePrototype
977 
997 static PROTOTYPE* MakeDegenerateProto( //this was MinSample
998  uint16_t N,
999  CLUSTER *Cluster,
1000  STATISTICS *Statistics,
1001  PROTOSTYLE Style,
1002  int32_t MinSamples) {
1003  PROTOTYPE *Proto = nullptr;
1004 
1005  if (MinSamples < MINSAMPLESNEEDED)
1006  MinSamples = MINSAMPLESNEEDED;
1007 
1008  if (Cluster->SampleCount < MinSamples) {
1009  switch (Style) {
1010  case spherical:
1011  Proto = NewSphericalProto (N, Cluster, Statistics);
1012  break;
1013  case elliptical:
1014  case automatic:
1015  Proto = NewEllipticalProto (N, Cluster, Statistics);
1016  break;
1017  case mixed:
1018  Proto = NewMixedProto (N, Cluster, Statistics);
1019  break;
1020  }
1021  Proto->Significant = false;
1022  }
1023  return (Proto);
1024 } // MakeDegenerateProto
1025 
1039 static PROTOTYPE* TestEllipticalProto(CLUSTERER* Clusterer,
1040  CLUSTERCONFIG *Config, CLUSTER* Cluster,
1041  STATISTICS* Statistics) {
1042  // Fraction of the number of samples used as a range around 1 within
1043  // which a cluster has the magic size that allows a boost to the
1044  // FTable by kFTableBoostMargin, thus allowing clusters near the
1045  // magic size (equal to the number of sample characters) to be more
1046  // likely to stay together.
1047  const double kMagicSampleMargin = 0.0625;
1048  const double kFTableBoostMargin = 2.0;
1049 
1050  int N = Clusterer->SampleSize;
1051  CLUSTER* Left = Cluster->Left;
1052  CLUSTER* Right = Cluster->Right;
1053  if (Left == nullptr || Right == nullptr)
1054  return nullptr;
1055  int TotalDims = Left->SampleCount + Right->SampleCount;
1056  if (TotalDims < N + 1 || TotalDims < 2)
1057  return nullptr;
1058  std::vector<float> Covariance(static_cast<size_t>(N) * N);
1059  std::vector<float> Inverse(static_cast<size_t>(N) * N);
1060  std::vector<float> Delta(N);
1061  // Compute a new covariance matrix that only uses essential features.
1062  for (int i = 0; i < N; ++i) {
1063  int row_offset = i * N;
1064  if (!Clusterer->ParamDesc[i].NonEssential) {
1065  for (int j = 0; j < N; ++j) {
1066  if (!Clusterer->ParamDesc[j].NonEssential)
1067  Covariance[j + row_offset] = Statistics->CoVariance[j + row_offset];
1068  else
1069  Covariance[j + row_offset] = 0.0f;
1070  }
1071  } else {
1072  for (int j = 0; j < N; ++j) {
1073  if (i == j)
1074  Covariance[j + row_offset] = 1.0f;
1075  else
1076  Covariance[j + row_offset] = 0.0f;
1077  }
1078  }
1079  }
1080  double err = InvertMatrix(&Covariance[0], N, &Inverse[0]);
1081  if (err > 1) {
1082  tprintf("Clustering error: Matrix inverse failed with error %g\n", err);
1083  }
1084  int EssentialN = 0;
1085  for (int dim = 0; dim < N; ++dim) {
1086  if (!Clusterer->ParamDesc[dim].NonEssential) {
1087  Delta[dim] = Left->Mean[dim] - Right->Mean[dim];
1088  ++EssentialN;
1089  } else {
1090  Delta[dim] = 0.0f;
1091  }
1092  }
1093  // Compute Hotelling's T-squared.
1094  double Tsq = 0.0;
1095  for (int x = 0; x < N; ++x) {
1096  double temp = 0.0;
1097  for (int y = 0; y < N; ++y) {
1098  temp += static_cast<double>(Inverse[y + N * x]) * Delta[y];
1099  }
1100  Tsq += Delta[x] * temp;
1101  }
1102  // Changed this function to match the formula in
1103  // Statistical Methods in Medical Research p 473
1104  // By Peter Armitage, Geoffrey Berry, J. N. S. Matthews.
1105  // Tsq *= Left->SampleCount * Right->SampleCount / TotalDims;
1106  double F = Tsq * (TotalDims - EssentialN - 1) / ((TotalDims - 2)*EssentialN);
1107  int Fx = EssentialN;
1108  if (Fx > FTABLE_X)
1109  Fx = FTABLE_X;
1110  --Fx;
1111  int Fy = TotalDims - EssentialN - 1;
1112  if (Fy > FTABLE_Y)
1113  Fy = FTABLE_Y;
1114  --Fy;
1115  double FTarget = FTable[Fy][Fx];
1116  if (Config->MagicSamples > 0 &&
1117  TotalDims >= Config->MagicSamples * (1.0 - kMagicSampleMargin) &&
1118  TotalDims <= Config->MagicSamples * (1.0 + kMagicSampleMargin)) {
1119  // Give magic-sized clusters a magic FTable boost.
1120  FTarget += kFTableBoostMargin;
1121  }
1122  if (F < FTarget) {
1123  return NewEllipticalProto (Clusterer->SampleSize, Cluster, Statistics);
1124  }
1125  return nullptr;
1126 }
1127 
1139 static PROTOTYPE* MakeSphericalProto(CLUSTERER* Clusterer,
1140  CLUSTER* Cluster, STATISTICS* Statistics,
1141  BUCKETS* Buckets) {
1142  PROTOTYPE *Proto = nullptr;
1143  int i;
1144 
1145  // check that each dimension is a normal distribution
1146  for (i = 0; i < Clusterer->SampleSize; i++) {
1147  if (Clusterer->ParamDesc[i].NonEssential)
1148  continue;
1149 
1150  FillBuckets (Buckets, Cluster, i, &(Clusterer->ParamDesc[i]),
1151  Cluster->Mean[i],
1152  sqrt (static_cast<double>(Statistics->AvgVariance)));
1153  if (!DistributionOK (Buckets))
1154  break;
1155  }
1156  // if all dimensions matched a normal distribution, make a proto
1157  if (i >= Clusterer->SampleSize)
1158  Proto = NewSphericalProto (Clusterer->SampleSize, Cluster, Statistics);
1159  return (Proto);
1160 } // MakeSphericalProto
1161 
1173 static PROTOTYPE* MakeEllipticalProto(CLUSTERER* Clusterer,
1174  CLUSTER* Cluster, STATISTICS* Statistics,
1175  BUCKETS* Buckets) {
1176  PROTOTYPE *Proto = nullptr;
1177  int i;
1178 
1179  // check that each dimension is a normal distribution
1180  for (i = 0; i < Clusterer->SampleSize; i++) {
1181  if (Clusterer->ParamDesc[i].NonEssential)
1182  continue;
1183 
1184  FillBuckets (Buckets, Cluster, i, &(Clusterer->ParamDesc[i]),
1185  Cluster->Mean[i],
1186  sqrt (static_cast<double>(Statistics->
1187  CoVariance[i * (Clusterer->SampleSize + 1)])));
1188  if (!DistributionOK (Buckets))
1189  break;
1190  }
1191  // if all dimensions matched a normal distribution, make a proto
1192  if (i >= Clusterer->SampleSize)
1193  Proto = NewEllipticalProto (Clusterer->SampleSize, Cluster, Statistics);
1194  return (Proto);
1195 } // MakeEllipticalProto
1196 
1212 static PROTOTYPE* MakeMixedProto(CLUSTERER* Clusterer,
1213  CLUSTER* Cluster, STATISTICS* Statistics,
1214  BUCKETS* NormalBuckets, double Confidence) {
1215  PROTOTYPE *Proto;
1216  int i;
1217  BUCKETS *UniformBuckets = nullptr;
1218  BUCKETS *RandomBuckets = nullptr;
1219 
1220  // create a mixed proto to work on - initially assume all dimensions normal*/
1221  Proto = NewMixedProto (Clusterer->SampleSize, Cluster, Statistics);
1222 
1223  // find the proper distribution for each dimension
1224  for (i = 0; i < Clusterer->SampleSize; i++) {
1225  if (Clusterer->ParamDesc[i].NonEssential)
1226  continue;
1227 
1228  FillBuckets (NormalBuckets, Cluster, i, &(Clusterer->ParamDesc[i]),
1229  Proto->Mean[i],
1230  sqrt (static_cast<double>(Proto->Variance.Elliptical[i])));
1231  if (DistributionOK (NormalBuckets))
1232  continue;
1233 
1234  if (RandomBuckets == nullptr)
1235  RandomBuckets =
1236  GetBuckets(Clusterer, D_random, Cluster->SampleCount, Confidence);
1237  MakeDimRandom (i, Proto, &(Clusterer->ParamDesc[i]));
1238  FillBuckets (RandomBuckets, Cluster, i, &(Clusterer->ParamDesc[i]),
1239  Proto->Mean[i], Proto->Variance.Elliptical[i]);
1240  if (DistributionOK (RandomBuckets))
1241  continue;
1242 
1243  if (UniformBuckets == nullptr)
1244  UniformBuckets =
1245  GetBuckets(Clusterer, uniform, Cluster->SampleCount, Confidence);
1246  MakeDimUniform(i, Proto, Statistics);
1247  FillBuckets (UniformBuckets, Cluster, i, &(Clusterer->ParamDesc[i]),
1248  Proto->Mean[i], Proto->Variance.Elliptical[i]);
1249  if (DistributionOK (UniformBuckets))
1250  continue;
1251  break;
1252  }
1253  // if any dimension failed to match a distribution, discard the proto
1254  if (i < Clusterer->SampleSize) {
1255  FreePrototype(Proto);
1256  Proto = nullptr;
1257  }
1258  return (Proto);
1259 } // MakeMixedProto
1260 
1268 static void MakeDimRandom(uint16_t i, PROTOTYPE* Proto, PARAM_DESC* ParamDesc) {
1269  Proto->Distrib[i] = D_random;
1270  Proto->Mean[i] = ParamDesc->MidRange;
1271  Proto->Variance.Elliptical[i] = ParamDesc->HalfRange;
1272 
1273  // subtract out the previous magnitude of this dimension from the total
1274  Proto->TotalMagnitude /= Proto->Magnitude.Elliptical[i];
1275  Proto->Magnitude.Elliptical[i] = 1.0 / ParamDesc->Range;
1276  Proto->TotalMagnitude *= Proto->Magnitude.Elliptical[i];
1277  Proto->LogMagnitude = log (static_cast<double>(Proto->TotalMagnitude));
1278 
1279  // note that the proto Weight is irrelevant for D_random protos
1280 } // MakeDimRandom
1281 
1289 static void MakeDimUniform(uint16_t i, PROTOTYPE* Proto, STATISTICS* Statistics) {
1290  Proto->Distrib[i] = uniform;
1291  Proto->Mean[i] = Proto->Cluster->Mean[i] +
1292  (Statistics->Min[i] + Statistics->Max[i]) / 2;
1293  Proto->Variance.Elliptical[i] =
1294  (Statistics->Max[i] - Statistics->Min[i]) / 2;
1295  if (Proto->Variance.Elliptical[i] < MINVARIANCE)
1296  Proto->Variance.Elliptical[i] = MINVARIANCE;
1297 
1298  // subtract out the previous magnitude of this dimension from the total
1299  Proto->TotalMagnitude /= Proto->Magnitude.Elliptical[i];
1300  Proto->Magnitude.Elliptical[i] =
1301  1.0 / (2.0 * Proto->Variance.Elliptical[i]);
1302  Proto->TotalMagnitude *= Proto->Magnitude.Elliptical[i];
1303  Proto->LogMagnitude = log (static_cast<double>(Proto->TotalMagnitude));
1304 
1305  // note that the proto Weight is irrelevant for uniform protos
1306 } // MakeDimUniform
1307 
1322 static STATISTICS*
1323 ComputeStatistics (int16_t N, PARAM_DESC ParamDesc[], CLUSTER * Cluster) {
1324  STATISTICS *Statistics;
1325  int i, j;
1326  float *CoVariance;
1327  float *Distance;
1328  LIST SearchState;
1329  SAMPLE *Sample;
1330  uint32_t SampleCountAdjustedForBias;
1331 
1332  // allocate memory to hold the statistics results
1333  Statistics = static_cast<STATISTICS *>(Emalloc (sizeof (STATISTICS)));
1334  Statistics->CoVariance = static_cast<float *>(Emalloc(sizeof(float) * N * N));
1335  Statistics->Min = static_cast<float *>(Emalloc (N * sizeof (float)));
1336  Statistics->Max = static_cast<float *>(Emalloc (N * sizeof (float)));
1337 
1338  // allocate temporary memory to hold the sample to mean distances
1339  Distance = static_cast<float *>(Emalloc (N * sizeof (float)));
1340 
1341  // initialize the statistics
1342  Statistics->AvgVariance = 1.0;
1343  CoVariance = Statistics->CoVariance;
1344  for (i = 0; i < N; i++) {
1345  Statistics->Min[i] = 0.0;
1346  Statistics->Max[i] = 0.0;
1347  for (j = 0; j < N; j++, CoVariance++)
1348  *CoVariance = 0;
1349  }
1350  // find each sample in the cluster and merge it into the statistics
1351  InitSampleSearch(SearchState, Cluster);
1352  while ((Sample = NextSample (&SearchState)) != nullptr) {
1353  for (i = 0; i < N; i++) {
1354  Distance[i] = Sample->Mean[i] - Cluster->Mean[i];
1355  if (ParamDesc[i].Circular) {
1356  if (Distance[i] > ParamDesc[i].HalfRange)
1357  Distance[i] -= ParamDesc[i].Range;
1358  if (Distance[i] < -ParamDesc[i].HalfRange)
1359  Distance[i] += ParamDesc[i].Range;
1360  }
1361  if (Distance[i] < Statistics->Min[i])
1362  Statistics->Min[i] = Distance[i];
1363  if (Distance[i] > Statistics->Max[i])
1364  Statistics->Max[i] = Distance[i];
1365  }
1366  CoVariance = Statistics->CoVariance;
1367  for (i = 0; i < N; i++)
1368  for (j = 0; j < N; j++, CoVariance++)
1369  *CoVariance += Distance[i] * Distance[j];
1370  }
1371  // normalize the variances by the total number of samples
1372  // use SampleCount-1 instead of SampleCount to get an unbiased estimate
1373  // also compute the geometic mean of the diagonal variances
1374  // ensure that clusters with only 1 sample are handled correctly
1375  if (Cluster->SampleCount > 1)
1376  SampleCountAdjustedForBias = Cluster->SampleCount - 1;
1377  else
1378  SampleCountAdjustedForBias = 1;
1379  CoVariance = Statistics->CoVariance;
1380  for (i = 0; i < N; i++)
1381  for (j = 0; j < N; j++, CoVariance++) {
1382  *CoVariance /= SampleCountAdjustedForBias;
1383  if (j == i) {
1384  if (*CoVariance < MINVARIANCE)
1385  *CoVariance = MINVARIANCE;
1386  Statistics->AvgVariance *= *CoVariance;
1387  }
1388  }
1389  Statistics->AvgVariance = static_cast<float>(pow(static_cast<double>(Statistics->AvgVariance),
1390  1.0 / N));
1391 
1392  // release temporary memory and return
1393  free(Distance);
1394  return (Statistics);
1395 } // ComputeStatistics
1396 
1408 static PROTOTYPE* NewSphericalProto(uint16_t N, CLUSTER* Cluster,
1409  STATISTICS* Statistics) {
1410  PROTOTYPE *Proto;
1411 
1412  Proto = NewSimpleProto (N, Cluster);
1413 
1414  Proto->Variance.Spherical = Statistics->AvgVariance;
1415  if (Proto->Variance.Spherical < MINVARIANCE)
1416  Proto->Variance.Spherical = MINVARIANCE;
1417 
1418  Proto->Magnitude.Spherical =
1419  1.0 / sqrt(2.0 * M_PI * Proto->Variance.Spherical);
1420  Proto->TotalMagnitude = static_cast<float>(pow(static_cast<double>(Proto->Magnitude.Spherical),
1421  static_cast<double>(N)));
1422  Proto->Weight.Spherical = 1.0 / Proto->Variance.Spherical;
1423  Proto->LogMagnitude = log (static_cast<double>(Proto->TotalMagnitude));
1424 
1425  return (Proto);
1426 } // NewSphericalProto
1427 
1438 static PROTOTYPE* NewEllipticalProto(int16_t N, CLUSTER* Cluster,
1439  STATISTICS* Statistics) {
1440  PROTOTYPE *Proto;
1441  float *CoVariance;
1442  int i;
1443 
1444  Proto = NewSimpleProto (N, Cluster);
1445  Proto->Variance.Elliptical = static_cast<float *>(Emalloc (N * sizeof (float)));
1446  Proto->Magnitude.Elliptical = static_cast<float *>(Emalloc (N * sizeof (float)));
1447  Proto->Weight.Elliptical = static_cast<float *>(Emalloc (N * sizeof (float)));
1448 
1449  CoVariance = Statistics->CoVariance;
1450  Proto->TotalMagnitude = 1.0;
1451  for (i = 0; i < N; i++, CoVariance += N + 1) {
1452  Proto->Variance.Elliptical[i] = *CoVariance;
1453  if (Proto->Variance.Elliptical[i] < MINVARIANCE)
1454  Proto->Variance.Elliptical[i] = MINVARIANCE;
1455 
1456  Proto->Magnitude.Elliptical[i] =
1457  1.0 / sqrt(2.0 * M_PI * Proto->Variance.Elliptical[i]);
1458  Proto->Weight.Elliptical[i] = 1.0 / Proto->Variance.Elliptical[i];
1459  Proto->TotalMagnitude *= Proto->Magnitude.Elliptical[i];
1460  }
1461  Proto->LogMagnitude = log (static_cast<double>(Proto->TotalMagnitude));
1462  Proto->Style = elliptical;
1463  return (Proto);
1464 } // NewEllipticalProto
1465 
1479 static PROTOTYPE* NewMixedProto(int16_t N, CLUSTER* Cluster,
1480  STATISTICS* Statistics) {
1481  PROTOTYPE *Proto;
1482  int i;
1483 
1484  Proto = NewEllipticalProto (N, Cluster, Statistics);
1485  Proto->Distrib = static_cast<DISTRIBUTION *>(Emalloc (N * sizeof (DISTRIBUTION)));
1486 
1487  for (i = 0; i < N; i++) {
1488  Proto->Distrib[i] = normal;
1489  }
1490  Proto->Style = mixed;
1491  return (Proto);
1492 } // NewMixedProto
1493 
1502 static PROTOTYPE *NewSimpleProto(int16_t N, CLUSTER *Cluster) {
1503  PROTOTYPE *Proto;
1504  int i;
1505 
1506  Proto = static_cast<PROTOTYPE *>(Emalloc (sizeof (PROTOTYPE)));
1507  Proto->Mean = static_cast<float *>(Emalloc (N * sizeof (float)));
1508 
1509  for (i = 0; i < N; i++)
1510  Proto->Mean[i] = Cluster->Mean[i];
1511  Proto->Distrib = nullptr;
1512 
1513  Proto->Significant = true;
1514  Proto->Merged = false;
1515  Proto->Style = spherical;
1516  Proto->NumSamples = Cluster->SampleCount;
1517  Proto->Cluster = Cluster;
1518  Proto->Cluster->Prototype = true;
1519  return (Proto);
1520 } // NewSimpleProto
1521 
1540 static bool
1541 Independent(PARAM_DESC* ParamDesc,
1542  int16_t N, float* CoVariance, float Independence) {
1543  int i, j;
1544  float *VARii; // points to ith on-diagonal element
1545  float *VARjj; // points to jth on-diagonal element
1546  float CorrelationCoeff;
1547 
1548  VARii = CoVariance;
1549  for (i = 0; i < N; i++, VARii += N + 1) {
1550  if (ParamDesc[i].NonEssential)
1551  continue;
1552 
1553  VARjj = VARii + N + 1;
1554  CoVariance = VARii + 1;
1555  for (j = i + 1; j < N; j++, CoVariance++, VARjj += N + 1) {
1556  if (ParamDesc[j].NonEssential)
1557  continue;
1558 
1559  if ((*VARii == 0.0) || (*VARjj == 0.0))
1560  CorrelationCoeff = 0.0;
1561  else
1562  CorrelationCoeff =
1563  sqrt (sqrt (*CoVariance * *CoVariance / (*VARii * *VARjj)));
1564  if (CorrelationCoeff > Independence)
1565  return false;
1566  }
1567  }
1568  return true;
1569 } // Independent
1570 
1586 static BUCKETS *GetBuckets(CLUSTERER* clusterer,
1587  DISTRIBUTION Distribution,
1588  uint32_t SampleCount,
1589  double Confidence) {
1590  // Get an old bucket structure with the same number of buckets.
1591  uint16_t NumberOfBuckets = OptimumNumberOfBuckets(SampleCount);
1592  BUCKETS *Buckets =
1593  clusterer->bucket_cache[Distribution][NumberOfBuckets - MINBUCKETS];
1594 
1595  // If a matching bucket structure is not found, make one and save it.
1596  if (Buckets == nullptr) {
1597  Buckets = MakeBuckets(Distribution, SampleCount, Confidence);
1598  clusterer->bucket_cache[Distribution][NumberOfBuckets - MINBUCKETS] =
1599  Buckets;
1600  } else {
1601  // Just adjust the existing buckets.
1602  if (SampleCount != Buckets->SampleCount)
1603  AdjustBuckets(Buckets, SampleCount);
1604  if (Confidence != Buckets->Confidence) {
1605  Buckets->Confidence = Confidence;
1606  Buckets->ChiSquared = ComputeChiSquared(
1607  DegreesOfFreedom(Distribution, Buckets->NumberOfBuckets),
1608  Confidence);
1609  }
1610  InitBuckets(Buckets);
1611  }
1612  return Buckets;
1613 } // GetBuckets
1614 
1631 static BUCKETS *MakeBuckets(DISTRIBUTION Distribution,
1632  uint32_t SampleCount,
1633  double Confidence) {
1634  const DENSITYFUNC DensityFunction[] =
1635  { NormalDensity, UniformDensity, UniformDensity };
1636  int i, j;
1637  BUCKETS *Buckets;
1638  double BucketProbability;
1639  double NextBucketBoundary;
1640  double Probability;
1641  double ProbabilityDelta;
1642  double LastProbDensity;
1643  double ProbDensity;
1644  uint16_t CurrentBucket;
1645  bool Symmetrical;
1646 
1647  // allocate memory needed for data structure
1648  Buckets = static_cast<BUCKETS *>(Emalloc(sizeof(BUCKETS)));
1649  Buckets->NumberOfBuckets = OptimumNumberOfBuckets(SampleCount);
1650  Buckets->SampleCount = SampleCount;
1651  Buckets->Confidence = Confidence;
1652  Buckets->Count =
1653  static_cast<uint32_t *>(Emalloc(Buckets->NumberOfBuckets * sizeof(uint32_t)));
1654  Buckets->ExpectedCount = static_cast<float *>(
1655  Emalloc(Buckets->NumberOfBuckets * sizeof(float)));
1656 
1657  // initialize simple fields
1658  Buckets->Distribution = Distribution;
1659  for (i = 0; i < Buckets->NumberOfBuckets; i++) {
1660  Buckets->Count[i] = 0;
1661  Buckets->ExpectedCount[i] = 0.0;
1662  }
1663 
1664  // all currently defined distributions are symmetrical
1665  Symmetrical = true;
1666  Buckets->ChiSquared = ComputeChiSquared(
1667  DegreesOfFreedom(Distribution, Buckets->NumberOfBuckets), Confidence);
1668 
1669  if (Symmetrical) {
1670  // allocate buckets so that all have approx. equal probability
1671  BucketProbability = 1.0 / static_cast<double>(Buckets->NumberOfBuckets);
1672 
1673  // distribution is symmetric so fill in upper half then copy
1674  CurrentBucket = Buckets->NumberOfBuckets / 2;
1675  if (Odd (Buckets->NumberOfBuckets))
1676  NextBucketBoundary = BucketProbability / 2;
1677  else
1678  NextBucketBoundary = BucketProbability;
1679 
1680  Probability = 0.0;
1681  LastProbDensity =
1682  (*DensityFunction[static_cast<int>(Distribution)]) (BUCKETTABLESIZE / 2);
1683  for (i = BUCKETTABLESIZE / 2; i < BUCKETTABLESIZE; i++) {
1684  ProbDensity = (*DensityFunction[static_cast<int>(Distribution)]) (i + 1);
1685  ProbabilityDelta = Integral (LastProbDensity, ProbDensity, 1.0);
1686  Probability += ProbabilityDelta;
1687  if (Probability > NextBucketBoundary) {
1688  if (CurrentBucket < Buckets->NumberOfBuckets - 1)
1689  CurrentBucket++;
1690  NextBucketBoundary += BucketProbability;
1691  }
1692  Buckets->Bucket[i] = CurrentBucket;
1693  Buckets->ExpectedCount[CurrentBucket] +=
1694  static_cast<float>(ProbabilityDelta * SampleCount);
1695  LastProbDensity = ProbDensity;
1696  }
1697  // place any leftover probability into the last bucket
1698  Buckets->ExpectedCount[CurrentBucket] +=
1699  static_cast<float>((0.5 - Probability) * SampleCount);
1700 
1701  // copy upper half of distribution to lower half
1702  for (i = 0, j = BUCKETTABLESIZE - 1; i < j; i++, j--)
1703  Buckets->Bucket[i] =
1704  Mirror(Buckets->Bucket[j], Buckets->NumberOfBuckets);
1705 
1706  // copy upper half of expected counts to lower half
1707  for (i = 0, j = Buckets->NumberOfBuckets - 1; i <= j; i++, j--)
1708  Buckets->ExpectedCount[i] += Buckets->ExpectedCount[j];
1709  }
1710  return Buckets;
1711 } // MakeBuckets
1712 
1726 static uint16_t OptimumNumberOfBuckets(uint32_t SampleCount) {
1727  uint8_t Last, Next;
1728  float Slope;
1729 
1730  if (SampleCount < kCountTable[0])
1731  return kBucketsTable[0];
1732 
1733  for (Last = 0, Next = 1; Next < LOOKUPTABLESIZE; Last++, Next++) {
1734  if (SampleCount <= kCountTable[Next]) {
1735  Slope = static_cast<float>(kBucketsTable[Next] - kBucketsTable[Last]) /
1736  static_cast<float>(kCountTable[Next] - kCountTable[Last]);
1737  return (static_cast<uint16_t>(kBucketsTable[Last] +
1738  Slope * (SampleCount - kCountTable[Last])));
1739  }
1740  }
1741  return kBucketsTable[Last];
1742 } // OptimumNumberOfBuckets
1743 
1760 static double
1761 ComputeChiSquared (uint16_t DegreesOfFreedom, double Alpha)
1762 #define CHIACCURACY 0.01
1763 #define MINALPHA (1e-200)
1764 {
1765  static LIST ChiWith[MAXDEGREESOFFREEDOM + 1];
1766 
1767  CHISTRUCT *OldChiSquared;
1768  CHISTRUCT SearchKey;
1769 
1770  // limit the minimum alpha that can be used - if alpha is too small
1771  // it may not be possible to compute chi-squared.
1772  Alpha = ClipToRange(Alpha, MINALPHA, 1.0);
1773  if (Odd (DegreesOfFreedom))
1774  DegreesOfFreedom++;
1775 
1776  /* find the list of chi-squared values which have already been computed
1777  for the specified number of degrees of freedom. Search the list for
1778  the desired chi-squared. */
1779  SearchKey.Alpha = Alpha;
1780  OldChiSquared = reinterpret_cast<CHISTRUCT *>first_node (search (ChiWith[DegreesOfFreedom],
1781  &SearchKey, AlphaMatch));
1782 
1783  if (OldChiSquared == nullptr) {
1784  OldChiSquared = NewChiStruct (DegreesOfFreedom, Alpha);
1785  OldChiSquared->ChiSquared = Solve (ChiArea, OldChiSquared,
1786  static_cast<double>(DegreesOfFreedom),
1787  CHIACCURACY);
1788  ChiWith[DegreesOfFreedom] = push (ChiWith[DegreesOfFreedom],
1789  OldChiSquared);
1790  }
1791  else {
1792  // further optimization might move OldChiSquared to front of list
1793  }
1794 
1795  return (OldChiSquared->ChiSquared);
1796 
1797 } // ComputeChiSquared
1798 
1812 static double NormalDensity(int32_t x) {
1813  double Distance;
1814 
1815  Distance = x - kNormalMean;
1816  return kNormalMagnitude * exp(-0.5 * Distance * Distance / kNormalVariance);
1817 } // NormalDensity
1818 
1826 static double UniformDensity(int32_t x) {
1827  constexpr auto UniformDistributionDensity = 1.0 / BUCKETTABLESIZE;
1828 
1829  if ((x >= 0) && (x <= BUCKETTABLESIZE)) {
1830  return UniformDistributionDensity;
1831  } else {
1832  return 0.0;
1833  }
1834 } // UniformDensity
1835 
1844 static double Integral(double f1, double f2, double Dx) {
1845  return (f1 + f2) * Dx / 2.0;
1846 } // Integral
1847 
1868 static void FillBuckets(BUCKETS *Buckets,
1869  CLUSTER *Cluster,
1870  uint16_t Dim,
1871  PARAM_DESC *ParamDesc,
1872  float Mean,
1873  float StdDev) {
1874  uint16_t BucketID;
1875  int i;
1876  LIST SearchState;
1877  SAMPLE *Sample;
1878 
1879  // initialize the histogram bucket counts to 0
1880  for (i = 0; i < Buckets->NumberOfBuckets; i++)
1881  Buckets->Count[i] = 0;
1882 
1883  if (StdDev == 0.0) {
1884  /* if the standard deviation is zero, then we can't statistically
1885  analyze the cluster. Use a pseudo-analysis: samples exactly on
1886  the mean are distributed evenly across all buckets. Samples greater
1887  than the mean are placed in the last bucket; samples less than the
1888  mean are placed in the first bucket. */
1889 
1890  InitSampleSearch(SearchState, Cluster);
1891  i = 0;
1892  while ((Sample = NextSample (&SearchState)) != nullptr) {
1893  if (Sample->Mean[Dim] > Mean)
1894  BucketID = Buckets->NumberOfBuckets - 1;
1895  else if (Sample->Mean[Dim] < Mean)
1896  BucketID = 0;
1897  else
1898  BucketID = i;
1899  Buckets->Count[BucketID] += 1;
1900  i++;
1901  if (i >= Buckets->NumberOfBuckets)
1902  i = 0;
1903  }
1904  }
1905  else {
1906  // search for all samples in the cluster and add to histogram buckets
1907  InitSampleSearch(SearchState, Cluster);
1908  while ((Sample = NextSample (&SearchState)) != nullptr) {
1909  switch (Buckets->Distribution) {
1910  case normal:
1911  BucketID = NormalBucket (ParamDesc, Sample->Mean[Dim],
1912  Mean, StdDev);
1913  break;
1914  case D_random:
1915  case uniform:
1916  BucketID = UniformBucket (ParamDesc, Sample->Mean[Dim],
1917  Mean, StdDev);
1918  break;
1919  default:
1920  BucketID = 0;
1921  }
1922  Buckets->Count[Buckets->Bucket[BucketID]] += 1;
1923  }
1924  }
1925 } // FillBuckets
1926 
1938 static uint16_t NormalBucket(PARAM_DESC *ParamDesc,
1939  float x,
1940  float Mean,
1941  float StdDev) {
1942  float X;
1943 
1944  // wraparound circular parameters if necessary
1945  if (ParamDesc->Circular) {
1946  if (x - Mean > ParamDesc->HalfRange)
1947  x -= ParamDesc->Range;
1948  else if (x - Mean < -ParamDesc->HalfRange)
1949  x += ParamDesc->Range;
1950  }
1951 
1952  X = ((x - Mean) / StdDev) * kNormalStdDev + kNormalMean;
1953  if (X < 0)
1954  return 0;
1955  if (X > BUCKETTABLESIZE - 1)
1956  return (static_cast<uint16_t>(BUCKETTABLESIZE - 1));
1957  return static_cast<uint16_t>(floor(static_cast<double>(X)));
1958 } // NormalBucket
1959 
1971 static uint16_t UniformBucket(PARAM_DESC *ParamDesc,
1972  float x,
1973  float Mean,
1974  float StdDev) {
1975  float X;
1976 
1977  // wraparound circular parameters if necessary
1978  if (ParamDesc->Circular) {
1979  if (x - Mean > ParamDesc->HalfRange)
1980  x -= ParamDesc->Range;
1981  else if (x - Mean < -ParamDesc->HalfRange)
1982  x += ParamDesc->Range;
1983  }
1984 
1985  X = ((x - Mean) / (2 * StdDev) * BUCKETTABLESIZE + BUCKETTABLESIZE / 2.0);
1986  if (X < 0)
1987  return 0;
1988  if (X > BUCKETTABLESIZE - 1)
1989  return static_cast<uint16_t>(BUCKETTABLESIZE - 1);
1990  return static_cast<uint16_t>(floor(static_cast<double>(X)));
1991 } // UniformBucket
1992 
2003 static bool DistributionOK(BUCKETS* Buckets) {
2004  float FrequencyDifference;
2005  float TotalDifference;
2006  int i;
2007 
2008  // compute how well the histogram matches the expected histogram
2009  TotalDifference = 0.0;
2010  for (i = 0; i < Buckets->NumberOfBuckets; i++) {
2011  FrequencyDifference = Buckets->Count[i] - Buckets->ExpectedCount[i];
2012  TotalDifference += (FrequencyDifference * FrequencyDifference) /
2013  Buckets->ExpectedCount[i];
2014  }
2015 
2016  // test to see if the difference is more than expected
2017  if (TotalDifference > Buckets->ChiSquared)
2018  return false;
2019  else
2020  return true;
2021 } // DistributionOK
2022 
2028 static void FreeStatistics(STATISTICS *Statistics) {
2029  free(Statistics->CoVariance);
2030  free(Statistics->Min);
2031  free(Statistics->Max);
2032  free(Statistics);
2033 } // FreeStatistics
2034 
2040 static void FreeBuckets(BUCKETS *buckets) {
2041  Efree(buckets->Count);
2042  Efree(buckets->ExpectedCount);
2043  Efree(buckets);
2044 } // FreeBuckets
2045 
2053 static void FreeCluster(CLUSTER *Cluster) {
2054  if (Cluster != nullptr) {
2055  FreeCluster (Cluster->Left);
2056  FreeCluster (Cluster->Right);
2057  free(Cluster);
2058  }
2059 } // FreeCluster
2060 
2073 static uint16_t DegreesOfFreedom(DISTRIBUTION Distribution, uint16_t HistogramBuckets) {
2074  static uint8_t DegreeOffsets[] = { 3, 3, 1 };
2075 
2076  uint16_t AdjustedNumBuckets;
2077 
2078  AdjustedNumBuckets = HistogramBuckets - DegreeOffsets[static_cast<int>(Distribution)];
2079  if (Odd (AdjustedNumBuckets))
2080  AdjustedNumBuckets++;
2081  return (AdjustedNumBuckets);
2082 
2083 } // DegreesOfFreedom
2084 
2092 static void AdjustBuckets(BUCKETS *Buckets, uint32_t NewSampleCount) {
2093  int i;
2094  double AdjustFactor;
2095 
2096  AdjustFactor = ((static_cast<double>(NewSampleCount)) /
2097  (static_cast<double>(Buckets->SampleCount)));
2098 
2099  for (i = 0; i < Buckets->NumberOfBuckets; i++) {
2100  Buckets->ExpectedCount[i] *= AdjustFactor;
2101  }
2102 
2103  Buckets->SampleCount = NewSampleCount;
2104 
2105 } // AdjustBuckets
2106 
2112 static void InitBuckets(BUCKETS *Buckets) {
2113  int i;
2114 
2115  for (i = 0; i < Buckets->NumberOfBuckets; i++) {
2116  Buckets->Count[i] = 0;
2117  }
2118 
2119 } // InitBuckets
2120 
2133 static int AlphaMatch(void *arg1, //CHISTRUCT *ChiStruct,
2134  void *arg2) { //CHISTRUCT *SearchKey)
2135  auto *ChiStruct = static_cast<CHISTRUCT *>(arg1);
2136  auto *SearchKey = static_cast<CHISTRUCT *>(arg2);
2137 
2138  return (ChiStruct->Alpha == SearchKey->Alpha);
2139 
2140 } // AlphaMatch
2141 
2151 static CHISTRUCT *NewChiStruct(uint16_t DegreesOfFreedom, double Alpha) {
2152  CHISTRUCT *NewChiStruct;
2153 
2154  NewChiStruct = static_cast<CHISTRUCT *>(Emalloc (sizeof (CHISTRUCT)));
2155  NewChiStruct->DegreesOfFreedom = DegreesOfFreedom;
2156  NewChiStruct->Alpha = Alpha;
2157  return (NewChiStruct);
2158 
2159 } // NewChiStruct
2160 
2174 static double
2175 Solve (SOLVEFUNC Function,
2176 void *FunctionParams, double InitialGuess, double Accuracy)
2177 #define INITIALDELTA 0.1
2178 #define DELTARATIO 0.1
2179 {
2180  double x;
2181  double f;
2182  double Slope;
2183  double Delta;
2184  double NewDelta;
2185  double xDelta;
2186  double LastPosX, LastNegX;
2187 
2188  x = InitialGuess;
2189  Delta = INITIALDELTA;
2190  LastPosX = FLT_MAX;
2191  LastNegX = -FLT_MAX;
2192  f = (*Function) (static_cast<CHISTRUCT *>(FunctionParams), x);
2193  while (Abs (LastPosX - LastNegX) > Accuracy) {
2194  // keep track of outer bounds of current estimate
2195  if (f < 0)
2196  LastNegX = x;
2197  else
2198  LastPosX = x;
2199 
2200  // compute the approx. slope of f(x) at the current point
2201  Slope =
2202  ((*Function) (static_cast<CHISTRUCT *>(FunctionParams), x + Delta) - f) / Delta;
2203 
2204  // compute the next solution guess */
2205  xDelta = f / Slope;
2206  x -= xDelta;
2207 
2208  // reduce the delta used for computing slope to be a fraction of
2209  //the amount moved to get to the new guess
2210  NewDelta = Abs (xDelta) * DELTARATIO;
2211  if (NewDelta < Delta)
2212  Delta = NewDelta;
2213 
2214  // compute the value of the function at the new guess
2215  f = (*Function) (static_cast<CHISTRUCT *>(FunctionParams), x);
2216  }
2217  return (x);
2218 
2219 } // Solve
2220 
2239 static double ChiArea(CHISTRUCT *ChiParams, double x) {
2240  int i, N;
2241  double SeriesTotal;
2242  double Denominator;
2243  double PowerOfx;
2244 
2245  N = ChiParams->DegreesOfFreedom / 2 - 1;
2246  SeriesTotal = 1;
2247  Denominator = 1;
2248  PowerOfx = 1;
2249  for (i = 1; i <= N; i++) {
2250  Denominator *= 2 * i;
2251  PowerOfx *= x;
2252  SeriesTotal += PowerOfx / Denominator;
2253  }
2254  return ((SeriesTotal * exp (-0.5 * x)) - ChiParams->Alpha);
2255 
2256 } // ChiArea
2257 
2281 static bool
2282 MultipleCharSamples(CLUSTERER* Clusterer,
2283  CLUSTER* Cluster, float MaxIllegal)
2284 #define ILLEGAL_CHAR 2
2285 {
2286  static std::vector<uint8_t> CharFlags;
2287  LIST SearchState;
2288  SAMPLE *Sample;
2289  int32_t CharID;
2290  int32_t NumCharInCluster;
2291  int32_t NumIllegalInCluster;
2292  float PercentIllegal;
2293 
2294  // initial estimate assumes that no illegal chars exist in the cluster
2295  NumCharInCluster = Cluster->SampleCount;
2296  NumIllegalInCluster = 0;
2297 
2298  if (Clusterer->NumChar > CharFlags.size()) {
2299  CharFlags.resize(Clusterer->NumChar);
2300  }
2301 
2302  for (auto& CharFlag : CharFlags)
2303  CharFlag = false;
2304 
2305  // find each sample in the cluster and check if we have seen it before
2306  InitSampleSearch(SearchState, Cluster);
2307  while ((Sample = NextSample (&SearchState)) != nullptr) {
2308  CharID = Sample->CharID;
2309  if (CharFlags[CharID] == false) {
2310  CharFlags[CharID] = true;
2311  }
2312  else {
2313  if (CharFlags[CharID] == true) {
2314  NumIllegalInCluster++;
2315  CharFlags[CharID] = ILLEGAL_CHAR;
2316  }
2317  NumCharInCluster--;
2318  PercentIllegal = static_cast<float>(NumIllegalInCluster) / NumCharInCluster;
2319  if (PercentIllegal > MaxIllegal) {
2320  destroy(SearchState);
2321  return true;
2322  }
2323  }
2324  }
2325  return false;
2326 
2327 } // MultipleCharSamples
2328 
2334 static double InvertMatrix(const float* input, int size, float* inv) {
2335  // Allocate memory for the 2D arrays.
2336  GENERIC_2D_ARRAY<double> U(size, size, 0.0);
2337  GENERIC_2D_ARRAY<double> U_inv(size, size, 0.0);
2338  GENERIC_2D_ARRAY<double> L(size, size, 0.0);
2339 
2340  // Initialize the working matrices. U starts as input, L as I and U_inv as O.
2341  int row;
2342  int col;
2343  for (row = 0; row < size; row++) {
2344  for (col = 0; col < size; col++) {
2345  U[row][col] = input[row*size + col];
2346  L[row][col] = row == col ? 1.0 : 0.0;
2347  U_inv[row][col] = 0.0;
2348  }
2349  }
2350 
2351  // Compute forward matrix by inversion by LU decomposition of input.
2352  for (col = 0; col < size; ++col) {
2353  // Find best pivot
2354  int best_row = 0;
2355  double best_pivot = -1.0;
2356  for (row = col; row < size; ++row) {
2357  if (Abs(U[row][col]) > best_pivot) {
2358  best_pivot = Abs(U[row][col]);
2359  best_row = row;
2360  }
2361  }
2362  // Exchange pivot rows.
2363  if (best_row != col) {
2364  for (int k = 0; k < size; ++k) {
2365  double tmp = U[best_row][k];
2366  U[best_row][k] = U[col][k];
2367  U[col][k] = tmp;
2368  tmp = L[best_row][k];
2369  L[best_row][k] = L[col][k];
2370  L[col][k] = tmp;
2371  }
2372  }
2373  // Now do the pivot itself.
2374  for (row = col + 1; row < size; ++row) {
2375  double ratio = -U[row][col] / U[col][col];
2376  for (int j = col; j < size; ++j) {
2377  U[row][j] += U[col][j] * ratio;
2378  }
2379  for (int k = 0; k < size; ++k) {
2380  L[row][k] += L[col][k] * ratio;
2381  }
2382  }
2383  }
2384  // Next invert U.
2385  for (col = 0; col < size; ++col) {
2386  U_inv[col][col] = 1.0 / U[col][col];
2387  for (row = col - 1; row >= 0; --row) {
2388  double total = 0.0;
2389  for (int k = col; k > row; --k) {
2390  total += U[row][k] * U_inv[k][col];
2391  }
2392  U_inv[row][col] = -total / U[row][row];
2393  }
2394  }
2395  // Now the answer is U_inv.L.
2396  for (row = 0; row < size; row++) {
2397  for (col = 0; col < size; col++) {
2398  double sum = 0.0;
2399  for (int k = row; k < size; ++k) {
2400  sum += U_inv[row][k] * L[k][col];
2401  }
2402  inv[row*size + col] = sum;
2403  }
2404  }
2405  // Check matrix product.
2406  double error_sum = 0.0;
2407  for (row = 0; row < size; row++) {
2408  for (col = 0; col < size; col++) {
2409  double sum = 0.0;
2410  for (int k = 0; k < size; ++k) {
2411  sum += static_cast<double>(input[row * size + k]) * inv[k * size + col];
2412  }
2413  if (row != col) {
2414  error_sum += Abs(sum);
2415  }
2416  }
2417  }
2418  return error_sum;
2419 }
FTABLE_X
#define FTABLE_X
Definition: cluster.cpp:32
tesseract::GenericHeap
Definition: genericheap.h:58
BUCKETS
Definition: cluster.cpp:179
ClipToRange
T ClipToRange(const T &x, const T &lower_bound, const T &upper_bound)
Definition: helpers.h:106
emalloc.h
ClusteringContext::heap
ClusterHeap * heap
Definition: cluster.cpp:198
tesseract::GenericHeap::Pop
bool Pop(Pair *entry)
Definition: genericheap.h:118
CLUSTERCONFIG
Definition: cluster.h:45
PROTOTYPE::TotalMagnitude
float TotalMagnitude
Definition: cluster.h:74
CLUSTERCONFIG::MagicSamples
int MagicSamples
Definition: cluster.h:52
CLUSTERER::NumChar
int32_t NumChar
Definition: cluster.h:88
Mean
float Mean(PROTOTYPE *Proto, uint16_t Dimension)
Definition: cluster.cpp:602
CHIACCURACY
#define CHIACCURACY
FreeClusterer
void FreeClusterer(CLUSTERER *Clusterer)
Definition: cluster.cpp:514
sample::SampleCount
unsigned SampleCount
Definition: cluster.h:34
BUCKETS::ChiSquared
double ChiSquared
Definition: cluster.cpp:183
first_node
#define first_node(l)
Definition: oldlist.h:84
BUCKETS::Confidence
double Confidence
Definition: cluster.cpp:182
PARAM_DESC::Circular
bool Circular
Definition: ocrfeatures.h:42
Emalloc
void * Emalloc(int Size)
Definition: emalloc.cpp:31
destroy_nodes
void destroy_nodes(LIST list, void_dest destructor)
Definition: oldlist.cpp:138
elliptical
Definition: cluster.h:43
ClusteringContext
Definition: cluster.cpp:197
STATISTICS::AvgVariance
float AvgVariance
Definition: cluster.cpp:173
MakeKDTree
KDTREE * MakeKDTree(int16_t KeySize, const PARAM_DESC KeyDesc[])
Definition: kdtree.cpp:179
CLUSTERCONFIG::Independence
float Independence
Definition: cluster.h:50
PROTOTYPE::LogMagnitude
float LogMagnitude
Definition: cluster.h:75
ASSERT_HOST
#define ASSERT_HOST(x)
Definition: errcode.h:87
KDDelete
void KDDelete(KDTREE *Tree, float Key[], void *Data)
Definition: kdtree.cpp:252
BUCKETS::Distribution
DISTRIBUTION Distribution
Definition: cluster.cpp:180
list_rec
Definition: oldlist.h:73
STATISTICS
Definition: cluster.cpp:172
sample::Clustered
bool Clustered
Definition: cluster.h:32
SqrtOf2Pi
#define SqrtOf2Pi
Definition: cluster.cpp:219
language_specific.log
log
Definition: language_specific.py:25
sample::Mean
float Mean[1]
Definition: cluster.h:38
FreeProtoList
void FreeProtoList(LIST *ProtoList)
Definition: cluster.cpp:538
STATISTICS::Max
float * Max
Definition: cluster.cpp:176
PROTOTYPE::Magnitude
FLOATUNION Magnitude
Definition: cluster.h:77
normal
Definition: cluster.h:55
sample::Prototype
bool Prototype
Definition: cluster.h:33
SOLVEFUNC
double(*)(CHISTRUCT *, double) SOLVEFUNC
Definition: cluster.cpp:205
Config
CLUSTERCONFIG Config
Definition: commontraining.cpp:88
DENSITYFUNC
double(*)(int32_t) DENSITYFUNC
Definition: cluster.cpp:204
KDWalk
void KDWalk(KDTREE *Tree, void_proc action, void *context)
Definition: kdtree.cpp:314
CLUSTERER::SampleSize
int16_t SampleSize
Definition: cluster.h:82
CLUSTERCONFIG::ProtoStyle
PROTOSTYLE ProtoStyle
Definition: cluster.h:46
StandardDeviation
float StandardDeviation(PROTOTYPE *Proto, uint16_t Dimension)
Definition: cluster.cpp:613
INITIALDELTA
#define INITIALDELTA
KDNearestNeighborSearch
void KDNearestNeighborSearch(KDTREE *Tree, float Query[], int QuerySize, float MaxDistance, int *NumberOfResults, void **NBuffer, float DBuffer[])
Definition: kdtree.cpp:304
Odd
#define Odd(N)
Definition: cluster.cpp:207
PARAM_DESC::Min
float Min
Definition: ocrfeatures.h:44
CLUSTERER::Root
CLUSTER * Root
Definition: cluster.h:86
MAXBUCKETS
#define MAXBUCKETS
Definition: cluster.h:26
PARAM_DESC::Range
float Range
Definition: ocrfeatures.h:46
MAXDISTANCE
#define MAXDISTANCE
NIL_LIST
#define NIL_LIST
Definition: oldlist.h:68
GENERIC_2D_ARRAY< double >
CLUSTERCONFIG::MaxIllegal
float MaxIllegal
Definition: cluster.h:48
CHISTRUCT::Alpha
double Alpha
Definition: cluster.cpp:192
tesseract::KDPair::key
Key key
Definition: kdpair.h:46
KDTREE
Definition: kdtree.h:47
PROTOTYPE
Definition: cluster.h:62
PARAM_DESC::MidRange
float MidRange
Definition: ocrfeatures.h:48
sample::Right
struct sample * Right
Definition: cluster.h:36
ClusterSamples
LIST ClusterSamples(CLUSTERER *Clusterer, CLUSTERCONFIG *Config)
Definition: cluster.cpp:483
CLUSTERCONFIG::Confidence
double Confidence
Definition: cluster.h:51
CLUSTERER::ParamDesc
PARAM_DESC * ParamDesc
Definition: cluster.h:83
PROTOTYPE::Merged
bool Merged
Definition: cluster.h:64
Abs
#define Abs(N)
Definition: cluster.cpp:209
CLUSTERER::NumberOfSamples
int32_t NumberOfSamples
Definition: cluster.h:84
BUCKETS::NumberOfBuckets
uint16_t NumberOfBuckets
Definition: cluster.cpp:184
ClusteringContext::tree
KDTREE * tree
Definition: cluster.cpp:200
uniform
Definition: cluster.h:55
BUCKETS::SampleCount
uint32_t SampleCount
Definition: cluster.cpp:181
FLOATUNION::Elliptical
float * Elliptical
Definition: cluster.h:59
DELTARATIO
#define DELTARATIO
BUCKETS::Bucket
uint16_t Bucket[BUCKETTABLESIZE]
Definition: cluster.cpp:185
CHISTRUCT::ChiSquared
double ChiSquared
Definition: cluster.cpp:193
NextSample
CLUSTER * NextSample(LIST *SearchState)
Definition: cluster.cpp:580
FTABLE_Y
#define FTABLE_Y
Definition: cluster.cpp:33
PROTOTYPE::Weight
FLOATUNION Weight
Definition: cluster.h:78
matrix.h
tesseract::KDPair::data
Data data
Definition: kdpair.h:45
FreeKDTree
void FreeKDTree(KDTREE *Tree)
Definition: kdtree.cpp:330
KDStore
void KDStore(KDTREE *Tree, float *Key, void *Data)
Definition: kdtree.cpp:211
DISTRIBUTION
DISTRIBUTION
Definition: cluster.h:55
MergeClusters
int32_t MergeClusters(int16_t N, PARAM_DESC ParamDesc[], int32_t n1, int32_t n2, float m[], float m1[], float m2[])
Definition: cluster.cpp:824
ILLEGAL_CHAR
#define ILLEGAL_CHAR
mixed
Definition: cluster.h:43
MAXNEIGHBORS
#define MAXNEIGHBORS
CHISTRUCT::DegreesOfFreedom
uint16_t DegreesOfFreedom
Definition: cluster.cpp:191
MAXDEGREESOFFREEDOM
#define MAXDEGREESOFFREEDOM
Definition: cluster.cpp:230
MakeClusterer
CLUSTERER * MakeClusterer(int16_t SampleSize, const PARAM_DESC ParamDesc[])
Definition: cluster.cpp:376
PARAM_DESC::Max
float Max
Definition: ocrfeatures.h:45
RootOf
#define RootOf(T)
Definition: kdtree.h:56
helpers.h
kdpair.h
PROTOSTYLE
PROTOSTYLE
Definition: cluster.h:43
ClusterHeap
tesseract::GenericHeap< ClusterPair > ClusterHeap
Definition: cluster.cpp:170
FTable
const double FTable[FTABLE_Y][FTABLE_X]
Definition: cluster.cpp:36
push
LIST push(LIST list, void *element)
Definition: oldlist.cpp:172
ClusteringContext::next
int32_t next
Definition: cluster.cpp:201
STATISTICS::Min
float * Min
Definition: cluster.cpp:175
MINBUCKETS
#define MINBUCKETS
Definition: cluster.h:25
sample::CharID
int32_t CharID
Definition: cluster.h:37
tprintf.h
MINSAMPLESNEEDED
#define MINSAMPLESNEEDED
Definition: cluster.cpp:153
InitSampleSearch
#define InitSampleSearch(S, C)
Definition: cluster.h:100
sample
Definition: cluster.h:31
MINSAMPLES
#define MINSAMPLES
Definition: cluster.cpp:152
tesseract::KDPairInc
Definition: kdpair.h:51
PARAM_DESC
Definition: ocrfeatures.h:41
tesseract::GenericHeap::Push
void Push(Pair *entry)
Definition: genericheap.h:95
cluster.h
PROTOTYPE::Significant
bool Significant
Definition: cluster.h:63
PARAM_DESC::NonEssential
bool NonEssential
Definition: ocrfeatures.h:43
CHISTRUCT
Definition: cluster.cpp:190
PROTOTYPE::Mean
float * Mean
Definition: cluster.h:73
TEMPCLUSTER
Definition: cluster.cpp:164
FLOATUNION::Spherical
float Spherical
Definition: cluster.h:58
CLUSTERER
Definition: cluster.h:81
PARAM_DESC::HalfRange
float HalfRange
Definition: ocrfeatures.h:47
Efree
void Efree(void *ptr)
Definition: emalloc.cpp:45
iterate
#define iterate(l)
Definition: oldlist.h:92
FreePrototype
void FreePrototype(void *arg)
Definition: cluster.cpp:549
spherical
Definition: cluster.h:43
destroy
LIST destroy(LIST list)
Definition: oldlist.cpp:123
PROTOTYPE::Style
unsigned Style
Definition: cluster.h:69
pop
LIST pop(LIST list)
Definition: oldlist.cpp:161
tprintf
DLLSYM void tprintf(const char *format,...)
Definition: tprintf.cpp:34
TEMPCLUSTER::Neighbor
CLUSTER * Neighbor
Definition: cluster.cpp:166
Mirror
#define Mirror(N, R)
Definition: cluster.cpp:208
MINVARIANCE
#define MINVARIANCE
Definition: cluster.cpp:143
PROTOTYPE::Variance
FLOATUNION Variance
Definition: cluster.h:76
STATISTICS::CoVariance
float * CoVariance
Definition: cluster.cpp:174
BUCKETS::ExpectedCount
float * ExpectedCount
Definition: cluster.cpp:187
automatic
Definition: cluster.h:43
TEMPCLUSTER::Cluster
CLUSTER * Cluster
Definition: cluster.cpp:165
HOTELLING
#define HOTELLING
Definition: cluster.cpp:31
CLUSTERER::bucket_cache
BUCKETS * bucket_cache[DISTRIBUTION_COUNT][MAXBUCKETS+1 - MINBUCKETS]
Definition: cluster.h:90
sample::Left
struct sample * Left
Definition: cluster.h:35
NORMALEXTENT
#define NORMALEXTENT
Definition: cluster.cpp:162
CLUSTERCONFIG::MinSamples
float MinSamples
Definition: cluster.h:47
BUCKETTABLESIZE
#define BUCKETTABLESIZE
Definition: cluster.cpp:161
PROTOTYPE::NumSamples
unsigned NumSamples
Definition: cluster.h:70
DISTRIBUTION_COUNT
Definition: cluster.h:55
LOOKUPTABLESIZE
#define LOOKUPTABLESIZE
Definition: cluster.cpp:229
PROTOTYPE::Cluster
CLUSTER * Cluster
Definition: cluster.h:71
search
LIST search(LIST list, void *key, int_compare is_equal)
Definition: oldlist.cpp:202
CLUSTERER::ProtoList
LIST ProtoList
Definition: cluster.h:87
MakeSample
SAMPLE * MakeSample(CLUSTERER *Clusterer, const float *Feature, int32_t CharID)
Definition: cluster.cpp:429
genericheap.h
D_random
Definition: cluster.h:55
BUCKETS::Count
uint32_t * Count
Definition: cluster.cpp:186
CLUSTERER::KDTree
KDTREE * KDTree
Definition: cluster.h:85
PROTOTYPE::Distrib
DISTRIBUTION * Distrib
Definition: cluster.h:72
MINALPHA
#define MINALPHA
ClusteringContext::candidates
TEMPCLUSTER * candidates
Definition: cluster.cpp:199