test6.cpp

Go to the documentation of this file.

#include "indii/fmri/hemodynamic/NeuralBalloonModel.hpp"
#include "indii/fmri/hemodynamic/FlowBalloonModel.hpp"
#include "indii/fmri/hemodynamic/BOLDCalculator.hpp"
#include "indii/ml/ode/AdaptiveRungeKutta.hpp"

#include <iostream>
#include <string.h>
#include <math.h>

using namespace std;
using namespace indii::fmri::hemodynamic;
using namespace indii::ml::ode;

/**
 * @file test6.cpp
 *
 * Test of NeuralBalloonModel in comparison to Friston et
 * al. (2000). Results are output on stdout, tab delimited, with
 * columns representing:
 *
 * \li \f$t\f$; time
 * \li \f$u(t)\f$
 * \li \f$y\f$; BOLD response
 *
 * Results are as follows:
 *
 * \image html test6.png "Results, c.f. Friston et al. (2000) Figure 1."
 * \image latex test6.eps "Results, c.f. Friston et al. (2000) Figure 1."
 *
 * @remarks The results given in Friston et al. (2000) are not from
 * experimentation with the Balloon model. This does serve as a useful
 * comparison of the Balloon model against empirical results, however,
 * if nothing else demonstrating the need for a sustained stimulus to
 * get meaningful results.
 */


/**
 * Duration of simulation (s).
 */
static const double END = 35.0;

/**
 * Impulse times
 */
static const double IMPULSES[] = { 5.0, 6.0 };
static const double IMPULSES_LENGTH = 2;

/**
 * \f$u(t)\f$; impulse of 1 at times 5s and 6s
 */
double U(double t, const double y[], void* o) {
  int i;
  for (i = 0; i < IMPULSES_LENGTH; i++) {
    if (t == IMPULSES[i]) {
      return 1.0;
    }
  }
  return 0.0;
}

/**
 * Run tests.
 */
int main(int argc, const char* argv[]) {
  NeuralBalloonModel balloonModel;
  BOLDCalculator boldCalculator(&balloonModel);
  AdaptiveRungeKutta ode(&balloonModel, balloonModel.suggestInitialState());

  /* set up balloon model */
  balloonModel.setFunction(NeuralBalloonModel::U, U);

  /* simulate in such a way that impulse times are not skipped */
  double t = 0.0;
  int i;
  for (i = 0; i < IMPULSES_LENGTH; i++) {
    while (t < IMPULSES[i]) {
      cout << t << "\t";
      cout << balloonModel.getFunction(NeuralBalloonModel::U) << "\t";
      cout << boldCalculator.calculate(ode.getState()) << "\t";
      cout << endl;
      t = ode.step(IMPULSES[i]);
    }
    ode.setDiscontinuity(); // discontinuity
  }
  while (t < END) {
    cout << t << "\t";
    cout << balloonModel.getFunction(NeuralBalloonModel::U) << "\t";
    cout << boldCalculator.calculate(ode.getState()) << "\t";
    cout << endl;
    t = ode.step(END);
  }

  return 0;
}

---