C++ To::get_mean方法代码示例

//! Return the signal to noise ratio
float Pulsar::AdaptiveSNR::get_snr (const Profile* profile)
{
  Reference::To<PhaseWeight> baseline;

  if (baseline_estimator)
    baseline = baseline_estimator->operate (profile);
  else
    baseline = profile->baseline();

  Estimate<double> mean = baseline->get_mean();
  Estimate<double> variance = baseline->get_variance();

  unsigned nbin = profile->get_nbin();
  unsigned off_pulse = baseline->get_nonzero_weight_count();
  unsigned on_pulse = nbin - off_pulse;

  double energy = profile->sum() - mean.val * nbin;
  double snr = energy / sqrt(variance.val*on_pulse);

  if (Profile::verbose)
  {
    double rms = sqrt (variance.val);

    cerr << "Pulsar::AdaptiveSNR::get_snr " << off_pulse << " out of " << nbin
	 << " bins in baseline\n  mean=" << mean << " rms=" << rms 
	 << " energy=" << energy << " S/N=" << snr << endl;
  }

  if (energy < 0)
    return 0.0;

  return snr;
}    

void Pulsar::LastHarmonic::build ()
{
  if (!profile)
    throw Error (InvalidState, "Pulsar::LastHarmonic::build",
		 "Profile not set");

  if (!baseline_estimator)
    throw Error (InvalidState, "Pulsar::LastHarmonic::build",
		 "Baseline estimator not set");

  Reference::To<PhaseWeight> baseline = baseline_estimator->baseline (profile);

  Estimate<double> mean = baseline->get_mean ();
  Estimate<double> var = baseline->get_variance ();
  Estimate<double> rms = sqrt(var);

#ifdef _DEBUG
  cerr << "Pulsar::LastHarmonic::build mean=" << mean
       << " rms=" << rms << endl;
#endif

  significant.find (profile, rms.get_value());

  // skip the DC term
  bin_rise = 1;
  bin_fall = significant.get();
}

void Pulsar::RemoveBaseline::Each::transform (Archive* archive)
{
  const unsigned nsub = archive->get_nsubint();
  const unsigned nchan = archive->get_nchan();
  const unsigned npol = archive->get_npol();

  bool pscrunch = (archive->get_state() == Signal::Coherence ||
		   archive->get_state() == Signal::PPQQ);

  for (unsigned isub=0; isub < nsub; isub++)
  {
    Integration* subint = archive->get_Integration (isub);
    for (unsigned ichan=0; ichan < nchan; ichan++)
    {
      Reference::To<Profile> profile = subint->get_Profile (0,ichan);
      if (pscrunch)
      {
	profile = profile->clone();
	profile->sum (subint->get_Profile (1,ichan));
      }

      Reference::To<PhaseWeight> baseline = profile->baseline();

      for (unsigned ipol=0; ipol < npol; ipol++)
      {
	Profile* p = subint->get_Profile(ipol, ichan);
	baseline->set_Profile (p);
	p->offset (-baseline->get_mean().val);
      }
    }
  }
};

void Pulsar::ComponentModel::fit (const Profile *profile)
{
  build ();

  unsigned nbin = profile->get_nbin(); 
  vector<float> data (nbin);

  float mean = 0.0;
  float variance = 1.0;

  if (fit_derivative)
  {
    vector<complex<float> > spec(nbin/2+2);

    // FFT
    FTransform::frc1d(nbin, (float*)&spec[0], profile->get_amps());

    unsigned ic, nc=nbin/2+1;
    float scale = 1.0/nbin;
    for (ic=0; ic < nc; ic++)
      spec[ic] *= complex<float>(0.0f, ic) // scale by i*frequency
	* scale; // correction for DFT normalisation

    // IFFT
    FTransform::bcr1d (nbin, &data[0], (float*)&spec[0]);
  }
  else // normal profile
  {
    Reference::To<PhaseWeight> baseline = profile->baseline();
    variance = baseline->get_variance().get_value();
    mean = baseline->get_mean().get_value();

    if (verbose)
      cerr << "Pulsar::ComponentModel::fit mean=" << mean
	   << " variance=" << variance << endl;

    // copy data
    for (unsigned i=0; i < nbin; i++)
      data[i] = profile->get_amps()[i] - mean;
  }

  Axis<double> argument; 
  model->set_argument (0, &argument);

  vector< Axis<double>::Value > xval;  
  vector< Estimate<double> > yval;

  unsigned nfree = 0;

  for (unsigned i=0; i < nbin; i++)
  {
    double radians = (i+0.5)/nbin * 2*M_PI;

    xval.push_back( argument.get_Value(radians) );
    yval.push_back( Estimate<double>(data[i], variance) );

    nfree ++;
  }
  
  for (unsigned i=0; i<model->get_nparam(); i++)
    if (model->get_infit(i))
    {
      if (nfree == 0)
	throw Error (InvalidState, "Pulsar::ComponentModel::fit",
		     "more free parameters than data");
      nfree --;
    }

  // fit

  LevenbergMarquardt<double> fit;
  fit.verbose = Function::verbose;
    
  chisq = fit.init (xval, yval, *model);

  // fit.singular_threshold = 1e-15; // dodge covariance problems
  unsigned iter = 1;
  unsigned not_improving = 0;

  string result = "RMS";
  if (variance)
    result = "chisq";

  unsigned max_iter = 100;

  while (not_improving < 25 && iter < max_iter)
  { 
    if (verbose)
      cerr << "iteration " << iter;

    float nchisq = fit.iter (xval, yval, *model);

    if (verbose)
      cerr << "     " << result << " = " << nchisq << endl;
    
    if (nchisq < chisq)
    {
      float diffchisq = chisq - nchisq;
      chisq = nchisq;
      not_improving = 0;
//.........这里部分代码省略.........