derivative of the covariance matrix with respect to the parameters.
329850
def update_gradients_full(self, dL_dK, X, X2=None):
                """"""
                X,slices = X[:,:-1],index_to_slices(X[:,-1])
                if X2 is None:
                    X2,slices2 = X,slices
                    K = np.zeros((X.shape[0], X.shape[0]))
                else:
                    X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])


                vyt = self.variance_Yt

                lyt = 1./(2*self.lengthscale_Yt)

                tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
                ttdist = (X[:,0][:,None] - X2[:,0][None,:])
                #rdist = [tdist,xdist]
                
                rd=tdist.shape[0]

                dka = np.zeros([rd,rd])
                dkc = np.zeros([rd,rd])
                dkYdvart = np.zeros([rd,rd])
                dkYdlent = np.zeros([rd,rd])

                dkdubias = np.zeros([rd,rd])

                kyy = lambda tdist: np.exp(-lyt*(tdist))
                dkyydlyt = lambda tdist: kyy(tdist)*(-tdist)

                k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )

                k4 = lambda ttdist: 2*lyt*(ttdist)

                dk1dlyt = lambda tdist: 2. - 4*2.*lyt*tdist

                dk4dlyt = lambda ttdist: 2*(ttdist)

                for i, s1 in enumerate(slices):
                    for j, s2 in enumerate(slices2):
                        for ss1 in s1:
                            for ss2 in s2:
                                if i==0 and j==0:
                                    dkYdvart[ss1,ss2] = kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 0
                                elif i==0 and j==1:
                                    dkYdvart[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    #dkYdvart[ss1,ss2] = ((2*lyt*ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (k4(ttdist[ss1,ss2])+1.)+\
                                    vyt*kyy(tdist[ss1,ss2])*(dk4dlyt(ttdist[ss1,ss2]))
                                    #dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (2*lyt*(ttdist[ss1,ss2])+1.)+\
                                    #vyt*kyy(tdist[ss1,ss2])*(2*ttdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 0
                                elif i==1 and j==1:
                                    dkYdvart[ss1,ss2] = (k1(tdist[ss1,ss2]) + 1. )* kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])*( k1(tdist[ss1,ss2]) + 1. ) +\
                          			vyt*kyy(tdist[ss1,ss2])*dk1dlyt(tdist[ss1,ss2])
                                    dkdubias[ss1,ss2] = 1
                                else:
                                    dkYdvart[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    #dkYdvart[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
                                    dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-k4(ttdist[ss1,ss2])+1.)+\
                                    vyt*kyy(tdist[ss1,ss2])*(-dk4dlyt(ttdist[ss1,ss2]) )
                                    dkdubias[ss1,ss2] = 0
                                    #dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-2*lyt*(ttdist[ss1,ss2])+1.)+\
                                    #vyt*kyy(tdist[ss1,ss2])*(-2)*(ttdist[ss1,ss2])
   

                self.variance_Yt.gradient = np.sum(dkYdvart * dL_dK)

                self.lengthscale_Yt.gradient =  np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK)

                self.ubias.gradient = np.sum(dkdubias * dL_dK)
1179869
def setCovariance(self,cov):
        
        self.setParams(sp.log(sp.diagonal(cov)))
1195563
def fgrad_y_psi(self, y, return_covar_chain=False):
        
        mpsi = self.psi

        w, s, r, d = self.fgrad_y(y, return_precalc=True)
        gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
        for i in range(len(mpsi)):
            a,b,c  = mpsi[i]
            gradients[:, :, i, 0] = (b * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 1] = a * (d[i] - 2.0 * s[i] * r[i] * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 2] = (-2.0 * a * (b ** 2) * r[i] * ((1.0 / np.cosh(s[i])) ** 2)).T
        gradients[:, :, 0, 3] = 1.0

        if return_covar_chain:
            covar_grad_chain = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
            for i in range(len(mpsi)):
                a,b,c = mpsi[i]
                covar_grad_chain[:, :, i, 0] = (r[i]).T
                covar_grad_chain[:, :, i, 1] = (a * (y + c) * ((1.0 / np.cosh(s[i])) ** 2).T)
                covar_grad_chain[:, :, i, 2] = a * b * ((1.0 / np.cosh(s[i])) ** 2).T
            covar_grad_chain[:, :, 0, 3] = y
            return gradients, covar_grad_chain

        return gradients
1138276
def derivatives(self, x, y, sigma0, Rs, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        e = abs(1 - q)
        x_ = (cos_phi*x_shift+sin_phi*y_shift)*np.sqrt(1 - e)
        y_ = (-sin_phi*x_shift+cos_phi*y_shift)*np.sqrt(1 + e)

        f_x_prim, f_y_prim = self.spherical.derivatives(x_, y_, sigma0, Rs)
        f_x_prim *= np.sqrt(1 - e)
        f_y_prim *= np.sqrt(1 + e)
        f_x = cos_phi*f_x_prim-sin_phi*f_y_prim
        f_y = sin_phi*f_x_prim+cos_phi*f_y_prim
        return f_x, f_y
1364948
def hessian(self, x, y, theta_E, s, q):
        
        alpha_ra, alpha_dec = self.derivatives(x, y, theta_E, s, q)
        diff = self._diff
        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, theta_E, s, q)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, theta_E, s, q)

        f_xx = (alpha_ra_dx - alpha_ra) / diff
        f_xy = (alpha_ra_dy - alpha_ra) / diff
        # f_yx = (alpha_dec_dx - alpha_dec)/diff
        f_yy = (alpha_dec_dy - alpha_dec) / diff
        return f_xx, f_yy, f_xy
956087
def gradient_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        _, h_init, _ = self.info_init_params
        _, _, _, h1, h2, _ = self.info_dynamics_params
        H_diag, H_upper_diag = self.sparse_J_prior

        # Compute the gradient from the prior
        g = -1 * symm_block_tridiag_matmul(H_diag, H_upper_diag, x)
        g[0] += h_init
        g[:-1] += h1
        g[1:] += h2

        # Compute gradient from the likelihood terms
        g += self.grad_local_log_likelihood(x)

        return g
729968
def _mdens_deriv(self,m):
        """"""
        return -2.*self.a2**2*(self.a+2.*m)/m**3/(self.a+m)**3
1208542
def derivativeW(self,*args):
        '''
        
        '''
        if self.n_dims >= 4:
            j = 0
        else:
            assert False, "Derivative with respect to W can't be called when n_dims < 4!"
        if self.i_dim == j:
            return np.ones_like(*args[0])
        else:
            return np.zeros_like(*args[0])
411808
def obfn_dfd(self):
        
        """

        return 0.5*np.linalg.norm((self.D.dot(self.obfn_fvar()) - self.S))**2
414411
def _partial_diff(var, *params):
    
    if isinstance(var, sympy.Derivative):
        return sympy.Derivative(var.expr, *(var.variables + params))
    else:
        return D(var, *params)
1089007
def _derY(self,x,y):
        '''
        
        '''
        x_pos, y_pos = self.findSector(x,y)
        alpha, beta = self.findCoords(x,y,x_pos,y_pos)

        # Get four corners data for each point
        xA = self.x_values[x_pos,y_pos]
        xB = self.x_values[x_pos+1,y_pos]
        xC = self.x_values[x_pos,y_pos+1]
        xD = self.x_values[x_pos+1,y_pos+1]
        yA = self.y_values[x_pos,y_pos]
        yB = self.y_values[x_pos+1,y_pos]
        yC = self.y_values[x_pos,y_pos+1]
        yD = self.y_values[x_pos+1,y_pos+1]
        fA = self.f_values[x_pos,y_pos]
        fB = self.f_values[x_pos+1,y_pos]
        fC = self.f_values[x_pos,y_pos+1]
        fD = self.f_values[x_pos+1,y_pos+1]

        # Calculate components of the alpha,beta --> x,y delta translation matrix
        alpha_x = (1-beta)*(xB-xA) + beta*(xD-xC)
        alpha_y = (1-beta)*(yB-yA) + beta*(yD-yC)
        beta_x  = (1-alpha)*(xC-xA) + alpha*(xD-xB)
        beta_y  = (1-alpha)*(yC-yA) + alpha*(yD-yB)

        # Invert the delta translation matrix into x,y --> alpha,beta
        det = alpha_x*beta_y - beta_x*alpha_y
        y_alpha = -beta_x/det
        y_beta  = alpha_x/det

        # Calculate the derivative of f w.r.t. alpha and beta
        dfda = (1-beta)*(fB-fA) + beta*(fD-fC)
        dfdb = (1-alpha)*(fC-fA) + alpha*(fD-fB)

        # Calculate the derivative with respect to x (and return it)
        dfdy = y_alpha*dfda + y_beta*dfdb
        return dfdy
341650
def _maximization(self, X):
        """"""

        # Iterate through clusters and recalculate mean and covariance
        for i in range(self.k):
            resp = np.expand_dims(self.responsibility[:, i], axis=1)
            mean = (resp * X).sum(axis=0) / resp.sum()
            covariance = (X - mean).T.dot((X - mean) * resp) / resp.sum()
            self.parameters[i]["mean"], self.parameters[i]["cov"] = mean, covariance

        # Update weights
        n_samples = np.shape(X)[0]
        self.priors = self.responsibility.sum(axis=0) / n_samples
515973
def _param_grad_helper(self, dL_dK, X, X2, target):
        """"""
        self._K_computations(X, X2)
        base = self.variance*self.degree*self._K_poly_arg**(self.degree-1)
        base_cov_grad = base*dL_dK


            
        target[0] += np.sum(self._K_dvar*dL_dK)
        target[1] += (self._K_inner_prod*base_cov_grad).sum()
        target[2] += base_cov_grad.sum()
642700
def derivatives(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):
        

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        if Rs < 0.0000001:
            Rs = 0.0000001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        dx, dy = self.coreBurkAlpha(R, Rs, rho0, r_core, x_, y_)

        return dx, dy
218427
def parameters_changed(self):
        

        #np.set_printoptions(16)
        #print(self.param_array)
        

        # Get the model matrices from the kernel
        (F,L,Qc,H,P_inf, P0, dFt,dQct,dP_inft, dP0t) = self.kern.sde()
        
        
        # necessary parameters
        measurement_dim = self.output_dim
        grad_params_no = dFt.shape[2]+1 # we also add measurement noise as a parameter

        # add measurement noise as a parameter and get the gradient matrices
        dF    = np.zeros([dFt.shape[0],dFt.shape[1],grad_params_no])
        dQc   = np.zeros([dQct.shape[0],dQct.shape[1],grad_params_no])
        dP_inf = np.zeros([dP_inft.shape[0],dP_inft.shape[1],grad_params_no])
        dP0 = np.zeros([dP0t.shape[0],dP0t.shape[1],grad_params_no])

        # Assign the values for the kernel function
        dF[:,:,:-1] = dFt
        dQc[:,:,:-1] = dQct
        dP_inf[:,:,:-1] = dP_inft
        dP0[:,:,:-1] = dP0t

        # The sigma2 derivative
        dR = np.zeros([measurement_dim,measurement_dim,grad_params_no])
        dR[:,:,-1] = np.eye(measurement_dim)

        # Balancing
        if self.balance:
            (F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf,dP0) = ssm.balance_ss_model(F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf, dP0)
            print("SSM parameters_changed balancing!")
        # Use the Kalman filter to evaluate the likelihood
        grad_calc_params = {}
        grad_calc_params['dP_inf'] = dP_inf
        grad_calc_params['dF'] = dF
        grad_calc_params['dQc'] = dQc
        grad_calc_params['dR'] = dR
        grad_calc_params['dP_init'] = dP0

        kalman_filter_type = self.kalman_filter_type

        # The following code is required because sometimes the shapes of self.Y
        # becomes 3D even though is must be 2D. The reason is undiscovered.
        Y = self.Y
        if self.ts_number is None:
            Y.shape = (self.num_data,1)
        else:
            Y.shape = (self.num_data,1,self.ts_number)

        (filter_means, filter_covs, log_likelihood,
         grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(F,L,Qc,H,
                                      float(self.Gaussian_noise.variance),P_inf,self.X,Y,m_init=None,
                                      P_init=P0, p_kalman_filter_type = kalman_filter_type, calc_log_likelihood=True,
                                      calc_grad_log_likelihood=True,
                                      grad_params_no=grad_params_no,
                                      grad_calc_params=grad_calc_params)

        if np.any( np.isfinite(log_likelihood) == False):
            #import pdb; pdb.set_trace()
            print("State-Space: NaN valkues in the log_likelihood")

        if np.any( np.isfinite(grad_log_likelihood) == False):
            #import pdb; pdb.set_trace()
            print("State-Space: NaN values in the grad_log_likelihood")
        #print(grad_log_likelihood)

        grad_log_likelihood_sum = np.sum(grad_log_likelihood,axis=1)
        grad_log_likelihood_sum.shape = (grad_log_likelihood_sum.shape[0],1)
        self._log_marginal_likelihood = np.sum( log_likelihood,axis=1 )
        self.likelihood.update_gradients(grad_log_likelihood_sum[-1,0])

        self.kern.sde_update_gradient_full(grad_log_likelihood_sum[:-1,0])
68859
def conditional(mean, covar, dims_in, dims_out, covariance_type='full'):
    
    in_in = covar[ix_(dims_in, dims_in)]
    in_out = covar[ix_(dims_in, dims_out)]
    out_in = covar[ix_(dims_out, dims_in)]
    out_out = covar[ix_(dims_out, dims_out)]
    in_in_inv = inv(in_in)
    out_in_dot_in_in_inv = out_in.dot(in_in_inv)

    cond_covar = out_out - out_in_dot_in_in_inv.dot(in_out)
    cond_mean = lambda x: mean[dims_out] + out_in_dot_in_in_inv.dot(x - mean[dims_in])
    return lambda x: [cond_mean(x), cond_covar]
865166
def deriv(self, t, y, data=None):
        " "
        theta, dtheta_dt = y
        return np.array([dtheta_dt, - self.g_l * gv.sin(theta)])