PLearn 0.1
CorrelationKernel.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // CorrelationKernel.cc
00004 //
00005 // Copyright (C) 2005 Olivier Delalleau 
00006 // 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 /* *******************************************************      
00036  * $Id: CorrelationKernel.cc 4462 2005-11-07 17:27:21Z tihocan $ 
00037  ******************************************************* */
00038 
00039 // Authors: Olivier Delalleau
00040 
00044 #include "CorrelationKernel.h"
00045 #include <plearn/math/TMat_sort.h>          
00046 #include <plearn/vmat/TransposeVMatrix.h>
00047 #include <plearn/vmat/VMat_basic_stats.h>
00048 
00049 namespace PLearn {
00050 using namespace std;
00051 
00053 // CorrelationKernel //
00055 CorrelationKernel::CorrelationKernel() 
00056     : correlation("linear"),
00057       transform(""),
00058       var_threshold(0)
00059 {
00060     transform_prg = new VMatLanguage();
00061 }
00062 
00063 PLEARN_IMPLEMENT_OBJECT(CorrelationKernel,
00064                         "Compute a similarity measure between two variables by their correlation.",
00065                         "Here, the two examples x and y in K(x,y) are understood as being samples\n"
00066                         "of a two random variables.\n"
00067                         "If the 'correlation' option is set to 'mutual_information_norm', the formula used\n"
00068                         "is the following:\n"
00069                         " M = - ln |sin(u,v)| where u and v are the centered variables\n"
00070                         " H_u = 1/2 + ln(sqrt(2Pi) * (stddev(u) + epsilon)) with epsilon = e^{1/2} / sqrt(2Pi)\n"
00071                         " M_norm = M / (sqrt(H_u) sqrt(H_v))\n"
00072     );
00073 
00075 // declareOptions //
00077 void CorrelationKernel::declareOptions(OptionList& ol)
00078 {
00079     // ### Declare all of this object's options here
00080     // ### For the "flags" of each option, you should typically specify  
00081     // ### one of OptionBase::buildoption, OptionBase::learntoption or 
00082     // ### OptionBase::tuningoption. Another possible flag to be combined with
00083     // ### is OptionBase::nosave
00084 
00085     declareOption(ol, "correlation", &CorrelationKernel::correlation, OptionBase::buildoption,
00086                   "The correlation method used to compute the similarity, among:\n"
00087                   "- 'linear'                  : linear correlation\n"
00088                   "- 'mutual_information_norm' : normalized mutual information\n");
00089 
00090     declareOption(ol, "transform", &CorrelationKernel::transform, OptionBase::buildoption,
00091                   "An additional transformation applied on the similarity, in VPL language.\n");
00092 
00093     declareOption(ol, "var_threshold", &CorrelationKernel::var_threshold, OptionBase::buildoption,
00094                   "If set to a value > 0, denote by x_i the training point whose variance is the\n"
00095                   "'var_threshold' quantile of all training variances. If v_i is its variance, then\n"
00096                   "all pairs of points whose product of variances is less than v_i^2 will be given a\n"
00097                   "similarity (correlation) of 1e-6.\n");
00098 
00099     // Now call the parent class' declareOptions
00100     inherited::declareOptions(ol);
00101 }
00102 
00104 // build //
00106 void CorrelationKernel::build()
00107 {
00108     inherited::build();
00109     build_();
00110 }
00111 
00113 // build_ //
00115 void CorrelationKernel::build_()
00116 {
00117     // ### This method should do the real building of the object,
00118     // ### according to set 'options', in *any* situation. 
00119     // ### Typical situations include:
00120     // ###  - Initial building of an object from a few user-specified options
00121     // ###  - Building of a "reloaded" object: i.e. from the complete set of all serialised options.
00122     // ###  - Updating or "re-building" of an object after a few "tuning" options have been modified.
00123     // ### You should assume that the parent class' build_() has already been called.
00124     transform_prg->setSourceFieldNames(TVec<string>(1)); // Dummy fieldnames.
00125     transform_prg->compileString(transform, transform_prg_fields);
00126     result_vec.resize(1);
00127     result_transformed_vec.resize(1);
00128 }
00129 
00131 // evaluate //
00133 real CorrelationKernel::evaluate(const Vec& x1, const Vec& x2) const {
00134     real result = 0;
00135     int n = x1.length();
00136 #ifdef BOUNDCHECK
00137     if (x2.length() != n)
00138         PLERROR("In CorrelationKernel::evaluate - x1 and x2 must have same size");
00139 #endif
00140     if (correlation == "linear") {
00141         VMat x1_(x1.toMat(n, 1));
00142         VMat x2_(x2.toMat(n, 1));
00143         x1_->defineSizes(n, 0, 0);
00144         x2_->defineSizes(n, 0, 0);
00145         correlations(x1_, x2_, correl, pvalues, true);
00146         result = correl(0,0);
00147     } else if (correlation == "mutual_information_norm") {
00148         Vec x1_c(n);
00149         Vec x2_c(n);
00150         x1_c << x1;
00151         x2_c << x2;
00152         x1_c -= mean(x1);
00153         x2_c -= mean(x2);
00154         real cos12 = dot(x1_c, x2_c) / (norm(x1_c) * norm(x2_c));
00155         if (cos12 > 1)
00156             cos12 = 1;
00157         real sin12 = sqrt(1 - cos12 * cos12);   // |sin(x1,x2)|
00158         if (fast_is_equal(sin12, 0.0, 1.0, 1e-5))
00159             // The variables are considered equal.
00160             result = 1;
00161         else {
00162             real m = - pl_log(sin12);
00163             real epsilon = exp(0.5) / sqrt(2 * Pi);
00164             real sigma1 = sqrt(variance(x1_c.toMat(1,n), 0.0));
00165             real sigma2 = sqrt(variance(x2_c.toMat(1,n), 0.0));
00166             real h1 = 0.5 + pl_log(sqrt(2 * Pi) * (sigma1 + epsilon));
00167             real h2 = 0.5 + pl_log(sqrt(2 * Pi) * (sigma2 + epsilon));
00168             if (h1 <= 0 || h2 <= 0)
00169                 PLERROR("In CorrelationKernel::evaluate - Entropy should always be > 0");
00170             result = m / (sqrt(h1) * sqrt(h2));
00171         }
00172     } else
00173         PLERROR("In CorrelationKernel::evaluate - Unknown value for 'correlation': "
00174                 "%s", correlation.c_str());
00175     if (var_threshold > 0) {
00176         real v_1 = variance(x1, mean(x1));
00177         real v_2 = variance(x2, mean(x2));
00178         if (v_1 * v_2 < min_product_var)
00179             result = 1e-6;
00180     }
00181     result_vec[0] = result;
00182     transform_prg->run(result_vec, result_transformed_vec);
00183     return result_transformed_vec[0];
00184     /*
00185       if (transform.empty())
00186       return result;
00187       else if (transform == "minus_log")
00188       return -log(result);
00189       else if (transform == "minus_log_abs")
00190       return -log(fabs(result));
00191       else {
00192       PLERROR("In CorrelationKernel::evaluate - Unknown value for 'transform': "
00193       "%s", transform.c_str());
00194       return 0; // To make the compiler happy.
00195       }
00196     */
00197 }
00198 
00199 /* ### This method will very often be overridden.
00201 // evaluate_i_j //
00203 real CorrelationKernel::evaluate_i_j(int i, int j) const {
00204 // ### Evaluate the kernel on a pair of training points.
00205 }
00206 */
00207 
00209 // makeDeepCopyFromShallowCopy //
00211 void CorrelationKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00212 {
00213     inherited::makeDeepCopyFromShallowCopy(copies);
00214     deepCopyField(correl, copies);
00215     deepCopyField(pvalues, copies);
00216     deepCopyField(mean_vec, copies);
00217     deepCopyField(var_vec, copies);
00218     deepCopyField(transform_prg, copies);
00219     deepCopyField(transform_prg_fields, copies);
00220     deepCopyField(result_vec, copies);
00221     deepCopyField(result_transformed_vec, copies);
00222 
00223 }
00224 
00226 // setDataForKernelMatrix //
00228 void CorrelationKernel::setDataForKernelMatrix(VMat the_data) {
00229     inherited::setDataForKernelMatrix(the_data);
00230     if (var_threshold > 0) {
00231         // Compute variances.
00232         VMat transp = new TransposeVMatrix(the_data);
00233         computeMeanAndVariance(transp, mean_vec, var_vec);
00234         // Compute quantile.
00235         sortElements(var_vec);
00236         int q = (int) floor(var_threshold * (var_vec.length() - 1));
00237         min_product_var = var_vec[q] * var_vec[q];
00238     }
00239 }
00240 
00241 } // end of namespace PLearn
00242 
00243 
00244 /*
00245   Local Variables:
00246   mode:c++
00247   c-basic-offset:4
00248   c-file-style:"stroustrup"
00249   c-file-offsets:((innamespace . 0)(inline-open . 0))
00250   indent-tabs-mode:nil
00251   fill-column:79
00252   End:
00253 */
00254 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines