diff --git a/Main.cpp b/Main.cpp index 9095b5e..d38e8c3 100644 --- a/Main.cpp +++ b/Main.cpp @@ -28,7 +28,6 @@ #include #include #include "Thalamic_Column.h" -#include "ODE.h" /****************************************************************************************************/ /* Fixed simulation settings */ @@ -44,6 +43,22 @@ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ +/* Constants for SRK4 iteration */ +/****************************************************************************************************/ +extern const vector B1 = {0, + 0.626708569400000081728308032325, + 1.7296310295000001389098542858846, + 1.2703689705000000831347506391467}; +extern const vector B2 = {0, + 0.78000033203198970710445792065002, + 1.28727807507536762265942797967, + 0.44477273249350995909523476257164}; +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + + +/****************************************************************************************************/ /* Main simulation routine */ /****************************************************************************************************/ int main(void) { @@ -56,7 +71,7 @@ int main(void) { /* Simulation */ for (int t=0; t< T*res; ++t) { - ODE (Thalamus); + Thalamus.iterate_ODE(); } time (&end); diff --git a/ODE.h b/ODE.h deleted file mode 100644 index cb2165f..0000000 --- a/ODE.h +++ /dev/null @@ -1,59 +0,0 @@ -/* -* Copyright (c) 2014 Michael Schellenberger Costa -* -* 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. -*/ - -/****************************************************************************************************/ -/* Implementation of the ODE solver */ -/****************************************************************************************************/ -#pragma once -#include "Thalamic_Column.h" - -/****************************************************************************************************/ -/* Evaluation of SRK4 */ -/****************************************************************************************************/ -inline void ODE(Thalamic_Column& Thalamus) { - /* First calculating every ith RK moment. Has to be in order, 1th moment first */ - for (int i=1; i<=4; ++i) { - Thalamus.set_RK(i); - } - - /* Add moments */ - Thalamus.add_RK(); -} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Parameters for SRK4 integration */ -/****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 94593af..1ca2007 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -303,3 +303,18 @@ void Thalamic_Column::add_RK(void) { /****************************************************************************************************/ /* end */ /****************************************************************************************************/ + + +/****************************************************************************************************/ +/* Evaluation of SRK4 */ +/****************************************************************************************************/ +void Thalamic_Column::iterate_ODE(void) { + /* First calculating every ith RK moment. Has to be in order, 1th moment first */ + for (int i=1; i<=4; ++i) { + set_RK(i); + } + add_RK(); +} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/Thalamic_Column.h b/Thalamic_Column.h index f42bfc5..02b2371 100644 --- a/Thalamic_Column.h +++ b/Thalamic_Column.h @@ -101,6 +101,7 @@ class Thalamic_Column { /* ODE functions */ void set_RK (int); void add_RK (void); + void iterate_ODE (void); /* Data storage access */ friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 2)); @@ -144,8 +145,8 @@ class Thalamic_Column { const double theta_r = -58.6; /* Sigmoid gain in mV */ - const double sigma_t = 3; - const double sigma_r = 3; + const double sigma_t = 4; + const double sigma_r = 4; /* Scaling parameter for sigmoidal mapping (dimensionless) */ const double C1 = (3.14159265/sqrt(3)); @@ -207,9 +208,9 @@ class Thalamic_Column { double input = 0.0; /* Connectivities (dimensionless) */ - const double N_tr = 5; + const double N_tr = 6; const double N_rt = 5; - const double N_rr = 50; + const double N_rr = 100; friend class Stim; }; diff --git a/Thalamus.cpp b/Thalamus.cpp index 0d9d672..797e17b 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -27,13 +27,11 @@ /* The Simulation requires the following boost libraries: Preprocessor */ /* Random */ /****************************************************************************************************/ -#include #include "mex.h" #include "matrix.h" #include "Thalamic_Column.h" #include "Stimulation.h" #include "saves.h" -#include "ODE.h" /****************************************************************************************************/ /* Fixed simulation settings */ @@ -49,6 +47,22 @@ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ +/* Constants for SRK4 iteration */ +/****************************************************************************************************/ +extern const vector B1 = {0, + 0.626708569400000081728308032325, + 1.7296310295000001389098542858846, + 1.2703689705000000831347506391467}; +extern const vector B2 = {0, + 0.78000033203198970710445792065002, + 1.28727807507536762265942797967, + 0.44477273249350995909523476257164}; +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + + +/****************************************************************************************************/ /* Simulation routine */ /* lhs defines outputs */ /* rhs defines inputs */ @@ -80,7 +94,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Simulation */ int count = 0; for (int t=0; t=onset*res && t%red==0){ get_data(count, Thalamus, Pr_Vt, Pr_Vr);