/*! \file ECModels.h \brief Contains declaration of classes that implements models for long time-constant eddy currents. \author Jesper Andersson \version 1.0b, Oct., 2022. */ // Declarations of classes that implements a hirearchy // of models for long-time constant fields from eddy // currents induced by diffusion gradients. // // LongECModels.h // // Jesper Andersson, FMRIB Image Analysis Group // // Copyright (C) 2022 University of Oxford // #ifndef LongECModels_h #define LongECModels_h #include #include #include #include #include #include "armawrap/newmat.h" #include "newimage/newimageall.h" #include "miscmaths/miscmaths.h" #include "EddyHelperClasses.h" // Cumbersome forward declaration of nlohmann::json to avoid exposing it to the nvcc compiler // I _hope_ the macro below doesn't change in json.hpp #ifndef INCLUDE_NLOHMANN_JSON_HPP_ namespace nlohmann { class json; } #endif namespace EDDY { /****************************************************************//** * * \brief Translates an EC time constant and intercept to weights * * Given TR, an EDDY::MultiBandGroup object and parameters * consisting of time-constant (half-life in seconds), and * scale factor/intercept it will return a vector of weights * for the current and the following volume eddy current-induced * field. * This is a helper class for two of the models below. * ********************************************************************/ class EC_TimeAndIntercept2Weights { public: EC_TimeAndIntercept2Weights(double TR, const MultiBandGroups& mbg) EddyTry : _TR(TR), _mbg(mbg) { } EddyCatch /// Get weights for EC-field from previous volume, arma::colvec GetWeightsForPrevious(double Thl, // Half-life in seconds double Intcpt) const; // Intercept /// Get weights for EC field from current volume. arma::colvec GetWeightsForCurrent(double Thl, // Half-life in seconds double Intcpt) const; // Intercept /// Return the number of MB-groups unsigned int NGroups() const EddyTry { return(_mbg.NGroups()); } EddyCatch private: const double _TR; /// Repetition time const MultiBandGroups _mbg; /// Multi-band information }; /****************************************************************//** * * \brief This class translates an EC time constant to weights * * Given TR, an EDDY::MultiBandGroup object and a time-constant * (half-life in seconds) of the eddy currents it will return * a vector of weights for the previous and the current volume * eddy current-induced field. This is a helper class for two of * the models below. * * It is currently retired in favour of the model above. * ********************************************************************/ class EC_TimeConst2Weights { public: EC_TimeConst2Weights(double TR, const MultiBandGroups& mbg) EddyTry : _TR(TR), _mbg(mbg) { } EddyCatch /// Get weights for EC field from previous volume. Thl: Half life in seconds arma::colvec GetWeightsForPrevious(double Thl) const; /// Get weights for EC field from current volume. Thl: Half life in seconds arma::colvec GetWeightsForCurrent(double Thl) const; /// Return the number of MB-groups unsigned int NGroups() const EddyTry { return(_mbg.NGroups()); } EddyCatch private: const double _TR; /// Repetition time const MultiBandGroups _mbg; /// Multi-band information void calculate_weights(double Thl, // Half life in seconds arma::colvec& prev, // Weights for EC from previous volume arma::colvec& curr) const; // Weights for EC from current volume }; /****************************************************************//** * * \brief This class manages the parameters for the long time-constant * (time varying) eddy currents. * * This class manages the parameters for the long time-constant * (time varying) eddy currents. I have chosen to _not_ use the * strategy of a virtual base class, as I do with ECModel below, * and instead have if-else blocks in the methods. This is because * I want some of the methods to be friends of ECScan, and I don't * want to have to friend declare every new subclass. * * All the Long EC models are based on the same assumption that * the first few MB-groups of a volume experiences residual EC * from the previous volume, and also that those same first few * volumes does not yet experince the full EC-field from the * diffusion gradient of the current volume. * There are four concrete subclasses of the pure virtual * parent class. They differ in if they model the weights directly * or through a model with a time-constant. They also differ in * if they model the different volumes individually, or if they * assume that the effects are the same across all volumes. * ********************************************************************/ class ECScan; // Forward declaration to allow us to declare ECScan as parameter to functions class ECScanManager; // Forward declaration to allow us to declare ECScanManager as parameter to functions class LongECModel { public: LongECModel() EddyTry {} EddyCatch ~LongECModel() {} /// Returns the model virtual LongECModelType Model() const = 0; /// Returns true for models that estimate individual parameters for each volume virtual bool IsIndividual() const = 0; /// Returns true for models that estimate joint parameters for all volumes virtual bool IsJoint() const = 0; /// Returns true if model estimates weights directly virtual bool EstimatesWeights() const = 0; /// Returns true if model estimates time constant/constants virtual bool EstimatesTimeConst() const = 0; /// Returns the number of MB-groups virtual unsigned int NGroups() const = 0; /// Return number of derivatives for a single scan virtual int NDerivs() const = 0; /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. virtual double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const = 0; /// Gets suitable scale for numerical derivative for parameter i virtual double GetDerivScale(unsigned int i) const = 0; /// Sets the value of a single parameter (which has a derivative) for scan sc /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. virtual void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const = 0; /// Get weight for diagonal of Hessian virtual double GetLevenbergLambda(unsigned int s) const = 0; /// Set weight for diagonal of Hessian virtual void SetLevenbergLambda(unsigned int s, double lambda) = 0; /// Returns all the parameters of the model virtual std::vector GetParameters(const ECScanManager& sm) const = 0; /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any virtual arma::colvec GetParameters(const ECScan& sc, unsigned int s) const = 0; /// Updates all the parameters of the model. This is the method that should be used for setting. virtual void UpdateParameters(ECScanManager& sm, const std::vector& update) = 0; /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. virtual void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update) = 0; /// Return a json-struct report and (optionally) write a plain ascii text file virtual nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const = 0; }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class IndividualWeightsModel // // This class models the weights directly on a volume-by-volume // basis. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class IndividualWeightsModel : public LongECModel { public: IndividualWeightsModel(ECScanManager& sm); ~IndividualWeightsModel() {} /// Returns the model LongECModelType Model() const EddyTry { return(LongECModelType::Individual); } EddyCatch /// Returns true for models that estimate individual parameters for each volume bool IsIndividual() const EddyTry { return(true); } EddyCatch /// Returns true for models that estimate joint parameters for all volumes bool IsJoint() const EddyTry { return(!IsIndividual()); } EddyCatch /// Returns true if model estimates weights directly bool EstimatesWeights() const EddyTry { return(true); } EddyCatch /// Returns true if model estimates time constant/constants bool EstimatesTimeConst() const EddyTry { return(!EstimatesWeights()); } EddyCatch /// Returns the number of MB-groups unsigned int NGroups() const EddyTry { return(_ngroups); } EddyCatch /// Return number of derivatives for a single scan int NDerivs() const EddyTry { return(2*_ngroups); } EddyCatch /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const; /// Gets suitable scale for numerical derivative for parameter i double GetDerivScale(unsigned int i) const; /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. /// Sets the value of a single parameter (which has a derivative) for scan sc void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const; /// Get weight for diagonal of Hessian double GetLevenbergLambda(unsigned int s) const; /// Set weight for diagonal of Hessian void SetLevenbergLambda(unsigned int s, double lambda); /// Returns all the parameters of the model std::vector GetParameters(const ECScanManager& sm) const; /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any arma::colvec GetParameters(const ECScan& sc, unsigned int s) const; /// Updates all the parameters of the model. N.B. that "update" means that the update is added to the current parameters. void UpdateParameters(ECScanManager& sm, const std::vector& update); /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. /// N.B. that "update" means that the update is added to the current parameters. /// N.B. that the index s implicitly assumes ScanType st==ScanType::Any void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update); /// Return a json-struct report and (optionally) write a plain ascii text file nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const; private: unsigned int _nscans; unsigned int _ngroups; std::vector _prev; std::vector _curr; std::vector _ll; // Levenberg lambda bool model_and_scan_agree(const ECScan& sc, unsigned int s) const; }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class IndividualTimeConstantsModel // // This class models the time constants (half-life in seconds) // and intercept on a volume-by-volume basis. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class IndividualTimeConstantsModel : public LongECModel { public: IndividualTimeConstantsModel(ECScanManager& sm); ~IndividualTimeConstantsModel() {} /// Returns the model LongECModelType Model() const EddyTry { return(LongECModelType::IndividualTimeConstant); } EddyCatch /// Returns true for models that estimate individual parameters for each volume bool IsIndividual() const EddyTry { return(true); } EddyCatch /// Returns true for models that estimate joint parameters for all volumes bool IsJoint() const EddyTry { return(!IsIndividual()); } EddyCatch /// Returns true if model estimates weights directly bool EstimatesWeights() const EddyTry { return(false); } EddyCatch /// Returns true if model estimates time constant/constants bool EstimatesTimeConst() const EddyTry { return(!EstimatesWeights()); } EddyCatch /// Returns the number of MB-groups unsigned int NGroups() const EddyTry { return(_tai2w.NGroups()); } EddyCatch /// Return number of derivatives for a single scan int NDerivs() const EddyTry { return(4); } EddyCatch /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const; /// Gets suitable scale for numerical derivative for parameter i double GetDerivScale(unsigned int i) const; /// Sets the value of a single parameter (which has a derivative) for scan sc /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const; /// Get weight for diagonal of Hessian double GetLevenbergLambda(unsigned int s) const; /// Set weight for diagonal of Hessian void SetLevenbergLambda(unsigned int s, double lambda); /// Returns all the parameters of the model std::vector GetParameters(const ECScanManager& sm) const; /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any arma::colvec GetParameters(const ECScan& sc, unsigned int s) const; /// Updates all the parameters of the model. N.B. that "update" means that the update is added to the current parameters. void UpdateParameters(ECScanManager& sm, const std::vector& update); /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. /// N.B. that "update" means that the update is added to the current parameters. /// N.B. that the index s implicitly assumes ScanType st==ScanType::Any void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update); /// Return a json-struct report and (optionally) write a plain ascii text file nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const; private: unsigned int _nscans; EC_TimeAndIntercept2Weights _tai2w; std::vector _par; // Vector of [previous[TC Intercept] current[TC Intercept]]] std::vector _ll; // Levenberg lambda bool model_and_scan_agree(const ECScan& sc, unsigned int s) const; }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class JointWeightsModel // // This class models the weights jointly across all volumes. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class JointWeightsModel : public LongECModel { public: JointWeightsModel(ECScanManager& sm); ~JointWeightsModel() {} /// Returns the model LongECModelType Model() const EddyTry { return(LongECModelType::Joint); } EddyCatch /// Returns true for models that estimate individual parameters for each volume bool IsIndividual() const EddyTry { return(false); } EddyCatch /// Returns true for models that estimate joint parameters for all volumes bool IsJoint() const EddyTry { return(!IsIndividual()); } EddyCatch /// Returns true if model estimates weights directly bool EstimatesWeights() const EddyTry { return(true); } EddyCatch /// Returns true if model estimates time constant/constants bool EstimatesTimeConst() const EddyTry { return(!EstimatesWeights()); } EddyCatch /// Returns the number of MB-groups unsigned int NGroups() const EddyTry { return(_ngroups); } EddyCatch /// Return number of derivatives for a single scan int NDerivs() const EddyTry { return(2*_ngroups); } EddyCatch /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const; /// Gets suitable scale for numerical derivative for parameter i double GetDerivScale(unsigned int i) const; /// Sets the value of a single parameter (which has a derivative) for scan sc /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const; /// Get weight for diagonal of Hessian double GetLevenbergLambda(unsigned int s) const; /// Set weight for diagonal of Hessian void SetLevenbergLambda(unsigned int s, double lambda); /// Returns all the parameters of the model std::vector GetParameters(const ECScanManager& sm) const; /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any arma::colvec GetParameters(const ECScan& sc, unsigned int s) const; /// Updates all the parameters of the model. N.B. that "update" means that the update is added to the current parameters. void UpdateParameters(ECScanManager& sm, const std::vector& update); /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. /// N.B. that "update" means that the update is added to the current parameters. /// N.B. that the index s implicitly assumes ScanType st==ScanType::Any void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update); /// Return a json-struct report and (optionally) write a plain ascii text file nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const; private: unsigned int _nscans; unsigned int _ngroups; arma::colvec _prev; arma::colvec _curr; double _ll; // Levenberg lambda bool model_and_scan_agree(const ECScan& sc, unsigned int s) const; }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class JointTimeConstantModel // // This class models the time-constant (half-life in seconds) // jointly across all volumes. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class JointTimeConstantModel : public LongECModel { public: JointTimeConstantModel(ECScanManager& sm); ~JointTimeConstantModel() {} /// Returns the model LongECModelType Model() const EddyTry { return(LongECModelType::JointTimeConstant); } EddyCatch /// Returns true for models that estimate individual parameters for each volume bool IsIndividual() const EddyTry { return(false); } EddyCatch /// Returns true for models that estimate joint parameters for all volumes bool IsJoint() const EddyTry { return(!IsIndividual()); } EddyCatch /// Returns true if model estimates weights directly bool EstimatesWeights() const EddyTry { return(false); } EddyCatch /// Returns true if model estimates time constant/constants bool EstimatesTimeConst() const EddyTry { return(!EstimatesWeights()); } EddyCatch /// Returns the number of MB-groups unsigned int NGroups() const EddyTry { return(_tai2w.NGroups()); } EddyCatch /// Return number of derivatives for a single scan int NDerivs() const EddyTry { return(1); } EddyCatch /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const; /// Gets suitable scale for numerical derivative for parameter i double GetDerivScale(unsigned int i) const; /// Sets the value of a single parameter (which has a derivative) for scan sc /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const; /// Get weight for diagonal of Hessian double GetLevenbergLambda(unsigned int s) const; /// Set weight for diagonal of Hessian void SetLevenbergLambda(unsigned int s, double lambda); /// Returns all the parameters of the model std::vector GetParameters(const ECScanManager& sm) const; /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any arma::colvec GetParameters(const ECScan& sc, unsigned int s) const; /// Updates all the parameters of the model. N.B. that "update" means that the update is added to the current parameters. void UpdateParameters(ECScanManager& sm, const std::vector& update); /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. /// N.B. that "update" means that the update is added to the current parameters. /// N.B. that the index s implicitly assumes ScanType st==ScanType::Any void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update); /// Return a json-struct report and (optionally) write a plain ascii text file nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const; private: unsigned int _nscans; EC_TimeAndIntercept2Weights _tai2w; arma::colvec _par; // previous[TC Intercept] current[TC Intercept] double _ll; // Levenberg lambda bool model_and_scan_agree(const ECScan& sc, unsigned int s) const; }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class NoLongECModel // // This is a "null-model" for the cases where we don't want to // model long eddy currents. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class NoLongECModel : public LongECModel { public: NoLongECModel() {} ~NoLongECModel() {} /// Returns the model LongECModelType Model() const EddyTry { return(LongECModelType::None); } EddyCatch /// Returns true for models that estimate individual parameters for each volume bool IsIndividual() const EddyTry { return(false); } EddyCatch /// Returns true for models that estimate joint parameters for all volumes bool IsJoint() const EddyTry { return(false); } EddyCatch /// Returns true if model estimates weights directly bool EstimatesWeights() const EddyTry { return(false); } EddyCatch /// Returns true if model estimates time constant/constants bool EstimatesTimeConst() const EddyTry { return(false); } EddyCatch /// Returns the number of MB-groups unsigned int NGroups() const EddyTry { throw EddyException("NoLongECModel::NGroups: unknown number of groups"); } EddyCatch /// Return number of derivatives for a single scan int NDerivs() const EddyTry { return(0); } EddyCatch /// Gets the value of a single parameter (which has derivative) for scan sc. /// It is intended for use on temporary copies of the scan sc, so this function /// does not check for agreement between the scan and the model. double GetDerivParam(const EDDY::ECScan& sc, unsigned int i, unsigned int s) const EddyTry { throw EddyException("NoLongECModel::GetDerivParam: model has no parameters"); } EddyCatch /// Gets suitable scale for numerical derivative for parameter i double GetDerivScale(unsigned int i) const EddyTry { throw EddyException("NoLongECModel::GetDerivScale: model has no parameters"); } EddyCatch /// Sets the value of a single parameter (which has a derivative) for scan sc /// It is intended for use on temporary copies of the scan sc, so this function /// does not set the corresponding parameter in the model, only in the scan. void SetDerivParam(ECScan& sc, unsigned int i, unsigned int s, double val) const EddyTry { throw EddyException("NoLongECModel::SetDerivParam: model has no parameters"); } EddyCatch /// Get weight for diagonal of Hessian double GetLevenbergLambda(unsigned int s) const EddyTry { throw EddyException("NoLongECModel::GetLevenbergLambda: model has no lambda"); } EddyCatch /// Set weight for diagonal of Hessian virtual void SetLevenbergLambda(unsigned int s, double lambda) EddyTry { throw EddyException("NoLongECModel::SetLevenbergLambda: model has no lambda"); } EddyCatch /// Returns all the parameters of the model std::vector GetParameters(const ECScanManager& sm) const EddyTry { throw EddyException("NoLongECModel::GetParameters(sm): model has no parameters"); } EddyCatch /// Returns the parameters for a single volume /// N.B. that the index i implicitly assumes ScanType st==ScanType::Any arma::colvec GetParameters(const ECScan& sc, unsigned int s) const EddyTry { throw EddyException("NoLongECModel::GetParameters(sc,s): model has no parameters"); } EddyCatch /// Updates all the parameters of the model. N.B. that "update" means that the update is added to the current parameters. void UpdateParameters(ECScanManager& sm, const std::vector& update) EddyTry { throw EddyException("NoLongECModel::UpdateParameters(sm,update): model has no parameters"); } EddyCatch /// Updates the parameters for a single volume. This is the method to use when rolling back updates on individual volumes. /// N.B. that "update" means that the update is added to the current parameters. /// N.B. that the index s implicitly assumes ScanType st==ScanType::Any void UpdateParameters(ECScan& sc, unsigned int s, const arma::colvec& update) EddyTry { throw EddyException("NoLongECModel::UpdateParameters(sc,s,update): model has no parameters"); } EddyCatch nlohmann::json WriteReport(const ECScanManager& sm, const std::string& fname, bool write_old_style_file) const; }; } // End namespace EDDY #endif // End #ifndef LongECModels_h