tesseract 4.1.1
Loading...
Searching...
No Matches
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 "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%.
36const 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
167};
168
171
175 float *Min; // largest negative distance from the mean
176 float *Max; // largest positive distance from the mean
177};
178
179struct 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
192 double Alpha;
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
204using DENSITYFUNC = double (*)(int32_t);
205using 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
220static const double kNormalStdDev = BUCKETTABLESIZE / (2.0 * NORMALEXTENT);
221static const double kNormalVariance =
223static const double kNormalMagnitude =
225static const double kNormalMean = BUCKETTABLESIZE / 2;
226
229#define LOOKUPTABLESIZE 8
230#define MAXDEGREESOFFREEDOM MAXBUCKETS
231
232static const uint32_t kCountTable[LOOKUPTABLESIZE] = {
233 MINSAMPLES, 200, 400, 600, 800, 1000, 1500, 2000
234}; // number of samples
235
236static 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--------------------------------------------------------------------------*/
243static void CreateClusterTree(CLUSTERER* Clusterer);
244
245static void MakePotentialClusters(ClusteringContext* context, CLUSTER* Cluster,
246 int32_t Level);
247
248static CLUSTER* FindNearestNeighbor(KDTREE*Tree, CLUSTER* Cluster,
249 float* Distance);
250
251static CLUSTER* MakeNewCluster(CLUSTERER* Clusterer, TEMPCLUSTER* TempCluster);
252
253static void ComputePrototypes(CLUSTERER* Clusterer, CLUSTERCONFIG* Config);
254
255static PROTOTYPE* MakePrototype(CLUSTERER* Clusterer, CLUSTERCONFIG* Config,
256 CLUSTER* Cluster);
257
258static PROTOTYPE* MakeDegenerateProto(uint16_t N,
259 CLUSTER* Cluster, STATISTICS* Statistics,
260 PROTOSTYLE Style, int32_t MinSamples);
261
262static PROTOTYPE* TestEllipticalProto(CLUSTERER* Clusterer,
263 CLUSTERCONFIG* Config, CLUSTER* Cluster,
264 STATISTICS* Statistics);
265
266static PROTOTYPE* MakeSphericalProto(CLUSTERER* Clusterer,
267 CLUSTER* Cluster, STATISTICS* Statistics,
268 BUCKETS* Buckets);
269
270static PROTOTYPE* MakeEllipticalProto(CLUSTERER* Clusterer,
271 CLUSTER* Cluster, STATISTICS* Statistics,
272 BUCKETS* Buckets);
273
274static PROTOTYPE* MakeMixedProto(CLUSTERER* Clusterer,
275 CLUSTER* Cluster, STATISTICS* Statistics,
276 BUCKETS* NormalBuckets, double Confidence);
277
278static void MakeDimRandom(uint16_t i, PROTOTYPE* Proto, PARAM_DESC* ParamDesc);
279
280static void MakeDimUniform(uint16_t i, PROTOTYPE* Proto, STATISTICS* Statistics);
281
282static STATISTICS* ComputeStatistics(int16_t N, PARAM_DESC ParamDesc[],
283 CLUSTER* Cluster);
284
285static PROTOTYPE* NewSphericalProto(uint16_t N, CLUSTER* Cluster,
286 STATISTICS* Statistics);
287
288static PROTOTYPE* NewEllipticalProto(int16_t N, CLUSTER* Cluster,
289 STATISTICS* Statistics);
290
291static PROTOTYPE* NewMixedProto(int16_t N, CLUSTER *Cluster, STATISTICS *Statistics);
292
293static PROTOTYPE* NewSimpleProto(int16_t N, CLUSTER *Cluster);
294
295static bool Independent(PARAM_DESC* ParamDesc,
296 int16_t N, float* CoVariance, float Independence);
297
298static BUCKETS *GetBuckets(CLUSTERER* clusterer,
299 DISTRIBUTION Distribution,
300 uint32_t SampleCount,
301 double Confidence);
302
303static BUCKETS *MakeBuckets(DISTRIBUTION Distribution,
304 uint32_t SampleCount,
305 double Confidence);
306
307static uint16_t OptimumNumberOfBuckets(uint32_t SampleCount);
308
309static double ComputeChiSquared(uint16_t DegreesOfFreedom, double Alpha);
310
311static double NormalDensity(int32_t x);
312
313static double UniformDensity(int32_t x);
314
315static double Integral(double f1, double f2, double Dx);
316
317static void FillBuckets(BUCKETS *Buckets,
318 CLUSTER *Cluster,
319 uint16_t Dim,
320 PARAM_DESC *ParamDesc,
321 float Mean,
322 float StdDev);
323
324static uint16_t NormalBucket(PARAM_DESC *ParamDesc,
325 float x,
326 float Mean,
327 float StdDev);
328
329static uint16_t UniformBucket(PARAM_DESC *ParamDesc,
330 float x,
331 float Mean,
332 float StdDev);
333
334static bool DistributionOK(BUCKETS* Buckets);
335
336static void FreeStatistics(STATISTICS *Statistics);
337
338static void FreeBuckets(BUCKETS *Buckets);
339
340static void FreeCluster(CLUSTER *Cluster);
341
342static uint16_t DegreesOfFreedom(DISTRIBUTION Distribution, uint16_t HistogramBuckets);
343
344static void AdjustBuckets(BUCKETS *Buckets, uint32_t NewSampleCount);
345
346static void InitBuckets(BUCKETS *Buckets);
347
348static int AlphaMatch(void *arg1, // CHISTRUCT *ChiStruct,
349 void *arg2); // CHISTRUCT *SearchKey);
350
351static CHISTRUCT *NewChiStruct(uint16_t DegreesOfFreedom, double Alpha);
352
353static double Solve(SOLVEFUNC Function,
354 void *FunctionParams,
355 double InitialGuess,
356 double Accuracy);
357
358static double ChiArea(CHISTRUCT *ChiParams, double x);
359
360static bool MultipleCharSamples(CLUSTERER* Clusterer,
361 CLUSTER* Cluster,
362 float MaxIllegal);
363
364static double InvertMatrix(const float* input, int size, float* inv);
365
366//--------------------------Public Code--------------------------------------
375CLUSTERER *
376MakeClusterer (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
429SAMPLE* 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
514void 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
538void FreeProtoList(LIST *ProtoList) {
539 destroy_nodes(*ProtoList, FreePrototype);
540} // FreeProtoList
541
549void 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
580CLUSTER *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
602float Mean(PROTOTYPE *Proto, uint16_t Dimension) {
603 return (Proto->Mean[Dimension]);
604} // Mean
605
613float 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]);
627 ASSERT_HOST(!"Distribution count not allowed!");
628 }
629 }
630 return 0.0f;
631} // StandardDeviation
632
633
634/*---------------------------------------------------------------------------
635 Private Code
636----------------------------------------------------------------------------*/
650static 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
716static 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
745static CLUSTER*
746FindNearestNeighbor(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
781static 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
824int32_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
865static 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
908static 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
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,
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,
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,
972 break;
973 }
974 FreeStatistics(Statistics);
975 return Proto;
976} // MakePrototype
977
997static 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
1039static 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
1139static 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
1173static 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
1212static 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
1268static 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
1289static 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
1322static STATISTICS*
1323ComputeStatistics (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
1408static 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)
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
1438static 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
1479static 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
1502static 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
1540static bool
1541Independent(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
1586static 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
1631static 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
1726static 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
1760static double
1761ComputeChiSquared (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
1812static 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
1826static 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
1844static double Integral(double f1, double f2, double Dx) {
1845 return (f1 + f2) * Dx / 2.0;
1846} // Integral
1847
1868static 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
1938static 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
1971static 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
2003static 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
2028static void FreeStatistics(STATISTICS *Statistics) {
2029 free(Statistics->CoVariance);
2030 free(Statistics->Min);
2031 free(Statistics->Max);
2032 free(Statistics);
2033} // FreeStatistics
2034
2040static void FreeBuckets(BUCKETS *buckets) {
2041 Efree(buckets->Count);
2042 Efree(buckets->ExpectedCount);
2043 Efree(buckets);
2044} // FreeBuckets
2045
2053static void FreeCluster(CLUSTER *Cluster) {
2054 if (Cluster != nullptr) {
2055 FreeCluster (Cluster->Left);
2056 FreeCluster (Cluster->Right);
2057 free(Cluster);
2058 }
2059} // FreeCluster
2060
2073static 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
2092static 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
2112static 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
2133static 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
2151static 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
2174static double
2175Solve (SOLVEFUNC Function,
2176void *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
2239static 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
2281static bool
2282MultipleCharSamples(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
2334static 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}
#define ASSERT_HOST(x)
Definition: errcode.h:88
T ClipToRange(const T &x, const T &lower_bound, const T &upper_bound)
Definition: helpers.h:108
DLLSYM void tprintf(const char *format,...)
Definition: tprintf.cpp:35
#define MAXDISTANCE
#define MINALPHA
#define Abs(N)
Definition: cluster.cpp:209
#define Odd(N)
Definition: cluster.cpp:207
#define BUCKETTABLESIZE
Definition: cluster.cpp:161
#define MAXDEGREESOFFREEDOM
Definition: cluster.cpp:230
const double FTable[FTABLE_Y][FTABLE_X]
Definition: cluster.cpp:36
#define FTABLE_X
Definition: cluster.cpp:32
tesseract::GenericHeap< ClusterPair > ClusterHeap
Definition: cluster.cpp:170
#define INITIALDELTA
float Mean(PROTOTYPE *Proto, uint16_t Dimension)
Definition: cluster.cpp:602
#define Mirror(N, R)
Definition: cluster.cpp:208
double(*)(CHISTRUCT *, double) SOLVEFUNC
Definition: cluster.cpp:205
#define LOOKUPTABLESIZE
Definition: cluster.cpp:229
#define MINSAMPLES
Definition: cluster.cpp:152
#define DELTARATIO
#define SqrtOf2Pi
Definition: cluster.cpp:219
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
float StandardDeviation(PROTOTYPE *Proto, uint16_t Dimension)
Definition: cluster.cpp:613
#define MINSAMPLESNEEDED
Definition: cluster.cpp:153
#define HOTELLING
Definition: cluster.cpp:31
void FreeProtoList(LIST *ProtoList)
Definition: cluster.cpp:538
#define MAXNEIGHBORS
#define ILLEGAL_CHAR
void FreeClusterer(CLUSTERER *Clusterer)
Definition: cluster.cpp:514
void FreePrototype(void *arg)
Definition: cluster.cpp:549
#define CHIACCURACY
#define MINVARIANCE
Definition: cluster.cpp:143
LIST ClusterSamples(CLUSTERER *Clusterer, CLUSTERCONFIG *Config)
Definition: cluster.cpp:483
CLUSTER * NextSample(LIST *SearchState)
Definition: cluster.cpp:580
#define FTABLE_Y
Definition: cluster.cpp:33
#define NORMALEXTENT
Definition: cluster.cpp:162
double(*)(int32_t) DENSITYFUNC
Definition: cluster.cpp:204
SAMPLE * MakeSample(CLUSTERER *Clusterer, const float *Feature, int32_t CharID)
Definition: cluster.cpp:429
CLUSTERER * MakeClusterer(int16_t SampleSize, const PARAM_DESC ParamDesc[])
Definition: cluster.cpp:376
#define MAXBUCKETS
Definition: cluster.h:27
PROTOSTYLE
Definition: cluster.h:44
@ elliptical
Definition: cluster.h:44
@ spherical
Definition: cluster.h:44
@ automatic
Definition: cluster.h:44
@ mixed
Definition: cluster.h:44
#define MINBUCKETS
Definition: cluster.h:26
DISTRIBUTION
Definition: cluster.h:56
@ DISTRIBUTION_COUNT
Definition: cluster.h:56
@ D_random
Definition: cluster.h:56
@ uniform
Definition: cluster.h:56
@ normal
Definition: cluster.h:56
#define InitSampleSearch(S, C)
Definition: cluster.h:101
void KDStore(KDTREE *Tree, float *Key, void *Data)
Definition: kdtree.cpp:212
void FreeKDTree(KDTREE *Tree)
Definition: kdtree.cpp:331
void KDNearestNeighborSearch(KDTREE *Tree, float Query[], int QuerySize, float MaxDistance, int *NumberOfResults, void **NBuffer, float DBuffer[])
Definition: kdtree.cpp:305
void KDDelete(KDTREE *Tree, float Key[], void *Data)
Definition: kdtree.cpp:253
void KDWalk(KDTREE *Tree, void_proc action, void *context)
Definition: kdtree.cpp:315
KDTREE * MakeKDTree(int16_t KeySize, const PARAM_DESC KeyDesc[])
Definition: kdtree.cpp:180
void(*)(...) void_proc
Definition: kdtree.h:26
#define RootOf(T)
Definition: kdtree.h:57
void * Emalloc(int Size)
Definition: emalloc.cpp:31
void Efree(void *ptr)
Definition: emalloc.cpp:45
void destroy_nodes(LIST list, void_dest destructor)
Definition: oldlist.cpp:157
LIST search(LIST list, void *key, int_compare is_equal)
Definition: oldlist.cpp:258
LIST destroy(LIST list)
Definition: oldlist.cpp:141
LIST pop(LIST list)
Definition: oldlist.cpp:201
LIST push(LIST list, void *element)
Definition: oldlist.cpp:213
#define iterate(l)
Definition: oldlist.h:101
#define first_node(l)
Definition: oldlist.h:92
#define NIL_LIST
Definition: oldlist.h:76
CLUSTERCONFIG Config
bool Pop(Pair *entry)
Definition: genericheap.h:118
void Push(Pair *entry)
Definition: genericheap.h:95
CLUSTER * Neighbor
Definition: cluster.cpp:166
CLUSTER * Cluster
Definition: cluster.cpp:165
float * Max
Definition: cluster.cpp:176
float AvgVariance
Definition: cluster.cpp:173
float * Min
Definition: cluster.cpp:175
float * CoVariance
Definition: cluster.cpp:174
double ChiSquared
Definition: cluster.cpp:183
float * ExpectedCount
Definition: cluster.cpp:187
uint16_t NumberOfBuckets
Definition: cluster.cpp:184
uint32_t SampleCount
Definition: cluster.cpp:181
double Confidence
Definition: cluster.cpp:182
uint32_t * Count
Definition: cluster.cpp:186
DISTRIBUTION Distribution
Definition: cluster.cpp:180
uint16_t Bucket[BUCKETTABLESIZE]
Definition: cluster.cpp:185
double Alpha
Definition: cluster.cpp:192
double ChiSquared
Definition: cluster.cpp:193
uint16_t DegreesOfFreedom
Definition: cluster.cpp:191
ClusterHeap * heap
Definition: cluster.cpp:198
TEMPCLUSTER * candidates
Definition: cluster.cpp:199
Definition: cluster.h:32
struct sample * Right
Definition: cluster.h:37
int32_t CharID
Definition: cluster.h:38
bool Prototype
Definition: cluster.h:34
unsigned SampleCount
Definition: cluster.h:35
bool Clustered
Definition: cluster.h:33
float Mean[1]
Definition: cluster.h:39
struct sample * Left
Definition: cluster.h:36
int MagicSamples
Definition: cluster.h:53
float Independence
Definition: cluster.h:51
PROTOSTYLE ProtoStyle
Definition: cluster.h:47
float MinSamples
Definition: cluster.h:48
float MaxIllegal
Definition: cluster.h:49
double Confidence
Definition: cluster.h:52
float Spherical
Definition: cluster.h:59
float * Elliptical
Definition: cluster.h:60
FLOATUNION Magnitude
Definition: cluster.h:78
unsigned NumSamples
Definition: cluster.h:71
FLOATUNION Variance
Definition: cluster.h:77
unsigned Style
Definition: cluster.h:70
float * Mean
Definition: cluster.h:74
float LogMagnitude
Definition: cluster.h:76
bool Significant
Definition: cluster.h:64
bool Merged
Definition: cluster.h:65
float TotalMagnitude
Definition: cluster.h:75
DISTRIBUTION * Distrib
Definition: cluster.h:73
FLOATUNION Weight
Definition: cluster.h:79
CLUSTER * Cluster
Definition: cluster.h:72
CLUSTER * Root
Definition: cluster.h:87
int32_t NumberOfSamples
Definition: cluster.h:85
int16_t SampleSize
Definition: cluster.h:83
BUCKETS * bucket_cache[DISTRIBUTION_COUNT][MAXBUCKETS+1 - MINBUCKETS]
Definition: cluster.h:91
int32_t NumChar
Definition: cluster.h:89
KDTREE * KDTree
Definition: cluster.h:86
PARAM_DESC * ParamDesc
Definition: cluster.h:84
LIST ProtoList
Definition: cluster.h:88
Definition: kdtree.h:48
float HalfRange
Definition: ocrfeatures.h:48
float Range
Definition: ocrfeatures.h:47
bool Circular
Definition: ocrfeatures.h:43
float Max
Definition: ocrfeatures.h:46
float MidRange
Definition: ocrfeatures.h:49
bool NonEssential
Definition: ocrfeatures.h:44
float Min
Definition: ocrfeatures.h:45