diff --git a/fastlib/u/niche/entropy/kee/kde_risk.m b/fastlib/u/niche/entropy/kee/kde_risk.m
index 60fa14f69e..440725b36e 100644
--- a/fastlib/u/niche/entropy/kee/kde_risk.m
+++ b/fastlib/u/niche/entropy/kee/kde_risk.m
@@ -2,18 +2,21 @@
 % estimate for the given data and kernel bandwidth
 function E_risk = kde_risk(h, X, N);
 
-h
+d = size(X,1);
 
 h_star = sqrt(2) *h;
 two_h_squared =      2 * (h      ^ 2);
 two_h_star_squared = 2 * (h_star ^ 2);
-h_sqrt_2_pi =       h      * sqrt(2*pi);
-h_star_sqrt_2_pi =  h_star * sqrt(2*pi);
+h_sqrt_2_pi =       (h^d)      * sqrt(2*pi);
+h_star_sqrt_2_pi =  sqrt(2) * (h^d) * sqrt(2*pi);
 
 sum = 0;
 for i=1:N
+  v = X(:,i);
+%  disp(i);
   for j=1:N
-    sum = sum + exp(-((X(j) - X(i))^2) / two_h_star_squared);
+    dist_sq = norm(X(:,j) - v)^2;
+    sum = sum + exp(-dist_sq / two_h_star_squared);
   end
 end
 
@@ -25,7 +28,8 @@ int_f_hat_squared = sum / (h_star_sqrt_2_pi * (N^2));
 sum = 0;
 for i=1:N
   for j=1:N
-    sum = sum + exp(-((X(i) - X(j))^2) / two_h_squared);
+    dist_sq = norm(X(j,:) - X(i,:))^2;
+    sum = sum + exp(-dist_sq / two_h_squared);
   end
 end
 
diff --git a/fastlib/u/niche/entropy/kee/kde_risk_cv.m b/fastlib/u/niche/entropy/kee/kde_risk_cv.m
index a46cfdf601..f1283b9cd9 100644
--- a/fastlib/u/niche/entropy/kee/kde_risk_cv.m
+++ b/fastlib/u/niche/entropy/kee/kde_risk_cv.m
@@ -1,8 +1,10 @@
 % kde_risk_cv() - cross validate to select the optimal bandwidth
 % according to risk minimization
-function [risk] = kde_risk_cv(N, h_array);
+function [risk] = kde_risk_cv(x, h_array);
 
-x = normrnd(zeros(N,1), 1);
+%x = load('/home/niche/scaled_refined_astroset_20k.ds');
+
+%x = normrnd(zeros(N,1), 1);
 
 num_h = length(h_array);
 
@@ -10,5 +12,5 @@ risk = zeros(num_h,1);
 
 for i=1:num_h
 %  h_array(i)
-  risk(i) = kde_risk(x, h_array(i));
+  risk(i) = kde_risk(h_array(i), x, size(x,2));
 end
diff --git a/fastlib/u/niche/entropy/tree/bi_fdkde.cc b/fastlib/u/niche/entropy/tree/bi_fdkde.cc
index 9c594bddd4..85f8ea86f5 100644
--- a/fastlib/u/niche/entropy/tree/bi_fdkde.cc
+++ b/fastlib/u/niche/entropy/tree/bi_fdkde.cc
@@ -713,44 +713,45 @@ public:
 	}
 	break;
       case 3:
-
-	double right_diff = right_bandwidth - mid_bandwidth;
-	double left_diff = mid_bandwidth - left_bandwidth;
-
-	if((right_diff < EPSILON) && (left_diff < EPSILON)) {
-	  search = false;
-	}
-	else {
-	  if(right_diff > left_diff) {
-	    choose = 1; //choose right
-	  }
-	  else if(left_diff > right_diff) {
-	    choose = 0; // choose left
+	{
+	  double right_diff = right_bandwidth - mid_bandwidth;
+	  double left_diff = mid_bandwidth - left_bandwidth;
+	  
+	  if((right_diff < EPSILON) && (left_diff < EPSILON)) {
+	    search = false;
 	  }
 	  else {
-	    if(left_score < right_score) {
+	    if(right_diff > left_diff) {
+	      choose = 1; //choose right
+	    }
+	    else if(left_diff > right_diff) {
 	      choose = 0; // choose left
 	    }
 	    else {
-	      choose = 1; // choose right
+	      if(left_score < right_score) {
+		choose = 0; // choose left
+	      }
+	      else {
+		choose = 1; // choose right
+	      }
 	    }
+	    
+	    if(choose == 0) {
+	      param->bandwidth = (left_bandwidth + mid_bandwidth) / 2;
+	      printf("left = %f, mid = %f\n", left_bandwidth, mid_bandwidth);
+	      printf("splitting on left and mid: ");
+	    }
+	    else {
+	      param->bandwidth = (right_bandwidth + mid_bandwidth) / 2;
+	      printf("right = %f, mid = %f\n", right_bandwidth, mid_bandwidth);
+	      printf("splitting on right and mid: ");
+	    }
+	    printf("%f\n", param->bandwidth);
 	  }
-
-	  if(choose == 0) {
-	    param->bandwidth = (left_bandwidth + mid_bandwidth) / 2;
-	    printf("left = %f, mid = %f\n", left_bandwidth, mid_bandwidth);
-	    printf("splitting on left and mid: ");
-	  }
-	  else {
-	    param->bandwidth = (right_bandwidth + mid_bandwidth) / 2;
-	    printf("right = %f, mid = %f\n", right_bandwidth, mid_bandwidth);
-	    printf("splitting on right and mid: ");
-	  }
-	  printf("%f\n", param->bandwidth);
 	}
 	break;
       default:
-	;//impossible!
+	printf("impossible\n"); //impossible!
       }
 
 
@@ -1055,40 +1056,41 @@ public:
 	}
 	break;
       case 3:
-
-	double right_diff = right_bandwidth - mid_bandwidth;
-	double left_diff = mid_bandwidth - left_bandwidth;
-
-	if((right_diff < EPSILON) && (left_diff < EPSILON)) {
-	  search = false;
-	}
-	else {
-	  if(right_diff > left_diff) {
-	    choose = 1; //choose right
-	  }
-	  else if(left_diff > right_diff) {
-	    choose = 0; // choose left
+	{
+	  double right_diff = right_bandwidth - mid_bandwidth;
+	  double left_diff = mid_bandwidth - left_bandwidth;
+	  
+	  if((right_diff < EPSILON) && (left_diff < EPSILON)) {
+	    search = false;
 	  }
 	  else {
-	    if(left_score < right_score) {
+	    if(right_diff > left_diff) {
+	      choose = 1; //choose right
+	    }
+	    else if(left_diff > right_diff) {
 	      choose = 0; // choose left
 	    }
 	    else {
-	      choose = 1; // choose right
+	      if(left_score < right_score) {
+		choose = 0; // choose left
+	      }
+	      else {
+		choose = 1; // choose right
+	      }
 	    }
+	    
+	    if(choose == 0) {
+	      param->bandwidth = (left_bandwidth + mid_bandwidth) / 2;
+	      printf("left = %f, mid = %f\n", left_bandwidth, mid_bandwidth);
+	      printf("splitting on left and mid: ");
+	    }
+	    else {
+	      param->bandwidth = (right_bandwidth + mid_bandwidth) / 2;
+	      printf("right = %f, mid = %f\n", right_bandwidth, mid_bandwidth);
+	      printf("splitting on right and mid: ");
+	    }
+	    printf("%f\n", param->bandwidth);
 	  }
-
-	  if(choose == 0) {
-	    param->bandwidth = (left_bandwidth + mid_bandwidth) / 2;
-	    printf("left = %f, mid = %f\n", left_bandwidth, mid_bandwidth);
-	    printf("splitting on left and mid: ");
-	  }
-	  else {
-	    param->bandwidth = (right_bandwidth + mid_bandwidth) / 2;
-	    printf("right = %f, mid = %f\n", right_bandwidth, mid_bandwidth);
-	    printf("splitting on right and mid: ");
-	  }
-	  printf("%f\n", param->bandwidth);
 	}
 	break;
       default:
@@ -1200,16 +1202,16 @@ public:
 	}
 
 	/*
-	Matrix ref_densities_matrix;
-	ref_densities_matrix.AliasRowVector(ref_densities);
+	  Matrix ref_densities_matrix;
+	  ref_densities_matrix.AliasRowVector(ref_densities);
 
-	Matrix linspace_densities_matrix;
-	linspace_densities_matrix.AliasRowVector(linspace_densities);
+	  Matrix linspace_densities_matrix;
+	  linspace_densities_matrix.AliasRowVector(linspace_densities);
 	
-	data::Save(fx_param_str(module, "ref_densities", "ref_densities.csv"),
-		   ref_densities_matrix);
-	data::Save(fx_param_str(module, "linspace_densities", "linspace_densities.csv"),
-		   linspace_densities_matrix);
+	  data::Save(fx_param_str(module, "ref_densities", "ref_densities.csv"),
+	  ref_densities_matrix);
+	  data::Save(fx_param_str(module, "linspace_densities", "linspace_densities.csv"),
+	  linspace_densities_matrix);
 	*/
 
 
@@ -1220,9 +1222,9 @@ public:
 	  ((double)(param->r_count));
 	
 	/*
-	double second_term =
+	  double second_term =
 	  la::Dot(h_KL_log_densities, ref_densities) / 
-	 ((double)(param->r_count));
+	  ((double)(param->r_count));
 	*/
 
 
diff --git a/fastlib/u/niche/functional/p300_funcica.m b/fastlib/u/niche/functional/p300_funcica.m
index e12b94f501..6802123d98 100644
--- a/fastlib/u/niche/functional/p300_funcica.m
+++ b/fastlib/u/niche/functional/p300_funcica.m
@@ -133,7 +133,7 @@ pca_results = pca_fd(centered_myfd, p, myfdPar);
 %  ic_scores(i,:) = scale_up_factor * ic_scores(i,:);
 %end
 
-save p300_filtered_lambda3Eneg6_correct_results_120.mat;
+%save p300_filtered_lambda3Eneg6_correct_results_120.mat;
 
 
 % given a set of curves, identify component curves of variation
diff --git a/fastlib/u/niche/pica/pica.m b/fastlib/u/niche/pica/pica.m
index a34ef58331..81c87f60b4 100644
--- a/fastlib/u/niche/pica/pica.m
+++ b/fastlib/u/niche/pica/pica.m
@@ -18,6 +18,7 @@ for i = 1:D
   S(i,:) = laplacinv(rand(N, 1), mu, b);
 end
 
+length(find(S < 0))
 while 1
   neg_indices = find(S < 0);
 
@@ -30,6 +31,12 @@ while 1
   end
 end
 
+% impose unit variance on each row of S
+
+for i = 1:D
+  S(i,:) = S(i,:) / std(S(i,:));
+end
+
 
 % X = A S
 % set mixing matrix A
diff --git a/fastlib/u/niche/pica/pica.tex b/fastlib/u/niche/pica/pica.tex
index c63bafdb76..db21c4fbd7 100644
--- a/fastlib/u/niche/pica/pica.tex
+++ b/fastlib/u/niche/pica/pica.tex
@@ -1,6 +1,6 @@
 \documentclass{article}
 \usepackage{amsmath,amsthm,amsfonts}
-\title{Postitive ICA}
+\title{Positive ICA}
 
 \begin{document}
 
@@ -9,5 +9,21 @@ $X \in \mathbf{R}^N$
 
 model: $X = A S$
 
+$\tilde{X} = Z X$ for whitening matrix $Z$ (note that $\tilde{X}$ is not centered).
+
+We want to find some
+$Y = W \tilde{X} = W Z X = V X$ ($\Rightarrow X = V^{-1} Y$)
+
+such that $Y$ has minimum sum marginal entropy under constraints that $V^{-1}$ and $Y$ are positive (nonnegative)
+
+
+
+Plumbley's nonnegative ICA is globally converges onto the independent components if they are well-grounded. What if they are not well-grounded?
+
+Can we quantify the number of nonnegative matrix factorization candidate solutions? How many rotations satisfy the positivity constraint
+
+
+
+
 
 \end{document}
\ No newline at end of file