test9.cpp

Go to the documentation of this file.

#include "indii/fmri/hemodynamic/LogFlowBalloonModel.hpp"
#include "indii/fmri/hemodynamic/LogBOLDCalculator.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 test9.cpp
 *
 * Basic test of LogFlowBalloonModel.
 *
 * This test is a repeat of test2.cpp, but using the
 * LogFlowBalloonModel. Results are output on stdout, tab delimited,
 * with columns representing:
 *
 * \li \f$t\f$; time
 * \li \f$f_{in}(t)\f$
 * \li \f$f_{out}(t)\f$
 * \li \f$CMRO_2(t) = \frac{E(t)}{E_0}f_{in}(t)\f$ (Buxton et al. 2004).
 * \li \f$\ln(q)\f$
 * \li \f$\ln(v)\f$
 * \li \f$y\f$; BOLD response
 *
 * Results are as follows:
 *
 * \image html test9.png "Results, c.f. test2.cpp"
 * \image latex test9.eps "Results, c.f. test2.cpp"
 */

/**
 * \f$\tau_0\f$
 */
static const double TAU_0 = 2.0;

/**
 * \f$E_0\f$
 */
static const double E_0 = 0.4;

/**
 * \f$V_0\f$
 */
static const double V_0 = 0.01;

/**
 * \f$\alpha\f$
 */
static const double ALPHA = 0.5;

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

/**
 * \f$f_{in}(t)\f$; Trapezoidal function with a rise time of 4 s,
 * duration of 4 s and fall time of 4 s.
 */
double F_IN(double t, const double y[], void* o) {
  double f;
  if (t < 4) {
    f = (0.7/4.0)*t + 1;
  } else if (t < 8) {
    f = 1.7;
  } else if (t < 12) {
    f = (-0.7/4.0)*(t - 8) + 1.7;
  } else {
    f = 1.0;
  }
  return f;
}

double F_OUT(double t, const double y[], void* o) {
  double v = exp(y[LogFlowBalloonModel::LN_V]);

  return 0.7/8.0 * pow((v - 1.0)*2.0/0.3, 3.0) + 1.0;
}

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

  /* set up balloon model */
  balloonModel.setFunction(LogFlowBalloonModel::F_IN, F_IN);
  balloonModel.setFunction(LogFlowBalloonModel::F_OUT, F_OUT);
  balloonModel.setParameter(LogFlowBalloonModel::TAU_0, TAU_0);
  balloonModel.setParameter(LogFlowBalloonModel::E_0, E_0);
  balloonModel.setParameter(LogFlowBalloonModel::V_0, V_0);
  balloonModel.setParameter(LogFlowBalloonModel::ALPHA, ALPHA);

  /* set up BOLD calculator */
  double cmro2;

  double t = 0.0;
  while (t < END) {
    cmro2 = balloonModel.getFunction(LogFlowBalloonModel::E) /
        balloonModel.getParameter(LogFlowBalloonModel::E_0) *
        balloonModel.getFunction(LogFlowBalloonModel::F_IN);
    cout << t << "\t";
    cout << balloonModel.getFunction(LogFlowBalloonModel::F_IN) << "\t";
    cout << balloonModel.getFunction(LogFlowBalloonModel::F_OUT) << "\t";
    cout << cmro2 << "\t";
    cout << ode.getVariable(LogFlowBalloonModel::LN_Q) << "\t";
    cout << ode.getVariable(LogFlowBalloonModel::LN_V) << "\t";
    cout << boldCalculator.calculate(ode.getState()) * 100.0 << "\t";
    cout << endl;
    t = ode.step(END);
  }

  return 0;
}

---