From 2e1e4bef89a4102759470a27d19ab952fb2a7148 Mon Sep 17 00:00:00 2001
From: Michael T Smith <lionfishy@gmail.com>
Date: Thu, 10 Nov 2016 18:16:24 +0000
Subject: [PATCH] Fixed CCD now

---
 GPy/core/model.py       | 17 ++++++++++++-----
 GPy/testing/ccd_test.py | 38 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 5 deletions(-)
 create mode 100644 GPy/testing/ccd_test.py

diff --git a/GPy/core/model.py b/GPy/core/model.py
index a5dd5795..3ab8094b 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -199,23 +199,30 @@ class Model(ParamzModel, Priorizable):
         import paramz
 
         transformed_points = param_points.copy()
+        print transformed_points.shape
         #alan's original code to transform those parameters, to the true space of parameters again
         #mike's change: some parameters have no transform, and thus won't be called by any of the
         #iterations through m2.constraints.items(). To handle these we keep track of those parameters
         #not included, and then add them (untransformed) at the end.
         f = np.ones(self.size).astype(bool)
         f[self.constraints[paramz.transformations.__fixed__]] = paramz.transformations.FIXED
+        
         new_t_points = [] 
-        todo = range(0,self.size)
+        todo = range(0,sum(f))
         new_t_points = np.zeros_like(transformed_points)
+        print new_t_points.shape
         for c, ind in self.constraints.items():
             print c,ind
             if c != paramz.transformations.__fixed__:
-                new_t_points[:,ind] = (c.f(transformed_points[:, ind[f[ind]]]))
-                todo.remove(ind)
+                for i in ind[f[ind]]:
+                    new_t_points[:,i] = (c.f(transformed_points[:, i]))
+                    print i
+                    todo.remove(i)
 
-        for param_ind in todo:
-            new_t_points[:,param_ind] = transformed_points[:, param_ind]
+        print todo
+        for i in todo:
+            print i
+            new_t_points[:,i] = transformed_points[:, i]
             
             
         return new_t_points, point_densities
diff --git a/GPy/testing/ccd_test.py b/GPy/testing/ccd_test.py
new file mode 100644
index 00000000..48f2574c
--- /dev/null
+++ b/GPy/testing/ccd_test.py
@@ -0,0 +1,38 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import GPy
+%matplotlib inline
+np.set_printoptions(suppress=True, precision=10)
+from paramz.transformations import Logexp
+
+class SimpleModel(GPy.Model):
+    def __init__(self, name):
+        super(SimpleModel, self).__init__(name)
+        self.params = []
+        self.peak = np.array([1,2])
+        for i,pos in enumerate(self.peak):
+            p = GPy.Param('param%d' % i, 1.0)
+            self.params.append(p)
+            self.link_parameter(p)
+    
+    def log_likelihood(self):
+        like = 0
+        for i,pos in enumerate(self.peak):
+            like -= ((self.params[i])-pos)**2
+        return like
+    
+    def parameters_changed(self):
+        for i,pos in enumerate(self.peak):
+            self.params[i].gradient = -2*((self.params[i])-pos)
+    
+m2 = SimpleModel('simple')
+m2.optimize()
+assert np.all(np.isclose(m2.numerical_parameter_hessian(),np.array([[2,0],[0,2]]))), "Numerical approximation to Hessian doesn't match. Error in numerical_parameter_hessian()."
+assert np.isclose(m2.param0,1,atol=0.01), "Failed to find likelihood maximum of test model's param0 while testing CCD"
+assert np.isclose(m2.param1,2,atol=0.01), "Failed to find likelihood maximum of test model's param1 while testing CCD"
+
+ccdpos,ccdres = m2.CCD()
+m2.optimize()
+
+assert np.all(np.isclose(np.sum((ccdpos-np.array([1,2]))**2,1)[1:],1.21,atol=0.01)), "CCD placement error - should be 1.21 from mode"
+s