C++ Gaussian类代码示例

  /** @note Because of the virtual public inheritance, we create the Pdf
      baseclass ourselves
  */
  CondGausAddNoise::CondGausAddNoise(const Gaussian& additiveNoise,
				     int num_conditional_arguments)
    : AnalyticConditionalGaussian(additiveNoise.DimensionGet(),
			  num_conditional_arguments)
    , _additiveNoise_Mu   (additiveNoise.ExpectedValueGet())
    , _additiveNoise_Sigma(additiveNoise.CovarianceGet())
  {}

void SumFactor::_internal_update(
    Variable* var,
    std::vector<Gaussian*> y,
    std::vector<Gaussian*> fy,
    std::vector<double>* a) {

    double sum_pi = 0.0, sum_tau = 0.0, new_pi, new_tau, da;
    unsigned int i = 0, size = a->size();
    Gaussian gy, gfy;

    for (i = 0; i < size; ++i) {
        da = (*a)[i];
        gy = *y[i];
        gfy = *fy[i];

        sum_pi = sum_pi + ((da * da) / (gy.pi - gfy.pi));
        sum_tau = sum_tau + (da * (gy.tau - gfy.tau) / (gy.pi - gfy.pi));
    }

    new_pi = 1.0 / sum_pi;
    new_tau = new_pi * sum_tau;

    Gaussian* gaussian = new Gaussian();
    gaussian->init_pi_tau(new_pi, new_tau);
    var->update_message(this, gaussian);
}

//----------------------------------------------------------------------------------------------------------------------  
void SMCEnvironmentViz::drawGaussian(const Gaussian& g) const
{
  const int numPoints = 50;
  ci::Vec2f start, end;
  
  // Three sigmas cover approx. 99% of the gaussian
  float a = -3.0 * g.getW();
  float inc = 6.0 * g.getW() / numPoints;
  
  ci::Path2d path;
  path.moveTo(g.getPointAt(a));
  for(int i = 0; i < numPoints; ++i)
  {
    end = g.getPointAt((a+inc));
    path.lineTo(end);
    a += inc;
  }
  path.close();
  
  glColor4f(0.5,0.5,0.5,0.25);
  glDisable(GL_CULL_FACE);
  glDisable(GL_DEPTH_TEST);
  ci::gl::drawSolid(path);
  glEnable(GL_CULL_FACE);
  glEnable(GL_DEPTH_TEST);
}

    bool is_approx(const Gaussian<OtherVariate>& other)
    {


        return fl::are_similar(mean(), other.mean()) &&
               fl::are_similar(covariance(), other.covariance());
    }

// load the HMM from a file
void HMMState::load(FileInput &file, unsigned char iEstimationMethod) {

	// phonetic symbol
	char strPhone[MAX_PHONETIC_SYMBOL_LENGTH+1];	
	IOBase::readBytes(file.getStream(),reinterpret_cast<char*>(strPhone),MAX_PHONETIC_SYMBOL_LENGTH+1);
	m_iPhone = m_phoneSet->getPhoneIndex(strPhone);
	assert(m_iPhone != UCHAR_MAX);
	
	// state
	IOBase::read(file.getStream(),&m_iState);
	assert(m_iState < NUMBER_HMM_STATES);
	
	// within word position (DEPRECATED)
	IOBase::read(file.getStream(),&m_iPosition);
	
	// Gaussian components
	int iGaussianComponents = -1;
	IOBase::read(file.getStream(),&iGaussianComponents);
	for(int iGaussian = 0 ; iGaussian < iGaussianComponents ; ++iGaussian) {
	
		Gaussian *gaussian = new Gaussian(m_iDim,m_iCovarianceModeling);	
		
		IOBase::read(file.getStream(),&gaussian->weight());
		gaussian->mean().readData(file.getStream());		
		if (m_iCovarianceModeling == COVARIANCE_MODELLING_TYPE_DIAGONAL) {
			gaussian->covarianceDiag().readData(file.getStream());
		} else {
			gaussian->covarianceFull().readData(file.getStream());
		}	
		m_gaussianMixture->addGaussianComponent(gaussian);
	}
}

// store the HMM into a file
void HMMState::store(FileOutput &file) {
	
	// phonetic symbol
	char strPhone[MAX_PHONETIC_SYMBOL_LENGTH+1];
	memset(strPhone,0,MAX_PHONETIC_SYMBOL_LENGTH+1);
	strcpy(strPhone,m_phoneSet->getStrPhone(m_iPhone));	
	IOBase::writeBytes(file.getStream(),reinterpret_cast<char*>(strPhone),MAX_PHONETIC_SYMBOL_LENGTH+1);
	
	// state
	IOBase::write(file.getStream(),m_iState);
	
	// within word position (DEPRECATED)
	IOBase::write(file.getStream(),m_iPosition);

	// Gaussian components
	int iGaussianComponents = m_gaussianMixture->getNumberComponents();
	IOBase::write(file.getStream(),iGaussianComponents);
	for(int iGaussian = 0 ; iGaussian < iGaussianComponents ; ++iGaussian) {
		Gaussian *gaussian = (*m_gaussianMixture)(iGaussian); 
		IOBase::write(file.getStream(),gaussian->weight());
		gaussian->mean().writeData(file.getStream());
		if (m_iCovarianceModeling == COVARIANCE_MODELLING_TYPE_DIAGONAL) {
			gaussian->covarianceDiag().writeData(file.getStream());
		} else {
			gaussian->covarianceFull().writeData(file.getStream());
		}
	}
}

void OpenCLExecuter::ocl_filterGaussian_init(float sigma, float size)
{
	Gaussian gauss;
	gauss.sigma = sigma;
	gauss.size = size;
	filter_width = (size-1)/2.0;
	filter_kernel = gauss.gaussianmask1Df(gauss.sigma, gauss.size);
}

void testGaussianNumerics(){
  
  std::cout.precision(16);
  std::cout << std::scientific;
    
  Gaussian g = Gaussian(5, 1);
  for(unsigned i = 0; i < 100; i++){
    std::cout << i << ": L: " << g.likelihood((fracfloat_t)i) << " -LL: " << (-INV_LN_2 * log(g.likelihood(i))) << ", S: " << g.surprisal((fracfloat_t)i) << std::endl;
  }
}

void GMMfitData::defaultAonly(Gaussian &g)
{
    Vector2d mu;
    mu<< 1.0, 0.1;
    g.setMu(mu);

    Matrix2d sigma;
    sigma<< 0.020,   -0.001,
           -0.001,    0.005;
    g.setSigma(sigma);
}

void GMMfitData::defaultDonly(Gaussian &g)
{
    Vector2d mu;
    mu<< 0.0, 0.98;
    g.setMu(mu);

    Matrix2d sigma;
    sigma<< 0.0012,    0.0000,
            0.0000,    0.0015;
    g.setSigma(sigma);
}

void GMMfitData::defaultLowFRET(Gaussian &g)
{
    Vector2d mu;
    mu<< 0.15, 0.5;
    g.setMu(mu);

    Matrix2d sigma;
    sigma<< 0.050,    0.000,
            0.000,    0.020;
    g.setSigma(sigma);
}

void vwcgauss
(const RealArray& energy, const RealArray& parameter, 
 /*@[email protected]*/ int spectrum, RealArray& flux, /*@[email protected]*/ RealArray& fluxError,
 /*@[email protected]*/ const string& init)
{
  fluxError.resize (0);
  Real LineEnergy, Sigma;
  vwcgaussProcessParameter (parameter, LineEnergy, Sigma);
  Gaussian G (energy, LineEnergy, Sigma);
  G.getFlux (flux);
  return;
}

void Test::RunCannyTest(const RGBImage* in) {
	Gaussian smoother;
	Sobel sobel;
	Canny canny;
	StudentPreProcessing converter;

	double sigma = 1.6;
	int kernelSize = 5;

	// Convert RGB to Intensity
	IntensityImage* grey = ImageFactory::newIntensityImage(in->getWidth(), in->getHeight());
	grey = converter.stepToIntensityImage(*in);

	BaseTimer timer;
	int total = 0;
	for (int i = 0; i < 10; ++i) {
		timer.start();

		/* Start of Canny */

		// 1. Smooth Intensity with Gaussian
		IntensityImage* gaussian = ImageFactory::newIntensityImage(in->getWidth() - kernelSize, in->getHeight() - kernelSize);
		smoother.smoothImage(grey, gaussian, sigma, kernelSize);

		// 2. Find edges with Sobel (convolve X and Y separately for Canny)
		int sobelKernelSize = 3;
		IntensityImage* sobelX = ImageFactory::newIntensityImage(gaussian->getWidth() - sobelKernelSize, gaussian->getHeight() - sobelKernelSize);
		IntensityImage* sobelY = ImageFactory::newIntensityImage(gaussian->getWidth() - sobelKernelSize, gaussian->getHeight() - sobelKernelSize);
		sobel.filterXY(gaussian, sobelX, sobelY);

		// 3. Apply non-maximum suppresion
		IntensityImage* nonMaxSuppression = ImageFactory::newIntensityImage(sobelX->getWidth(), sobelX->getHeight());
		canny.nonMaximumSurpression(sobelX, sobelY, nonMaxSuppression);

		// 4. Apply hysteresis threshold
		IntensityImage* hysteresisThreshold = ImageFactory::newIntensityImage(sobelX->getWidth(), sobelX->getHeight());
		canny.threshold(nonMaxSuppression, hysteresisThreshold, 60, 70);

		timer.stop();
		std::cout << timer.elapsedMilliSeconds() << " ms" << std::endl;
		total += timer.elapsedMilliSeconds();
		timer.reset();

		if (i == 9)
			ImageIO::saveIntensityImage(*hysteresisThreshold, ImageIO::getDebugFileName("Canny.png"));
	}
	std::cout << "Average time per Canny edge detection: " << (total / 10) << " ms" << std::endl;
	std::cout << "Press the X in the console window to exit program" << std::endl;
}

/**
 * Create a time series of input data.
 *
 * Description:\n
 * 	Creates a time series of complex floats representing a signal at
 * 	a given frequency.
 */
ComplexFloat32 *
createTdData(float64_t fMHz, float64_t driftHz, int32_t samples,
		float64_t noise)
{
	size_t size = samples * sizeof(ComplexFloat32);
	ComplexFloat32 *td = static_cast<ComplexFloat32 *> (fftwf_malloc(size));
	float64_t dfMHz = fMHz - centerFreqMHz;
	float64_t widthMHz = subchannels * subchannelWidthMHz;
	Gaussian gen;
	// no noise, just signal
	gen.setup(0, widthMHz, noise);
	gen.addCwSignal(dfMHz, driftHz, 1);
	gen.getSamples(td, samples);
	return (td);
}

  // constructor
  NonLinearAnalyticConditionalGaussian_Ginac::NonLinearAnalyticConditionalGaussian_Ginac
  (const GiNaC::matrix& func,
   const vector<GiNaC::symbol>& u,
   const vector<GiNaC::symbol>& x,
   const Gaussian& additiveNoise,
   const vector<GiNaC::symbol>& cond )
    :AnalyticConditionalGaussianAdditiveNoise(additiveNoise,3),
     func_sym        (func),
     cond_sym        (cond),
     u_sym           (u),
     x_sym           (x),
     cond_size       (cond_sym.size()),
     u_size          (u_sym.size()),
     x_size          (x_sym.size()),
     func_size       (func_sym.rows()),
     dfunc_dcond     (cond_size),
     dfunc_dx        (x_size)
  {
    // test for consistent input
    assert (func_sym.cols() == 1);
    assert (additiveNoise.DimensionGet() == cond_size);

    // derive func to cond
    for (unsigned int i=0; i < cond_size; i++)
      dfunc_dcond[i] = func_sym.diff(cond_sym[i]);

    // derive func to x
    for (unsigned int i=0; i < x_size; i++)
      dfunc_dx[i] = func_sym.diff(x_sym[i]);
  }