General handling of quadrature expectations for Gaussians and DiagonalGaussians
    Fallback method for missing analytic expectations
540337
def _quadrature_expectation(p, obj1, feature1, obj2, feature2, num_gauss_hermite_points):
    
    num_gauss_hermite_points = 100 if num_gauss_hermite_points is None else num_gauss_hermite_points

    if obj2 is None:
        eval_func = lambda x: get_eval_func(obj1, feature1)(x)
    elif obj1 is None:
        raise NotImplementedError("First object cannot be None.")
    else:
        eval_func = lambda x: (get_eval_func(obj1, feature1, np.s_[:, :, None])(x) *
                               get_eval_func(obj2, feature2, np.s_[:, None, :])(x))

    if isinstance(p, DiagonalGaussian):
        if isinstance(obj1, kernels.Kernel) and isinstance(obj2, kernels.Kernel) \
                and obj1.on_separate_dims(obj2):  # no joint expectations required

            eKxz1 = quadrature_expectation(p, (obj1, feature1),
                                           num_gauss_hermite_points=num_gauss_hermite_points)
            eKxz2 = quadrature_expectation(p, (obj2, feature2),
                                           num_gauss_hermite_points=num_gauss_hermite_points)
            return eKxz1[:, :, None] * eKxz2[:, None, :]

        else:
            cov = tf.matrix_diag(p.cov)
    else:
        cov = p.cov
    return mvnquad(eval_func, p.mu, cov, num_gauss_hermite_points)
372028
def solve_gamlasso(self, lam):
        ''''''
        weights = lam / (1 + self.gamma * np.abs(self.beta[self.trails[::2]] - self.beta[self.trails[1::2]]))
        s = self.solve_gfl(u)
        self.steps.append(s)
        return self.beta
1120144
def _gauss_funct(p, fjac=None, x=None, y=None, err=None,
                 weights=None):

    
    if p[2] != 0.0:
        Z = (x - p[1]) / p[2]
        model = p[0] * np.e ** (-Z ** 2 / 2.0)
    else:
        model = np.zeros(np.size(x))

    status = 0
    if weights is not None:
        if err is not None:
            print("Warning: Ignoring errors and using weights.\n")

        return [status, (y - model) * weights]

    elif err is not None:
        return [status, (y - model) / err]

    else:
        return [status, y - model]
1106108
def _expectation(p, kern, feat, mean, none, nghp=None):
    
    Xmu, Xcov = p.mu, p.cov

    with tf.control_dependencies([tf.assert_equal(
            tf.shape(Xmu)[1], tf.constant(kern.input_dim, settings.int_type),
            message="Currently cannot handle slicing in exKxz.")]):
        Xmu = tf.identity(Xmu)

    with params_as_tensors_for(kern, feat):
        N = tf.shape(Xmu)[0] - 1
        var_Z = kern.variance * feat.Z  # MxD
        tiled_Z = tf.tile(tf.expand_dims(var_Z, 0), (N, 1, 1))  # NxMxD
        eXX = Xcov[1, :-1] + (Xmu[:-1][..., None] * Xmu[1:][:, None, :])  # NxDxD
        return tf.matmul(tiled_Z, eXX)
1345552
def kernel(x1, x2, method='gaussian', sigma=1, **kwargs):
    """"""
    if method.lower() in ['gaussian', 'gauss', 'rbf']:
        K = np.exp(-dist(x1, x2) / (2 * sigma**2))
    return K
542295
def _expectation(p, kern, none1, none2, none3, nghp=None):
    
    return kern.Kdiag(p.mu)
41025
def _compute_errors(self):
        

        try:
            coeffs = fit_first_and_second_harmonics(self.sample.values[0],
                                                    self.sample.values[2])
            covariance = coeffs[1]
            coeffs = coeffs[0]
            model = first_and_second_harmonic_function(self.sample.values[0],
                                                       coeffs)
            residual_rms = np.std(self.sample.values[2] - model)
            errors = np.diagonal(covariance) * residual_rms

            eps = self.sample.geometry.eps
            pa = self.sample.geometry.pa

            # parameter errors result from direct projection of
            # coefficient errors.  These showed to be the error estimators
            # that best convey the errors measured in Monte Carlo
            # experiments (see Busko 1996; ASPC 101, 139).
            ea = abs(errors[2] / self.grad)
            eb = abs(errors[1] * (1. - eps) / self.grad)
            self.x0_err = np.sqrt((ea * np.cos(pa))**2 + (eb * np.sin(pa))**2)
            self.y0_err = np.sqrt((ea * np.sin(pa))**2 + (eb * np.cos(pa))**2)
            self.ellip_err = (abs(2. * errors[4] * (1. - eps) / self.sma /
                                  self.grad))
            if (abs(eps) > np.finfo(float).resolution):
                self.pa_err = (abs(2. * errors[3] * (1. - eps) / self.sma /
                                   self.grad / (1. - (1. - eps)**2)))
            else:
                self.pa_err = 0.
        except Exception:    # we want to catch everything
            self.x0_err = self.y0_err = self.pa_err = self.ellip_err = 0.
1354308
def _preoptimize_model(self):
        
        gaussian_model = DynReg(formula=self.formula, data=self.data_original)
        gaussian_model.fit()

        for i in range(self.z_no-self.family_z_no):
            self.latent_variables.z_list[i].start = gaussian_model.latent_variables.get_z_values()[i+1]

        if self.model_name2 == 't':

            def temp_function(params):
                return -np.sum(ss.t.logpdf(x=self.data, df=np.exp(params[0]), 
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2])))

            p = optimize.minimize(temp_function,np.array([2.0,0.0,-1.0]),method='L-BFGS-B')
            self.latent_variables.z_list[-2].start = p.x[2]
            self.latent_variables.z_list[-1].start = p.x[0]

        elif self.model_name2 == 'Skewt':

            def temp_function(params):
                return -np.sum(fam.Skewt.logpdf_internal(x=self.data,df=np.exp(params[0]),
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2]),gamma=np.exp(params[3])))

            p = optimize.minimize(temp_function,np.array([2.0,0.0,-1.0,0.0]),method='L-BFGS-B')
            self.latent_variables.z_list[-3].start = p.x[3]
            self.latent_variables.z_list[-2].start = p.x[2]
            self.latent_variables.z_list[-1].start = p.x[0]

        return gaussian_model.latent_variables
977586
def _expectation(p, kern1, feat1, kern2, feat2, nghp=None):
    
    """
    crossexps = []

    if kern1 == kern2 and feat1 == feat2:  # avoid duplicate computation by using transposes
        for i, k1 in enumerate(kern1.kernels):
            crossexps.append(expectation(p, (k1, feat1), (k1, feat1), nghp=nghp))

            for k2 in kern1.kernels[:i]:
                eKK = expectation(p, (k1, feat1), (k2, feat2), nghp=nghp)
                eKK += tf.matrix_transpose(eKK)
                crossexps.append(eKK)
    else:
        for k1, k2 in it.product(kern1.kernels, kern2.kernels):
            crossexps.append(expectation(p, (k1, feat1), (k2, feat2), nghp=nghp))

    return functools.reduce(tf.add, crossexps)
748837
def _initialize_likelihood_prior(self, positions, log_likelihoods, log_priors):
        
        func = SimpleCLFunction.from_string('''
            void compute(global mot_float_type* chain_position,
                         global mot_float_type* log_likelihood,
                         global mot_float_type* log_prior,
                         local mot_float_type* x_tmp,
                         void* data){
                
                bool is_first_work_item = get_local_id(0) == 0;
    
                if(is_first_work_item){
                    for(uint i = 0; i < ''' + str(self._nmr_params) + '''; i++){
                        x_tmp[i] = chain_position[i];
                    }
                    *log_prior = _computeLogPrior(x_tmp, data);
                }
                barrier(CLK_LOCAL_MEM_FENCE);
                    
                *log_likelihood = _computeLogLikelihood(x_tmp, data);
            }
        ''', dependencies=[self._get_log_prior_cl_func(), self._get_log_likelihood_cl_func()])

        kernel_data = {
            'chain_position': Array(positions, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_likelihood': Array(log_likelihoods, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_prior': Array(log_priors, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'x_tmp': LocalMemory('mot_float_type', self._nmr_params),
            'data': self._data
        }

        func.evaluate(kernel_data, self._nmr_problems,
                      use_local_reduction=all(env.is_gpu for env in self._cl_runtime_info.cl_environments),
                      cl_runtime_info=self._cl_runtime_info)
913529
def make_gaussian_kernel(sigma, npix=501, cdelt=0.01, xpix=None, ypix=None):
    

    sigma /= cdelt

    def fn(t, s): return 1. / (2 * np.pi * s ** 2) * np.exp(
        -t ** 2 / (s ** 2 * 2.0))
    dxy = make_pixel_distance(npix, xpix, ypix)
    k = fn(dxy, sigma)
    k /= (np.sum(k) * np.radians(cdelt) ** 2)

    return k
63206
public int updateMixture( float[] pixelValue , float[] dataRow , int modelIndex ) {

		// see which gaussian is the best fit based on Mahalanobis distance
		int index = modelIndex;
		float bestDistance = maxDistance*numBands;
		int bestIndex=-1;

		int ng; // number of gaussians in use
		for (ng = 0; ng < maxGaussians; ng++, index += gaussianStride) {
			float variance = dataRow[index+1];
			if( variance <= 0 ) {
				break;
			}

			float mahalanobis = 0;
			for (int i = 0; i < numBands; i++) {
				float mean = dataRow[index+2+i];
				float delta = pixelValue[i]-mean;
				mahalanobis += delta*delta/variance;
			}

			if( mahalanobis < bestDistance ) {
				bestDistance = mahalanobis;
				bestIndex = index;
			}
		}

		// Update the model for the best gaussian
		if( bestIndex != -1 ) {
			// If there is a good fit update the model
			float weight = dataRow[bestIndex];
			float variance = dataRow[bestIndex+1];

			weight += learningRate*(1f-weight);
			dataRow[bestIndex]   = 1; // set to one so that it can't possible go negative

			float sumDeltaSq = 0;
			for (int i = 0; i < numBands; i++) {
				float mean = dataRow[bestIndex+2+i];
				float delta = pixelValue[i]-mean;
				dataRow[bestIndex+2+i] = mean + delta*learningRate/weight;
				sumDeltaSq += delta*delta;
			}
			sumDeltaSq /= numBands;
			dataRow[bestIndex+1] = variance + (learningRate /weight)*(sumDeltaSq*1.2F - variance);
			// 1.2f is a fudge factor. Empirical testing shows that the above equation is biased. Can't be bothered
			// to verify the derivation and see if there's a mistake

			// Update Gaussian weights and prune models
			updateWeightAndPrune(dataRow, modelIndex, ng, bestIndex, weight);

			return weight >= significantWeight ? 0 : 1;
		} else if( ng < maxGaussians ) {
			// if there is no good fit then create a new model, if there is room

			bestIndex = modelIndex + ng*gaussianStride;
			dataRow[bestIndex]   = 1; // weight is changed later or it's the only model
			dataRow[bestIndex+1] = initialVariance;
			for (int i = 0; i < numBands; i++) {
				dataRow[bestIndex+2+i] = pixelValue[i];
			}
			// There are no models. Return unknown
			if( ng == 0 )
				return unknownValue;

			updateWeightAndPrune(dataRow, modelIndex, ng+1, bestIndex, learningRate);
			return 1;
		} else {
			// didn't match any models and can't create a new model
			return 1;
		}
	}
1226488
def gaussian(self, x, p):
        
        y = (1.0 / (p[1] * np.sqrt(2 * np.pi)) *
             np.exp(-(x - p[0]) ** 2 / (2 * p[1] ** 2))) * p[2]
        return y
258587
@Override
  public FittingFunctionResult eval(double x, double[] params) {
    final int len = params.length;
    // We always need triples: (mean, stddev, scaling)
    assert (len % 3) == 0;

    double y = 0.0;
    double[] gradients = new double[len];

    // Loosely based on the book:
    // Numerical Recipes in C: The Art of Scientific Computing
    // Due to their license, we cannot use their code, but we have to implement
    // the mathematics ourselves. We hope the loss in precision is not too big.
    for(int i = 2; i < params.length; i += 3) {
      // Standardized Gaussian parameter (centered, scaled by stddev)
      double stdpar = (x - params[i - 2]) / params[i - 1];
      double e = FastMath.exp(-.5 * stdpar * stdpar);
      double localy = params[i] / (params[i - 1] * MathUtil.SQRTTWOPI) * e;

      y += localy;
      // mean gradient
      gradients[i - 2] = localy * stdpar;
      // stddev gradient
      gradients[i - 1] = (stdpar * stdpar - 1.0) * localy;
      // amplitude gradient
      gradients[i] = e / (params[i - 1] * MathUtil.SQRTTWOPI);
    }

    return new FittingFunctionResult(y, gradients);
  }
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
946167
def laplacian_of_gaussian(image, mask, size, sigma):
    '''
    '''
    half_size = size//2
    i,j = np.mgrid[-half_size:half_size+1,
                   -half_size:half_size+1].astype(float) / float(sigma)
    distance = (i**2 + j**2)/2
    gaussian = np.exp(-distance)
    #
    # Normalize the Gaussian
    #
    gaussian = gaussian / np.sum(gaussian)

    log = (distance - 1) * gaussian
    #
    # Normalize the kernel to have a sum of zero
    #
    log = log - np.mean(log)
    if mask is None:
        mask = np.ones(image.shape[:2], bool)
    masked_image = image.copy()
    masked_image[~mask] = 0
    output = convolve(masked_image, log, mode='constant', cval=0)
    #
    # Do the LoG of the inverse of the mask. This finds the magnitude of the
    # contribution of the masked pixels. We then fudge by multiplying by the
    # value at the pixel of interest - this effectively sets the value at a
    # masked pixel to that of the pixel of interest.
    #
    # It underestimates the LoG, that's not a terrible thing.
    #
    correction = convolve((~ mask).astype(float), log, mode='constant', cval = 1)
    output += correction * image
    output[~ mask] = image[~ mask]
    return output
763014
def Kdiag(self, X):
        """"""
        ret = np.empty(X.shape[0])
        ret[:] = self.variance
        return ret