fastlib2/contrib/niche/kernel_vector_mult.h

@@ -131,9 +131,9 @@ class KernelVectorMult {
 
 
   // kernel function
-  EpanKernel epan_kernel;
+  GaussianKernel kernel;
 
-  double epan_cutoff_dist_;
+  double cutoff_dist_;
 
 
 
@@ -276,7 +276,7 @@ class KernelVectorMult {
             la::DistanceSqEuclidean(query_point, reference_point);
 
 	weighted_sums_[query_index] +=
-	  weights_[reference_index] * epan_kernel.EvalUnnormOnSq(distance);
+	  weights_[reference_index] * kernel.EvalUnnormOnSq(distance);
 
       } /* for reference_index */
 
@@ -305,7 +305,7 @@ class KernelVectorMult {
     DEBUG_SAME_DOUBLE(lower_bound_distance,
         MinNodeDistSq_(query_node, reference_node));
 
-    if (lower_bound_distance > epan_cutoff_dist_) {
+    if (0) {//lower_bound_distance > cutoff_dist_) {
 
       // execute prune
 
@@ -465,11 +465,11 @@ class KernelVectorMult {
     DEBUG_ASSERT(bandwidth > 0);
 
     // kernel function
-    epan_kernel.Init(bandwidth);
+    kernel.Init(bandwidth);
     
-    epan_cutoff_dist_ = bandwidth*bandwidth;
+    cutoff_dist_ = bandwidth*bandwidth;
     
-    DEBUG_ONLY(printf("epan_cutoff_dist = %f\n", epan_cutoff_dist_));
+    DEBUG_ONLY(printf("cutoff_dist = %f\n", cutoff_dist_));
   
 
   } /* Init */
@@ -502,7 +502,7 @@ class KernelVectorMult {
     GNPRecursion_(query_tree_, reference_tree_,
         MinNodeDistSq_(query_tree_, reference_tree_));
     //printf("queries_.n_rows() = %d\n", queries_.n_rows());
-    la::Scale(1 / epan_kernel.CalcNormConstant(queries_.n_rows()), &weighted_sums_);
+    la::Scale(1 / kernel.CalcNormConstant(queries_.n_rows()), &weighted_sums_);
 
     fx_timer_stop(module_, "dual_tree_computation");
 

fastlib2/contrib/niche/solvelinsys.h

@@ -30,6 +30,7 @@ private:
 
   // used for a simple linear transformation operator data_^T * data  
   Matrix data_;
+  float** K_;
   //Matrix K_;
 
   //EpanKernel epan_kernel_;
@@ -55,6 +56,7 @@ public:
 
 
 
+    
 
 
 
@@ -73,9 +75,30 @@ public:
 
     sigma_squared_ = sigma_squared_in;
 
-    /*
+    
     // for debugging purposes, explicitly represent kernel matrix K
-    K_.Init(n_points_, n_points_);
+
+    if((K_ = (float**) malloc(40000 * sizeof(float*))) == NULL) {
+      printf("failed\n");
+      exit(1);
+    }
+
+    for(int i = 0; i < 40000; i++) {
+      if((K_[i] = (float*) malloc(40000 * sizeof(float))) == NULL) {
+	printf("failed on %d\n", i);
+	exit(1);
+      }
+    }
+
+    printf("start writing\n");
+    for(int i = 0; i < 40000; i++) {
+      for(int j = 0; j < 40000; j++) {
+	K_[i][j] = 0;
+      }
+    }
+    printf("done writing\n");
+
+
     for(int i = 0; i < n_points_; i++) {
       Vector v_i;
       data_.MakeColumnVector(i, &v_i);
@@ -84,22 +107,18 @@ public:
 	Vector v_j;
 	data_.MakeColumnVector(j, &v_j);
 	
-	//double dist = la::Dot(v_i, v_j);
-	//K_i[j] = dist;
-	
 	double dist = la::DistanceSqEuclidean(v_i, v_j);
-	K_.set(i, j, epan_kernel_.EvalUnnormOnSq(dist));
+	K_[i][j] = kernel_.EvalUnnormOnSq(dist) / norm_constant_;
       }
     }
 
-    la::Scale(1 / norm_constant_, &K_);
     for(int i = 0; i < n_points_; i++) {
-      K_.set(i, i, K_.get(i, i) + sigma_squared_);
+      K_[i][i] += sigma_squared_;
     }
     
 
-    const char *K_file_name = "K.txt";
-    data::Save(K_file_name, K_);
+    //const char *K_file_name = "K.txt";
+    //data::Save(K_file_name, K_);
 
     
 //     int errors = 0;
@@ -116,7 +135,7 @@ public:
     
 
     // end debugging explicit representation of K
-    */
+    
 
 
   }
@@ -344,8 +363,8 @@ void SolveLinearSystem(Matrix references, Vector rhs, double bandwidth, double s
   //iterative_solver.SetAztecOption(AZ_precond, AZ_Jacobi);
   
   // Use Conjugate Gradient
-  //iterative_solver.SetAztecOption(AZ_solver, AZ_cg);
-  iterative_solver.SetAztecOption(AZ_solver, AZ_gmres);
+  iterative_solver.SetAztecOption(AZ_solver, AZ_cg);
+  //iterative_solver.SetAztecOption(AZ_solver, AZ_gmres);
   //iterative_solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
   
   // Use modified Gram-Schmidt.