indii/fmri/model/MultiSeedModel.hpp

#ifndef INDII_FMRI_MODEL_MULTISEEDMODEL_HPP
#define INDII_FMRI_MODEL_MULTISEEDMODEL_HPP

#include "../neural/MultiSeedNeuralModel.hpp"
#include "../hemodynamic/LogNeuralBalloonModel.hpp"
#include "../hemodynamic/LogBOLDCalculator.hpp"

#include "indii/ml/aux/GaussianPdf.hpp"
#include "indii/ml/ode/DifferentialModel.hpp"
#include "indii/ml/filter/UnscentedKalmanFilterModel.hpp"
#include "indii/ml/filter/ParticleSmootherModel.hpp"
#include "indii/ml/ode/AdaptiveRungeKutta.hpp"

#include <vector>

namespace indii {
  namespace fmri {
    namespace model {

/**
 * Multiple seed time-series %model combining
 * indii::fmri::neural::MultiSeedNeuralModel for %neural interactions,
 * indii::fmri::hemodynamic::LogNeuralBalloonModel for independent
 * %hemodynamic activity at each seed and
 * indii::fmri::hemodynamic::LogBOLDCalculator for predicted
 * measurements.
 *
 * @author Lawrence Murray <lawrence@indii.org>
 * @version $Rev: 317 $
 * @date $Date: 2007-10-05 17:18:14 +0100 (Fri, 05 Oct 2007) $
 *
 * Parameters of the %model may be estimated using
 * indii::ml::filter::UnscentedKalmanFilter<double>,
 * indii::ml::filter::ParticleFilter<double> or
 * indii::ml::filter::ParticleSmoother<double>. When retrieving the
 * state at any time from one of these methods, it will be arranged as
 * follows (for \f$N\f$ seeds):
 *
 * @li \f$\big( A_{1,1},A_{1,2},\ldots,A_{N,N-1},A_{N,N} \big)^T\f$,
 * %neural efficacies between seeds, arranged column-wise.
 *
 * @li \f$\big( (\ln q_1,\ln v_1,\ln f_1,s_1,u_1,b_1),\ldots,(\ln
 * q_N,\ln v_N,\ln f_N,s_N,u_N,b_N) \big)^T\f$, for each seed, its
 * %neural activity, four %hemodynamic parameters and baseline BOLD
 * signal.
 */
class MultiSeedModel :
    public indii::ml::ode::DifferentialModel,
    virtual public indii::ml::filter::UnscentedKalmanFilterModel<double>,
    virtual public indii::ml::filter::ParticleSmootherModel<double> {
public:
  /**
   * Order of state variables within each seed block of the state vector.
   */
  enum SeedVariables {
    LN_Q,
    LN_V,
    LN_F,
    S,
    U,
    B,
    /**
     * Number of state variables for each seed.
     */
    SEED_SIZE
  };

  /**
   * Construct new %model.
   *
   * @param N Number of seeds.
   */
  MultiSeedModel(unsigned int N);

  /**
   * Construct new %model.
   *
   * @param N Number of seeds.
   * @param w System noise.
   * @param v Measurement noise.
   * @param r Resampling noise.
   */
  MultiSeedModel(unsigned int N, const indii::ml::aux::GaussianPdf& w,
      const indii::ml::aux::GaussianPdf& v,
      const indii::ml::aux::GaussianPdf& r);

  /**
   * Destructor.
   */
  virtual ~MultiSeedModel();

  /**
   * Get state size.
   */
  virtual unsigned int getStateSize();

  /**
   * Get measurement size.
   */
  virtual unsigned int getMeasurementSize();

  /**
   * Get the %neural %model. Settings of this %model may be adjusted.
   *
   * @return The %neural %model.
   */
  indii::fmri::neural::MultiSeedNeuralModel* getNeuralModel();

  /**
   * Get the %hemodynamic %model for a particular seed. Settings of this
   * %model may be adjusted.
   *
   * @param seed The seed.
   *
   * @return The %hemodynamic %model for the given seed.
   */
  indii::fmri::hemodynamic::LogNeuralBalloonModel* getHemodynamicModel(
      const unsigned int seed);

  /**
   * Get the BOLD calculator for a particular seed. Settings of this
   * may be adjusted.
   *
   * @param seed The seed.
   *
   * @return The BOLD calculator for the given seed.
   */
  indii::fmri::hemodynamic::LogBOLDCalculator* getBOLDCalculator(
      const unsigned int seed);

  /**
   * Set system noise structure.
   *
   * @param w System noise.
   */
  void setSystemNoise(const indii::ml::aux::GaussianPdf& w);

  /**
   * Set measurement noise structure.
   *
   * @param v Measurement noise.
   */
  void setMeasurementNoise(const indii::ml::aux::GaussianPdf& v);

  /**
   * Set resampling noise structure.
   *
   * @param r Resampling noise.
   */
  void setResamplingNoise(const indii::ml::aux::GaussianPdf& r);

  /**
   * Calculate a suggested prior over the initial state of the system.
   * 
   * @return Suggested prior over the initial state of the system.
   */
  indii::ml::aux::GaussianPdf suggestPrior();

  /**
   * Calculate a suggested system noise structure.
   * 
   * @return Suggested system noise.
   */
  indii::ml::aux::GaussianPdf suggestSystemNoise();

  /**
   * Calculate a suggested measurement noise structure.
   * 
   * @return Suggested measurement noise.
   */
  indii::ml::aux::GaussianPdf suggestMeasurementNoise();

  /**
   * Calculate a suggested resampling noise structure.
   *
   * @return Suggested resampling noise.
   */
  indii::ml::aux::GaussianPdf suggestResamplingNoise();

  /**
   * @see indii::ml::model::UnscentedKalmanFilterModel
   */
  virtual indii::ml::aux::vector transition(const indii::ml::aux::vector &x,
      const indii::ml::aux::vector &w, double delta);

  /**
   * @see indii::ml::model::ParticleFilterModel
   */
  virtual indii::ml::aux::vector transition(const indii::ml::aux::vector& x,
      const double start, const double delta);

  /**
   * @see indii::ml::model::UnscentedKalmanFilterModel
   */
  virtual indii::ml::aux::vector measure(const indii::ml::aux::vector &x,
      const indii::ml::aux::vector &v);

  /**
   * @see indii::ml::model::ParticleFilterModel
   */
  virtual indii::ml::aux::vector measure(const indii::ml::aux::vector &x);

  /**
   * @see indii::ml::model::ParticleFilterModel
   */
  virtual double weight(const indii::ml::aux::vector& x,
      const indii::ml::aux::vector& y);

  /**
   * @see indii::ml::model::ParticleFilterModel
   */
  virtual indii::ml::aux::vector resample(const indii::ml::aux::vector& x);

  /**
   * @see indii::ml::model::ParticleSmootherModel
   */
  virtual indii::ml::aux::sparse_matrix alpha(
      const indii::ml::aux::DiracMixturePdf& p_xtn_ytn,
      const indii::ml::aux::DiracMixturePdf& p_xtnp1_ytnp1,
      const double start, const double delta);

  /**
   * @see indii::ml::ode::DifferentialModel
   */
  virtual void calculateDerivatives(double t, const double y[], double dydt[]);

private:
  /**
   * Number of seeds.
   */
  unsigned int N;

  /**
   * System noise.
   */
  indii::ml::aux::GaussianPdf w;

  /**
   * Measurement noise.
   */
  indii::ml::aux::GaussianPdf v;

  /**
   * Resampling noise.
   */
  indii::ml::aux::GaussianPdf r;

  /**
   * Neural model.
   */
  indii::fmri::neural::MultiSeedNeuralModel* neural;

  /**
   * Hemodynamic models.
   */
  std::vector<indii::fmri::hemodynamic::LogNeuralBalloonModel*> hemodynamics;

  /**
   * BOLD calculators.
   */
  std::vector<indii::fmri::hemodynamic::LogBOLDCalculator*> bolds;

  /**
   * Runge--Kutta.
   */
  indii::ml::ode::AdaptiveRungeKutta ode;

  /**
   * Initialise models.
   */
  void init();

  /**
   * Neural function for passing values from neural model to
   * hemodynamic models.
   */
  static double U_PASSER(double t, const double y[], void *o);

};

    }
  }
}

#endif

---