[mlpack-svn] r10179 - mlpack/trunk/src/contrib/nslagle/myKDE

[fastlab-svn at coffeetalk-1.cc.gatech.edu]

Author: nslagle
Date: 2011-11-08 11:23:17 -0500 (Tue, 08 Nov 2011)
New Revision: 10179

Modified:
   mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree.hpp
   mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp
Log:
mlpack/contrib/nslagle: the algorithm almost works...

Modified: mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree.hpp
===================================================================
--- mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree.hpp	2011-11-08 04:28:44 UTC (rev 10178)
+++ mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree.hpp	2011-11-08 16:23:17 UTC (rev 10179)
@@ -85,6 +85,7 @@
 
   void SetDefaults();
   double Priority(TTree* Q, TTree* T);
+  size_t NaiveLogLikelihood();
   size_t MultiBandwidthDualTree();
   void MultiBandwidthDualTreeBase(TTree* Q,
                                   TTree* T, size_t QIndex,

Modified: mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp
===================================================================
--- mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp	2011-11-08 04:28:44 UTC (rev 10178)
+++ mlpack/trunk/src/contrib/nslagle/myKDE/kde_dual_tree_impl.hpp	2011-11-08 16:23:17 UTC (rev 10179)
@@ -94,20 +94,25 @@
                            (log(DBL_EPSILON)/* + log(inverseBandwidths[bIndex])*/));
   }
 
-  upperBoundQPointByBandwidth.zeros(queryRoot->count(),bandwidthCount);
+  upperBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount);
+  lowerBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount);
   for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex)
   {
-    arma::vec col = upperBoundQPointByBandwidth.unsafe_col(bIndex);
+    arma::vec col = upperBoundQNodeByBandwidth.unsafe_col(bIndex);
     col.fill(referenceRoot->count() * inverseBandwidths[bIndex]);
+    arma::vec col2 = lowerBoundQNodeByBandwidth.unsafe_col(bIndex);
+    col2.fill(referenceRoot->count() * DBL_EPSILON);
   }
+
+  upperBoundQPointByBandwidth.zeros(queryRoot->count(),bandwidthCount);
   lowerBoundQPointByBandwidth.zeros(queryRoot->count(),bandwidthCount);
-  upperBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount);
   for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex)
   {
-    arma::vec col = upperBoundQNodeByBandwidth.unsafe_col(bIndex);
+    arma::vec col = upperBoundQPointByBandwidth.unsafe_col(bIndex);
     col.fill(referenceRoot->count() * inverseBandwidths[bIndex]);
+    arma::vec col2 = lowerBoundQPointByBandwidth.unsafe_col(bIndex);
+    col2.fill(referenceRoot->count() * DBL_EPSILON);
   }
-  lowerBoundQNodeByBandwidth.zeros(queryTreeSize,bandwidthCount);
 
   arma::vec dl;
   arma::vec du;
@@ -120,6 +125,7 @@
   ++nextAvailableNodeIndex;
   nodePriorityQueue.push(firstNode);
   size_t finalLevel = MultiBandwidthDualTree();
+  //size_t finalLevel = NaiveLogLikelihood();
   std::cout << "the best level is " << finalLevel << "\n";
 
   size_t maxIndex = -1;
@@ -164,6 +170,30 @@
 }
 
 template<typename TKernel, typename TTree>
+size_t KdeDualTree<TKernel, TTree>::NaiveLogLikelihood()
+{
+  for (size_t bIndex = 0; bIndex < bandwidthCount; ++bIndex)
+  {
+    bestLevelByBandwidth[bIndex] = 0;
+    double inverseBandwidth = inverseBandwidths[bIndex];
+    for (size_t q = queryRoot->begin(); q < queryRoot->end(); ++q)
+    {
+      double total = 0.0;
+      for (size_t t = referenceRoot->begin(); t < referenceRoot->end(); ++t)
+      {
+        arma::vec diff = queryData.unsafe_col(q) - referenceData.unsafe_col(t);
+        double dist = pow(arma::dot(diff,diff),0.5);
+        total += inverseBandwidth * kernel.Evaluate(inverseBandwidth * dist);
+      }
+      lowerBoundLevelByBandwidth(0,bIndex) -= log(DBL_EPSILON);
+      upperBoundLevelByBandwidth(0,bIndex) -= log(inverseBandwidth);
+      lowerBoundLevelByBandwidth(0,bIndex) += log(total / referenceRoot->count());
+      upperBoundLevelByBandwidth(0,bIndex) += log(total / referenceRoot->count());
+    }
+  }
+  return 0;
+}
+template<typename TKernel, typename TTree>
 size_t KdeDualTree<TKernel, TTree>::MultiBandwidthDualTree()
 {
   /* current level */
@@ -245,7 +275,6 @@
          *   to the current Q */
         if (fabs((du - dl)/(lowerBoundQNodeByBandwidth(QIndex, bIndex) + dl)) < delta)
         {
-          std::cout << "do we ever meet delta?\n";
           for (size_t q = Q->begin(); q < Q->end(); ++q)
           {
             lowerBoundQPointByBandwidth(q,bIndex) += deltaLower(bIndex);
@@ -616,7 +645,7 @@
   {
     /* subtract out the current log-likelihood amount for this Q node so we can readjust
      *   the Q node bounds by the current bandwidth */
-    std::cout << "BEFORE========\n" << lowerBoundLevelByBandwidth << std::endl << upperBoundLevelByBandwidth << std::endl;
+    //std::cout << "BEFORE========\n" << lowerBoundLevelByBandwidth << std::endl << upperBoundLevelByBandwidth << std::endl;
     upperBoundLevelByBandwidth(levelOfQ, bIndex) -=
         sizeOfQNode * log(upperBoundQNodeByBandwidth(QIndex, bIndex));
     double value = lowerBoundQNodeByBandwidth(QIndex, bIndex);
@@ -642,6 +671,19 @@
       {
         maximumUpper = upperBoundQPointByBandwidth(q,bIndex);
       }
+      upperBoundLevelByBandwidth(levelOfQ, bIndex) +=
+          log(upperBoundQPointByBandwidth(q, bIndex));
+      value = lowerBoundQPointByBandwidth(q, bIndex);
+      if (value > DBL_EPSILON)
+      {
+        lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
+            log(value);
+      }
+      else
+      {
+        lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
+            log(DBL_EPSILON);
+      }
     }
     /* adjust Q node bounds, then add the new quantities to the level by bandwidth
      *   log-likelihood bounds */
@@ -649,20 +691,20 @@
     //std::cout << "maximum upper - TNodeSize = " << maximumUpper - sizeOfTNode << "\n";
     upperBoundQNodeByBandwidth(QIndex, bIndex) = maximumUpper;/* -
                                    inverseBandwidths[bIndex] * sizeOfTNode;*/
-    upperBoundLevelByBandwidth(levelOfQ, bIndex) +=
-        sizeOfQNode * log(upperBoundQNodeByBandwidth(QIndex, bIndex));
-    value = lowerBoundQNodeByBandwidth(QIndex, bIndex);
-    if (value > DBL_EPSILON)
-    {
-      lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
-          sizeOfQNode * log(value);
-    }
-    else
-    {
-      lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
-          sizeOfQNode * log(DBL_EPSILON);
-    }
-    std::cout << "AFTER========\n" << lowerBoundLevelByBandwidth << std::endl << upperBoundLevelByBandwidth << std::endl;
+//    upperBoundLevelByBandwidth(levelOfQ, bIndex) +=
+//        sizeOfQNode * log(upperBoundQNodeByBandwidth(QIndex, bIndex));
+//    value = lowerBoundQNodeByBandwidth(QIndex, bIndex);
+//    if (value > DBL_EPSILON)
+//    {
+//      lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
+//          sizeOfQNode * log(value);
+//    }
+//    else
+//    {
+//      lowerBoundLevelByBandwidth(levelOfQ, bIndex) +=
+//          sizeOfQNode * log(DBL_EPSILON);
+//    }
+//    std::cout << "AFTER========\n" << lowerBoundLevelByBandwidth << std::endl << upperBoundLevelByBandwidth << std::endl;
   }
 }
 template<typename TKernel, typename TTree>