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