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

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 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 );
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

  /* Process any args */
  int opt=0;
  int verb=0;

  string expression;

  while ((opt=getopt(argc,argv,"hve:E:"))!=-1) {
    switch (opt) {

      case 'v':
        verb++;
        Archive::set_verbosity(verb);
        break;
        
      case 'e':
        output_ext.assign(optarg);
        break;

      case 'E':
	expression = optarg;
	break;

      case 'h':
      default:
        usage();
        usage_interactive();
        exit(0);
        break;

    }
  }

  if (optind==argc) {
    usage();
    cerr << PROG ": No filename given" << endl;
    exit(-1);
  }

  /* Load file */
  string filename = argv[optind];
  Reference::To<Archive> orig_arch = Archive::load(filename);
  Reference::To<Archive> arch = orig_arch->clone();
  arch->dedisperse();
  arch->remove_baseline();
  double bw = arch->get_bandwidth();
  string output_filename = replace_extension(filename, output_ext);
  string psrsh_filename = replace_extension(filename,"psh");

  // Create profile plots
  Reference::To<Archive> tot_arch=NULL,pf_arch=NULL,pt_arch=NULL;
  ProfilePlot *totplot=NULL;
  PhaseVsFrequency *pfplot=NULL;
  PhaseVsTime *ptplot=NULL;
  int totplot_id;

  totplot = new ProfilePlot;
  pfplot = new PhaseVsFrequency;
  ptplot= new PhaseVsTime;

  // Create window and subdivide
  totplot_id = cpgopen("/xs");
  cpgpap(0.0,1.5);
  cpgsubp(1,3);

  // Create Dynamic Spectrum Plot
  DynamicBaselineSpectrumPlot *dsplot = new DynamicBaselineSpectrumPlot;
  if (!expression.empty())
    dsplot->configure("exp="+expression);
  else
    dsplot->configure("var=1");

  dsplot->set_reuse_baseline();
  int dsplot_id = cpgopen("/xs");
  if (dsplot_id<=0) {
    cerr << PROG ": PGPLOT xwindows device open failed, exiting." << endl;
    exit(1);
  }
  cpgask(0);

  /* Input loop */
  char ch='\0';
  enum cursor_type curs=both_cursor;
  float x0=0.0, y0=0.0, x1, y1;
  int click=0, mode=0;
  int pol=0,izap=0;
  bool redraw=true, var=true, log=false, resum=true,remove_baseline=false,method=false;
  struct zap_range zap;
  vector<struct zap_range> zap_list;
  do {
    /* Redraw the plot if necessary */
    if (redraw) {

      char conf[256];
      sprintf(conf, "above:c=$file\\n%s %s, %s scale, %s baseline, pol %d.", 
	      method ? "total" : "off-pulse" ,
	      var ? "variance" : "mean",
	      log ? "log" : "linear",
	      remove_baseline ? "variable" : "constant",
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  for (unsigned j = 0; j < arch->get_nchan(); j++) {
	    labels[j] = arch->get_Integration(i)->get_centre_frequency(j);
	  }
	  for (unsigned j = 0; j < arch->get_nchan(); j++) {
	    double new_frequency = labels[labels.size()-1-j];
	    arch->get_Integration(i)->set_centre_frequency(j,new_frequency);
	  }
	}
	arch->set_bandwidth(-1.0 * arch->get_bandwidth());
      }

      if (flip_freq) {
        for (unsigned isub = 0; isub < arch->get_nsubint(); isub++) {
          Reference::To<Pulsar::Integration> 
            subint = arch->get_Integration(isub);
          for (unsigned ichan = 0; ichan < arch->get_nchan(); ichan++) {
            double new_freq = flip_freq_mhz 
              - (subint->get_centre_frequency(ichan) - flip_freq_mhz);
            subint->set_centre_frequency(ichan, new_freq);
          }
        }
        arch->set_bandwidth(-1.0 * arch->get_bandwidth());
      }

      if( reverse_freqs ) {
	// Of course it would be nice to do this with pointers.... but oh well I guess copying will have to do HSK 27/8/04

	unsigned nchan = arch->get_nchan();

	for( unsigned isub=0; isub<arch->get_nsubint(); isub++){
	  for( unsigned ipol =0; ipol<arch->get_npol(); ipol++){
	    for( unsigned ichan=0; ichan<nchan/2; ichan++){
	      Reference::To<Pulsar::Profile> lo = arch->get_Profile(isub,ipol,ichan);	      
	      Reference::To<Pulsar::Profile> tmp = lo->clone();

	      Reference::To<Pulsar::Profile> hi = arch->get_Profile(isub,ipol,nchan-1-ichan);

	      lo->operator=(*hi);
	      hi->operator=(*tmp);
	    }
	  }
	}
	arch->set_bandwidth( -1.0 * arch->get_bandwidth() );
      }

      if (reset_weights) {
	arch->uniform_weight(new_weight);
	if (verbose)
	  cout << "All profile weights set to " << new_weight << endl;
      }
      
      if (rotate)
	arch->rotate_phase (rphase);

      if (scattered_power_correction) {

	Pulsar::ScatteredPowerCorrection spc;
	if (arch->get_state() == Signal::Stokes)
	  arch->convert_state(Signal::Coherence);

	spc.correct (arch);

      }

      if (newdm)
      {

//.........这里部分代码省略.........
  for (ic=0; ic < nc; ic++)
  {
    if (sscanf(new_components[ic].c_str(), 
	       "%lf %lf %lf %n", &centre, &concentration, &height,
	       &name_start)!=3)
    {
      cerr << "Could not parse component " << ic << endl;
      return -1;
    }

    model.add_component(centre, concentration, height, 
			new_components[ic].c_str()+name_start);
  }

  bool iterate = true;

  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);
    }

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 main(int argc, char* argv[]) try
{
  if (argc < 2) {
    usage();
    return EXIT_SUCCESS;
  }
  
  int gotc = 0;
  while ((gotc = getopt(argc, argv, "hvV")) != -1) {
    switch (gotc) {
    case 'h':
      usage();
      return EXIT_SUCCESS;
    case 'V':
      Pulsar::Archive::set_verbosity(3);
      break;
    case 'v':
      Pulsar::Archive::set_verbosity(2);
      break;
    }
  }
  
  if (optind >= argc)
    {
      cerr << "pazi: please specify filename" << endl;
      return -1;
    }
  
  string filename = argv[optind];
  string extension = filename.substr(filename.length() - 2, 2);
  
  if (extension == "rf")
    extension = "rz";
  else if (extension == "cf")
    extension = "cz";
  else
    extension = "pazi";
  
  string write_filename = filename + ".";
  write_filename += extension;
  
  cerr << "pazi: loading data" << endl;
  base_archive = Archive::load(filename);
  
  if (base_archive->get_npol() == 4)
  {
    original_state = base_archive->get_state();
    base_archive->convert_state( Signal::Stokes );
  }
  
  backup_archive = base_archive->clone();
  
  cerr << "pazi: making fscrunched clone" << endl;
  mod_archive = base_archive->clone();
  mod_archive->pscrunch();
  mod_archive->remove_baseline();
  mod_archive->dedisperse();
  mod_archive->fscrunch();

  scrunched_archive = mod_archive->clone();
  scrunched_archive->tscrunch();

  ranges.second = get_max_value(base_archive, plot_type);
  positive_direction = base_archive->get_bandwidth() < 0.0;

  time_orig_plot = factory.construct("time");
  time_mod_plot = factory.construct("time");
  time_fui = time_mod_plot->get_frame_interface();

  freq_orig_plot = factory.construct("freq");
  freq_mod_plot = factory.construct("freq");
  freq_fui = freq_mod_plot->get_frame_interface();

  total_plot = factory.construct("flux");
  total_plot->configure("info=1");

  subint_orig_plot = new ProfilePlot;
  subint_mod_plot = new ProfilePlot;
  subint_fui = subint_mod_plot->get_frame_interface();

  subint_orig_plot->configure("info=1");
  subint_mod_plot->configure("info=1");

  unsigned window = 0;
  char device [8];

  for (unsigned i=0; i<2; i++) do {
    window ++;
    snprintf (device, 8, "%u/XS", window);
  }
  while ( cpgopen (device) < 0 );

  cpgask(0);

  cerr << endl << "Total S/N = " <<
    scrunched_archive->get_Profile(0,0,0)->snr() << endl << endl;

  total_plot->plot(scrunched_archive);
  cpgslct(1);
  time_orig_plot->plot(mod_archive);
//.........这里部分代码省略.........