//! 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;
}    

// Calculate SNR to fill in to plot array
void Pulsar::DynamicSNSpectrumPlot::get_plot_array( const Archive *data, 
    float *array )
{

  StandardSNR standard_snr;
  Reference::To<Archive> data_std = data->total();
  standard_snr.set_standard( data_std->get_Profile(0,0,0) );

  std::pair<unsigned,unsigned> srange = get_subint_range (data);
  std::pair<unsigned,unsigned> crange = get_chan_range (data);

  int ii = 0;

  for( int ichan = crange.first; ichan < crange.second; ichan++) 
  {
    for( int isub = srange.first; isub < srange.second; isub++ )
    {
      Reference::To<const Profile> prof = 
        data->get_Profile ( isub, pol, ichan);

      float snr;

      if (prof->get_weight()==0.0)
        snr = 0.0;
      else
        snr = standard_snr.get_snr(prof);

      array[ii] = snr;
      ii++;

    }
  }
}

void paz::setup ()
{
  if (!killfile.empty ())
  {
    FILE *fptr = fopen (killfile.c_str (), "r");

    if (!fptr)
      throw Error (FailedSys, "paz::setup", "fopen " + killfile);

    char *buffer = new char[4096];
    while (fgets (buffer, 4096, fptr))
    {
      char* key = strtok (buffer, whitespace);
      while (key)
      {
	chans_to_zero.push_back (fromstring<unsigned>(key));
	key = strtok (NULL, whitespace);
      }
    }
  }

  zero_channels = chans_to_zero.size()
    || freqs_to_zero.size() || freq_ranges_to_zero.size();

  bool zero_subints = subs_to_zero.size ();

  if (zero_intersection && !(zero_channels && zero_subints))
    throw Error (InvalidState, "paz::setup",
		 "must use -I with both -s|S|k *and* -w|W");
  
  if (mow_all_subints || subints_to_mow.size())
    mower = new Pulsar::LawnMower;
  
  if (median_zap || median_zap_bybin || median_zap_window)
  {
    median_zapper = new Pulsar::ChannelZapMedian;
    median_zapper->set_bybin (median_zap_bybin);
    median_zapper->set_paz_report (true);

    if (median_zap_window)
      median_zapper->set_window_size (median_zap_window);
  }

  if (edge_zap_percent < 0.0 || edge_zap_percent >= 100.0)
    throw Error (InvalidState, "paz::setup",
		 "invalid parameter to option -E");

  if (!std_filename.empty())
  {
    Reference::To<Pulsar::Archive> data = Pulsar::Archive::load (std_filename);
    data->pscrunch ();
    data->fscrunch ();
    data->tscrunch ();
    thestd = data->get_Profile (0, 0, 0);
    
    standard_snr.set_standard (thestd);
    Pulsar::Profile::snr_strategy.get_value().set (&standard_snr,
						   &Pulsar::StandardSNR::get_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();
}

std::string Pulsar::ArrivalTime::get_value (const std::string& key,
                                            const Tempo::toa& toa) try
{
  if(key == "parang") return get_parang( observation );
  else if(key == "tsub") return get_tsub( observation );
  else if(key == "observer") return get_observer( observation );
  else if(key == "projid") return get_projid( observation );
  else if(key == "rcvr") return get_rcvr( observation );
  else if(key == "backend") return get_backend( observation );
  else if(key == "period") return get_period_as_string( observation );
  else if(key == "be_delay") return get_be_delay( observation );
  else if(key == "subint") return tostring(toa_subint);
  else if(key == "chan") return tostring(toa_chan);
  else if(key == "gof") return tostring( toa.get_reduced_chisq() );

  else
  {
    Reference::To<TextInterface::Parser> interface;
    interface = standard_interface (const_cast<Archive*>(observation.get()));

    interface->set_prefix_name (false);
    interface->set_indentation ("");

    interface->process( "subint=" + tostring(toa_subint) );
    interface->process( "chan=" + tostring(toa_chan) );

    return process( interface, key );
  }
} 
catch (Error& e)
{
  return "*error*";
}

void Pulsar::Archive::set_model (const Predictor* new_model, bool apply)
{
  if (!new_model)
  {
    model = 0;
    return;
  }

  if (!good_model (new_model))
    throw Error (InvalidParam, "Pulsar::Archive::set_model",
                 "supplied model does not span Integrations");

  // swap the old with the new
  Reference::To<Predictor> oldmodel = model;
  model = new_model->clone();

  if (!apply)
    return;

  if (verbose == 3)
    cerr << "Pulsar::Archive::set_model apply the new model" << endl;

  // correct Integrations
  for (unsigned isub = 0; isub < get_nsubint(); isub++)
    apply_model (get_Integration(isub), oldmodel.ptr());
}

  void setup (Integration* subint, unsigned ichan)
  {
    Pulsar::MoreProfiles* more = NULL;

    unsigned npol = subint->get_npol();
    if (auxiliary)
    {
      more = subint->get_Profile (0,ichan)->get<Pulsar::MoreProfiles>();
      npol = more->get_size();
    }

    try {
      vector<unsigned> pols = ::indeces (npol, indeces);
      profiles.resize (pols.size());

      for (unsigned ipol=0; ipol < pols.size(); ipol++)
      {
	if (more)
	  profiles[ipol] = more->get_Profile( pols[ipol] );
	else
	  profiles[ipol] = subint->get_Profile ( pols[ipol],ichan);
      }
    } 
    catch (...)
    {
      // if parsing indeces fails, then assume that letter codes are used
      Reference::To<const PolnProfile> profile = new_Stokes (subint, ichan);
      profiles.resize (indeces.length());
      for (unsigned ipol=0; ipol < indeces.length(); ipol++)
	profiles[ipol] = new_Profile (profile, indeces[ipol]);
    } 
  }

//! Initialize and resize the output before calling unpack
void dsp::OneBitCorrection::transformation ()
{
  if (input->get_nbit() != 1)
    throw Error(InvalidState,"dsp::OneBitCorrection::transformation()",
		"input not 1-bit digitized");

  if (verbose)
    cerr << "Inside dsp::OneBitCorrection::transformation" << endl;

#if FIX_THIS
  if( input->has<PMDAQ_Extension>() ){
    set_first_chan( input->get<PMDAQ_Extension>()->get_chan_begin() );
    set_end_chan( input->get<PMDAQ_Extension>()->get_chan_end() );
  }
#endif

  unsigned end_channel = min(input->get_nchan(),end_chan);
  unsigned required_nchan = end_channel - first_chan;

  // Make sure output nchan doesn't change if it's already been set:
  {
    Reference::To<Observation> dummy = new Observation(*input);
    dummy->set_nchan( required_nchan );
    output->Observation::operator=(*dummy);
  }

  // output will contain floating point values
  output->set_nbit (8 * sizeof(float));
  
  // resize the output 
  output->resize (input->get_ndat());

  // unpack the data
  unpack ();
}

//! Return the system response as determined by the SingleAxis
::MEAL::Complex2*
Pulsar::SingleAxisCalibrator::solve (const vector<Estimate<double> >& source,
				     const vector<Estimate<double> >& sky)
{
  Reference::To<Calibration::SingleAxis> model = new Calibration::SingleAxis;

  if ( !source_set || source.size() != 4 )
  {
    if (verbose > 2)
      cerr << "Pulsar::SingleAxisCalibrator::solve" << endl;
    model->solve (source);
    return model.release();
  }

  if (verbose > 2)
    cerr << "Pulsar::SingleAxisCalibrator::solve reference source=" 
         << reference_source << endl;

  // Convert the coherency vectors into Stokes parameters.
  Stokes< Estimate<double> > stokes_source = coherency( convert (source) );

  if (!solver)
    solver = new Calibration::SingleAxisSolver;

  solver->set_input (reference_source);
  solver->set_output (stokes_source);
  solver->solve (model);

  return model.release();
}

void psredit::process (Pulsar::Archive* archive)
{
  Reference::To<TextInterface::Parser> interface = archive->get_interface();

  if (commands.size() == 0)
  {
    cout << interface->help (true) << endl;;
    return;
  }

  // so that a space precedes each parameter processed
  interface->set_indentation (" ");
  interface->set_prefix_name (prefix_name);

  if (output_filename)
    cout << archive->get_filename();

  for (unsigned j = 0; j < commands.size(); j++)
  {
    if (verbose)
      cerr << "psredit: processing command '" << commands[j] << "'" << endl;

    cout << interface->process (commands[j]);
  }

  cout << endl;
}

void DigitiserCountsPlot::CheckCounts( const Archive *const_data )
{
  Archive *data = const_cast<Archive*>( const_data );

  Reference::To<DigitiserCounts> counts = data->get<DigitiserCounts>();

  if( !counts )
  {
    if( verbose > 1 )
      cerr << "Attempted to plot DigitiserCounts on an archive without a DigitiserCounts countsension" << endl;
    return;
  }

  int npthist = counts->get_npthist();
  int ndigr = counts->get_ndigr();
  int nlev = counts->get_nlev();

  if( npthist == 0 || nlev == 0 || ndigr == 0 )
  {
    if( verbose > 1 )
      cerr << "Attempted to plot DigitiserCounts on an archive without parameter (npthist,ndigr,nlev)" << endl;
    return;
  }
  
  valid_data = true;
}

void Pulsar::foreach (Integration* subint, const Integration* operand,
		      Reference::To< Combination<Profile> > xform)
{
  for (unsigned ipol = 0; ipol < subint->get_npol(); ipol++)
    for (unsigned ichan = 0; ichan < subint->get_nchan(); ichan++) {
      xform->set_operand( operand->get_Profile (ipol, ichan) );
      xform->transform( subint->get_Profile (ipol, ichan) );
    }
}

double Pulsar::FITSArchive::get_offs_sub( unsigned int isub ) const
{
  Reference::To<const Integration> integ = get_Integration( isub );
  
  if( !integ )
    return -1;
  
  return (integ->get_epoch() - reference_epoch).in_seconds();
}

void insert (Pulsar::Archive* archive, unsigned isub, const MJD& epoch)
{
  Reference::To<Pulsar::Integration> empty;
  empty = archive->get_Integration(0)->clone();
  empty->zero();
  empty->set_epoch( epoch );

  archive->expert()->insert (isub, empty);
}

void Pulsar::CalibratorSpectrum::prepare (const Archive* data)
{
  Reference::To<Archive> clone;

  if (plot_Ip && data->get_state() != Signal::Stokes)
  {
    clone = data -> clone();
    clone->convert_state(Signal::Stokes);
    data = clone;
  }

  vector< vector< Estimate<double> > > hi;
  vector< vector< Estimate<double> > > lo;

  unsigned nchan = data->get_nchan();
  unsigned npol = data->get_npol();

  Reference::To<const Integration> subint = get_Integration(data, isubint);
  ReferenceCalibrator::get_levels (subint, nchan, hi, lo);

  assert (hi.size() == npol);

  unsigned ipol, ipt, npt = hi[0].size();

  if (!plot_total)
    for (ipol=0; ipol<npol; ipol++)
      for (ipt=0; ipt<npt; ipt++)
	hi[ipol][ipt] -= lo[ipol][ipt];

  if (plot_low)
    for (ipol=0; ipol<npol; ipol++)
      for (ipt=0; ipt<npt; ipt++)
	hi[ipol][ipt] = lo[ipol][ipt];

  if (plot_Ip)
  {
    for (ipt=0; ipt<npt; ipt++)
      if (hi[0][ipt].get_variance() != 0)
	hi[1][ipt] = sqrt( sqr(hi[1][ipt]) +
			   sqr(hi[2][ipt]) +
			   sqr(hi[3][ipt]) );
    npol = 2;
  }

  double cfreq = data->get_centre_frequency();
  double bw = data->get_bandwidth();

  plotter.clear ();
  plotter.set_xrange (cfreq-0.5*bw, cfreq+0.5*bw);

  for (ipol=0; ipol<npol; ipol++)
    plotter.add_plot (hi[ipol]);

  get_frame()->get_y_scale()->set_minmax (plotter.get_y_min(), 
					  plotter.get_y_max());
}

Pulsar::Predictor*
load_polyco (fitsfile* fptr, double* pred_phs, bool verbose)
{
  if (verbose)
    cerr << "load_polyco entered" << endl;

  // Load the polyco from the FITS file

  int status = 0;
  fits_movnam_hdu (fptr, BINARY_TBL, "POLYCO", 0, &status);

  if (status) {
    if (verbose)
      cerr << "load_polyco no POLYCO HDU" << endl;
    return 0;
  }

  // ask for the number of rows in the binary table
  long nrows = 0;
  fits_get_num_rows (fptr, &nrows, &status);
  if (status != 0)
    throw FITSError (status, "load_polyco", "fits_get_num_rows");

  if (nrows == 0) {
    if (verbose)
      cerr << "load_polyco no rows in POLYCO HDU" << endl;
    return 0;
  }

  Reference::To<polyco> model = new polyco;
  load (fptr, model);

  if (verbose)
    cerr << "load_polyco loaded\n" << *model << endl;

  if (pred_phs) {

    // ask for the number of rows in the binary table
    long nrows = 0;
    fits_get_num_rows (fptr, &nrows, &status);
    if (status != 0)
      throw FITSError (status, "load_polyco", "fits_get_num_rows");
    
    long firstelem = 1;
    long onelement = 1;
    int colnum = 0;
    int anynul = 0;
    fits_get_colnum (fptr, CASEINSEN, "PRED_PHS", &colnum, &status);
    fits_read_col (fptr, TDOUBLE, colnum, nrows, firstelem, onelement,
		   NULL, pred_phs, &anynul, &status);

  }

  return model.release();
}

void psradd::append (Pulsar::Archive* archive) try
{
  if (phase_align)
  {
    Reference::To<Pulsar::Archive> standard;
    standard = total->total();
    Pulsar::Profile* std = standard->get_Profile(0,0,0);
    
    Reference::To<Pulsar::Archive> observation;
    observation = archive->total();
    Pulsar::Profile* obs = observation->get_Profile(0,0,0);
    
    archive->rotate_phase( obs->shift(std).get_value() );
  }

  if (archive->get_state() != total->get_state())
  {
    if (verbose)
      cerr << "psradd: converting state" 
	   << " from " << archive->get_state()
	   << " to " << total->get_state() << endl;
    
    archive->convert_state( total->get_state() );
  }
  
  if (patch)
  {
    if (verbose)
      cerr << "psradd: patching any missing sub-integrations" << endl;
    patch->operate (total, archive);
  }
  
  if (verbose)
    cerr << "psradd: appending " << archive->get_filename() << endl;

  if (time_direction)
    time.append (total, archive);
  else
    frequency.append (total, archive);

  if (log_file)
    fprintf (log_file, " %s", archive->get_filename().c_str());

  if (very_verbose)
    cerr << "psradd: after append, instance count = " 
	 << Reference::Able::get_instance_count() << endl;
}
catch (Error& error)
{
  cerr << "psradd: Archive::append exception:\n" << error << endl;
  if (auto_add)
    reset_total = true;
}
catch (...)
{
  cerr << "psradd: Archive::append exception thrown" << endl;
  if (auto_add)
    reset_total = true;
}

void psrflux::set_standard(Archive *arch)
{
  // Convert
  stdarch = arch->total();
  stdarch->convert_state(Signal::Intensity);

  // Set up DS calculation
  Reference::To<StandardFlux> flux = new StandardFlux;
  flux->set_fit_shift(true); // always init with true, choose later
  flux->set_standard(stdarch->get_Profile(0,0,0));
  ds.set_flux_method(flux);
}

//
// Function to select a standard profile from a set of filenames
//
Pulsar::Profile* Pulsar::find_standard (const Archive* data, 
					const string& path) {
  
  // Store the current working directory
  char* directory = new char[1024];
  getcwd(directory, 1024);
  
  if (chdir(path.c_str()) != 0) {
    chdir(directory);
    cerr << "No suitable standard profile found." << endl;
    return 0;
  }
  
  vector <string> the_stds;
  
  char temp[128];
  FILE* fptr = popen("ls -1 *.std", "r");
  if (ferror(fptr)==0) {
    while(fscanf(fptr, "%s\n", temp) == 1) {
      the_stds.push_back(temp);
    }
  }
  
  Reference::To <Archive> candidate = 0;
  
  for (unsigned i = 0; i < the_stds.size(); i++) {
    candidate = Archive::load("./" + the_stds[i]);
    candidate -> centre();
    if ((data -> get_source()) != (candidate -> get_source())) {
      candidate = 0;
      continue;
    }
    if ((data -> get_centre_frequency() != 
	 candidate -> get_centre_frequency())) {
      candidate = 0;
      continue;
    }
    candidate -> tscrunch();
    candidate -> pscrunch();
    candidate -> fscrunch();
    
    chdir(directory);
    
    cerr << "Selected standard profile: " << the_stds[i] << endl;
    
    return (candidate.release() -> get_Profile(0,0,0));
  }
  
  chdir(directory);
  cerr << "No standard profile available." << endl;
  return 0;
  
}

//! Create new Input based on command line options
dsp::Input* dsp::SingleThread::Config::open (int argc, char** argv)
{
  vector<string> filenames;

  if (command_line_header)
  {
    CommandLineHeader clh;
    filenames.push_back ( clh.convert(argc,argv) );
  }

  else 
  {
    for (int ai=optind; ai<argc; ai++)
      dirglob (&filenames, argv[ai]);
  }

  unsigned nfile = filenames.size();

  if (nfile == 0)
  {
    std::cerr << "please specify filename[s] (or -h for help)" << endl;
    exit (-1);
  }

  if (Operation::verbose)
  {
    if (nfile > 1)
    {
      std::cerr << "opening contiguous data files: " << endl;
      for (unsigned ii=0; ii < filenames.size(); ii++)
        std::cerr << "  " << filenames[ii] << endl;
    }
    else
      std::cerr << "opening data file " << filenames[0] << endl;
  }

  Reference::To<File> file;

  if (nfile == 1)
    file = dsp::File::create( filenames[0] );
  else
  {
    dsp::MultiFile* multi = new dsp::MultiFile;
    file = multi;

    if (force_contiguity)
      multi->force_contiguity();

    multi->open (filenames);
  }

  return file.release();
}