Construct the coupling matrix (and adjacancy matrix) from predefined models
            or via sampling.
14207
def set_coupl(self, Coupl=None):
        
        self.varNames = collections.OrderedDict([(str(i), i) for i in range(self.dim)])
        if (self.model not in self.availModels.keys()
            and Coupl is None):
            self.read_model()
        elif 'var' in self.model:
            # vector auto regressive process
            self.Coupl = Coupl
            self.boolRules = collections.OrderedDict(
                              [(s, '') for s in self.varNames.keys()])
            names = list(self.varNames.keys())
            for gp in range(self.dim):
                pas = []
                for g in range(self.dim):
                    if np.abs(self.Coupl[gp,g] > 1e-10):
                        pas.append(names[g])
                self.boolRules[
                    names[gp]] = ''.join(pas[:1]
                                         + [' or ' + pa for pa in pas[1:]])
                self.Adj_signed = np.sign(Coupl)
        elif self.model in ['6','7','8','9','10']:
            self.Adj_signed = np.zeros((self.dim,self.dim))
            n_sinknodes = 2
#             sinknodes = np.random.choice(np.arange(0,self.dim),
#                                              size=n_sinknodes,replace=False)
            sinknodes = np.array([0,1])
            # assume sinknodes have feeback
            self.Adj_signed[sinknodes,sinknodes] = np.ones(n_sinknodes)
#             # allow negative feedback
#             if self.model == 10:
#                 plus_minus = (np.random.randint(0,2,n_sinknodes) - 0.5)*2
#                 self.Adj_signed[sinknodes,sinknodes] = plus_minus
            leafnodes = np.array(sinknodes)
            availnodes = np.array([i for i in range(self.dim) if i not in sinknodes])
#             settings.m(0,leafnodes,availnodes)
            while len(availnodes) != 0:
                # parent
                parent_idx = np.random.choice(np.arange(0,len(leafnodes)),
                                                    size=1,replace=False)
                parent = leafnodes[parent_idx]
                # children
                children_ids = np.random.choice(np.arange(0,len(availnodes)),
                                                       size=2,replace=False)
                children = availnodes[children_ids]
                settings.m(0,parent,children)
                self.Adj_signed[children,parent] = np.ones(2)
                if self.model == 8:
                    self.Adj_signed[children[0],children[1]] = -1
                if self.model in [9,10]:
                    self.Adj_signed[children[0],children[1]] = -1
                    self.Adj_signed[children[1],children[0]] = -1
                # update leafnodes
                leafnodes = np.delete(leafnodes,parent_idx)
                leafnodes = np.append(leafnodes,children)
                # update availnodes
                availnodes = np.delete(availnodes,children_ids)
#                 settings.m(0,availnodes)
#                 settings.m(0,leafnodes)
#                 settings.m(0,self.Adj)
#                 settings.m(0,'-')
        else:
            self.Adj = np.zeros((self.dim,self.dim))
            for i in range(self.dim):
                indep = np.random.binomial(1,self.p_indep)
                if indep == 0:
                    # this number includes parents (other variables)
                    # and the variable itself, therefore its
                    # self.maxnpar+2 in the following line
                    nr = np.random.randint(1,self.maxnpar+2)
                    j_par = np.random.choice(np.arange(0,self.dim),
                                             size=nr,replace=False)
                    self.Adj[i,j_par] = 1
                else:
                    self.Adj[i,i] = 1
        #
        self.Adj = np.abs(np.array(self.Adj_signed))
531785
def _sample_mvn(n_samples, cov, prng):
    
    n_features, _ = cov.shape
    return prng.multivariate_normal(np.zeros(n_features), cov, size=n_samples)
1159939
def generate_synthObs(self, bases_wave, bases_flux, basesCoeff, Av_star, z_star, sigma_star, resample_range = None, resample_int = 1):
        
        ''''''
        nbases              = basesCoeff.shape[0]        
        bases_wave_resam    = arange(int(resample_range[0]), int(resample_range[-1]), resample_int, dtype=float)
        npix_resample       = len(bases_wave_resam)
        
        #Resampling the range
        bases_flux_resam = empty((nbases, npix_resample))
        for i in range(nbases):
#             print bases_wave[i][0], bases_wave[i][-1]
#             print bases_wave_resam[0], bases_wave_resam[-1]
            bases_flux_resam[i,:] = interp1d(bases_wave[i], bases_flux[i], bounds_error=True)(bases_wave_resam)            
        
        #Display physical parameters
        synth_wave                  = bases_wave_resam * (1 + z_star)
        
        Av_vector                   = Av_star * ones(nbases)
        Xx_redd                     = CCM89_Bal07(3.4, bases_wave_resam)   
        r_sigma                     = sigma_star/(synth_wave[1] - synth_wave[0])

        #Defining empty kernel
        box                         = int(3 * r_sigma) if int(3 * r_sigma) < 3 else 3
        kernel_len                  = 2 * box + 1
        kernel                      = zeros((1, kernel_len)) 
        kernel_range                = arange(0, 2 * box + 1)
        
        #Generating the kernel with sigma (the norm factor is the sum of the gaussian)
        kernel[0,:]                 = exp(-0.5 * ((square(kernel_range-box)/r_sigma)))
        norm                        = np_sum(kernel[0,:])        
        kernel                      = kernel / norm

        #Convove bases with respect to kernel for dispersion velocity calculation
        bases_grid_convolve         = convolve2d(bases_flux_resam, kernel, mode='same', boundary='symm')  
        
        #Interpolate bases to wavelength range
        interBases_matrix           = (interp1d(bases_wave_resam, bases_grid_convolve, axis=1, bounds_error=True)(bases_wave_resam)).T       

        #Generate final flux model including dust        
        dust_attenuation            = power(10, -0.4 * outer(Xx_redd, Av_vector))
        bases_grid_model            = interBases_matrix * dust_attenuation
                    
        #Generate combined flux
        synth_flux                  = np_sum(basesCoeff.T * bases_grid_model, axis=1)
                                      
        return synth_wave, synth_flux
398767
def _get_kriging_matrix(self, n, n_withdrifts):
        """"""

        xy = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
                             self.Y_ADJUSTED[:, np.newaxis]), axis=1)
        d = cdist(xy, xy, 'euclidean')
        if self.UNBIAS:
            a = np.zeros((n_withdrifts+1, n_withdrifts+1))
        else:
            a = np.zeros((n_withdrifts, n_withdrifts))
        a[:n, :n] = - self.variogram_function(self.variogram_model_parameters, d)
        np.fill_diagonal(a, 0.)

        i = n
        if self.regional_linear_drift:
            a[:n, i] = self.X_ADJUSTED
            a[i, :n] = self.X_ADJUSTED
            i += 1
            a[:n, i] = self.Y_ADJUSTED
            a[i, :n] = self.Y_ADJUSTED
            i += 1
        if self.point_log_drift:
            for well_no in range(self.point_log_array.shape[0]):
                log_dist = np.log(np.sqrt((self.X_ADJUSTED - self.point_log_array[well_no, 0])**2 +
                                          (self.Y_ADJUSTED - self.point_log_array[well_no, 1])**2))
                if np.any(np.isinf(log_dist)):
                    log_dist[np.isinf(log_dist)] = -100.0
                a[:n, i] = - self.point_log_array[well_no, 2] * log_dist
                a[i, :n] = - self.point_log_array[well_no, 2] * log_dist
                i += 1
        if self.external_Z_drift:
            a[:n, i] = self.z_scalars
            a[i, :n] = self.z_scalars
            i += 1
        if self.specified_drift:
            for arr in self.specified_drift_data_arrays:
                a[:n, i] = arr
                a[i, :n] = arr
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                a[:n, i] = func(self.X_ADJUSTED, self.Y_ADJUSTED)
                a[i, :n] = func(self.X_ADJUSTED, self.Y_ADJUSTED)
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in creating kriging matrix. Kriging may fail.",
                          RuntimeWarning)
        if self.UNBIAS:
            a[n_withdrifts, :n] = 1.0
            a[:n, n_withdrifts] = 1.0
            a[n:n_withdrifts + 1, n:n_withdrifts + 1] = 0.0

        return a
1150301
void constructMatrix(DMatrixRMaj L) {
		L.reshape(cameras.size*eqs,10);

		// Known aspect ratio constraint makes more sense as a square
		double RR = this.aspectRatio *this.aspectRatio;

		// F = P*Q*P'
		// F is an arbitrary variable name and not fundamental matrix
		int index = 0;
		for (int i = 0; i < cameras.size; i++) {
			Projective P = cameras.get(i);
			DMatrix3x3 A = P.A;
			DMatrix3 B = P.a;

			// hard to tell if this helped or hurt. Keeping commented out for future investigation on proper scaling
//			double scale = NormOps_DDF3.normF(P.A);
//			CommonOps_DDF3.divide(P.A,scale);
//			CommonOps_DDF3.divide(P.a,scale);

			// constraint for a zero principle point

			// row for F[0,2] == 0
			L.data[index++] = A.a11*A.a31;
			L.data[index++] = (A.a12*A.a31 + A.a11*A.a32);
			L.data[index++] = (A.a13*A.a31 + A.a11*A.a33);
			L.data[index++] = (B.a1*A.a31 + A.a11*B.a3);
			L.data[index++] = A.a12*A.a32;
			L.data[index++] = (A.a13*A.a32 + A.a12*A.a33);
			L.data[index++] = (B.a1*A.a32 + A.a12*B.a3);
			L.data[index++] = A.a13*A.a33;
			L.data[index++] = (B.a1*A.a33 + A.a13*B.a3);
			L.data[index++] = B.a1*B.a3;

			// row for F[1,2] == 0
			L.data[index++] = A.a21*A.a31;
			L.data[index++] = (A.a22*A.a31 + A.a21*A.a32);
			L.data[index++] = (A.a23*A.a31 + A.a21*A.a33);
			L.data[index++] = (B.a2*A.a31 + A.a21*B.a3);
			L.data[index++] = A.a22*A.a32;
			L.data[index++] = (A.a23*A.a32 + A.a22*A.a33);
			L.data[index++] = (B.a2*A.a32 + A.a22*B.a3);
			L.data[index++] = A.a23*A.a33;
			L.data[index++] = (B.a2*A.a33 + A.a23*B.a3);
			L.data[index++] = B.a2*B.a3;

			if( zeroSkew ) {
				// row for F[0,1] == 0
				L.data[index++] = A.a11*A.a21;
				L.data[index++] = (A.a12*A.a21 + A.a11*A.a22);
				L.data[index++] = (A.a13*A.a21 + A.a11*A.a23);
				L.data[index++] = (B.a1*A.a21 + A.a11*B.a2);
				L.data[index++] = A.a12*A.a22;
				L.data[index++] = (A.a13*A.a22 + A.a12*A.a23);
				L.data[index++] = (B.a1*A.a22 + A.a12*B.a2);
				L.data[index++] = A.a13*A.a23;
				L.data[index++] = (B.a1*A.a23 + A.a13*B.a2);
				L.data[index++] = B.a1*B.a2;
			}

			if( knownAspect ) {
				// aspect^2*F[0,0]-F[1,1]
				L.data[index++] = A.a11*A.a11*RR - A.a21*A.a21;
				L.data[index++] = 2*(A.a11*A.a12*RR - A.a21*A.a22);
				L.data[index++] = 2*(A.a11*A.a13*RR - A.a21*A.a23);
				L.data[index++] = 2*(A.a11*B.a1*RR - A.a21*B.a2);
				L.data[index++] = A.a12*A.a12*RR - A.a22*A.a22;
				L.data[index++] = 2*(A.a12*A.a13*RR - A.a22*A.a23);
				L.data[index++] = 2*(A.a12*B.a1*RR - A.a22*B.a2);
				L.data[index++] = A.a13*A.a13*RR - A.a23*A.a23;
				L.data[index++] = 2*(A.a13*B.a1*RR - A.a23*B.a2);
				L.data[index++] = B.a1*B.a1*RR - B.a2*B.a2;
			}
		}
	}
611029
def decompose_covariance(c):
    

    # make it a kickass copy of the original
    c = _n.array(c)

    # first get the error vector
    e = []
    for n in range(0, len(c[0])): e.append(_n.sqrt(c[n][n]))

    # now cycle through the matrix, dividing by e[1]*e[2]
    for n in range(0, len(c[0])):
        for m in range(0, len(c[0])):
            c[n][m] = c[n][m] / (e[n]*e[m])

    return [_n.array(e), _n.array(c)]
1091866
def _parallel_build_estimators(n_estimators, ensemble, X, y, cost_mat,
                               seeds, verbose):
    """"""
    # Retrieve settings
    n_samples, n_features = X.shape
    max_samples = ensemble.max_samples
    max_features = ensemble.max_features

    if (not isinstance(max_samples, (numbers.Integral, np.integer)) and
            (0.0 < max_samples <= 1.0)):
        max_samples = int(max_samples * n_samples)

    if (not isinstance(max_features, (numbers.Integral, np.integer)) and
            (0.0 < max_features <= 1.0)):
        max_features = int(max_features * n_features)

    bootstrap = ensemble.bootstrap
    bootstrap_features = ensemble.bootstrap_features

    # Build estimators
    estimators = []
    estimators_samples = []
    estimators_features = []

    for i in range(n_estimators):
        if verbose > 1:
            print(("building estimator %d of %d" % (i + 1, n_estimators)))

        random_state = check_random_state(seeds[i])
        seed = check_random_state(random_state.randint(MAX_INT))
        estimator = ensemble._make_estimator(append=False)

        try:  # Not all estimator accept a random_state
            estimator.set_params(random_state=seed)
        except ValueError:
            pass

        # Draw features
        if bootstrap_features:
            features = random_state.randint(0, n_features, max_features)
        else:
            features = sample_without_replacement(n_features,
                                                  max_features,
                                                  random_state=random_state)

        # Draw samples, using a mask, and then fit
        if bootstrap:
            indices = random_state.randint(0, n_samples, max_samples)
        else:
            indices = sample_without_replacement(n_samples,
                                                 max_samples,
                                                 random_state=random_state)

        sample_counts = np.bincount(indices, minlength=n_samples)

        estimator.fit((X[indices])[:, features], y[indices], cost_mat[indices, :])
        samples = sample_counts > 0.

        estimators.append(estimator)
        estimators_samples.append(samples)
        estimators_features.append(features)

    return estimators, estimators_samples, estimators_features
17299
def spm(x, model, y=None, n_samples=None, dx_min=-0.1,
        dx_max=0.1, n_dxs=5, dy_min=-0.1, dy_max=0.1, n_dys=5,
        angle_min=-30, angle_max=30, n_angles=31, black_border_size=0):
  
  if y is None:
    preds = model.get_probs(x)
    # Using model predictions as ground truth to avoid label leaking
    preds_max = reduce_max(preds, 1, keepdims=True)
    y = tf.to_float(tf.equal(preds, preds_max))
    y = tf.stop_gradient(y)
    del preds
  y = y / reduce_sum(y, 1, keepdims=True)

  # Define the range of transformations
  dxs = np.linspace(dx_min, dx_max, n_dxs)
  dys = np.linspace(dy_min, dy_max, n_dys)
  angles = np.linspace(angle_min, angle_max, n_angles)

  if n_samples is None:
    import itertools
    transforms = list(itertools.product(*[dxs, dys, angles]))
  else:
    sampled_dxs = np.random.choice(dxs, n_samples)
    sampled_dys = np.random.choice(dys, n_samples)
    sampled_angles = np.random.choice(angles, n_samples)
    transforms = zip(sampled_dxs, sampled_dys, sampled_angles)
  transformed_ims = parallel_apply_transformations(
      x, transforms, black_border_size)

  def _compute_xent(x):
    preds = model.get_logits(x)
    return tf.nn.softmax_cross_entropy_with_logits_v2(
        labels=y, logits=preds)

  all_xents = tf.map_fn(
      _compute_xent,
      transformed_ims,
      parallel_iterations=1)  # Must be 1 to avoid keras race conditions

  # Return the adv_x with worst accuracy

  # all_xents is n_total_samples x batch_size (SB)
  all_xents = tf.stack(all_xents)  # SB

  # We want the worst case sample, with the largest xent_loss
  worst_sample_idx = tf.argmax(all_xents, axis=0)  # B

  batch_size = tf.shape(x)[0]
  keys = tf.stack([
      tf.range(batch_size, dtype=tf.int32),
      tf.cast(worst_sample_idx, tf.int32)
  ], axis=1)
  transformed_ims_bshwc = tf.einsum('sbhwc->bshwc', transformed_ims)
  after_lookup = tf.gather_nd(transformed_ims_bshwc, keys)  # BHWC
  return after_lookup
333910
def _Moran_BV_Matrix_array(variables, w, permutations=0, varnames=None):
    
    if varnames is None:
        varnames = ['x{}'.format(i) for i in range(k)]

    k = len(variables)
    rk = list(range(0, k - 1))
    results = {}
    for i in rk:
        for j in range(i + 1, k):
            y1 = variables[i]
            y2 = variables[j]
            results[i, j] = Moran_BV(y1, y2, w, permutations=permutations)
            results[j, i] = Moran_BV(y2, y1, w, permutations=permutations)
            results[i, j].varnames = {'x': varnames[i], 'y': varnames[j]}
            results[j, i].varnames = {'x': varnames[j], 'y': varnames[i]}
    return results
682373
def _build_matrix_non_uniform(p, q, coords, k):
    """"""
    A = [[1] * (p+q+1)]
    for i in range(1, p + q + 1):
        line = [(coords[k+j] - coords[k])**i for j in range(-p, q+1)]
        A.append(line)
    return np.array(A)
737479
def convert_to_adjacency_matrix(matrix):
    
    for i in range(matrix.shape[0]):
        
        if isspmatrix(matrix):
            col = find(matrix[:,i])[2]
        else:
            col = matrix[:,i].T.tolist()[0]

        coeff = max( Fraction(c).limit_denominator().denominator for c in col )
        matrix[:,i] *= coeff

    return matrix
188936
def fitlin(imgarr,refarr):
    
    # Initialize variables
    _mat = np.zeros((3,3),dtype=np.float64)
    _xorg = imgarr[0][0]
    _yorg = imgarr[0][1]
    _xoorg = refarr[0][0]
    _yoorg = refarr[0][1]
    _sigxox = 0.
    _sigxoy = 0.
    _sigxo = 0.
    _sigyox = 0.
    _sigyoy = 0.
    _sigyo = 0.

    _npos = len(imgarr)
    # Populate matrices
    for i in range(_npos):
        _mat[0][0] += np.power((imgarr[i][0] - _xorg),2)
        _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg)
        _mat[0][2] += (imgarr[i][0] - _xorg)
        _mat[1][1] += np.power((imgarr[i][1] - _yorg),2)
        _mat[1][2] += imgarr[i][1] - _yorg

        _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg)
        _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg)
        _sigxo += refarr[i][0] - _xoorg
        _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg)
        _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg)
        _sigyo += refarr[i][1] - _yoorg

    _mat[2][2] = _npos
    _mat[1][0] = _mat[0][1]
    _mat[2][0] = _mat[0][2]
    _mat[2][1] = _mat[1][2]

    # Now invert this matrix
    _mat = linalg.inv(_mat)

    _a  = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2]
    _b  = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2])
    #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2]

    _c  = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2]
    _d  = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2]
    #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2]

    _xt = _xoorg - _a*_xorg+_b*_yorg
    _yt = _yoorg - _d*_xorg-_c*_yorg

    return [_a,_b,_xt],[_c,_d,_yt]
1364793
def build_model_from_xy(self, x_values, y_values):
        """"""
        self.init_ace(x_values, y_values)
        self.run_ace()
        self.build_interpolators()
311851
def get_matrix(self, indices):
        """"""
        new = numpy.empty(self.samples1.shape)
        for idx in range(len(indices)):
            if indices[idx]:
                new[idx] = self.samples1[idx]
            else:
                new[idx] = self.samples2[idx]

        if self.poly:
            new = self.poly(*new)
        return new
623725
def matrix(self, _):
        """"""
        x, y, z = self.coord
        return np.array([
            [1., 0., 0., x],
            [0., 1., 0., y],
            [0., 0., 1., z],
            [0., 0., 0., 1.]
        ])
591925
def build_model(self, n_features, regtype='none'):
        
        self._create_placeholders(n_features)
        self._create_variables(n_features)
        self.encode = self.sample_hidden_from_visible(self.input_data)[0]
        self.reconstruction = self.sample_visible_from_hidden(
            self.encode, n_features)

        hprob0, hstate0, vprob, hprob1, hstate1 = self.gibbs_sampling_step(
            self.input_data, n_features)
        positive = self.compute_positive_association(self.input_data,
                                                     hprob0, hstate0)

        nn_input = vprob

        for step in range(self.gibbs_sampling_steps - 1):
            hprob, hstate, vprob, hprob1, hstate1 = self.gibbs_sampling_step(
                nn_input, n_features)
            nn_input = vprob

        negative = tf.matmul(tf.transpose(vprob), hprob1)

        self.w_upd8 = self.W.assign_add(
            self.learning_rate * (positive - negative) / self.batch_size)

        self.bh_upd8 = self.bh_.assign_add(tf.multiply(self.learning_rate, tf.reduce_mean(
            tf.subtract(hprob0, hprob1), 0)))

        self.bv_upd8 = self.bv_.assign_add(tf.multiply(self.learning_rate, tf.reduce_mean(
            tf.subtract(self.input_data, vprob), 0)))

        variables = [self.W, self.bh_, self.bv_]
        regterm = Layers.regularization(variables, self.regtype, self.regcoef)

        self.cost = self.loss.compile(vprob, self.input_data, regterm=regterm)
619129
def inference(self, kern, X, likelihood, Y, Y_metadata=None):

        
        N = X.shape[0] #number of training points
        D = X.shape[1] #number of dimensions

        Kds = np.zeros(D, dtype=object) #vector for holding covariance per dimension
        Qs = np.zeros(D, dtype=object) #vector for holding eigenvectors of covariance per dimension
        QTs = np.zeros(D, dtype=object) #vector for holding transposed eigenvectors of covariance per dimension
        V_kron = 1 # kronecker product of eigenvalues

        # retrieve the one-dimensional variation of the designated kernel
        oneDkernel = kern.get_one_dimensional_kernel(D)

        for d in range(D):
            xg = list(set(X[:,d])) #extract unique values for a dimension
            xg = np.reshape(xg, (len(xg), 1))
            oneDkernel.lengthscale = kern.lengthscale[d]
            Kds[d] = oneDkernel.K(xg)
            [V, Q] = np.linalg.eig(Kds[d])
            V_kron = np.kron(V_kron, V)
            Qs[d] = Q
            QTs[d] = Q.T

        noise = likelihood.variance + 1e-8

        alpha_kron = self.kron_mvprod(QTs, Y)
        V_kron = V_kron.reshape(-1, 1)
        alpha_kron = alpha_kron / (V_kron + noise)
        alpha_kron = self.kron_mvprod(Qs, alpha_kron)

        log_likelihood = -0.5 * (np.dot(Y.T, alpha_kron) + np.sum((np.log(V_kron + noise))) + N*log_2_pi)

        # compute derivatives wrt parameters Thete
        derivs = np.zeros(D+2, dtype='object')
        for t in range(len(derivs)):
            dKd_dTheta = np.zeros(D, dtype='object')
            gamma = np.zeros(D, dtype='object')
            gam = 1
            for d in range(D):
                xg = list(set(X[:,d]))
                xg = np.reshape(xg, (len(xg), 1))
                oneDkernel.lengthscale = kern.lengthscale[d]
                if t < D:
                    dKd_dTheta[d] = oneDkernel.dKd_dLen(xg, (t==d), lengthscale=kern.lengthscale[t]) #derivative wrt lengthscale
                elif (t == D):
                    dKd_dTheta[d] = oneDkernel.dKd_dVar(xg) #derivative wrt variance
                else:
                    dKd_dTheta[d] = np.identity(len(xg)) #derivative wrt noise
                gamma[d] = np.diag(np.dot(np.dot(QTs[d], dKd_dTheta[d].T), Qs[d]))
                gam = np.kron(gam, gamma[d])
            
            gam = gam.reshape(-1,1)
            kappa = self.kron_mvprod(dKd_dTheta, alpha_kron)
            derivs[t] = 0.5*np.dot(alpha_kron.T,kappa) - 0.5*np.sum(gam / (V_kron + noise))

        # separate derivatives
        dL_dLen = derivs[:D]
        dL_dVar = derivs[D]
        dL_dThetaL = derivs[D+1]

        return GridPosterior(alpha_kron=alpha_kron, QTs=QTs, Qs=Qs, V_kron=V_kron), \
                log_likelihood, {'dL_dLen':dL_dLen, 'dL_dVar':dL_dVar, 'dL_dthetaL':dL_dThetaL}