001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.ml.kernel;
031
032import org.openimaj.citation.annotation.Reference;
033import org.openimaj.citation.annotation.ReferenceType;
034import org.openimaj.citation.annotation.References;
035import org.openimaj.feature.DoubleFV;
036import org.openimaj.feature.FeatureExtractor;
037import org.openimaj.feature.FeatureVector;
038import org.openimaj.math.util.MathUtils;
039import org.openimaj.math.util.MathUtils.ExponentAndMantissa;
040
041/**
042 * Implementation of the Homogeneous Kernel Map. The Homogeneous Kernel Map
043 * transforms data into a compact linear representation such that applying a
044 * linear SVM can approximate to a high degree of accuracy the application of a
045 * non-linear SVM to the original data. Additive kernels including Chi2,
046 * intersection, and Jensen-Shannon are supported.
047 * <p>
048 * This implementation is based directly on the VLFeat implementation written by
049 * Andrea Verdaldi, although it has been refactored to better fit with Java
050 * conventions.
051 *
052 * @see "http://www.vlfeat.org/api/homkermap.html"
053 * @see "http://www.robots.ox.ac.uk/~vgg/software/homkermap/"
054 *
055 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
056 * @author Based on code originally written by Andrea Verdaldi
057 */
058@References(
059                references = {
060                                @Reference(
061                                                type = ReferenceType.Article,
062                                                author = { "Vedaldi, A.", "Zisserman, A." },
063                                                title = "Efficient Additive Kernels via Explicit Feature Maps",
064                                                year = "2012",
065                                                journal = "Pattern Analysis and Machine Intelligence, IEEE Transactions on",
066                                                pages = { "480", "492" },
067                                                number = "3",
068                                                volume = "34",
069                                                customData = {
070                                                                "keywords",
071                                                                "approximation theory;computer vision;data handling;feature extraction;learning (artificial intelligence);spectral analysis;support vector machines;Nystrom approximation;additive homogeneous kernels;approximate finite-dimensional feature maps;approximation error;computer vision;data dependency;explicit feature maps;exponential decay;large scale nonlinear support vector machines;linear SVM;spectral analysis;Additives;Approximation methods;Histograms;Kernel;Measurement;Support vector machines;Training;Kernel methods;feature map;large scale learning;object detection.;object recognition",
072                                                                "doi", "10.1109/TPAMI.2011.153",
073                                                                "ISSN", "0162-8828"
074                                                }),
075                                @Reference(
076                                                type = ReferenceType.Inproceedings,
077                                                author = { "A. Vedaldi", "A. Zisserman" },
078                                                title = "Efficient Additive Kernels via Explicit Feature Maps",
079                                                year = "2010",
080                                                booktitle = "Proceedings of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)")
081                })
082public class HomogeneousKernelMap {
083        /**
084         * Types of supported kernel for the {@link HomogeneousKernelMap}
085         *
086         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
087         *
088         */
089        public enum KernelType {
090                /**
091                 * Intersection kernel
092                 *
093                 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
094                 *
095                 */
096                Intersection {
097                        @Override
098                        protected double getSpectrum(double omega) {
099                                return (2.0 / Math.PI) / (1 + 4 * omega * omega);
100                        }
101                },
102                /**
103                 * Chi^2 kernel
104                 *
105                 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
106                 *
107                 */
108                Chi2 {
109                        @Override
110                        protected double getSpectrum(double omega) {
111                                return 2.0 / (Math.exp(Math.PI * omega) + Math.exp(-Math.PI * omega));
112                        }
113                },
114                /**
115                 * Jenson-Shannon Kernel
116                 *
117                 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
118                 *
119                 */
120                JensonShannon {
121                        @Override
122                        protected double getSpectrum(double omega) {
123                                return (2.0 / Math.log(4.0)) * 2.0 / (Math.exp(Math.PI * omega) + Math.exp(-Math.PI * omega))
124                                                / (1 + 4 * omega * omega);
125                        }
126                };
127
128                protected abstract double getSpectrum(double omega);
129        }
130
131        /**
132         * Types of window supported by the {@link HomogeneousKernelMap}.
133         *
134         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
135         *
136         */
137        public enum WindowType {
138                /**
139                 * Uniform window
140                 *
141                 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
142                 *
143                 */
144                Uniform {
145                        @Override
146                        double computeKappaHat(KernelType type, double omega, HomogeneousKernelMap map) {
147                                return type.getSpectrum(omega);
148                        }
149                },
150                /**
151                 * Rectangular window
152                 *
153                 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
154                 *
155                 */
156                Rectangular {
157                        @Override
158                        double computeKappaHat(KernelType type, double omega, HomogeneousKernelMap map) {
159                                double kappa_hat = 0;
160                                final double epsilon = 1e-2;
161                                final double omegaRange = 2.0 / (map.period * epsilon);
162                                final double domega = 2 * omegaRange / (2 * 1024.0 + 1);
163
164                                for (double omegap = -omegaRange; omegap <= omegaRange; omegap += domega) {
165                                        double win = sinc((map.period / 2.0) * omegap);
166                                        win *= (map.period / (2.0 * Math.PI));
167                                        kappa_hat += win * type.getSpectrum(omegap + omega);
168                                }
169
170                                kappa_hat *= domega;
171
172                                // project on the positive orthant (see PAMI)
173                                kappa_hat = Math.max(kappa_hat, 0.0);
174
175                                return kappa_hat;
176                        }
177
178                        private double sinc(double x) {
179                                if (x == 0.0)
180                                        return 1.0;
181                                return Math.sin(x) / x;
182                        }
183                };
184
185                abstract double computeKappaHat(KernelType type, double omega, HomogeneousKernelMap map);
186
187                protected double getSmoothSpectrum(double omega, HomogeneousKernelMap map) {
188                        return computeKappaHat(map.kernelType, omega, map);
189                }
190        }
191
192        /**
193         * Helper implementation of a {@link FeatureExtractor} that wraps another
194         * {@link FeatureExtractor} and then applies the {@link HomogeneousKernelMap} to
195         * the output before returning the vector.
196         *
197         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
198         *
199         * @param <T>
200         *            Type of object that features can be extracted from
201         */
202        public static class ExtractorWrapper<T> implements FeatureExtractor<DoubleFV, T> {
203                private FeatureExtractor<? extends FeatureVector, T> inner;
204                private HomogeneousKernelMap map;
205
206                /**
207                 * Construct with the given internal extractor and homogeneous kernel map.
208                 *
209                 * @param inner
210                 *            the internal extractor
211                 * @param map
212                 *            the homogeneous kernel map
213                 */
214                public ExtractorWrapper(FeatureExtractor<? extends FeatureVector, T> inner, HomogeneousKernelMap map) {
215                        this.inner = inner;
216                        this.map = map;
217                }
218
219                @Override
220                public DoubleFV extractFeature(T object) {
221                        return map.evaluate(inner.extractFeature(object).asDoubleFV());
222                }
223        }
224
225        private KernelType kernelType;
226        private double period;
227        private double gamma;
228        private int order;
229        private int numSubdivisions;
230        private double subdivision;
231        private int minExponent;
232        private int maxExponent;
233        private double[] table;
234
235        /**
236         * Construct with the given kernel and window. The Gamma and order values are
237         * set at their defaults of 1. The period is computed automatically.
238         *
239         * @param kernelType
240         *            the type of kernel
241         * @param windowType
242         *            the type of window (use {@link WindowType#Rectangular} if unsure)
243         */
244        public HomogeneousKernelMap(KernelType kernelType, WindowType windowType) {
245                this(kernelType, 1, 1, -1, windowType);
246        }
247
248        /**
249         * Construct with the given kernel, gamma and window. The period is computed
250         * automatically and the approximation order is set to 1.
251         *
252         * @param kernelType
253         *            the type of kernel
254         * @param gamma
255         *            the gamma value. the standard kernels are 1-homogeneous, but
256         *            smaller values can work better in practice.
257         * @param windowType
258         *            the type of window (use {@link WindowType#Rectangular} if unsure)
259         */
260        public HomogeneousKernelMap(KernelType kernelType,
261                        double gamma,
262                        WindowType windowType)
263        {
264                this(kernelType, gamma, 1, -1, windowType);
265        }
266
267        /**
268         * Construct with the given kernel, gamma, order and window. The period is
269         * computed automatically.
270         *
271         * @param kernelType
272         *            the type of kernel
273         * @param gamma
274         *            the gamma value. the standard kernels are 1-homogeneous, but
275         *            smaller values can work better in practice.
276         * @param order
277         *            the approximation order (usually 1 is enough)
278         * @param windowType
279         *            the type of window (use {@link WindowType#Rectangular} if unsure)
280         */
281        public HomogeneousKernelMap(KernelType kernelType,
282                        double gamma,
283                        int order,
284                        WindowType windowType)
285        {
286                this(kernelType, gamma, order, -1, windowType);
287        }
288
289        /**
290         * Construct with the given kernel, gamma, order, period and window. If the
291         * period is negative, it will be replaced by the default.
292         *
293         * @param kernelType
294         *            the type of kernel
295         * @param gamma
296         *            the gamma value. the standard kernels are 1-homogeneous, but
297         *            smaller values can work better in practice.
298         * @param order
299         *            the approximation order (usually 1 is enough)
300         * @param period
301         *            the periodicity of the kernel spectrum
302         * @param windowType
303         *            the type of window (use {@link WindowType#Rectangular} if unsure)
304         */
305        public HomogeneousKernelMap(KernelType kernelType,
306                        double gamma,
307                        int order,
308                        double period,
309                        WindowType windowType)
310        {
311                if (gamma <= 0)
312                        throw new IllegalArgumentException("Gamma must be > 0");
313                final int tableWidth, tableHeight;
314
315                this.kernelType = kernelType;
316                this.gamma = gamma;
317                this.order = order;
318
319                if (period < 0) {
320                        period = computeDefaultPeriod(windowType, kernelType);
321                }
322                this.period = period;
323
324                this.numSubdivisions = 8 + 8 * order;
325                this.subdivision = 1.0 / this.numSubdivisions;
326                this.minExponent = -20;
327                this.maxExponent = 8;
328
329                tableHeight = 2 * this.order + 1;
330                tableWidth = this.numSubdivisions * (this.maxExponent - this.minExponent + 1);
331                this.table = new double[tableHeight * tableWidth + 2 * (1 + this.order)];
332
333                int tableOffset = 0;
334                final int kappaOffset = tableHeight * tableWidth;
335                final int freqOffset = kappaOffset + 1 + order;
336                final double L = 2.0 * Math.PI / this.period;
337
338                /* precompute the sampled periodicized spectrum */
339                int j = 0;
340                int i = 0;
341                while (i <= this.order) {
342                        table[freqOffset + i] = j;
343                        table[kappaOffset + i] = windowType.getSmoothSpectrum(j * L, this);
344                        j++;
345                        if (table[kappaOffset + i] > 0 || j >= 3 * i)
346                                i++;
347                }
348
349                /* fill table */
350                for (int exponent = minExponent; exponent <= maxExponent; exponent++) {
351                        double x, Lxgamma, Llogx, xgamma;
352                        double sqrt2kappaLxgamma;
353                        double mantissa = 1.0;
354
355                        for (i = 0; i < numSubdivisions; i++, mantissa += subdivision) {
356                                x = mantissa * Math.pow(2, exponent);
357                                xgamma = Math.pow(x, this.gamma);
358                                Lxgamma = L * xgamma;
359                                Llogx = L * Math.log(x);
360
361                                table[tableOffset++] = Math.sqrt(Lxgamma * table[kappaOffset]);
362                                for (j = 1; j <= this.order; ++j) {
363                                        sqrt2kappaLxgamma = Math.sqrt(2.0 * Lxgamma * table[kappaOffset + j]);
364                                        table[tableOffset++] = sqrt2kappaLxgamma * Math.cos(table[freqOffset + j] * Llogx);
365                                        table[tableOffset++] = sqrt2kappaLxgamma * Math.sin(table[freqOffset + j] * Llogx);
366                                }
367                        }
368                }
369        }
370
371        private double computeDefaultPeriod(WindowType windowType, KernelType kernelType) {
372                double period = 0;
373
374                // compute default period
375                switch (windowType) {
376                case Uniform:
377                        switch (kernelType) {
378                        case Chi2:
379                                period = 5.86 * Math.sqrt(order + 0) + 3.65;
380                                break;
381                        case JensonShannon:
382                                period = 6.64 * Math.sqrt(order + 0) + 7.24;
383                                break;
384                        case Intersection:
385                                period = 2.38 * Math.log(order + 0.8) + 5.6;
386                                break;
387                        }
388                        break;
389                case Rectangular:
390                        switch (kernelType) {
391                        case Chi2:
392                                period = 8.80 * Math.sqrt(order + 4.44) - 12.6;
393                                break;
394                        case JensonShannon:
395                                period = 9.63 * Math.sqrt(order + 1.00) - 2.93;
396                                break;
397                        case Intersection:
398                                period = 2.00 * Math.log(order + 0.99) + 3.52;
399                                break;
400                        }
401                        break;
402                }
403                return Math.max(period, 1.0);
404        }
405
406        /**
407         * Evaluate the kernel for the given <code>x</code> value. The output values
408         * will be written into the destination array at <code>offset + j*stride</code>
409         * intervals where <code>j</code> is between 0 and <code>2 * order + 1</code>.
410         *
411         * @param destination
412         *            the destination array
413         * @param stride
414         *            the stride
415         * @param offset
416         *            the offset
417         * @param x
418         *            the value to compute the kernel approximation for
419         */
420        public void evaluate(double[] destination, int stride, int offset, double x) {
421                final ExponentAndMantissa em = MathUtils.frexp(x);
422
423                double mantissa = em.mantissa;
424                int exponent = em.exponent;
425                final double sign = (mantissa >= 0.0) ? +1.0 : -1.0;
426                mantissa *= 2 * sign;
427                exponent--;
428
429                if (mantissa == 0 || exponent <= minExponent || exponent >= maxExponent) {
430                        for (int j = 0; j <= order; j++) {
431                                destination[offset + j * stride] = 0.0;
432                        }
433                        return;
434                }
435
436                final int featureDimension = 2 * order + 1;
437                int v1offset = (exponent - minExponent) * numSubdivisions * featureDimension;
438
439                mantissa -= 1.0;
440                while (mantissa >= subdivision) {
441                        mantissa -= subdivision;
442                        v1offset += featureDimension;
443                }
444
445                int v2offset = v1offset + featureDimension;
446                for (int j = 0; j < featureDimension; j++) {
447                        final double f1 = table[v1offset++];
448                        final double f2 = table[v2offset++];
449
450                        destination[offset + j * stride] = sign * ((f2 - f1) * (numSubdivisions * mantissa) + f1);
451                }
452        }
453
454        /**
455         * Compute the Homogeneous Kernel Map approximation of the given feature vector
456         *
457         * @param in
458         *            the feature vector
459         * @return the expanded feature vector
460         */
461        public DoubleFV evaluate(DoubleFV in) {
462                final int step = (2 * order + 1);
463                final DoubleFV out = new DoubleFV(step * in.length());
464
465                for (int i = 0; i < in.length(); i++) {
466                        evaluate(out.values, 1, i * step, in.values[i]);
467                }
468
469                return out;
470        }
471
472        /**
473         * Construct a new {@link ExtractorWrapper} that applies the map to features
474         * extracted by an internal extractor.
475         *
476         * @param inner
477         *            the internal extractor
478         * @return the wrapped {@link FeatureExtractor}
479         * @param <T>
480         *            Type of object that features can be extracted from
481         */
482        public <T> FeatureExtractor<DoubleFV, T> createWrappedExtractor(FeatureExtractor<? extends FeatureVector, T> inner) {
483                return new ExtractorWrapper<T>(inner, this);
484        }
485}