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

void psrspa::matched_finder ( const Archive* arch )
{
  string name = arch->get_filename ();
  for ( unsigned isub = 0; isub < arch->get_nsubint(); isub++ )
  {
    Reference::To<Profile> profile = arch->get_Profile(isub,0,0)->clone();

    while (profile->get_nbin() > 256)
    {
      finder->set_Profile( profile );  // might finder optimize on &profile?

      PhaseWeight weight;
      finder->get_weight( &weight );

      matched_report (name, isub, weight, *profile);

      profile->bscrunch (2);
    }
  }
}

void Pulsar::SpectrumPlot::get_spectra (const Archive* data)
{
  unsigned nchan = data->get_nchan();

  spectra.resize(1);
  spectra[0].resize(nchan);

  Reference::To<ProfileStats> stats;
  Reference::To<TextInterface::Parser> parser;

  if (!expression.empty())
  {
    stats = new ProfileStats;
    parser = stats->get_interface ();
  }

  Reference::To<const Integration> subint;
  subint = get_Integration (data, isubint);

  for (unsigned ichan=0; ichan<nchan; ichan++)
  {
    Reference::To<const Profile> profile;
    profile = get_Profile (subint, ipol, ichan);
    if (profile -> get_weight() == 0.0)
      spectra[0][ichan] = 0.0;
    else if (stats)
    {
      stats->set_Profile (profile);
      string value = process( parser, expression );
      spectra[0][ichan] = fromstring<float>( value );
    }
    else if (ibin.get_integrate())
      spectra[0][ichan] = profile->sum() / (float)profile->get_nbin();
    else
      spectra[0][ichan] = profile->get_amps()[ibin.get_value()];
  }
}

void psrpca::finalize ()
{
  arrival->set_observation ( total );
  arrival->get_toas(toas);

  if ( remove_std_baseline )
    std_archive -> remove_baseline ();
  Reference::To<Profile> std_prof = std_archive->get_Profile(0, 0, 0);
  float *s_amps = std_prof->get_amps ();
  const float nbin = std_prof->get_nbin ();

  double scale, offset, snr;

  if ( total_count < nbin )
    cerr << "WARNING: psrpca::finalize - not enough observations provided, "
	    "covariance matrix will not have full rank" << endl;

  //total->remove_baseline();

  gsl_matrix *profiles = gsl_matrix_alloc ( (unsigned)nbin, total_count );

  for (unsigned i_subint = 0; i_subint < total->get_nsubint(); i_subint++ )
  {
    Reference::To<Profile> prof = total->get_Profile ( i_subint, 0, 0 );
    if ( apply_shift )
      prof->rotate_phase ( toas[i_subint].get_phase_shift() );
    snr = prof->snr ();

    //calculate the scale
    float *p_amps = prof->get_amps ();
    scale = 0.0;

    for ( unsigned i_bin = 0; i_bin < nbin; i_bin++ )
    {
      scale += s_amps[i_bin] * p_amps[i_bin];
    }
    scale = (prof->get_nbin()* scale - prof->sum() * std_prof->sum()) /
	    (prof->get_nbin()* std_prof->sumsq() - std_prof->sum() * std_prof->sum());

    // calculate the baseline offset
    offset = (scale * std_prof->sum() - prof->sum()) / nbin;

    if ( prof_to_std )
    {
      //match the profile to standard and subtract the standard
      if ( apply_offset )
	prof->offset  ( offset );
      if ( apply_scale )
	prof->scale ( 1.0/scale );
      prof->diff ( std_prof );

      double* damps;
      damps = new double [ (unsigned)nbin ];
      transform( prof->get_amps(), prof->get_amps() + (unsigned)nbin, damps, CastToDouble() );

      gsl_vector_const_view view = gsl_vector_const_view_array( damps, nbin );
      gsl_matrix_set_col ( profiles, i_subint, &view.vector );
      t_cov->add_Profile ( prof, snr );
    }
    else
    {// prof_to_std is false
      Reference::To<Profile> diff = prof->clone ();
      diff->set_amps ( std_prof->get_amps () );
      if ( apply_offset ) 
	diff->offset( -offset );
      if ( apply_scale ) 
	diff->scale (scale);
      diff->diff ( prof );
      diff->scale (-1);

      double* damps;
      damps = new double [ (unsigned)nbin ];
      transform( diff->get_amps(), diff->get_amps() + (unsigned)nbin, damps, CastToDouble() );

      gsl_vector_const_view view = gsl_vector_const_view_array( damps, nbin );
      gsl_matrix_set_col ( profiles, i_subint, &view.vector );

      t_cov->add_Profile ( diff, snr );
      prof->set_amps ( diff->get_amps() );
    }
  }

  covariance = gsl_matrix_alloc ( (int) nbin, (int) nbin );
  t_cov->get_covariance_matrix_gsl ( covariance );

  // write the covariance matrix and difference profiles
  FILE *out;
  if ( save_covariance_matrix )
  {
    out = fopen ( (prefix+"_covariance.dat").c_str(), "w" );
    gsl_matrix_fprintf(out, covariance, "%g");
    fclose ( out );
  } // save covariance matrix

  if ( save_diffs ) 
    total->unload ( prefix+"_diffs.ar" );

  //solve the eigenproblem
  gsl_matrix_view m = gsl_matrix_submatrix ( covariance, 0, 0, (int)nbin, (int)nbin );
  gsl_vector *eval = gsl_vector_alloc ( (int)nbin );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  while (iterate)
  {
    iterate = false;

    // fit if specified
    if (fit)
    {
      // set fit flags
      model.set_infit(fit_flags.c_str());
      // fit
      model.fit(archive->get_Integration(0)->get_Profile(0,0));
    }
 
    // plot
    if (!pgdev.empty())
    {
      Reference::To<Pulsar::Archive> scrunched;
      scrunched = archive->clone();
      if (bscrunch > 1)
	scrunched->bscrunch(bscrunch);

      cpgpage();

      Profile *prof = scrunched->get_Integration(0)->get_Profile(0,0);
      float ymin = prof->min();
      float ymax = prof->max();
      float extra = 0.05*(ymax-ymin);

      ymin -= extra;
      ymax += extra;
      cpgswin(0.0, 1.0, ymin, ymax);
      cpgbox("bcnst", 0, 0, "bcnst", 0, 0);
      cpglab("Pulse phase", "Intensity", "");
      unsigned i, npts=prof->get_nbin();
      cpgsci(14);
      for (ic=0; ic < model.get_ncomponents(); ic++)
	plot(model, npts, true, ic);
      vector<float> xvals(npts);
      cpgsci(4);
      for (i=0; i < npts; i++)
	xvals[i] = i/((double)npts);
      if (line_plot)
	cpgline(npts, &xvals[0], prof->get_amps());
      else
	cpgbin(npts, &xvals[0], prof->get_amps(), 0);
      cpgsci(2);
      plot_difference(model, prof, line_plot);
      cpgsci(1);
      plot(model, npts, true);
    }
    
    
    while (interactive)
    {
      iterate = true;
      fit = false;

      cerr << "Enter command ('h' for help)" << endl;

      float curs_x=0, curs_y=0;
      char key = ' ';
      if (cpgband(6,0,0,0,&curs_x, &curs_y, &key) != 1 || key == 'q')
      {
	cerr << "Quitting ..." << endl;
	iterate = false;
	break;

void Pulsar::PolnSpectrumStats::build () try
{
  if (!profile)
    return;

  fourier = fourier_transform (profile, plan);

  // convert to Stokes parameters and drop the Nyquist bin
  fourier->convert_state (Signal::Stokes);
  fourier->resize( profile->get_nbin() );

  // form the power spectral density
  Reference::To<PolnProfile> psd = fourier_to_psd (fourier);

  // separate fourier into real and imaginary components
  Reference::To<PolnProfile> re = psd->clone();
  Reference::To<PolnProfile> im = psd->clone();

  unsigned npol = 4;
  unsigned nbin = psd -> get_nbin();

  for (unsigned ipol=0; ipol < npol; ipol++)
  {
    float* C_ptr = fourier->get_Profile(ipol)->get_amps();
    float* re_ptr = re->get_Profile(ipol)->get_amps();
    float* im_ptr = im->get_Profile(ipol)->get_amps();

    for (unsigned ibin=0; ibin < nbin; ibin++)
    {
      re_ptr[ibin] = C_ptr[ibin*2];
      im_ptr[ibin] = C_ptr[ibin*2+1];
    }
  }

  if (!regions_set)
  {
    LastHarmonic last;
    last.set_Profile( psd->get_Profile(0) );

    last.get_weight (&onpulse);
    last.get_baseline_estimator()->get_weight (&baseline);

    last_harmonic = last.get_last_harmonic();

#ifdef _DEBUG
    cerr << "Pulsar::PolnSpectrumStats::build last harmonic=" 
	 << last_harmonic << " nbin on=" << onpulse.get_weight_sum() << endl;
#endif

    real->set_regions (onpulse, baseline);
    imag->set_regions (onpulse, baseline);
  }

  if (onpulse.get_nbin () != re->get_nbin())
  {
    PhaseWeight on_temp = onpulse;
    PhaseWeight off_temp = baseline;

    on_temp.resize( re->get_nbin() );
    off_temp.resize( re->get_nbin() );

    if (re->get_nbin() > onpulse.get_nbin ())
    {
      copy_pad( on_temp, onpulse, 0 );
      copy_pad( off_temp, baseline, 1 );
    }

    real->set_regions (on_temp, off_temp);
    imag->set_regions (on_temp, off_temp);
  }

  real->set_profile (re.release());
  imag->set_profile (im.release());

}
catch (Error& error)
{
  throw error += "Pulsar::PolnSpectrumStats::build";
}

int scan_pulses(Reference::To<Pulsar::Archive> arch, vector<pulse>& data, 
		int method, float cphs, float dcyc)
{
  /* Find and store the flux and phase of each pulse in the file
     to the data vector. */
  
  pulse newentry;

  Reference::To<Pulsar::Profile> prof;
  
  newentry.file = arch->get_filename();

  double nm, nv, vm;
  
  int nbin = arch->get_nbin();
  int bwid = int(float(nbin) * dcyc);
  int cbin = 0;

  for (unsigned i = 0; i < arch->get_nsubint(); i++) {

    newentry.intg = i;
    
    prof = arch->get_Profile(i, 0, 0);
    
    prof->stats(prof->find_min_phase(), &nm, &nv, &vm);
    newentry.err = sqrt(nv);

    switch (method) {
      
    case 0: // Method of total flux
      
      newentry.flx = prof->sum();
      newentry.phs = prof->find_max_phase();
      break;
      
    case 1: // Method of maximum amplitude

      if (dcyc == 0.0) {
	newentry.flx = prof->max();
	newentry.phs  = prof->find_max_phase();
      }
      else {
	cbin = int(prof->find_max_phase(dcyc) * float(nbin));
	newentry.flx = prof->sum(cbin - bwid/2, cbin + bwid/2);
	newentry.phs  = prof->find_max_phase(dcyc);
      }
      break;

    case 2: // User-defined phase centre
      
      cbin = int(float(nbin) * cphs);

      if (dcyc == 0.0) {
	newentry.flx = (prof->get_amps())[cbin];
	newentry.phs = cphs;
      }
      else {
	newentry.flx = prof->sum(cbin - bwid/2, cbin + bwid/2);
	newentry.phs = cphs;
      }
      
      break;
      
    default:
      cerr << "No phase selection method chosen!" << endl;
    }
    
    data.push_back(newentry);
  }
  
  return data.size();
}