test7.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 "indii/ml/aux/vector.hpp"

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

namespace aux = indii::ml::aux;

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

/**
 * @file test7.cpp
 *
 * Test of jump features.
 *
 * This test uses a default NeuralBalloonModel with a driving neural activity
 * function \f$u(t)\f$ consisting of a 1 s long burst of activity followed by
 * no activity. The model is simulated for 60s, then reset to its state at the
 * 5s mark and simulated again for the remaining time. The idea is to ensure
 * that use of setState() and setTime() achieve consistent results with
 * running the model from time zero.
 *
 * Results are output into files \c results/test7.out and \c
 * results/test7_repeat.out, tab delimited, with columns representing:
 *
 * \li \f$t\f$; time
 * \li \f$u(t)\f$
 * \li \f$q\f$
 * \li \f$v\f$
 * \li \f$f\f$
 * \li \f$s\f$
 * \li \f$y\f$; BOLD response
 *
 * Results are as follows:
 *
 * \image html test7.png "Results"
 * \image latex test7.eps "Results"
 *
 * The original simulation is plotted with a dotted line in each
 * case. Starting from the 5s mark, the repeated simulation is plotted with a
 * solid line, which should exactly overlay the dotted line for the remainder
 * of the simulation.
 */

/**
 * Duration of simulation (s).
 */
static const double END = 60.0;
static const double JUMP = 5.0;

/**
 * \f$u(t)\f$; 1 for 1 s, then 0.
 */
double U(double t, const double y[], void* o) {
  double u;
  if (t < 1.0) {
    u = 1.0;
  } else {
    u = 0.0;
  }
  return u;
}

/**
 * Run tests.
 */
int main(int argc, const char* argv[]) {
  ofstream first("results/test7.out");
  ofstream repeat("results/test7_repeat.out");

  NeuralBalloonModel balloonModel;
  BOLDCalculator boldCalculator(&balloonModel);
  AdaptiveRungeKutta ode(&balloonModel, balloonModel.suggestInitialState());

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

  double t_jump;
  aux::vector y_jump(4);
  double t = 0.0;
  while (t < JUMP) {
    first << t << "\t";
    first << balloonModel.getFunction(NeuralBalloonModel::U) << "\t";
    first << ode.getVariable(FlowBalloonModel::Q) << "\t";
    first << ode.getVariable(FlowBalloonModel::V) << "\t";
    first << ode.getVariable(NeuralBalloonModel::F) << "\t";
    first << ode.getVariable(NeuralBalloonModel::S) << "\t";
    first << boldCalculator.calculate(ode.getState()) * 100.0 << "\t";
    first << endl;
    t = ode.step(JUMP);
  }
  t_jump = ode.getTime();
  y_jump = ode.getState();
  while (t < END) {
    first << t << "\t";
    first << balloonModel.getFunction(NeuralBalloonModel::U) << "\t";
    first << ode.getVariable(FlowBalloonModel::Q) << "\t";
    first << ode.getVariable(FlowBalloonModel::V) << "\t";
    first << ode.getVariable(NeuralBalloonModel::F) << "\t";
    first << ode.getVariable(NeuralBalloonModel::S) << "\t";
    first << boldCalculator.calculate(ode.getState()) * 100.0 << "\t";
    first << endl;
    t = ode.step(END);
  }
  ode.setTime(t_jump);
  ode.setState(y_jump);
  t = ode.getTime();
  while (t < END) {
    repeat << t << "\t";
    repeat << balloonModel.getFunction(NeuralBalloonModel::U) << "\t";
    repeat << ode.getVariable(FlowBalloonModel::Q) << "\t";
    repeat << ode.getVariable(FlowBalloonModel::V) << "\t";
    repeat << ode.getVariable(NeuralBalloonModel::F) << "\t";
    repeat << ode.getVariable(NeuralBalloonModel::S) << "\t";
    repeat << boldCalculator.calculate(ode.getState()) * 100.0 << "\t";
    repeat << endl;
    t = ode.step(END);
  }

  return 0;
}

---