diff --git a/.cproject b/.cproject
deleted file mode 100644
index dc875a9..0000000
--- a/.cproject
+++ /dev/null
@@ -1,187 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/.project b/.project
deleted file mode 100644
index b0dd08a..0000000
--- a/.project
+++ /dev/null
@@ -1,83 +0,0 @@
-
-
- Steyn-Ross
-
-
-
-
-
- org.eclipse.cdt.managedbuilder.core.genmakebuilder
- clean,full,incremental,
-
-
- ?name?
-
-
-
- org.eclipse.cdt.make.core.append_environment
- true
-
-
- org.eclipse.cdt.make.core.autoBuildTarget
- all
-
-
- org.eclipse.cdt.make.core.buildArguments
-
-
-
- org.eclipse.cdt.make.core.buildCommand
- make
-
-
- org.eclipse.cdt.make.core.buildLocation
- ${workspace_loc:/Steyn-Ross/Debug}
-
-
- org.eclipse.cdt.make.core.cleanBuildTarget
- clean
-
-
- org.eclipse.cdt.make.core.contents
- org.eclipse.cdt.make.core.activeConfigSettings
-
-
- org.eclipse.cdt.make.core.enableAutoBuild
- false
-
-
- org.eclipse.cdt.make.core.enableCleanBuild
- true
-
-
- org.eclipse.cdt.make.core.enableFullBuild
- true
-
-
- org.eclipse.cdt.make.core.fullBuildTarget
- all
-
-
- org.eclipse.cdt.make.core.stopOnError
- true
-
-
- org.eclipse.cdt.make.core.useDefaultBuildCmd
- true
-
-
-
-
- org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder
- full,incremental,
-
-
-
-
-
- org.eclipse.cdt.core.cnature
- org.eclipse.cdt.core.ccnature
- org.eclipse.cdt.managedbuilder.core.managedBuildNature
- org.eclipse.cdt.managedbuilder.core.ScannerConfigNature
-
-
diff --git a/saves.h b/Data_Storage.h
similarity index 100%
rename from saves.h
rename to Data_Storage.h
diff --git a/NM_Thalamus.pro b/NM_Thalamus.pro
new file mode 100644
index 0000000..c1f904d
--- /dev/null
+++ b/NM_Thalamus.pro
@@ -0,0 +1,20 @@
+TEMPLATE = app
+CONFIG += console
+CONFIG -= app_bundle
+CONFIG -= qt
+
+TARGET = Thalamus.cpp
+
+SOURCES += Thalamic_Column.cpp \
+ Thalamus_mex.cpp \
+ Thalamus.cpp
+
+HEADERS += Data_Storage.h \
+ ODE.h \
+ Random_Stream.h \
+ Stimulation.h \
+ Thalamic_Column.h
+
+QMAKE_CXXFLAGS += -std=c++11 -O3
+
+SOURCES -= Thalamus_mex.cpp
diff --git a/Random_Stream.h b/Random_Stream.h
new file mode 100644
index 0000000..abd6c10
--- /dev/null
+++ b/Random_Stream.h
@@ -0,0 +1,81 @@
+/*
+ * Copyright (c) 2015 University of Lübeck
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com
+ * Stefanie Gareis: gareis@inb.uni-luebeck.de
+ *
+ */
+
+/****************************************************************************************************/
+/* Random number streams */
+/****************************************************************************************************/
+#pragma once
+#include
+
+/****************************************************************************************************/
+/* Struct for normal distribution */
+/****************************************************************************************************/
+struct random_stream_normal
+{
+ /* Random number engine: Mersenne-Twister */
+ std::mt19937_64 mt;
+ /* Random number generator: Normal-distribution */
+ std::normal_distribution norm_dist;
+
+ /* Constructors */
+ random_stream_normal(){}
+ random_stream_normal(double mean, double stddev)
+ : mt(rand()) , norm_dist(mean, stddev)
+ {}
+
+ /* Overwrites the function-call operator "( )" */
+ double operator( )(void) {
+ return norm_dist(mt);
+ }
+};
+/****************************************************************************************************/
+/* end */
+/****************************************************************************************************/
+
+/****************************************************************************************************/
+/* Struct for uniform int distribution */
+/****************************************************************************************************/
+struct random_stream_uniform_int
+{
+ /* Random number engine: Mersenne-Twister */
+ std::mt19937_64 mt;
+ /* Random number generator: Uniform integer-distribution */
+ std::uniform_int_distribution<> uniform_dist;
+
+ /* Constructors */
+ random_stream_uniform_int(){}
+ random_stream_uniform_int(double lower_bound, double upper_bound)
+ : mt(rand()) , uniform_dist(lower_bound, upper_bound)
+ {}
+
+ /* Overwrites the function-call operator "( )" */
+ double operator( )(void) {
+ return uniform_dist(mt);
+ }
+};
+/****************************************************************************************************/
+/* end */
+/****************************************************************************************************/
diff --git a/T_model.pro b/T_model.pro
deleted file mode 100644
index 2c79a62..0000000
--- a/T_model.pro
+++ /dev/null
@@ -1,17 +0,0 @@
-TEMPLATE = app
-CONFIG += console
-CONFIG -= app_bundle
-CONFIG -= qt
-
-SOURCES += Main.cpp \
- Thalamic_Column.cpp \
- Thalamus.cpp
-
-HEADERS += ODE.h \
- saves.h \
- Thalamic_Column.h \
- Stimulation.h
-
-QMAKE_CXXFLAGS += -std=c++11 -O3
-
-SOURCES -= Thalamus.cpp
diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp
index 8b33f82..80371c0 100644
--- a/Thalamic_Column.cpp
+++ b/Thalamic_Column.cpp
@@ -46,8 +46,8 @@ void Thalamic_Column::set_RNG(void) {
/* Create RNG for each stream */
for (int i=0; i
#include
-#include
-#include
-#include
+#include "Random_Stream.h"
using std::vector;
/****************************************************************************************************/
-/* Typedefs for RNG */
-/****************************************************************************************************/
-typedef boost::mt11213b ENG; /* Mersenne Twister */
-typedef boost::normal_distribution DIST; /* Normal Distribution */
-typedef boost::variate_generator GEN; /* Variate generator */
-/****************************************************************************************************/
-/* end */
-/****************************************************************************************************/
-
-
-/****************************************************************************************************/
/* Macro for vector initialization */
/****************************************************************************************************/
#ifndef _INIT
@@ -80,9 +67,9 @@ class Thalamic_Column {
/* Synaptic currents */
double I_et (int) const;
- double I_it (int) const;
+ double I_rt (int) const;
double I_er (int) const;
- double I_ir (int) const;
+ double I_rr (int) const;
/* Activation functions */
double m_inf_T_t (int) const;
@@ -124,12 +111,12 @@ class Thalamic_Column {
vector Vt = _INIT(E_L_t), /* TC membrane voltage */
Vr = _INIT(E_L_r), /* RE membrane voltage */
Ca = _INIT(Ca_0), /* Calcium concentration of TC population */
- y_tt = _INIT(0.0), /* PostSP from TC population to TC population */
- y_tr = _INIT(0.0), /* PostSP from TC population to RE population */
+ y_et = _INIT(0.0), /* PostSP from TC population to TC population */
+ y_er = _INIT(0.0), /* PostSP from TC population to RE population */
y_rt = _INIT(0.0), /* PostSP from RE population to TC population */
y_rr = _INIT(0.0), /* PostSP from RE population to RE population */
- x_tt = _INIT(0.0), /* derivative of y_tt */
- x_tr = _INIT(0.0), /* derivative of y_tr */
+ x_et = _INIT(0.0), /* derivative of y_tt */
+ x_er = _INIT(0.0), /* derivative of y_tr */
x_rt = _INIT(0.0), /* derivative of y_rt */
x_rr = _INIT(0.0), /* derivative of y_rr */
h_T_t = _INIT(0.0), /* inactivation of T channel */
@@ -138,7 +125,7 @@ class Thalamic_Column {
m_h2 = _INIT(0.0); /* activation of h channel bound with protein */
/* Random number generators */
- vector MTRands;
+ vector MTRands;
/* Container for noise */
vector Rand_vars;
@@ -224,6 +211,10 @@ class Thalamic_Column {
const double N_rt = 5;
const double N_rr = 30;
+ /* Parameters for SRK iteration */
+ const vector A = {0.5, 0.5, 1.0, 1.0};
+ const vector B = {0.75, 0.75, 0.0, 0.0};
+
friend class Stim;
};
/****************************************************************************************************/
diff --git a/Thalamus.cpp b/Thalamus.cpp
index b096aaa..1b63c2b 100644
--- a/Thalamus.cpp
+++ b/Thalamus.cpp
@@ -21,78 +21,48 @@
*/
/****************************************************************************************************/
-/* Implementation of the simulation as MATLAB routine (mex compiler) */
-/* mex command is given by: */
-/* mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" Thalamus.cpp Thalamic_Column.cpp */
-/* The Simulation requires the following boost libraries: Preprocessor */
-/* Random */
+/* Main file for compilation tests */
/****************************************************************************************************/
-#include "mex.h"
-#include "matrix.h"
+#include
+#include
#include "Thalamic_Column.h"
-#include "Stimulation.h"
-#include "saves.h"
/****************************************************************************************************/
/* Fixed simulation settings */
/****************************************************************************************************/
-extern const int onset = 15; /* time until data is stored in s */
+typedef std::chrono::high_resolution_clock::time_point timer;
+extern const int T = 30; /* Simulation length s */
extern const int res = 1E4; /* number of iteration steps per s */
-extern const int red = res/100; /* number of iterations that is saved */
extern const double dt = 1E3/res; /* duration of a timestep in ms */
extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */
/****************************************************************************************************/
/* end */
/****************************************************************************************************/
+
/****************************************************************************************************/
-/* Simulation routine */
-/* lhs defines outputs */
-/* rhs defines inputs */
+/* Main simulation routine */
/****************************************************************************************************/
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- /* Initialize the seeder */
- srand(time(NULL));
-
- /* Fetch inputs */
- const int T = (int) (mxGetScalar(prhs[0])); /* Duration of simulation in s */
- const int Time = (T+onset)*res; /* Total number of iteration steps */
- double* Param_Thalamus = mxGetPr (prhs[1]); /* Parameters of cortical module */
- double* var_stim = mxGetPr (prhs[2]); /* Parameters of stimulation protocol */
-
+int main(void) {
/* Initialize the populations */
- Thalamic_Column Thalamus(Param_Thalamus);
+ Thalamic_Column Thalamus = Thalamic_Column();
- /* Initialize the stimulation protocol */
- Stim Stimulation(Thalamus, var_stim);
+ /* Take the time of the simulation */
+ timer start,end;
- /* Create data containers */
- mxArray* Vt = SetMexArray(1, T*res/red);
- mxArray* Vr = SetMexArray(1, T*res/red);
- mxArray* ah = SetMexArray(1, T*res/red);
+ /* Simulation */
+ start = std::chrono::high_resolution_clock::now();
+ for (int t=0; t< T*res; ++t) {
+ Thalamus.iterate_ODE();
+ }
+ end = std::chrono::high_resolution_clock::now();
- /* Pointer to the actual data block */
- double* Pr_Vt = mxGetPr(Vt);
- double* Pr_Vr = mxGetPr(Vr);
- double* Pr_ah = mxGetPr(ah);
-
- /* Simulation */
- int count = 0;
- for (int t=0; t=onset*res && t%red==0){
- get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_ah);
- ++count;
- }
- }
-
- /* Output of the simulation */
- plhs[0] = Vt;
- plhs[1] = Vr;
- plhs[2] = ah;
-return;
+ /* Time consumed by the simulation */
+ double dif = 1E-3*std::chrono::duration_cast( end - start ).count();
+ std::cout << "simulation done!\n";
+ std::cout << "took " << dif << " seconds" << "\n";
+ std::cout << "end\n";
}
/****************************************************************************************************/
-/* end */
+/* end */
/****************************************************************************************************/
diff --git a/Main.cpp b/Thalamus_mex.cpp
similarity index 57%
rename from Main.cpp
rename to Thalamus_mex.cpp
index 5346ed3..ed2bc1b 100644
--- a/Main.cpp
+++ b/Thalamus_mex.cpp
@@ -21,49 +21,75 @@
*/
/****************************************************************************************************/
-/* Main file for compilation tests */
-/* The Simulation requires the following boost libraries: Random */
+/* Implementation of the simulation as MATLAB routine (mex compiler) */
+/* mex command is given by: */
+/* mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" Thalamus_mew.cpp Thalamic_Column.cpp */
/****************************************************************************************************/
-#include
-#include
-#include "Thalamic_Column.h"
+#include "mex.h"
+#include "matrix.h"
+#include "Data_Storage.h"
+#include "Stimulation.h"
/****************************************************************************************************/
/* Fixed simulation settings */
/****************************************************************************************************/
-typedef std::chrono::high_resolution_clock::time_point timer;
-extern const int T = 30; /* Simulation length s */
+extern const int onset = 15; /* time until data is stored in s */
extern const int res = 1E4; /* number of iteration steps per s */
+extern const int red = res/100; /* number of iterations that is saved */
extern const double dt = 1E3/res; /* duration of a timestep in ms */
extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */
/****************************************************************************************************/
/* end */
/****************************************************************************************************/
-
/****************************************************************************************************/
-/* Main simulation routine */
+/* Simulation routine */
+/* lhs defines outputs */
+/* rhs defines inputs */
/****************************************************************************************************/
-int main(void) {
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+ /* Initialize the seeder */
+ srand(time(NULL));
+
+ /* Fetch inputs */
+ const int T = (int) (mxGetScalar(prhs[0])); /* Duration of simulation in s */
+ const int Time = (T+onset)*res; /* Total number of iteration steps */
+ double* Param_Thalamus = mxGetPr (prhs[1]); /* Parameters of cortical module */
+ double* var_stim = mxGetPr (prhs[2]); /* Parameters of stimulation protocol */
+
/* Initialize the populations */
- Thalamic_Column Thalamus = Thalamic_Column();
+ Thalamic_Column Thalamus(Param_Thalamus);
- /* Take the time of the simulation */
- timer start,end;
+ /* Initialize the stimulation protocol */
+ Stim Stimulation(Thalamus, var_stim);
- /* Simulation */
- start = std::chrono::high_resolution_clock::now();
- for (int t=0; t< T*res; ++t) {
- Thalamus.iterate_ODE();
- }
- end = std::chrono::high_resolution_clock::now();
+ /* Create data containers */
+ mxArray* Vt = SetMexArray(1, T*res/red);
+ mxArray* Vr = SetMexArray(1, T*res/red);
+ mxArray* ah = SetMexArray(1, T*res/red);
- /* Time consumed by the simulation */
- double dif = 1E-3*std::chrono::duration_cast( end - start ).count();
- std::cout << "simulation done!\n";
- std::cout << "took " << dif << " seconds" << "\n";
- std::cout << "end\n";
+ /* Pointer to the actual data block */
+ double* Pr_Vt = mxGetPr(Vt);
+ double* Pr_Vr = mxGetPr(Vr);
+ double* Pr_ah = mxGetPr(ah);
+
+ /* Simulation */
+ int count = 0;
+ for (int t=0; t=onset*res && t%red==0){
+ get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_ah);
+ ++count;
+ }
+ }
+
+ /* Output of the simulation */
+ plhs[0] = Vt;
+ plhs[1] = Vr;
+ plhs[2] = ah;
+return;
}
/****************************************************************************************************/
-/* end */
+/* end */
/****************************************************************************************************/