1 #include <string>
2 #include <vector>
3 #include <iostream>
4 #include "testlib/testlib_test.h"
5 #include "vpdl/vpdl_mixture_of.h"
6 #include "vpdl/vpdl_gaussian.h"
7 #include "vpdl/vpdl_gaussian_sphere.h"
8 #include "vpdl/vpdl_gaussian_indep.h"
9 #include "vnl/vnl_vector_fixed.h"
10 #include "vnl/vnl_vector.h"
11 #ifdef _MSC_VER
12 # include "vcl_msvc_warnings.h"
13 #endif
14
15
16 // function to sort by increasing weight
17 template <class dist_t>
18 bool
weight_less(const dist_t &,const typename vpdt_dist_traits<dist_t>::scalar_type & w1,const dist_t &,const typename vpdt_dist_traits<dist_t>::scalar_type & w2)19 weight_less(const dist_t & /*d1*/,
20 const typename vpdt_dist_traits<dist_t>::scalar_type & w1,
21 const dist_t & /*d2*/,
22 const typename vpdt_dist_traits<dist_t>::scalar_type & w2)
23 {
24 return w1 < w2;
25 }
26
27
28 template <class T>
29 void
test_mixture_of_type(T epsilon,const std::string & type_name)30 test_mixture_of_type(T epsilon, const std::string & type_name)
31 {
32 // test a bunch of instantiations to make sure they compile
33 {
34 vpdl_mixture_of<vpdl_gaussian_sphere<T, 3>> mix_gauss_s3;
35 vpdl_mixture_of<vpdl_gaussian_indep<T, 3>> mix_gauss_i3;
36 vpdl_mixture_of<vpdl_gaussian<T, 3>> mix_gauss_f3;
37 vpdl_mixture_of<vpdl_gaussian_sphere<T, 1>> mix_gauss_s1;
38 vpdl_mixture_of<vpdl_gaussian_indep<T, 1>> mix_gauss_i1;
39 vpdl_mixture_of<vpdl_gaussian<T, 1>> mix_gauss_f1;
40 vpdl_mixture_of<vpdl_gaussian_sphere<T>> mix_gauss_s;
41 vpdl_mixture_of<vpdl_gaussian_indep<T>> mix_gauss_i;
42 vpdl_mixture_of<vpdl_gaussian<T>> mix_gauss_f;
43 }
44
45 // create a few sample distributions
46 vnl_vector_fixed<T, 3> mean1(T(1), T(1), T(1));
47 vnl_vector_fixed<T, 3> var1(T(0.25), T(1), T(1));
48
49 vnl_vector_fixed<T, 3> mean2(T(2), T(0), T(-3));
50 vnl_vector_fixed<T, 3> var2(T(1), T(0.25), T(2));
51
52 vnl_vector_fixed<T, 3> mean3(T(1), T(0), T(0));
53 vnl_vector_fixed<T, 3> var3(T(0.5), T(1.5), T(0.25));
54
55 std::cout << "=================== fixed<3> =======================" << std::endl;
56 {
57 vpdl_gaussian_indep<T, 3> gauss1(mean1, var1);
58 vpdl_gaussian_indep<T, 3> gauss2(mean2, var2);
59 vpdl_gaussian_indep<T, 3> gauss3(mean3, var3);
60
61 typedef vpdl_gaussian_indep<T, 3> dist_t;
62 vpdl_mixture_of<dist_t> mixture;
63
64 TEST(("initial num_components <" + type_name + ">").c_str(), mixture.num_components(), 0);
65
66 mixture.insert(gauss1, T(0.2));
67 bool valid = mixture.num_components() == 1;
68 mixture.insert(gauss2, T(0.5));
69 valid = valid && mixture.num_components() == 2;
70 TEST(("num_components after insert <" + type_name + ">").c_str(), valid, true);
71
72 const dist_t & dist1 = mixture.distribution(0);
73 TEST(("distribution access <" + type_name + ">").c_str(),
74 dist1.mean() == gauss1.mean() && dist1.covariance() == gauss1.covariance(),
75 true);
76 TEST(("weight access <" + type_name + ">").c_str(), mixture.weight(1), T(0.5));
77
78 valid = mixture.remove_last();
79 TEST(("remove last <" + type_name + ">").c_str(), valid && mixture.num_components() == 1, true);
80
81 mixture.insert(gauss2, T(0.3));
82 mixture.insert(gauss3, T(0.15));
83 mixture.set_weight(0, T(0.05));
84 TEST(("set_weight <" + type_name + ">").c_str(), mixture.weight(0), T(0.05));
85
86 mixture.sort();
87 TEST(("default sort <" + type_name + ">").c_str(),
88 mixture.weight(0) == T(0.3) && mixture.weight(1) == T(0.15) && mixture.weight(2) == T(0.05) &&
89 &mixture.distribution(2) == &dist1,
90 true);
91
92 vpdl_mixture_of<dist_t> mixture2(mixture);
93 mixture.remove_last();
94 mixture.remove_last();
95 mixture.remove_last();
96 TEST(("copy constructor <" + type_name + ">").c_str(),
97 mixture2.num_components() == 3 && mixture2.weight(1) == T(0.15),
98 true);
99
100 vnl_vector_fixed<T, 3> pt(T(0), T(1.5), T(1));
101 T prob = T(0.1 * gauss1.prob_density(pt) + 0.6 * gauss2.prob_density(pt) + 0.3 * gauss3.prob_density(pt));
102 TEST_NEAR(("probability density <" + type_name + ">").c_str(), mixture2.prob_density(pt), prob, epsilon);
103
104 // test gradient virtual functions against numerical difference
105 vnl_vector_fixed<T, 3> g3;
106 T dp = std::sqrt(epsilon);
107 T den = mixture2.density(pt);
108 T den_x = mixture2.density(pt + vnl_vector_fixed<T, 3>(dp, 0, 0));
109 T den_y = mixture2.density(pt + vnl_vector_fixed<T, 3>(0, dp, 0));
110 T den_z = mixture2.density(pt + vnl_vector_fixed<T, 3>(0, 0, dp));
111 vnl_vector_fixed<T, 3> grad(den_x - den, den_y - den, den_z - den);
112 grad /= dp;
113 T density = mixture2.gradient_density(pt, g3);
114 TEST_NEAR(("gradient density <" + type_name + ">").c_str(), (g3 - grad).inf_norm(), 0, dp);
115 TEST_NEAR(("density <" + type_name + ">").c_str(), density, den, epsilon);
116
117 prob = T(0.1 * gauss1.cumulative_prob(pt) + 0.6 * gauss2.cumulative_prob(pt) + 0.3 * gauss3.cumulative_prob(pt));
118 TEST_NEAR(("cumulative probability <" + type_name + ">").c_str(), mixture2.cumulative_prob(pt), prob, epsilon);
119
120 vnl_vector_fixed<T, 3> pt2(T(10), T(5), T(8));
121 prob = T(0.1 * gauss1.box_prob(pt, pt2) + 0.6 * gauss2.box_prob(pt, pt2) + 0.3 * gauss3.box_prob(pt, pt2));
122 TEST_NEAR(("box probability <" + type_name + ">").c_str(), mixture2.box_prob(pt, pt2), prob, epsilon);
123
124 mixture2.sort(weight_less<dist_t>);
125 TEST(("sort by increasing weight <" + type_name + ">").c_str(),
126 mixture2.weight(0) == T(0.05) && mixture2.weight(1) == T(0.15) && mixture2.weight(2) == T(0.3),
127 true);
128
129 mixture2.normalize_weights();
130 TEST(("normalize <" + type_name + ">").c_str(),
131 mixture2.weight(0) == T(0.1) && mixture2.weight(1) == T(0.3) && mixture2.weight(2) == T(0.6),
132 true);
133 }
134
135 std::cout << "=================== scalar =======================" << std::endl;
136 {
137 vpdl_gaussian<T, 1> gauss1(mean1[0], var1[0]);
138 vpdl_gaussian<T, 1> gauss2(mean2[0], var2[0]);
139 vpdl_gaussian<T, 1> gauss3(mean3[0], var3[0]);
140
141 typedef vpdl_gaussian<T, 1> dist_t;
142 vpdl_mixture_of<dist_t> mixture;
143
144 TEST(("initial num_components <" + type_name + ">").c_str(), mixture.num_components(), 0);
145
146 mixture.insert(gauss1, T(0.2));
147 bool valid = mixture.num_components() == 1;
148 mixture.insert(gauss2, T(0.5));
149 valid = valid && mixture.num_components() == 2;
150 TEST(("num_components after insert <" + type_name + ">").c_str(), valid, true);
151
152 const dist_t & dist1 = mixture.distribution(0);
153 TEST(("distribution access <" + type_name + ">").c_str(),
154 dist1.mean() == gauss1.mean() && dist1.covariance() == gauss1.covariance(),
155 true);
156 TEST(("weight access <" + type_name + ">").c_str(), mixture.weight(1), T(0.5));
157
158 valid = mixture.remove_last();
159 TEST(("remove last <" + type_name + ">").c_str(), valid && mixture.num_components() == 1, true);
160
161 mixture.insert(gauss2, T(0.3));
162 mixture.insert(gauss3, T(0.15));
163 mixture.set_weight(0, T(0.05));
164 TEST(("set_weight <" + type_name + ">").c_str(), mixture.weight(0), T(0.05));
165
166 mixture.sort();
167 TEST(("default sort <" + type_name + ">").c_str(),
168 mixture.weight(0) == T(0.3) && mixture.weight(1) == T(0.15) && mixture.weight(2) == T(0.05) &&
169 &mixture.distribution(2) == &dist1,
170 true);
171
172 vpdl_mixture_of<dist_t> mixture2(mixture);
173 mixture.remove_last();
174 mixture.remove_last();
175 mixture.remove_last();
176 TEST(("copy constructor <" + type_name + ">").c_str(),
177 mixture2.num_components() == 3 && mixture2.weight(1) == T(0.15),
178 true);
179
180 T pt = T(1);
181 T prob = T(0.1 * gauss1.prob_density(pt) + 0.6 * gauss2.prob_density(pt) + 0.3 * gauss3.prob_density(pt));
182 TEST_NEAR(("probability density <" + type_name + ">").c_str(), mixture2.prob_density(pt), prob, epsilon);
183
184 // test gradient virtual functions against numerical difference
185 T g;
186 T dp = std::sqrt(epsilon);
187 T den = mixture2.density(pt);
188 T den_x = mixture2.density(pt + dp);
189 T grad = (den_x - den) / dp;
190 T density = mixture2.gradient_density(pt, g);
191 TEST_NEAR(("gradient density <" + type_name + ">").c_str(), g, grad, dp);
192 TEST_NEAR(("density <" + type_name + ">").c_str(), density, den, epsilon);
193
194 prob = T(0.1 * gauss1.cumulative_prob(pt) + 0.6 * gauss2.cumulative_prob(pt) + 0.3 * gauss3.cumulative_prob(pt));
195 TEST_NEAR(("cumulative probability <" + type_name + ">").c_str(), mixture2.cumulative_prob(pt), prob, epsilon);
196
197 T pt2 = T(10);
198 prob = T(0.1 * gauss1.box_prob(pt, pt2) + 0.6 * gauss2.box_prob(pt, pt2) + 0.3 * gauss3.box_prob(pt, pt2));
199 TEST_NEAR(("box probability <" + type_name + ">").c_str(), mixture2.box_prob(pt, pt2), prob, epsilon);
200
201 mixture2.sort(weight_less<dist_t>);
202 TEST(("sort by increasing weight <" + type_name + ">").c_str(),
203 mixture2.weight(0) == T(0.05) && mixture2.weight(1) == T(0.15) && mixture2.weight(2) == T(0.3),
204 true);
205
206 mixture2.normalize_weights();
207 TEST(("normalize <" + type_name + ">").c_str(),
208 mixture2.weight(0) == T(0.1) && mixture2.weight(1) == T(0.3) && mixture2.weight(2) == T(0.6),
209 true);
210
211 T mean = T(0.1 * gauss1.mean() + 0.6 * gauss2.mean() + 0.3 * gauss3.mean());
212 T cmp_mean;
213 mixture2.compute_mean(cmp_mean);
214 TEST_NEAR(("compute_mean <" + type_name + ">").c_str(), cmp_mean, mean, epsilon);
215
216 T var = T(0.1 * gauss1.covariance() + 0.6 * gauss2.covariance() + 0.3 * gauss3.covariance());
217 var += T(0.1 * gauss1.mean() * gauss1.mean() + 0.6 * gauss2.mean() * gauss2.mean() +
218 0.3 * gauss3.mean() * gauss3.mean());
219 var -= mean * mean;
220 T cmp_var;
221 mixture2.compute_covar(cmp_var);
222 TEST_NEAR(("compute_covar <" + type_name + ">").c_str(), cmp_var, var, epsilon);
223 }
224
225 std::cout << "=================== variable =======================" << std::endl;
226 {
227 vpdl_gaussian_indep<T> gauss1(mean1.as_ref(), var1.as_ref());
228 vpdl_gaussian_indep<T> gauss2(mean2.as_ref(), var2.as_ref());
229 vpdl_gaussian_indep<T> gauss3(mean3.as_ref(), var3.as_ref());
230
231 using dist_t = vpdl_gaussian_indep<T>;
232 vpdl_mixture_of<dist_t> mixture;
233
234 TEST(("initial num_components <" + type_name + ">").c_str(), mixture.num_components(), 0);
235 TEST(("initial dimension <" + type_name + ">").c_str(), mixture.dimension(), 0);
236
237 mixture.insert(gauss1, T(0.2));
238 bool valid = mixture.num_components() == 1;
239 mixture.insert(gauss2, T(0.5));
240 valid = valid && mixture.num_components() == 2;
241 TEST(("num_components after insert <" + type_name + ">").c_str(), valid, true);
242 TEST(("dimension after insert <" + type_name + ">").c_str(), mixture.dimension(), 3);
243
244 const dist_t & dist1 = mixture.distribution(0);
245 TEST(("distribution access <" + type_name + ">").c_str(),
246 dist1.mean() == gauss1.mean() && dist1.covariance() == gauss1.covariance(),
247 true);
248 TEST(("weight access <" + type_name + ">").c_str(), mixture.weight(1), T(0.5));
249
250 valid = mixture.remove_last();
251 TEST(("remove last <" + type_name + ">").c_str(), valid && mixture.num_components() == 1, true);
252
253 mixture.insert(gauss2, T(0.3));
254 mixture.insert(gauss3, T(0.15));
255 mixture.set_weight(0, T(0.05));
256 TEST(("set_weight <" + type_name + ">").c_str(), mixture.weight(0), T(0.05));
257
258 mixture.sort();
259 TEST(("default sort <" + type_name + ">").c_str(),
260 mixture.weight(0) == T(0.3) && mixture.weight(1) == T(0.15) && mixture.weight(2) == T(0.05) &&
261 &mixture.distribution(2) == &dist1,
262 true);
263
264 vpdl_mixture_of<dist_t> mixture2(mixture);
265 mixture.remove_last();
266 mixture.remove_last();
267 mixture.remove_last();
268 TEST(("copy constructor <" + type_name + ">").c_str(),
269 mixture2.num_components() == 3 && mixture2.weight(1) == T(0.15) && mixture2.dimension() == 3,
270 true);
271
272 vnl_vector_fixed<T, 3> pt(T(0), T(1.5), T(1));
273 T prob = T(0.1 * gauss1.prob_density(pt.as_ref()) + 0.6 * gauss2.prob_density(pt.as_ref()) +
274 0.3 * gauss3.prob_density(pt.as_ref()));
275 TEST_NEAR(("probability density <" + type_name + ">").c_str(), mixture2.prob_density(pt.as_ref()), prob, epsilon);
276
277 // test gradient virtual functions against numerical difference
278 vnl_vector<T> g;
279 T dp = std::sqrt(epsilon);
280 T den = mixture2.density(pt.as_ref());
281 T den_x = mixture2.density((pt + vnl_vector_fixed<T, 3>(dp, 0, 0)).as_ref());
282 T den_y = mixture2.density((pt + vnl_vector_fixed<T, 3>(0, dp, 0)).as_ref());
283 T den_z = mixture2.density((pt + vnl_vector_fixed<T, 3>(0, 0, dp)).as_ref());
284 vnl_vector_fixed<T, 3> grad(den_x - den, den_y - den, den_z - den);
285 grad /= dp;
286 T density = mixture2.gradient_density(pt.as_ref(), g);
287 TEST_NEAR(("gradient density <" + type_name + ">").c_str(), (g - grad).inf_norm(), 0, dp);
288 TEST_NEAR(("density <" + type_name + ">").c_str(), density, den, epsilon);
289
290 prob = T(0.1 * gauss1.cumulative_prob(pt.as_ref()) + 0.6 * gauss2.cumulative_prob(pt.as_ref()) +
291 0.3 * gauss3.cumulative_prob(pt.as_ref()));
292 TEST_NEAR(
293 ("cumulative probability <" + type_name + ">").c_str(), mixture2.cumulative_prob(pt.as_ref()), prob, epsilon);
294
295 vnl_vector_fixed<T, 3> pt2(T(10), T(5), T(8));
296 prob = T(0.1 * gauss1.box_prob(pt.as_ref(), pt2.as_ref()) + 0.6 * gauss2.box_prob(pt.as_ref(), pt2.as_ref()) +
297 0.3 * gauss3.box_prob(pt.as_ref(), pt2.as_ref()));
298 TEST_NEAR(
299 ("box probability <" + type_name + ">").c_str(), mixture2.box_prob(pt.as_ref(), pt2.as_ref()), prob, epsilon);
300
301 mixture2.sort(weight_less<dist_t>);
302 TEST(("sort by increasing weight <" + type_name + ">").c_str(),
303 mixture2.weight(0) == T(0.05) && mixture2.weight(1) == T(0.15) && mixture2.weight(2) == T(0.3),
304 true);
305
306 mixture2.normalize_weights();
307 TEST(("normalize <" + type_name + ">").c_str(),
308 mixture2.weight(0) == T(0.1) && mixture2.weight(1) == T(0.3) && mixture2.weight(2) == T(0.6),
309 true);
310 }
311 }
312
313
314 static void
test_mixture_of()315 test_mixture_of()
316 {
317 test_mixture_of_type(1e-5f, "float");
318 test_mixture_of_type(1e-13, "double");
319 }
320
321 TESTMAIN(test_mixture_of);
322