diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index a6478850..2efc2403 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -80,8 +80,6 @@ class SparseGP(GPBase): self.psi2 = None def _computations(self): - - # factor Kmm self.Lm = jitchol(self.Kmm) @@ -93,6 +91,8 @@ class SparseGP(GPBase): psi2_beta = self.psi2.sum(0) * self.likelihood.precision evals, evecs = linalg.eigh(psi2_beta) clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable + if not np.array_equal(evals, clipped_evals): + pass#print evals tmp = evecs * np.sqrt(clipped_evals) tmp = tmp.T else: @@ -114,7 +114,8 @@ class SparseGP(GPBase): # back substutue C into psi1Vf tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1Vf), lower=1, trans=0) self._LBi_Lmi_psi1Vf, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0) - tmp, info2 = dpotrs(self.LB, tmp, lower=1) + #tmp, info2 = dpotrs(self.LB, tmp, lower=1) + tmp, info2 = dtrtrs(self.LB, self._LBi_Lmi_psi1Vf, lower=1, trans=1) self.Cpsi1Vf, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1) # Compute dL_dKmm diff --git a/GPy/core/svigp.py b/GPy/core/svigp.py index 7d2bec36..8a2d889d 100644 --- a/GPy/core/svigp.py +++ b/GPy/core/svigp.py @@ -92,12 +92,56 @@ class SVIGP(GPBase): self._vb_steplength_trace = [] def getstate(self): - return GPBase.getstate(self) - + steplength_params = [self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength] + return GPBase.getstate(self) + \ + [self.get_vb_param(), + self.Z, + self.num_inducing, + self.has_uncertain_inputs, + self.X_variance, + self.X_batch, + self.X_variance_batch, + steplength_params, + self.batchcounter, + self.batchsize, + self.epochs, + self.momentum, + self.data_prop, + self._param_trace, + self._param_steplength_trace, + self._vb_steplength_trace, + self._ll_trace, + self._grad_trace, + self.Y, + self._permutation, + self.iterations + ] def setstate(self, state): - return GPBase.setstate(self, state) - + self.iterations = state.pop() + self._permutation = state.pop() + self.Y = state.pop() + self._grad_trace = state.pop() + self._ll_trace = state.pop() + self._vb_steplength_trace = state.pop() + self._param_steplength_trace = state.pop() + self._param_trace = state.pop() + self.data_prop = state.pop() + self.momentum = state.pop() + self.epochs = state.pop() + self.batchsize = state.pop() + self.batchcounter = state.pop() + steplength_params = state.pop() + (self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength) = steplength_params + self.X_variance_batch = state.pop() + self.X_batch = state.pop() + self.X_variance = state.pop() + self.has_uncertain_inputs = state.pop() + self.num_inducing = state.pop() + self.Z = state.pop() + vb_param = state.pop() + GPBase.setstate(self, state) + self.set_vb_param(vb_param) def _compute_kernel_matrices(self): # kernel computations, using BGPLVM notation diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py index d106cd4f..60e78461 100644 --- a/GPy/models/gplvm.py +++ b/GPy/models/gplvm.py @@ -63,6 +63,20 @@ class GPLVM(GP): return np.hstack((dL_dX.flatten(), GP._log_likelihood_gradients(self))) + def jacobian(self,X): + target = np.zeros((X.shape[0],X.shape[1],self.output_dim)) + for i in range(self.output_dim): + target[:,:,i]=self.kern.dK_dX(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X) + return target + + def magnification(self,X): + target=np.zeros(X.shape[0]) + J = np.zeros((X.shape[0],X.shape[1],self.output_dim)) + J=self.jacobian(X) + for i in range(X.shape[0]): + target[i]=np.sqrt(pb.det(np.dot(J[i,:,:],np.transpose(J[i,:,:])))) + return target + def plot(self): assert self.likelihood.Y.shape[1] == 2 pb.scatter(self.likelihood.Y[:, 0], self.likelihood.Y[:, 1], 40, self.X[:, 0].copy(), linewidth=0, cmap=pb.cm.jet) @@ -72,3 +86,6 @@ class GPLVM(GP): def plot_latent(self, *args, **kwargs): return util.plot_latent.plot_latent(self, *args, **kwargs) + + def plot_magnification(self, *args, **kwargs): + return util.plot_latent.plot_magnification(self, *args, **kwargs) diff --git a/GPy/util/datasets/connections.txt b/GPy/util/datasets/connections.txt new file mode 100644 index 00000000..e1886100 --- /dev/null +++ b/GPy/util/datasets/connections.txt @@ -0,0 +1,22 @@ +LFHD, RFHD +RFHD, RBHD +RBHD, LBHD +LBHD, LFHD +LELB, LWRB +LWRB, LFIN +LELB, LSHO +LSHO, RSHO +RSHO, STRN +LSHO, STRN +RSHO, RELB +RELB, RWRB +RWRB, RFIN +LSHO, LFWT +RSHO, RFWT +LFWT, RFWT +LFWT, LKNE +RFWT, RKNE +LKNE, LHEE +RKNE, RHEE +RMT5, RHEE +LMT5, LHEE diff --git a/GPy/util/plot_latent.py b/GPy/util/plot_latent.py index a5f81be4..dd21b2ea 100644 --- a/GPy/util/plot_latent.py +++ b/GPy/util/plot_latent.py @@ -98,3 +98,79 @@ def plot_latent(model, labels=None, which_indices=None, ax.figure.canvas.show() raw_input('Enter to continue') return ax + +def plot_magnification(model, labels=None, which_indices=None, + resolution=60, ax=None, marker='o', s=40, + fignum=None, plot_inducing=False, legend=True, + aspect='auto', updates=False): + """ + :param labels: a np.array of size model.num_data containing labels for the points (can be number, strings, etc) + :param resolution: the resolution of the grid on which to evaluate the predictive variance + """ + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) + util.plot.Tango.reset() + + if labels is None: + labels = np.ones(model.num_data) + + input_1, input_2 = most_significant_input_dimensions(model, which_indices) + + # first, plot the output variance as a function of the latent space + Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution) + Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) + def plot_function(x): + Xtest_full[:, [input_1, input_2]] = x + mf=model.magnification(Xtest_full) + return mf + view = ImshowController(ax, plot_function, tuple(xmin) + tuple(xmax), + resolution, aspect=aspect, interpolation='bilinear', + cmap=pb.cm.gray) + + # make sure labels are in order of input: + ulabels = [] + for lab in labels: + if not lab in ulabels: + ulabels.append(lab) + + marker = itertools.cycle(list(marker)) + + for i, ul in enumerate(ulabels): + if type(ul) is np.string_: + this_label = ul + elif type(ul) is np.int64: + this_label = 'class %i' % ul + else: + this_label = 'class %i' % i + m = marker.next() + + index = np.nonzero(labels == ul)[0] + if model.input_dim == 1: + x = model.X[index, input_1] + y = np.zeros(index.size) + else: + x = model.X[index, input_1] + y = model.X[index, input_2] + ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) + + ax.set_xlabel('latent dimension %i' % input_1) + ax.set_ylabel('latent dimension %i' % input_2) + + if not np.all(labels == 1.) and legend: + ax.legend(loc=0, numpoints=1) + + ax.set_xlim(xmin[0], xmax[0]) + ax.set_ylim(xmin[1], xmax[1]) + ax.grid(b=False) # remove the grid if present, it doesn't look good + ax.set_aspect('auto') # set a nice aspect ratio + + if plot_inducing: + ax.plot(model.Z[:, input_1], model.Z[:, input_2], '^w') + + if updates: + ax.figure.canvas.show() + raw_input('Enter to continue') + + pb.title('Magnification Factor') + return ax