#include <cstdlib>
#include <string>
#include <vector>
#include <limits>
#include "utils/options.h"
#include "newimage/newimage.h"
#include "EddyHelperClasses.h"

#ifndef EddyCommandLineOptions_h
#define EddyCommandLineOptions_h

namespace EDDY {

/****************************************************************//***
* \brief This is the exception that is being thrown by
*  EddyCommandLineOptions if there is a problem with the
*  user input.
*
*  This is the exception that is being thrown by
*  EddyCommandLineOptions if there is a problem with the
*  user input. It is different to EddyException because
*  I don't want it to generate a stack dump that might
*  confuse the user.
*
********************************************************************/
class EddyInputError: public std::exception
{
public:
  EddyInputError(const std::string& msg) noexcept : m_msg(std::string("EddyInputError:  ") + msg) {}
  ~EddyInputError() noexcept {}
  virtual const char * what() const noexcept { return(m_msg.c_str()); }
private:
  std::string m_msg;
};

class S2VParam {
public:
  S2VParam() EddyTry : _order(std::vector<int>(1,0)), _lambda(std::vector<float>(1,0.0)),
	       _fwhm(std::vector<float>(1,0.0)), _niter(std::vector<int>(1,0)) {} EddyCatch
  S2VParam(const std::vector<int>& order, const std::vector<float>& lambda,
	   const std::vector<float>& fwhm, const std::vector<int>& niter);
  unsigned int NOrder() const;
  unsigned int Order(unsigned int oi) const;
  std::vector<unsigned int> Order() const;
  unsigned int NIter(unsigned int oi) const;
  double Lambda(unsigned int oi) const;
  std::vector<double> Lambda() const;
  float FWHM(unsigned int oi, unsigned int iter) const;
  std::vector<float> FWHM(unsigned int oi) const;
private:
  std::vector<int>      _order;
  std::vector<float>    _lambda;
  std::vector<float>    _fwhm;
  std::vector<int>      _niter;
  unsigned int total_niter() const EddyTry {
    unsigned int rval=0;
    for (unsigned int i=0; i<_niter.size(); i++) rval += static_cast<unsigned int>(_niter[i]);
    return(rval);
  } EddyCatch
};

class DebugIndexClass {
public:
  DebugIndexClass(const std::string& in);
  DebugIndexClass(unsigned int frst, unsigned int lst) EddyTry { _indx.resize(lst-frst+1); for (unsigned int i=0; i<lst-frst+1; i++) _indx[i]=frst+i; } EddyCatch
  DebugIndexClass() : _indx(0) {}
  bool IsAmongIndicies(unsigned int indx) const EddyTry { for (unsigned int i=0; i<_indx.size(); i++) if (indx==_indx[i]) return(true); return(false); } EddyCatch
  unsigned int Min() const EddyTry {
    unsigned int rval=std::numeric_limits<unsigned int>::max(); for (unsigned int i=0; i<_indx.size(); i++) rval=min(rval,_indx[i]); return(rval);
  } EddyCatch
  unsigned int Max() const EddyTry {
    unsigned int rval=std::numeric_limits<unsigned int>::min(); for (unsigned int i=0; i<_indx.size(); i++) rval=max(rval,_indx[i]); return(rval);
  } EddyCatch
  const std::vector<unsigned int>& GetIndicies() const EddyTry { return(_indx); } EddyCatch
  void SetIndicies(const std::vector<unsigned int>& indx) EddyTry { _indx=indx; } EddyCatch
private:
  std::vector<unsigned int> _indx;

  std::vector<unsigned int> parse_commaseparated_numbers(const std::string& list) const;
  std::vector<std::string> parse_commaseparated_list(const std::string&  list) const;
  unsigned int min(unsigned int i, unsigned int j) const { if (i<=j) return(i); else return(j); }
  unsigned int max(unsigned int i, unsigned int j) const { if (i>=j) return(i); else return(j); }
};

class EddyCommandLineOptions {
public:
  EddyCommandLineOptions(int argc, char *argv[]);
  // These are utility functions used to enquire about the user input.
  // Below they are sorted  into functions relevant to both diffusion and fMRI,
  // only to diffusion or only to fMRI.
  //
  // Functions relevant to both diffusion and fMRI
  bool IsfMRI() const { return(_is_fmri); }
  bool IsDiffusion() const { return(!IsfMRI()); }
  // Names of input files
  std::string ImaFname() const EddyTry { return(_imain.value()); } EddyCatch
  std::string MaskFname() const EddyTry { return(_mask.value()); } EddyCatch
  std::string AcqpFname() const EddyTry { return(_acqp.value()); } EddyCatch
  std::string IndexFname() const EddyTry { return(_index.value()); } EddyCatch
  std::string TopupFname() const EddyTry { return(_topup.value()); } EddyCatch
  std::string FieldFname() const EddyTry { return(_field.value()); } EddyCatch
  // Names of output files
  std::string ParOutFname() const EddyTry { return(_out.value()+std::string(".eddy_parameters")); } EddyCatch
  std::string MovementOverTimeOutFname() const EddyTry { return(_out.value()+std::string(".eddy_movement_over_time")); } EddyCatch
  std::string IOutFname() const EddyTry { return(_out.value()); } EddyCatch
  std::string OutMaskFname() const EddyTry { return((_dont_mask_output.value()) ? _out.value()+std::string(".eddy_output_mask") : std::string(""));} EddyCatch
  std::string MoveBySuscFirstOrderFname() const EddyTry { return(_out.value()+std::string(".eddy_mbs_first_order_fields")); } EddyCatch
  std::string MoveBySuscSecondOrderFname() const EddyTry { return(_out.value()+std::string(".eddy_mbs_second_order_fields")); } EddyCatch
  std::string RMSOutFname() const EddyTry { return(_out.value()+std::string(".eddy_movement_rms")); } EddyCatch
  std::string DFieldOutFname() const EddyTry { return(_out.value()+std::string(".eddy_displacement_fields")); } EddyCatch
  std::string OLReportFname() const EddyTry { return(_out.value()+std::string(".eddy_outlier_report")); } EddyCatch
  std::string OLMapReportFname(bool problem=false) const EddyTry {
    if (problem) return(_out.value()+std::string(".eddy_outlier_map_written_when_all_volumes_in_a_shell_has_outliers"));
    else return(_out.value()+std::string(".eddy_outlier_map"));
  } EddyCatch
  std::string OLNStDevMapReportFname(bool problem=false) const EddyTry {
    if (problem) return(_out.value()+std::string(".eddy_outlier_n_stdev_map_written_when_all_volumes_in_a_shell_has_outliers"));
    else return(_out.value()+std::string(".eddy_outlier_n_stdev_map"));
  } EddyCatch
  std::string OLNSqrStDevMapReportFname(bool problem=false) const EddyTry {
    if (problem) return(_out.value()+std::string(".eddy_outlier_n_sqr_stdev_map_written_when_all_volumes_in_a_shell_has_outliers"));
    else return(_out.value()+std::string(".eddy_outlier_n_sqr_stdev_map"));
  } EddyCatch
  std::string OLFreeDataFname() const EddyTry { return(_out.value()+std::string(".eddy_outlier_free_data")); } EddyCatch
  std::string ResidualsOutFname() const EddyTry { return((_residuals.value()) ? _out.value()+std::string(".eddy_residuals") : std::string("")); } EddyCatch
  std::string ShellIndexOutFname() const EddyTry { return(_out.value()+std::string(".eddy_shell_indicies.json")); } EddyCatch
  std::string JsonOutFname() const EddyTry { return(_out.value()+std::string(".eddy.json")); } EddyCatch
  std::string AdditionalWithOutliersOutFname() const EddyTry { return(_out.value()+std::string(".eddy_results_with_outliers_retained")); } EddyCatch
  std::string PredictionsOutFname() const EddyTry { return(_out.value()+std::string(".eddy_predictions")); } EddyCatch
  std::string PredictionsInScanSpaceOutFname() const EddyTry { return(_out.value()+std::string(".eddy_predictions_in_scan_space")); } EddyCatch
  std::string LoggerFname() const EddyTry { return(_out.value()+std::string(".eddy_timing_log")); } EddyCatch
  std::string FieldMatFname() const EddyTry { return(_field_mat.value()); } EddyCatch
  bool SupressOldStyleTextFiles() const EddyTry { return(_no_text_files.value()); } EddyCatch;
  std::vector<unsigned int> Indicies() const EddyTry { return(_indvec); } EddyCatch
  std::vector<unsigned int> SessionIndicies() const EddyTry { return(_sessvec); } EddyCatch
  unsigned int Session(unsigned int i) const EddyTry { return((i < _sessvec.size()) ? _sessvec[i] : 0); } EddyCatch
  float RepetitionTime() const EddyTry { return(IsfMRI() ? _fmri._rep_time.value() : _diff._rep_time.value()); } EddyCatch
  bool RepetitionTimeIsSet() const EddyTry { return(IsfMRI() ? true : _diff._rep_time.set()); } EddyCatch
  float FWHM(unsigned int iter) const;
  std::vector<float> FWHM() const EddyTry { return(_fwhm); } EddyCatch
  unsigned int NIter() const { return(_niter); }
  void SetNIterAndFWHM(unsigned int niter, const std::vector<float>& fwhm);
  unsigned int S2V_NIter(unsigned int  oi) const EddyTry { return(_s2vparam.NIter(oi)); } EddyCatch
  void SetS2VParam(unsigned int order, float lambda, float fwhm, unsigned int niter);
  unsigned int Index(unsigned int i) const EddyTry { return(_indvec[i]); } EddyCatch
  EDDY::MultiBandGroups MultiBand() const;
  unsigned int NumOfNonZeroMovementModelOrder() const EddyTry { return(_s2vparam.NOrder()); } EddyCatch
  unsigned int MovementModelOrder(unsigned int oi) const EddyTry { return(_s2vparam.Order(oi)); } EddyCatch
  std::vector<unsigned int> MovementModelOrder() const EddyTry { return(_s2vparam.Order()); } EddyCatch
  double S2V_Lambda(unsigned int oi) const EddyTry { return(_s2vparam.Lambda(oi)); } EddyCatch
  std::vector<double> S2V_Lambda() const EddyTry { return(_s2vparam.Lambda()); } EddyCatch
  float S2V_FWHM(unsigned int oi, unsigned int iter) const EddyTry { return(_s2vparam.FWHM(oi,iter)); } EddyCatch
  std::vector<float> S2V_FWHM(unsigned int oi) const EddyTry { return(_s2vparam.FWHM(oi)); } EddyCatch
  bool IsSliceToVol() const EddyTry { return(NumOfNonZeroMovementModelOrder()>0); } EddyCatch
  bool EstimateMoveBySusc() const EddyTry { return(_estimate_mbs.value()); } EddyCatch
  unsigned int MoveBySuscNiter() const EddyTry { return(static_cast<unsigned int>(_mbs_niter.value())); } EddyCatch
  unsigned int N_MBS_Interleaves() const { return(3); }
  double MoveBySuscLambda() const EddyTry { return(static_cast<double>(_mbs_lambda.value())); } EddyCatch
  std::vector<unsigned int> MoveBySuscParam() const EddyTry { std::vector<unsigned int> rval(2); rval[0] = 3; rval[1] = 4; return(rval); } EddyCatch
  unsigned int MoveBySuscOrder() const { return(1); }
  bool MoveBySuscUseJacobian() const { return(true); }
  double MoveBySuscKsp() const EddyTry { return(static_cast<double>(_mbs_ksp.value())); } EddyCatch
  EDDY::ECModelType FirstLevelModel() const;
  EDDY::ECModelType b0_FirstLevelModel() const;
  bool ReplaceOutliers() const EddyTry { return(_rep_ol.value()); } EddyCatch
  bool WriteOutlierFreeData() const EddyTry { return(_rep_ol.value()); } EddyCatch
  const EDDY::OutlierDefinition& OLDef() const EddyTry { return(_oldef); } EddyCatch
  unsigned int OLErrorType() const EddyTry { return(_ol_ec.value()); } EddyCatch
  unsigned int RefScanNumber() const EddyTry { return(static_cast<unsigned int>(_ref_scan_no.value())); } EddyCatch
  EDDY::OLType OLType() const EddyTry { if (_ol_type.value()==std::string("sw")) return(OLType::SliceWise);
                                    else if (_ol_type.value()==std::string("gw")) return(OLType::GroupWise);
                                    else return(OLType::Both); } EddyCatch
  double OLUpperThresholdJacobianMask() const EddyTry { return(_ol_jacut.value()); } EddyCatch
  bool WriteMovementRMS() const { return(true); }
  bool WriteDisplacementFields() const EddyTry { return(_dfields.value()); } EddyCatch
  bool WriteResiduals() const EddyTry { return(_residuals.value()); } EddyCatch
  bool WriteAdditionalResultsWithOutliersRetained() const EddyTry { return(_with_outliers.value()); } EddyCatch
  bool WritePredictions() const EddyTry { return(_write_predictions.value()); } EddyCatch
  bool LogTimings() const EddyTry { return(_log_timings.value()); } EddyCatch
  bool History() const EddyTry { return(_history.value()); } EddyCatch
  bool FillEmptyPlanes() const EddyTry { return(_fep.value() ? true : !IsDiffusion() ? false : _diff._resamp.value() == std::string("lsr") ? true : false); } EddyCatch
  std::string InitFname() const EddyTry { return(_init.value()); } EddyCatch
  std::string InitS2VFname() const EddyTry { return(_init_s2v.value()); } EddyCatch
  std::string InitMBSFname() const EddyTry { return(_init_mbs.value()); } EddyCatch
  unsigned int NumberOfThreads() const EddyTry { return(static_cast<unsigned int>(_nthr.value())); } EddyCatch
  bool Verbose() const EddyTry { return(_verbose.value()); } EddyCatch
  bool VeryVerbose() const EddyTry { return(_very_verbose.value()); } EddyCatch
  bool WriteSliceStats() const EddyTry { return(_write_slice_stats.value()); } EddyCatch
  int DebugLevel() const { return(_debug); }
  const DebugIndexClass& DebugIndicies() const EddyTry { return(_dbg_indx); } EddyCatch
  void SetDebug(unsigned int level, const std::vector<unsigned int>& indx) EddyTry { _debug=level; _dbg_indx.SetIndicies(indx); } EddyCatch
  FinalResamplingType ResamplingMethod() const;
  CovarianceFunctionType CovarianceFunction() const;
  HyParCostFunctionType HyParCostFunction() const;
  bool MaskOutput() const EddyTry { return(!_dont_mask_output.value()); } EddyCatch
  int  InitRand() const { return(_init_rand); }
  unsigned int NVoxHp() const { return(_nvoxhp_internal); }
  void SetNVoxHp(unsigned int n) { _nvoxhp_internal = n; }
  double HyParFudgeFactor() const { return(_hypar_ff_internal); }
  void SetHyParFudgeFactor(double ff) { _hypar_ff_internal = ff; }
  bool HyperParFixed() const { return(_fixed_hpar); }
  void SetHyperParFixed(bool val=true) { _fixed_hpar = val; }
  NEWMAT::ColumnVector HyperParValues() const;
  void SetHyperParValues(const NEWMAT::ColumnVector& p);
  NEWIMAGE::interpolation InterpolationMethod() const;
  NEWIMAGE::extrapolation ExtrapolationMethod() const;
  bool ExtrapolationValidInPE() const EddyTry { return(_epvalid.value()); } EddyCatch
  NEWIMAGE::interpolation S2VInterpolationMethod() const;
  EDDY::PolationPara PolationParameters() const EddyTry {
    return(EDDY::PolationPara(InterpolationMethod(),ExtrapolationMethod(),ExtrapolationValidInPE(),S2VInterpolationMethod()));
  } EddyCatch
  bool ShapeReferencesSet() const EddyTry { return(_shape_references_set); } EddyCatch

  // Functions relevant only to diffusion
  std::string BVecsFname() const EddyTry { return(IsDiffusion() ? _diff._bvecs.value() : std::string("")); } EddyCatch
  std::string BValsFname() const EddyTry { return(IsDiffusion() ? _diff._bvals.value() : std::string("")); } EddyCatch
  std::string BDeltasFname() const EddyTry { return(IsDiffusion() ? _diff._bdeltas.value() : std::string("")); } EddyCatch
  std::string ECFOutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_fields") : std::string("")); } EddyCatch
  std::string LongECOutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_long_ec_parameters") : std::string("")); } EddyCatch
  std::string RotatedBVecsOutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_rotated_bvecs") : std::string("")); } EddyCatch
  std::string RotatedBVecsLSROutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_rotated_bvecs_for_SLR") : std::string("")); } EddyCatch
  std::string BVecsLSROutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_bvecs_for_SLR") : std::string("")); } EddyCatch
  std::string RestrictedRMSOutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_restricted_movement_rms") : std::string("")); } EddyCatch
  std::string DwiMssHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_dwi_mss_history") : std::string("")); } EddyCatch
  std::string DwiParHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_dwi_parameter_history") : std::string("")); } EddyCatch
  std::string B0MssHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_b0_mss_history") : std::string("")); } EddyCatch
  std::string B0ParHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_b0_parameter_history") : std::string("")); } EddyCatch
  std::string DwiMssS2VHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_slice_to_vol_dwi_mss_history") : std::string("")); } EddyCatch
  std::string DwiParS2VHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_slice_to_vol_dwi_parameter_history") : std::string("")); } EddyCatch
  std::string B0MssS2VHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_slice_to_vol_b0_mss_history") : std::string("")); } EddyCatch
  std::string B0ParS2VHistoryFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_slice_to_vol_b0_parameter_history") : std::string("")); } EddyCatch
  std::string PeasReportFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_post_eddy_shell_alignment_parameters") : std::string("")); } EddyCatch
  std::string PeasAlongPEReportFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_post_eddy_shell_PE_translation_parameters") : std::string("")); } EddyCatch
  std::string CNROutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_cnr_maps") : std::string("")); } EddyCatch
  std::string RangeCNROutFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_range_cnr_maps") : std::string("")); } EddyCatch
  std::string ScatterBrainPredictionsOutFname() const EddyTry { return(_out.value()+std::string(".eddy_scatter_brain_predictions")); } EddyCatch
  std::string VolumeScatterBrainPredictionsOutFname() const EddyTry { return(_out.value()+std::string(".eddy_volume_scatter_brain_predictions")); } EddyCatch
  std::string MIPrintFname() const EddyTry { return(IsDiffusion() ? _out.value()+std::string(".eddy_PEAS_MI_values") : std::string("")); } EddyCatch
  double BValueRange() const { return(IsDiffusion() ? static_cast<double>(_diff._brange.value()) : 0.0); }
  EDDY::SecondLevelECModelType SecondLevelModel() const { return(IsDiffusion() ? _diff._slm : SecondLevelECModelType::None); }
  EDDY::SecondLevelECModelType b0_SecondLevelModel() const { return(IsDiffusion() ? _diff._b0_slm : SecondLevelECModelType::None); }
  void SetSecondLevelModel(EDDY::SecondLevelECModelType slm) { if (IsDiffusion()) _diff._slm = slm; }
  void Set_b0_SecondLevelModel(EDDY::SecondLevelECModelType b0_slm) { if (IsDiffusion()) _diff._b0_slm = b0_slm; }
  bool HasSecondLevelModel() const EddyTry { return(SecondLevelModel() != EDDY::SecondLevelECModelType::None); } EddyCatch
  bool Has_b0_SecondLevelModel() const EddyTry { return(b0_SecondLevelModel() != EDDY::SecondLevelECModelType::None); } EddyCatch
  EDDY::LongECModelType LongECModel() const EddyTry { return(IsDiffusion() ? _diff._long_ec_model : EDDY::LongECModelType::None); } EddyCatch
  std::string LongECModelString() const;
  bool EstimateLongEC() const EddyTry { return(IsDiffusion() ? this->LongECModel() != EDDY::LongECModelType::None : false); } EddyCatch
  unsigned int LongECNIter() const EddyTry { return(IsDiffusion() ? 5 : 0); } EddyCatch
  bool ReestimateECWhenEstimatingLongEC() const EddyTry { return(IsDiffusion() ? !_diff._long_ec_dont_reest.value() : false); } EddyCatch 
  bool SeparateOffsetFromMovementWhenEstimatingLongEC() const EddyTry { return(IsDiffusion() ? _diff._long_ec_sep_offs_move.value() : false); } EddyCatch
  bool SeparateOffsetFromMovement() const EddyTry { return(IsDiffusion() ? !_diff._dont_sep_offs_move_internal : false); } EddyCatch
  void SetSeparateOffsetFromMovement(bool val=true) EddyTry { if (IsDiffusion()) _diff._dont_sep_offs_move_internal = val; } EddyCatch 
  EDDY::OffsetModelType OffsetModel() const;
  EDDY::OLSumStats OLSummaryStatistics() const;
  bool AddNoiseToReplacements() const EddyTry { return(IsDiffusion() ? _diff._rep_noise.value() : false); } EddyCatch
  bool RegisterDWI() const { return(IsDiffusion() ? _diff._rdwi : false); }
  bool Registerb0() const { return(IsDiffusion() ? _diff._rb0 : false); }
  bool WriteFields() const EddyTry { return(IsDiffusion() ? _diff._fields.value() : false); } EddyCatch
  bool WriteRotatedBVecs() const { return(IsDiffusion() ? true : false); }
  bool WriteCNRMaps() const EddyTry { return(IsDiffusion() ? _diff._write_cnr_maps.value() : false); } EddyCatch
  bool WriteRangeCNRMaps() const EddyTry { return(IsDiffusion() ? _diff._write_range_cnr_maps.value() : false); } EddyCatch
  bool WriteScatterBrainPredictions() const EddyTry { return(IsDiffusion() ? _diff._write_scatter_brain_predictions.value() : false); } EddyCatch
  bool WriteVolumeScatterBrainPredictions() const EddyTry { return(IsDiffusion() ? _diff._write_scatter_brain_predictions.value() : false); } EddyCatch
  bool RotateBVecsDuringEstimation() const EddyTry { return(IsDiffusion() ? _diff._rbvde_internal : false); } EddyCatch
  void SetRotateBVecsDuringEstimation(bool val=true) EddyTry { if (IsDiffusion()) _diff._rbvde_internal = val; } EddyCatch
  bool AlignShellsPostEddy() const EddyTry { return(IsDiffusion() ? !_diff._dont_peas.value() : false); } EddyCatch
  bool UseB0sToAlignShellsPostEddy() const EddyTry { return(IsDiffusion() ? _diff._use_b0s_for_peas.value() : false); } EddyCatch
  bool DontCheckShelling() const EddyTry { return(IsDiffusion() ? _diff._data_is_shelled.value() : true); } EddyCatch
  bool DoTestRot() const EddyTry { return(IsDiffusion() ? _diff._test_rot.set() : false); } EddyCatch
  bool PrintMIValues() const EddyTry { return(IsDiffusion() ? _diff._print_mi_values.value() : false); } EddyCatch
  bool PrintMIPlanes() const EddyTry { return(IsDiffusion() ? _diff._print_mi_planes.value() : false); } EddyCatch
  std::vector<float> TestRotAngles() const;
  double LSResamplingLambda() const EddyTry { return(IsDiffusion() ? _diff._lsr_lambda.value() : 0.0); } EddyCatch
  int B0ShapeReference() const EddyTry { return(IsDiffusion() ? _diff._b0_shape_reference : -1); } EddyCatch 
  std::pair<int,double> ShellShapeReference(unsigned int i) const;

  // Functions relevant only to fMRI
  std::string fMRIMssHistoryFname() const EddyTry { return(IsfMRI() ? _out.value()+std::string(".eddy_fmri_mss_history") : std::string("")); } EddyCatch
  std::string fMRIParHistoryFname() const EddyTry { return(IsfMRI() ? _out.value()+std::string(".eddy_fmri_parameter_history") : std::string("")); } EddyCatch
  std::string fMRIMssS2VHistoryFname() const EddyTry { return(IsfMRI() ? _out.value()+std::string(".eddy_slice_to_vol_fmri_mss_history") : std::string("")); } EddyCatch
  std::string fMRIParS2VHistoryFname() const EddyTry { return(IsfMRI() ? _out.value()+std::string(".eddy_slice_to_vol_fmri_parameter_history") : std::string("")); } EddyCatch
  float EchoTime() const EddyTry { return(IsfMRI() ? _fmri._echo_time.value() : 0.0f); } EddyCatch
  bool WriteSNRMaps() const EddyTry { return(IsfMRI() ? _fmri._write_snr_maps.value() : false); } EddyCatch
  int fMRIShapeReference() const EddyTry { return(IsfMRI() ? _fmri._shape_reference : -1); } EddyCatch
private:
  /// Struct that hold options specific to diffusion
  struct diff
  {
    diff();
    Utilities::Option<std::string>               _bvecs;
    Utilities::Option<std::string>               _bvals;
    Utilities::FmribOption<std::string>          _bdeltas;
    Utilities::Option<bool>                      _data_is_shelled;
    Utilities::Option<float>                     _brange;
    Utilities::Option<std::string>               _resamp;
    Utilities::FmribOption<float>                _lsr_lambda;
    Utilities::Option<std::string>               _flm;
    Utilities::Option<std::string>               _slm_str;
    Utilities::FmribOption<std::string>          _b0_flm;
    Utilities::FmribOption<std::string>          _b0_slm_str;
    Utilities::FmribOption<float>                _rep_time;
    Utilities::FmribOption<std::string>          _long_ec_str;
    Utilities::HiddenOption<bool>                _long_ec_dont_reest;
    Utilities::HiddenOption<bool>                _long_ec_sep_offs_move;
    Utilities::Option<bool>                      _rep_noise;
    Utilities::FmribOption<std::string>          _ol_ss;
    Utilities::Option<bool>                      _ol_pos;
    Utilities::Option<bool>                      _ol_sqr;
    Utilities::FmribOption<std::string>          _covfunc;
    Utilities::FmribOption<bool>                 _rbvde;
    Utilities::Option<bool>                      _dont_sep_offs_move;
    Utilities::FmribOption<std::string>          _offset_model;
    Utilities::Option<bool>                      _dont_peas;
    Utilities::FmribOption<bool>                 _use_b0s_for_peas;
    Utilities::Option<std::string>               _session;
    Utilities::FmribOption<bool>                 _fields;
    Utilities::Option<bool>                      _write_cnr_maps;
    Utilities::FmribOption<bool>                 _write_range_cnr_maps;
    Utilities::FmribOption<bool>                 _write_scatter_brain_predictions;
    Utilities::FmribOption<bool>                 _dwi_only;
    Utilities::FmribOption<bool>                 _b0_only;
    Utilities::HiddenOption<std::vector<float> > _test_rot;
    Utilities::HiddenOption<bool>                _print_mi_values;
    Utilities::HiddenOption<bool>                _print_mi_planes;
    Utilities::HiddenOption<int>                 _mb_offs;             // Deprecated
    Utilities::HiddenOption<bool>                _sep_offs_move;       // Deprecated
    Utilities::HiddenOption<bool>                _peas;                // Deprecated
    Utilities::HiddenOption<bool>                _rms;                 // Deprecated
    Utilities::HiddenOption<std::string>         _slorder;             // Deprecated
    EDDY::SecondLevelECModelType                 _slm;
    EDDY::SecondLevelECModelType                 _b0_slm;
    EDDY::LongECModelType                        _long_ec_model;
    bool                                         _rdwi;
    bool                                         _rb0;
    bool                                         _rbvde_internal;
    bool                                         _dont_sep_offs_move_internal;
    int                                          _b0_shape_reference;
    struct shell_shape_ref_struct {
      shell_shape_ref_struct(int indx, double b_val) : _indx(indx), _b_val(b_val) {} 
      int    _indx;
      double _b_val;
    };
    std::vector<shell_shape_ref_struct>          _shell_shape_references;
  };
  /// Struct that hold options specific to fmri
  struct fmri
  {
    fmri(); 
    Utilities::Option<std::string>               _session;
    Utilities::Option<float>                     _rep_time;
    Utilities::Option<float>                     _echo_time;
    Utilities::FmribOption<std::string>          _covfunc;
    Utilities::Option<bool>                      _write_snr_maps;
    Utilities::FmribOption<std::string>          _fast_map;
    int                                          _shape_reference;
  };
  std::string                                  _title;
  std::string                                  _examples;
  Utilities::Option<std::string>               _imain;
  Utilities::Option<std::string>               _mask;
  Utilities::Option<std::string>               _acqp;
  Utilities::Option<std::string>               _index;
  Utilities::Option<std::string>               _out;
  Utilities::Option<std::string>               _topup;
  Utilities::Option<std::string>               _field;
  Utilities::Option<std::string>               _field_mat;
  Utilities::Option<int>                       _niter_tmp;
  Utilities::Option<std::vector<float> >       _fwhm_tmp;
  Utilities::FmribOption<int>                  _ref_scan_no;
  Utilities::Option<std::vector<int> >         _shape_ref_scan_nos;
  Utilities::Option<std::string>               _interp;
  Utilities::FmribOption<std::string>          _extrap;
  Utilities::FmribOption<bool>                 _epvalid;
  Utilities::Option<bool>                      _rep_ol;
  Utilities::Option<float>                     _ol_nstd;
  Utilities::Option<int>                       _ol_nvox;
  Utilities::FmribOption<int>                  _ol_ec;
  Utilities::Option<std::string>               _ol_type;
  Utilities::FmribOption<float>                _ol_jacut;
  Utilities::Option<int>                       _mb;
  Utilities::Option<std::string>               _slspec;
  Utilities::Option<std::string>               _json;
  Utilities::Option<std::vector<int> >         _mp_order;
  Utilities::Option<std::vector<float> >       _s2v_lambda;
  Utilities::Option<std::vector<int> >         _s2v_niter;
  Utilities::FmribOption<std::vector<float> >  _s2v_fwhm;
  Utilities::Option<std::string>               _s2v_interp;
  Utilities::Option<bool>                      _estimate_mbs;
  Utilities::Option<int>                       _mbs_niter;
  Utilities::Option<float>                     _mbs_lambda;
  Utilities::Option<float>                     _mbs_ksp;
  Utilities::FmribOption<std::string>          _hyparcostfunc;
  Utilities::Option<int>                       _nvoxhp;
  Utilities::Option<int>                       _initrand;
  Utilities::Option<float>                     _hyparfudgefactor;
  Utilities::FmribOption<std::vector<float> >  _hypar;
  Utilities::Option<bool>                      _fep;
  Utilities::FmribOption<bool>                 _dont_mask_output;
  Utilities::FmribOption<std::string>          _init;
  Utilities::FmribOption<std::string>          _init_s2v;
  Utilities::FmribOption<std::string>          _init_mbs;
  Utilities::FmribOption<bool>                 _dfields;
  Utilities::Option<bool>                      _residuals;
  Utilities::FmribOption<bool>                 _with_outliers;
  Utilities::FmribOption<bool>                 _history;
  Utilities::FmribOption<bool>                 _write_slice_stats;
  Utilities::FmribOption<bool>                 _write_predictions;
  Utilities::Option<bool>                      _no_text_files;
  Utilities::FmribOption<int>                  _debug_tmp;
  Utilities::FmribOption<std::string>          _dbg_indx_str;
  Utilities::HiddenOption<bool>                _log_timings;
  Utilities::Option<int>                       _nthr;
  Utilities::FmribOption<bool>                 _very_verbose;
  Utilities::Option<bool>                      _verbose;
  Utilities::Option<bool>                      _help;
  diff                                         _diff; // Diffusion specific arguments
  fmri                                         _fmri; // fMRI specific arguments
  bool                                         _is_fmri;
  bool                                         _init_rand;
  bool                                         _fixed_hpar;
  double                                       _hypar_ff_internal;
  unsigned int                                 _nvoxhp_internal;
  std::vector<float>                           _hypar_internal;
  std::vector<unsigned int>                    _indvec;
  std::vector<unsigned int>                    _sessvec;
  unsigned int                                 _debug;
  DebugIndexClass                              _dbg_indx;
  EDDY::OutlierDefinition                      _oldef;
  unsigned int                                 _niter;
  std::vector<float>                           _fwhm;
  EDDY::S2VParam                               _s2vparam;
  bool                                         _shape_references_set;

  void do_initial_parsing(int argc, char *argv[]);
  bool indicies_kosher(NEWMAT::Matrix indx, NEWMAT::Matrix acqp);
  bool session_indicies_kosher(NEWMAT::Matrix indx);
  std::vector<unsigned int> get_slorder(const std::string& fname,
					unsigned int       ngrp) const;
  void set_shape_refs(const std::vector<int>&   shape_refs,
		      const std::vector<int>&   shells,
		      const std::vector<double> shell_b_vals);
};

} // End namespace EDDY

#endif // End #ifndef EddyCommandLineOptions_h