-
Notifications
You must be signed in to change notification settings - Fork 1
/
Bazhenov_mex.cpp
120 lines (106 loc) · 5.43 KB
/
Bazhenov_mex.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
/*
* Copyright (c) 2016 Michael Schellenberger Costa [email protected]
*
* 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 Bazhenov2002 model as MATLAB routine (mex compiler) */
/* mex command is given by: */
/* mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3 -fopenmp" */
/* Bazhenov.cpp Pyramidal_Neuron.cpp Inhibitory_Neuron.cpp */
/* The Simulation requires the following boost libraries: Random */
/****************************************************************************************************/
#include "mex.h"
#include "matrix.h"
#include <vector>
#include "Data_Storage.h"
#include "Initialize_Neurons.h"
#include "Iterate_ODE.h"
/****************************************************************************************************/
/* Fixed simulation settings */
/****************************************************************************************************/
typedef std::chrono::high_resolution_clock::time_point timer;
extern const int T = 1; /* Simulation length in s */
extern const int res = 5E4; /* Number of iteration steps per s */
extern const double dt = 1E3/res; /* Duration of a timestep in ms */
extern const std::vector<int> NumCells = {128, /* Number of pyramidal cells */
32, /* Number of inhibitory cells */
128, /* Number of thalamocortical cells */
32}; /* Number of reticular cells */
extern const int N_Cores= 7; /* Number of CPU cores */
/****************************************************************************************************/
/* end */
/****************************************************************************************************/
mxArray* SetMexArray(int N, int M);
/****************************************************************************************************/
/* Simulation routine */
/* lhs defines outputs */
/* rhs defines inputs */
/****************************************************************************************************/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
/* Seed the random number generator */
srand(time(NULL));
/* Initialize the populations */
std::vector<Pyramidal_Neuron> PY;
std::vector<Inhibitory_Neuron> IN;
std::vector<Thalamocortical_Neuron> TC;
std::vector<Reticular_Neuron> RE;
setupNetwork(PY, IN, TC, RE);
/* Data container in MATLAB format */
std::vector<mxArray*> Data;
Data.push_back(SetMexArray(NumCells[PYRAMIDAL], T*res/red)); // Ve
Data.push_back(SetMexArray(NumCells[INHIBITORY], T*res/red)); // Vi
Data.push_back(SetMexArray(NumCells[PYRAMIDAL], T*res/red)); // Ca
/* Pointer to the data blocks */
std::vector<double*> pData;
pDate.reserve(Data.size());
for(const auto& arrayptr : Data) {
pData.push_back(mxGetPr(arrayptr));
}
/* Simulation */
int count = 0;
for (int t = 0; t < T*res; ++t) {
Iterate_ODE(PY, IN, TC, RE);
if(t%red==0){
get_data(count++, PY, IN, TC, RE, pData);
}
}
/* Return the data containers */
nlhs = Data.size();
for(const auto& arrayptr : Data) {
plhs[&arrayptr - Data.data()] = arrayptr;
}
return;
}
/****************************************************************************************************/
/* end */
/****************************************************************************************************/
/****************************************************************************************************/
/* Create MATLAB data container */
/****************************************************************************************************/
mxArray* SetMexArray(int N, int M) {
mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL);
mxSetM(Array, N);
mxSetN(Array, M);
mxSetData(Array, mxMalloc(sizeof(double)*M*N));
return Array;
}
/****************************************************************************************************/
/* end */
/****************************************************************************************************/