indii/fmri/model/MultiSeedModelMk2.hpp

#ifndef INDII_FMRI_MODEL_MULTISEEDMODELMK2_HPP
#define INDII_FMRI_MODEL_MULTISEEDMODELMK2_HPP

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

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

#include <vector>

/**
 * Range of variables for a given seed.
 */
#define SEED_RANGE(i) ublas::range(N*N+N+i*SEED_SIZE,N*N+N+(i+1)*SEED_SIZE)

namespace indii {
  namespace fmri {
    namespace model {

/**
 * Multiple seed time-series %model combining
 * indii::fmri::neural::MultiSeedNeuralModelMk2 for %neural interactions,
 * indii::fmri::hemodynamic::LogNeuralBalloonModel for independent
 * %hemodynamic activity at each seed and
 * indii::fmri::hemodynamic::LogBOLDCalculator for predicted BOLD
 * signal 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::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( c_1,\ldots,c_N \big)^T\f$, efficacies of input on
 * seeds.
 * @li \f$\big( (\ln q_1,\ln v_1,\ln f_1,s_1,u_1),\ldots,(\ln q_N,\ln
 * v_N,\ln f_N,s_N,u_N) \big)^T\f$, for each seed, its %neural activity
 * and four %hemodynamic parameters.
 * @li \f$\big( b_1,\ldots,b_N \big)^T\f$, for each seed, its baseline
 * BOLD signal.
 */
class MultiSeedModelMk2 : public indii::ml::ode::DifferentialModel,
    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,
    /**
     * Number of state variables for each seed.
     */
    SEED_SIZE
  };

  /**
   * Construct new %model.
   *
   * @param N Number of seeds.
   * @param u Experimental stimulus, for example a
   * indii::fmri::experimental::BlockDesignModel.
   */
  MultiSeedModelMk2(const unsigned int N, indii::ml::ode::FunctionModel* u);

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

  /**
   * 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::MultiSeedNeuralModelMk2* 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);

  /**
   * Get system noise structure.
   *
   * @return System noise.
   */
  const indii::ml::aux::WienerProcess<double>& getSystemNoise();

  /**
   * Get measurement noise structure.
   *
   * @return Measurement noise.
   */
  const indii::ml::aux::GaussianPdf& getMeasurementNoise();

  /**
   * Get resampling noise structure.
   *
   * @return Resampling noise.
   */
  const indii::ml::aux::GaussianPdf& getResamplingNoise();

  /**
   * Set system noise structure.
   *
   * @param w System noise.
   */
  void setSystemNoise(const indii::ml::aux::WienerProcess<double>& 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::WienerProcess<double> 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::ParticleFilterModel
   */
  virtual indii::ml::aux::vector transition(const indii::ml::aux::vector& x,
      const double start, const double delta);

  /**
   * @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::model::ParticleSmootherModel
   */
  virtual void alphaPrecalculate(
      const indii::ml::aux::DiracMixturePdf& p_xtn_ytn, const double start,
      const double delta);

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

  /**
   * Range of \f$A\f$ in state vector.
   */
  const ublas::range A_RANGE;

  /**
   * Range of \f$\mathbf{c}\f$ in state vector.
   */
  const ublas::range C_RANGE;

  /**
   * Range of seed variables in state vector.
   */
  const ublas::range SEEDS_RANGE;

  /**
   * Range of baselines in state vector.
   */
  const ublas::range B_RANGE;

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

  /**
   * Size of state.
   */
  const unsigned int STATE_SIZE;

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

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

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

  /**
   * System noise.
   */
  indii::ml::aux::WienerProcess<double> w;

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

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

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

  /**
   * alpha() precalculations.
   */
  std::vector<indii::ml::aux::vector> transitions;

  /**
   * As transition(), but without adding noise.
   */
  virtual indii::ml::aux::vector noiselessTransition(
      const indii::ml::aux::vector& x, const double start, const double delta);
  /**
   * As measure(), but without adding noise.
   */
  virtual indii::ml::aux::vector noiselessMeasure(
      const indii::ml::aux::vector &x);

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

};

    }
  }
}

#endif

---