Merge branch 'Test_param'

Conflicts: Plots.m Thalamic_Column.cpp Thalamic_Column.h Thalamus.cpp saves.h
miscco · 3a9920a2c5d0be313ef12875589b1881085e6bca · miscco committed Feb 8, 2016 · Feb 8, 2016 · 3a9920a
2 parents 8a73b3f + a3bfdae
commit 3a9920a2c5d0be313ef12875589b1881085e6bca
diff --git a/.cproject b/.cproject
diff --git a/.project b/.project
@@ -1,83 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<projectDescription>
-	<name>Steyn-Ross</name>
-	<comment></comment>
-	<projects>
-	</projects>
-	<buildSpec>
-		<buildCommand>
-			<name>org.eclipse.cdt.managedbuilder.core.genmakebuilder</name>
-			<triggers>clean,full,incremental,</triggers>
-			<arguments>
-				<dictionary>
-					<key>?name?</key>
-					<value></value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.append_environment</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.autoBuildTarget</key>
-					<value>all</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.buildArguments</key>
-					<value></value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.buildCommand</key>
-					<value>make</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.buildLocation</key>
-					<value>${workspace_loc:/Steyn-Ross/Debug}</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
-					<value>clean</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.contents</key>
-					<value>org.eclipse.cdt.make.core.activeConfigSettings</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableAutoBuild</key>
-					<value>false</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableCleanBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableFullBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
-					<value>all</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.stopOnError</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
-					<value>true</value>
-				</dictionary>
-			</arguments>
-		</buildCommand>
-		<buildCommand>
-			<name>org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder</name>
-			<triggers>full,incremental,</triggers>
-			<arguments>
-			</arguments>
-		</buildCommand>
-	</buildSpec>
-	<natures>
-		<nature>org.eclipse.cdt.core.cnature</nature>
-		<nature>org.eclipse.cdt.core.ccnature</nature>
-		<nature>org.eclipse.cdt.managedbuilder.core.managedBuildNature</nature>
-		<nature>org.eclipse.cdt.managedbuilder.core.ScannerConfigNature</nature>
-	</natures>
-</projectDescription>
diff --git a/Data_Storage.h b/Data_Storage.h
@@ -0,0 +1,47 @@
+/*
+ *	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
+ *
+ *	Based on:	A thalamocortical neural mass model of the EEG during NREM sleep and its response
+ *				to auditory stimulation.
+ *				M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
+ *				JC Claussen.
+ *				PLoS Computational Biology In Review (in review).
+ */
+
+/****************************************************************************************************/
+/*									Functions for data storage										*/
+/****************************************************************************************************/
+#pragma once
+#include "Thalamic_Column.h"
+
+/****************************************************************************************************/
+/*											Save data												*/
+/****************************************************************************************************/
+void get_data(int counter, Thalamic_Column& Col,	double* Vt,	double* Vr,	double* ah) {
+	Vt 		[counter] = Col.Vt 		[0];
+	Vr 		[counter] = Col.Vr 		[0];
+	ah 		[counter] = Col.act_h 	();
+}
+/****************************************************************************************************/
+/*										 		end													*/
+/****************************************************************************************************/
diff --git a/Main.cpp b/Main.cpp
@@ -1,86 +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.
-*/
-
-/****************************************************************************************************/
-/*		Main file for compilation tests																*/
-/*		The Simulation requires the following boost libraries:	Preprocessor						*/
-/*																Random								*/
-/****************************************************************************************************/
-#include <iostream>
-#include <ctime>
-#include "Thalamic_Column.h"
-
-/****************************************************************************************************/
-/*										Fixed simulation settings									*/
-/****************************************************************************************************/
-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			 										*/
-/****************************************************************************************************/
-
-
-/****************************************************************************************************/
-/*									Constants for SRK4 iteration									*/
-/****************************************************************************************************/
-extern const vector<double> B1 = {0,
-								  0.626708569400000081728308032325,
-								  1.7296310295000001389098542858846,
-		 	 	 	 	 	 	  1.2703689705000000831347506391467};
-extern const vector<double> B2 = {0,
-								  0.78000033203198970710445792065002,
-								  1.28727807507536762265942797967,
-								  0.44477273249350995909523476257164};
-/****************************************************************************************************/
-/*										 		end													*/
-/****************************************************************************************************/
-
-
-/****************************************************************************************************/
-/*										Main simulation routine										*/
-/****************************************************************************************************/
-int main(void) {
-	/* Initialize the populations */
-	Thalamic_Column Thalamus;
-
-	/* Takes the time of the simulation */
-	time_t start,end;
-	time (&start);
-
-	/* Simulation */
-	for (int t=0; t< T*res; ++t) {
-		Thalamus.iterate_ODE();
-	}
-	time (&end);
-
-	/* Time consumed by the simulation */
-	double dif = difftime(end,start);
-	std::cout << "simulation done!\n";
-	std::cout << "took " << dif 	<< " seconds" << "\n";
-	std::cout << "end\n";
-}
-/****************************************************************************************************/
-/*										 		end			 										*/
-/****************************************************************************************************/
diff --git a/NM_Thalamus.pro b/NM_Thalamus.pro
@@ -0,0 +1,19 @@
+TEMPLATE = app
+CONFIG += console
+CONFIG -= app_bundle
+CONFIG -= qt
+
+TARGET = Thalamus.cpp
+
+SOURCES +=  Thalamic_Column.cpp \
+	    Thalamus.cpp	\
+	    Thalamus_mex.cpp
+
+HEADERS +=  ODE.h		\
+	    Data_Storage.h	\
+	    Random_Stream.h	\
+	    Thalamic_Column.h
+
+QMAKE_CXXFLAGS += -std=c++11 -O3
+
+SOURCES -= Thalamus_mex.cpp
diff --git a/Plots.m b/Plots.m
@@ -1,51 +1,24 @@
 % mex command is given by: 
-% mex CXXFLAGS="\$CXXFLAGS -std=c++11" Thalamus.cpp Thalamic_Column.cpp
+% mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" Thalamus_mex.cpp Thalamic_Column.cpp
 
 function Plots(T)
 
 if nargin == 0
-    Con     	= [ 6;      % sigma
-                    4;		% N_tr
-                    4;		% N_rt
-                    20];		% N_rr 
-               
+    Con     	= [ 0.01;       % g_h
+                    0.05];      % g_LK
 
-    % stimulation parameters
-    % first number is the mode of stimulation
-    % 0 == none
-    % 1 == periodic
-    % 2 == phase dependend up state
-    % 3 == phase dependend down state
-    
-    var_stim    = [ 0;           % mode of stimulation
-                    50;          % strength of the stimulus      in Hz (spikes per second)
-                    120;       	% duration of the stimulus      in ms
-                    7;          % time between stimuli          in s    
-                    1];         % time until stimuli after min 	in ms
     T       	= 30;  		% duration of the simulation
 end
-[Vt, Vr, I_t, I_h] = Thalamus(T, Con, var_stim);
+[Vt, Vr, ah] = Thalamus_mex(T, Con, var_stim);
 
 L        = max(size(Vt));
 timeaxis = linspace(0,T,L);
 
-Fs      = L/T;
-[Pxx,f] = pwelch(Vt-mean(Vt), hamming(4*Fs), 2*Fs, [], Fs,'ConfidenceLevel', 0.9, 'power');
-n       = find(f<=60, 1, 'last' );
-
+% Plot the data
 figure(1)
 subplot(311), plot(timeaxis,Vt)
-title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV')
+title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt [mV]')
 subplot(312), plot(timeaxis,Vr)
-title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
-subplot(313), plot(f(1:n),Pxx(1:n))
-title('Powerspectrum of Steyn-Ross model with pwelch'), xlabel('frequency in Hz'), ylabel('Power')
-%save('Thalamus.mat','Vt','Vr')
-% 
-% figure(2)
-% subplot(311), plot(timeaxis,Vt)
-% title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV')
-% subplot(312), plot(timeaxis,I_t)
-% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
-% subplot(313), plot(timeaxis,I_h)
-% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
+title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr [mV]')
+subplot(313), plot(timeaxis,ah)
+title('Thalamic relay I_h activation'), xlabel('time in s'), ylabel('m_h')
diff --git a/Random_Stream.h b/Random_Stream.h
@@ -0,0 +1,80 @@
+/*
+ *	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 <random>
+
+/****************************************************************************************************/
+/*									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<double>    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													*/
+/****************************************************************************************************/