diff --git a/Main.cpp b/Main.cpp index 076dfa4..11dea93 100644 --- a/Main.cpp +++ b/Main.cpp @@ -4,7 +4,6 @@ #include #include #include -#include "randoms.h" #include "Thalamic_Column.h" #include "ODE.h" @@ -17,12 +16,8 @@ extern const double h = sqrt(dt); // simulation of the thalamic model int main(void) { - // Initializing the mersenne twister. - MTRand mtrand; - - // creating the random input - vector u_t1 = rand_inp(mtrand, res, T, 0, 10, 1E3, phi_st, phi_st, phi_inp); - vector u_t2 = rand_inp(mtrand, res, T, 0, 10, 1E3, phi_st, phi_st, phi_inp); + // Initializing the seeder. + srand(time(0)); // Initializing the populations; Thalamic_Column Col; @@ -33,8 +28,7 @@ int main(void) { // simulation for (int t=0; t< T*res; ++t) { - ODE (Col, u_t1[t], u_t2[t]); - //ODE2(Col, u_t1[t]); + ODE (Col); } time (&end); diff --git a/MersenneTwister.h b/MersenneTwister.h deleted file mode 100644 index b212569..0000000 --- a/MersenneTwister.h +++ /dev/null @@ -1,462 +0,0 @@ -// MersenneTwister.h -// Mersenne Twister random number generator -- a C++ class MTRand -// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus -// Richard J. Wagner v1.1 28 September 2009 wagnerr@umich.edu - -// The Mersenne Twister is an algorithm for generating random numbers. It -// was designed with consideration of the flaws in various other generators. -// The period, 2^19937-1, and the order of equidistribution, 623 dimensions, -// are far greater. The generator is also fast; it avoids multiplication and -// division, and it benefits from caches and pipelines. For more information -// see the inventors' web page at -// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html - -// Reference -// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally -// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on -// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. - -// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, -// Copyright (C) 2000 - 2009, Richard J. Wagner -// All rights reserved. -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions -// are met: -// -// 1. Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. -// -// 2. Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in the -// documentation and/or other materials provided with the distribution. -// -// 3. The names of its contributors may not be used to endorse or promote -// products derived from this software without specific prior written -// permission. -// -// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE -// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. - -// The original code included the following notice: -// -// When you use this, send an email to: m-mat@math.sci.hiroshima-u.ac.jp -// with an appropriate reference to your work. -// -// It would be nice to CC: wagnerr@umich.edu and Cokus@math.washington.edu -// when you write. - -#pragma once -#ifndef MERSENNETWISTER_H -#define MERSENNETWISTER_H - -// Not thread safe (unless auto-initialization is avoided and each thread has -// its own MTRand object) - -#include -#include -#include -#include -#include - -class MTRand { -// Data -public: - typedef unsigned long uint32; // unsigned integer type, at least 32 bits - - enum { N = 624 }; // length of state vector - enum { SAVE = N + 1 }; // length of array for save() - -protected: - enum { M = 397 }; // period parameter - - uint32 state[N]; // internal state - uint32 *pNext; // next value to get from state - int left; // number of values left before reload needed - -// Methods -public: - MTRand( const uint32 oneSeed ); // initialize with a simple uint32 - MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or array - MTRand(); // auto-initialize with /dev/urandom or time() and clock() - MTRand( const MTRand& o ); // copy - - // Do NOT use for CRYPTOGRAPHY without securely hashing several returned - // values together, otherwise the generator state can be learned after - // reading 624 consecutive values. - - // Access to 32-bit random numbers - uint32 randInt(); // integer in [0,2^32-1] - uint32 randInt( const uint32 n ); // integer in [0,n] for n < 2^32 - double rand(); // real number in [0,1] - double rand( const double n ); // real number in [0,n] - double randExc(); // real number in [0,1) - double randExc( const double n ); // real number in [0,n) - double randDblExc(); // real number in (0,1) - double randDblExc( const double n ); // real number in (0,n) - double operator()(); // same as rand() - - // Access to 53-bit random numbers (capacity of IEEE double precision) - double rand53(); // real number in [0,1) - - // Access to nonuniform random number distributions - double randNorm( const double mean = 0.0, const double stddev = 1.0 ); - - // Re-seeding functions with same behavior as initializers - void seed( const uint32 oneSeed ); - void seed( uint32 *const bigSeed, const uint32 seedLength = N ); - void seed(); - - // Saving and loading generator state - void save( uint32* saveArray ) const; // to array of size SAVE - void load( uint32 *const loadArray ); // from such array - friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ); - friend std::istream& operator>>( std::istream& is, MTRand& mtrand ); - MTRand& operator=( const MTRand& o ); - -protected: - void initialize( const uint32 oneSeed ); - void reload(); - uint32 hiBit( const uint32 u ) const { return u & 0x80000000UL; } - uint32 loBit( const uint32 u ) const { return u & 0x00000001UL; } - uint32 loBits( const uint32 u ) const { return u & 0x7fffffffUL; } - uint32 mixBits( const uint32 u, const uint32 v ) const - { return hiBit(u) | loBits(v); } - uint32 magic( const uint32 u ) const - { return loBit(u) ? 0x9908b0dfUL : 0x0UL; } - uint32 twist( const uint32 m, const uint32 s0, const uint32 s1 ) const - { return m ^ (mixBits(s0,s1)>>1) ^ magic(s1); } - static uint32 hash( time_t t, clock_t c ); -}; - -// Functions are defined in order of usage to assist inlining - -inline MTRand::uint32 MTRand::hash( time_t t, clock_t c ) -{ - // Get a uint32 from t and c - // Better than uint32(x) in case x is floating point in [0,1] - // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk) - - static uint32 differ = 0; // guarantee time-based seeds will change - - uint32 h1 = 0; - unsigned char *p = (unsigned char *) &t; - for( size_t i = 0; i < sizeof(t); ++i ) - { - h1 *= UCHAR_MAX + 2U; - h1 += p[i]; - } - uint32 h2 = 0; - p = (unsigned char *) &c; - for( size_t j = 0; j < sizeof(c); ++j ) - { - h2 *= UCHAR_MAX + 2U; - h2 += p[j]; - } - return ( h1 + differ++ ) ^ h2; -} - -inline void MTRand::initialize( const uint32 seed ) -{ - // Initialize generator state with seed - // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. - // In previous versions, most significant bits (MSBs) of the seed affect - // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. - register uint32 *s = state; - register uint32 *r = state; - register int i = 1; - *s++ = seed & 0xffffffffUL; - for( ; i < N; ++i ) - { - *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL; - r++; - } -} - -inline void MTRand::reload() -{ - // Generate N new values in state - // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com) - static const int MmN = int(M) - int(N); // in case enums are unsigned - register uint32 *p = state; - register int i; - for( i = N - M; i--; ++p ) - *p = twist( p[M], p[0], p[1] ); - for( i = M; --i; ++p ) - *p = twist( p[MmN], p[0], p[1] ); - *p = twist( p[MmN], p[0], state[0] ); - - left = N, pNext = state; -} - -inline void MTRand::seed( const uint32 oneSeed ) -{ - // Seed the generator with a simple uint32 - initialize(oneSeed); - reload(); -} - -inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength ) -{ - // Seed the generator with an array of uint32's - // There are 2^19937-1 possible initial states. This function allows - // all of those to be accessed by providing at least 19937 bits (with a - // default seed length of N = 624 uint32's). Any bits above the lower 32 - // in each element are discarded. - // Just call seed() if you want to get array from /dev/urandom - initialize(19650218UL); - register int i = 1; - register uint32 j = 0; - register int k = ( N > seedLength ? N : seedLength ); - for( ; k; --k ) - { - state[i] = - state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL ); - state[i] += ( bigSeed[j] & 0xffffffffUL ) + j; - state[i] &= 0xffffffffUL; - ++i; ++j; - if( i >= N ) { state[0] = state[N-1]; i = 1; } - if( j >= seedLength ) j = 0; - } - for( k = N - 1; k; --k ) - { - state[i] = - state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL ); - state[i] -= i; - state[i] &= 0xffffffffUL; - ++i; - if( i >= N ) { state[0] = state[N-1]; i = 1; } - } - state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array - reload(); -} - -inline void MTRand::seed() -{ - // Seed the generator with an array from /dev/urandom if available - // Otherwise use a hash of time() and clock() values - - // First try getting an array from /dev/urandom - FILE* urandom = fopen( "/dev/urandom", "rb" ); - if( urandom ) - { - uint32 bigSeed[N]; - register uint32 *s = bigSeed; - register int i = N; - register bool success = true; - while( success && i-- ) - success = fread( s++, sizeof(uint32), 1, urandom ); - fclose(urandom); - if( success ) { seed( bigSeed, N ); return; } - } - - // Was not successful, so use time() and clock() instead - seed( hash( time(NULL), clock() ) ); -} - -inline MTRand::MTRand( const uint32 oneSeed ) - { seed(oneSeed); } - -inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength ) - { seed(bigSeed,seedLength); } - -inline MTRand::MTRand() - { seed(); } - -inline MTRand::MTRand( const MTRand& o ) -{ - register const uint32 *t = o.state; - register uint32 *s = state; - register int i = N; - for( ; i--; *s++ = *t++ ) {} - left = o.left; - pNext = &state[N-left]; -} - -inline MTRand::uint32 MTRand::randInt() -{ - // Pull a 32-bit integer from the generator state - // Every other access function simply transforms the numbers extracted here - - if( left == 0 ) reload(); - --left; - - register uint32 s1; - s1 = *pNext++; - s1 ^= (s1 >> 11); - s1 ^= (s1 << 7) & 0x9d2c5680UL; - s1 ^= (s1 << 15) & 0xefc60000UL; - return ( s1 ^ (s1 >> 18) ); -} - -inline MTRand::uint32 MTRand::randInt( const uint32 n ) -{ - // Find which bits are used in n - // Optimized by Magnus Jonsson (magnus@smartelectronix.com) - uint32 used = n; - used |= used >> 1; - used |= used >> 2; - used |= used >> 4; - used |= used >> 8; - used |= used >> 16; - - // Draw numbers until one is found in [0,n] - uint32 i; - do - i = randInt() & used; // toss unused bits to shorten search - while( i > n ); - return i; -} - -inline double MTRand::rand() - { return double(randInt()) * (1.0/4294967295.0); } - -inline double MTRand::rand( const double n ) - { return rand() * n; } - -inline double MTRand::randExc() - { return double(randInt()) * (1.0/4294967296.0); } - -inline double MTRand::randExc( const double n ) - { return randExc() * n; } - -inline double MTRand::randDblExc() - { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); } - -inline double MTRand::randDblExc( const double n ) - { return randDblExc() * n; } - -inline double MTRand::rand53() -{ - uint32 a = randInt() >> 5, b = randInt() >> 6; - return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada -} - -inline double MTRand::randNorm( const double mean, const double stddev ) -{ - // Return a real number from a normal (Gaussian) distribution with given - // mean and standard deviation by polar form of Box-Muller transformation - double x, y, r; - do - { - x = 2.0 * rand() - 1.0; - y = 2.0 * rand() - 1.0; - r = x * x + y * y; - } - while ( r >= 1.0 || r == 0.0 ); - double s = sqrt( -2.0 * log(r) / r ); - return mean + x * s * stddev; -} - -inline double MTRand::operator()() -{ - return rand(); -} - -inline void MTRand::save( uint32* saveArray ) const -{ - register const uint32 *s = state; - register uint32 *sa = saveArray; - register int i = N; - for( ; i--; *sa++ = *s++ ) {} - *sa = left; -} - -inline void MTRand::load( uint32 *const loadArray ) -{ - register uint32 *s = state; - register uint32 *la = loadArray; - register int i = N; - for( ; i--; *s++ = *la++ ) {} - left = *la; - pNext = &state[N-left]; -} - -inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ) -{ - register const MTRand::uint32 *s = mtrand.state; - register int i = mtrand.N; - for( ; i--; os << *s++ << "\t" ) {} - return os << mtrand.left; -} - -inline std::istream& operator>>( std::istream& is, MTRand& mtrand ) -{ - register MTRand::uint32 *s = mtrand.state; - register int i = mtrand.N; - for( ; i--; is >> *s++ ) {} - is >> mtrand.left; - mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left]; - return is; -} - -inline MTRand& MTRand::operator=( const MTRand& o ) -{ - if( this == &o ) return (*this); - register const uint32 *t = o.state; - register uint32 *s = state; - register int i = N; - for( ; i--; *s++ = *t++ ) {} - left = o.left; - pNext = &state[N-left]; - return (*this); -} - -#endif // MERSENNETWISTER_H - -// Change log: -// -// v0.1 - First release on 15 May 2000 -// - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus -// - Translated from C to C++ -// - Made completely ANSI compliant -// - Designed convenient interface for initialization, seeding, and -// obtaining numbers in default or user-defined ranges -// - Added automatic seeding from /dev/urandom or time() and clock() -// - Provided functions for saving and loading generator state -// -// v0.2 - Fixed bug which reloaded generator one step too late -// -// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew -// -// v0.4 - Removed trailing newline in saved generator format to be consistent -// with output format of built-in types -// -// v0.5 - Improved portability by replacing static const int's with enum's and -// clarifying return values in seed(); suggested by Eric Heimburg -// - Removed MAXINT constant; use 0xffffffffUL instead -// -// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits -// - Changed integer [0,n] generator to give better uniformity -// -// v0.7 - Fixed operator precedence ambiguity in reload() -// - Added access for real numbers in (0,1) and (0,n) -// -// v0.8 - Included time.h header to properly support time_t and clock_t -// -// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto -// - Allowed for seeding with arrays of any length -// - Added access for real numbers in [0,1) with 53-bit resolution -// - Added access for real numbers from normal (Gaussian) distributions -// - Increased overall speed by optimizing twist() -// - Doubled speed of integer [0,n] generation -// - Fixed out-of-range number generation on 64-bit machines -// - Improved portability by substituting literal constants for long enum's -// - Changed license from GNU LGPL to BSD -// -// v1.1 - Corrected parameter label in randNorm from "variance" to "stddev" -// - Changed randNorm algorithm from basic to polar form for efficiency -// - Updated includes from deprecated to standard forms -// - Cleaned declarations and definitions to please Intel compiler -// - Revised twist() operator to work on ones'-complement machines -// - Fixed reload() function to work when N and M are unsigned -// - Added copy constructor and copy operator from Salvador Espana diff --git a/ODE.h b/ODE.h index 9c51073..7ab66c4 100644 --- a/ODE.h +++ b/ODE.h @@ -4,10 +4,26 @@ #include #include "Thalamic_Column.h" +/*****************************************************************************************************/ +/***************************** parameters for integration of SRK4 **************************/ +/*****************************************************************************************************/ +// vectors needed for stochastic RK +extern const vector B1 = {0, + 0.626708569400000081728308032325, + 1.7296310295000001389098542858846, + 1.2703689705000000831347506391467}; +extern const vector B2 = {0, + 0.78000033203198970710445792065002, + 1.28727807507536762265942797967, + 0.44477273249350995909523476257164}; +/*****************************************************************************************************/ +/********************************** end **********************************/ +/*****************************************************************************************************/ + // function that evaluates ODE using stochastic Runge Kutta -inline void ODE(Thalamic_Column& Col, double ut1, double ut2) { +inline void ODE(Thalamic_Column& Col) { for (auto i=1; i<=4; ++i) { - Col.set_RK(i, ut1, ut2); + Col.set_RK(i); } - Col.add_RK(ut1); + Col.add_RK(); } diff --git a/Spindle.png b/Spindle.png deleted file mode 100644 index 51971cc..0000000 Binary files a/Spindle.png and /dev/null differ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 38642f1..eb33e1d 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -227,11 +227,11 @@ double Thalamic_Column::I_CAN (int N) const{ /*****************************************************************************************************/ /********************************** RK noise scaling **********************************/ /*****************************************************************************************************/ -// function that returns the noise to excitatory population for stochastic RK4 -double Thalamic_Column::noise_xRK(int N, double u_1, double u_2) const{ +// function that returns the noise to exitatory population for stochastic RK4 +double Cortical_Column::noise_xRK(int N) const{ extern const double h; extern const vector B1, B2; - double n = 1 / h * (B1[N-1] * u_1 + B2[N-1] * u_2); + double n = 1 / h * (B1[N-1] * Rand_vars[0] + B2[N-1] * Rand_vars[1]); return n; } /*****************************************************************************************************/ @@ -243,7 +243,7 @@ double Thalamic_Column::noise_xRK(int N, double u_1, double u_2) const{ /********************************** ODE functions **********************************/ /*****************************************************************************************************/ // function that calculates the Nth RK term -void Thalamic_Column::set_RK (int N, double u_t1, double u_t2) { +void Thalamic_Column::set_RK (int N) { extern const double dt; _SWITCH((Cat) (Car) (Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr) @@ -267,10 +267,10 @@ void Thalamic_Column::set_RK (int N, double u_t1, double u_t2) { Phi_tr [N] = dt*(var_x_tr); Phi_rt [N] = dt*(var_x_rt); Phi_rr [N] = dt*(var_x_rr); - x_tt [N] = dt*(pow(gamma_t, 2) * (noise_xRK(N, u_t1, u_t2) - var_Phi_tt) - 2 * gamma_t * var_x_tt); - x_tr [N] = dt*(pow(gamma_t, 2) * (N_tr * get_Qt(N) - var_Phi_tr) - 2 * gamma_t * var_x_tr); - x_rt [N] = dt*(pow(gamma_r, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_r * var_x_rt); - x_rr [N] = dt*(pow(gamma_r, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_r * var_x_rr); + x_tt [N] = dt*(pow(gamma_t, 2) * (noise_xRK(N) - var_Phi_tt) - 2 * gamma_t * var_x_tt); + x_tr [N] = dt*(pow(gamma_t, 2) * (N_tr * get_Qt(N) - var_Phi_tr) - 2 * gamma_t * var_x_tr); + x_rt [N] = dt*(pow(gamma_r, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_r * var_x_rt); + x_rr [N] = dt*(pow(gamma_r, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_r * var_x_rr); } // function that ads all the RK terms together @@ -284,7 +284,7 @@ void Thalamic_Column::add_RK(double u_t) { Phi_tr [0] += (Phi_tr [1] + Phi_tr [2] * 2 + Phi_tr [3] * 2 + Phi_tr [4])/6; Phi_rt [0] += (Phi_rt [1] + Phi_rt [2] * 2 + Phi_rt [3] * 2 + Phi_rt [4])/6; Phi_rr [0] += (Phi_rr [1] + Phi_rr [2] * 2 + Phi_rr [3] * 2 + Phi_rr [4])/6; - x_tt [0] += (x_tt [1] + x_tt [2] * 2 + x_tt [3] * 2 + x_tt [4])/6 + pow(gamma_t, 2) * h * u_t; + x_tt [0] += (x_tt [1] + x_tt [2] * 2 + x_tt [3] * 2 + x_tt [4])/6 + pow(gamma_t, 2) * h * Rand_vars[0]; x_tr [0] += (x_tr [1] + x_tr [2] * 2 + x_tr [3] * 2 + x_tr [4])/6; x_rt [0] += (x_rt [1] + x_rt [2] * 2 + x_rt [3] * 2 + x_rt [4])/6; x_rr [0] += (x_rr [1] + x_rr [2] * 2 + x_rr [3] * 2 + x_rr [4])/6; @@ -297,6 +297,10 @@ void Thalamic_Column::add_RK(double u_t) { m_h [0] += (m_h [1] + m_h [2] * 2 + m_h [3] * 2 + m_h [4])/6; m_h2 [0] += (m_h2 [1] + m_h2 [2] * 2 + m_h2 [3] * 2 + m_h2 [4])/6; P_h [0] += (P_h [1] + P_h [2] * 2 + P_h [3] * 2 + P_h [4])/6; + // generating the noise for the next iteration + for (unsigned i=0; i #include #include "macros.h" #include "parameters.h" +#include +#include +#include using std::vector; -// implementation of the cortical module after Zavaglia2006 +// typedefs for the RNG +typedef boost::mt11213b ENG; // Mersenne Twister +typedef boost::normal_distribution DIST; // Normal Distribution +typedef boost::variate_generator GEN; // Variate generator +// implementation of the thalamic module class Thalamic_Column { public: // Constructors @@ -21,7 +27,8 @@ class Thalamic_Column { m_KCa (_INIT(0.0)), m_CAN (_INIT(0.0)), m_h (_INIT(0.0)), m_h2 (_INIT(0.0)), P_h (_INIT(0.0)), N_tr (210), N_rt (410), N_rr (800) - {} + {MTRands = {{ENG(rand()), DIST (mphi_t, dphi_t)}, {ENG(rand()), DIST (mphi_t, dphi_t)}}; + Rand_vars = {MTRands[0](), MTRands[1]()};} // Constructors Thalamic_Column(double* Con) @@ -32,7 +39,8 @@ class Thalamic_Column { m_KCa (_INIT(0.0)), m_CAN (_INIT(0.0)), m_h (_INIT(0.0)), m_h2 (_INIT(0.0)), P_h (_INIT(0.0)), N_tr (Con[0]), N_rt (Con[1]), N_rr (Con[2]) - {} + {MTRands = {{ENG(rand()), DIST (mphi_t, dphi_t)}, {ENG(rand()), DIST (mphi_t, dphi_t)}}; + Rand_vars = {MTRands[0](), MTRands[1]()};} // firing rate functions double get_Qt (int) const; @@ -70,14 +78,14 @@ class Thalamic_Column { double I_CAN (int) const; // noise functions - double noise_xRK (int, double, double) const; + double noise_xRK (int) const; // ODE functions - void set_RK (int, double, double); - void add_RK (double); + void set_RK (int); + void add_RK (void); // function to extract the data - friend void get_data (int, Thalamic_Column&, _REPEAT(vector&, 11)); + friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 1)); private: // population variables @@ -107,5 +115,11 @@ class Thalamic_Column { double N_tr, // TC to RE loop N_rt, // RE to TC loop N_rr; // RE self loop + + // random number generators + vector MTRands; + + // container for random numbers + vector Rand_vars; }; diff --git a/Thalamus.cpp b/Thalamus.cpp index 32e5e99..f69f74e 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -4,7 +4,6 @@ #include #include "mex.h" #include "matrix.h" -#include "randoms.h" #include "Thalamic_Column.h" #include "saves.h" #include "ODE.h" @@ -17,8 +16,8 @@ extern const double h = sqrt(dt); // input arguments are a vector of length 8 with the connectivities and an integer setting the resolution of the grid void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - // Initializing the mersenne twister. - MTRand mtrand; + // Initializing the seeder. + srand(time(0)); // inputs double* Connectivity = mxGetPr (prhs[0]); @@ -28,51 +27,22 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { const int Time = (T+onset)*res; - // random input - vector u_t1 = rand_var(mtrand, (T+onset)*res, phi_st, phi_st); - vector u_t2 = rand_var(mtrand, (T+onset)*res, phi_st, phi_st); - - // random input with phase independent stimulation - //vector u_t1 = rand_inp(mtrand, Time, mphi_st, dphi_st, onset, var_stim); - //vector u_t2 = rand_inp(mtrand, Time, mphi_st, dphi_st, onset, var_stim); - - // Initializing the populations; Thalamic_Column Col(Connectivity); // setting up the data containers - vector Vt (T*res); - vector Vr (T*res); - vector Cat (T*res); - vector Car (T*res); - vector Phi_tr (T*res); - vector Phi_rt (T*res); - vector I_T_t (T*res); - vector I_T_r (T*res); - vector I_KCa (T*res); - vector I_CAN (T*res); - vector I_h (T*res); + mxArray* Vt = SetMexArray(1, T*red); int count = 0; for (int t=0; t=onset*res){ - get_data(count, Col, Vt, Vr, Cat, Car, Phi_tr, Phi_rt, I_T_t, I_T_r, I_KCa, I_CAN, I_h); + get_data(count, Col, Vt); ++count; } } // output of the simulation plhs[0] = getMexArray(Vt); - plhs[1] = getMexArray(Vr); - plhs[2] = getMexArray(Cat); - plhs[3] = getMexArray(Car); - plhs[4] = getMexArray(Phi_tr); - plhs[5] = getMexArray(Phi_rt); - plhs[6] = getMexArray(I_T_t); - plhs[7] = getMexArray(I_T_r); - plhs[8] = getMexArray(I_KCa); - plhs[9] = getMexArray(I_CAN); - plhs[10]= getMexArray(I_h); return; } diff --git a/Thalamus.mexa64 b/Thalamus.mexa64 deleted file mode 100755 index d52e854..0000000 Binary files a/Thalamus.mexa64 and /dev/null differ diff --git a/exportfig.m b/exportfig.m deleted file mode 100644 index 4af7c0e..0000000 --- a/exportfig.m +++ /dev/null @@ -1,1057 +0,0 @@ -function varargout = exportfig(varargin) -%EXPORTFIG Export a figure. -% EXPORTFIG(H, FILENAME) writes the figure H to FILENAME. H is -% a figure handle and FILENAME is a string that specifies the -% name of the output file. -% -% EXPORTFIG(H, FILENAME, OPTIONS) writes the figure H to FILENAME -% with options initially specified by the structure OPTIONS. The -% field names of OPTIONS must be legal parameters listed below -% and the field values must be legal values for the corresponding -% parameter. Default options can be set in releases prior to R12 -% by storing the OPTIONS structure in the root object's appdata -% with the command -% setappdata(0,'exportfigdefaults', OPTIONS) -% and for releases after R12 by setting the preference with the -% command -% setpref('exportfig', 'defaults', OPTIONS) -% -% EXPORTFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies -% parameters that control various characteristics of the output -% file. Any parameter value can be the string 'auto' which means -% the parameter uses the default factory behavior, overriding -% any other default for the parameter. -% -% Format Paramter: -% 'Format' a string -% specifies the output format. Defaults to 'eps'. For a -% list of export formats type 'help print'. -% 'Preview' one of the strings 'none', 'tiff' -% specifies a preview for EPS files. Defaults to 'none'. -% -% Size Parameters: -% 'Width' a positive scalar -% specifies the width in the figure's PaperUnits -% 'Height' a positive scalar -% specifies the height in the figure's PaperUnits -% 'Bounds' one of the strings 'tight', 'loose' -% specifies a tight or loose bounding box. Defaults to 'tight'. -% 'BoundsCode' one of the strings 'internal', 'mcode' -% specifies whether to use the internal EPS code or -% general m-code for computing bounding boxes. Defaults to -% 'internal'. -% 'Reference' an axes handle or a string -% specifies that the width and height parameters -% are relative to the given axes. If a string is -% specified then it must evaluate to an axes handle. -% -% Specifying only one dimension sets the other dimension -% so that the exported aspect ratio is the same as the -% figure's or reference axes' current aspect ratio. -% If neither dimension is specified the size defaults to -% the width and height from the figure's or reference -% axes' size. Tight bounding boxes are only computed for -% 2-D views and in that case the computed bounds enclose all -% text objects. -% -% Rendering Parameters: -% 'Color' one of the strings 'bw', 'gray', 'cmyk' -% 'bw' specifies that lines and text are exported in -% black and all other objects in grayscale -% 'gray' specifies that all objects are exported in grayscale -% 'rgb' specifies that all objects are exported in color -% using the RGB color space -% 'cmyk' specifies that all objects are exported in color -% using the CMYK color space -% 'Renderer' one of 'painters', 'zbuffer', 'opengl' -% specifies the renderer to use -% 'Resolution' a positive scalar -% specifies the resolution in dots-per-inch. -% 'LockAxes' one of 0 or 1 -% specifies that all axes limits and ticks should be fixed -% while exporting. -% -% The default color setting is 'bw'. -% -% Font Parameters: -% 'FontMode' one of the strings 'scaled', 'fixed' -% 'FontSize' a positive scalar -% in 'scaled' mode multiplies with the font size of each -% text object to obtain the exported font size -% in 'fixed' mode specifies the font size of all text -% objects in points -% 'DefaultFixedFontSize' a positive scalar -% in 'fixed' mode specified the default font size in -% points -% 'FontSizeMin' a positive scalar -% specifies the minimum font size allowed after scaling -% 'FontSizeMax' a positive scalar -% specifies the maximum font size allowed after scaling -% 'FontEncoding' one of the strings 'latin1', 'adobe' -% specifies the character encoding of the font -% 'SeparateText' one of 0 or 1 -% specifies that the text objects are stored in separate -% file as EPS with the base filename having '_t' appended. -% -% If FontMode is 'scaled' but FontSize is not specified then a -% scaling factor is computed from the ratio of the size of the -% exported figure to the size of the actual figure. -% -% The default 'FontMode' setting is 'scaled'. -% -% Line Width Parameters: -% 'LineMode' one of the strings 'scaled', 'fixed' -% 'LineWidth' a positive scalar -% 'DefaultFixedLineWidth' a positive scalar -% 'LineWidthMin' a positive scalar -% specifies the minimum line width allowed after scaling -% 'LineWidthMax' a positive scalar -% specifies the maximum line width allowed after scaling -% The semantics of 'Line' parameters are exactly the -% same as the corresponding 'Font' parameters, except that -% they apply to line widths instead of font sizes. -% -% Style Map Parameter: -% 'LineStyleMap' one of [], 'bw', or a function name or handle -% specifies how to map line colors to styles. An empty -% style map means styles are not changed. The style map -% 'bw' is a built-in mapping that maps lines with the same -% color to the same style and otherwise cycles through the -% available styles. A user-specified map is a function -% that takes as input a cell array of line objects and -% outputs a cell array of line style strings. The default -% map is []. -% -% Examples: -% exportfig(gcf,'fig1.eps','height',3); -% Exports the current figure to the file named 'fig1.eps' with -% a height of 3 inches (assuming the figure's PaperUnits is -% inches) and an aspect ratio the same as the figure's aspect -% ratio on screen. -% -% opts = struct('FontMode','fixed','FontSize',10,'height',3); -% exportfig(gcf, 'fig2.eps', opts, 'height', 5); -% Exports the current figure to 'fig2.eps' with all -% text in 10 point fonts and with height 5 inches. -% -% See also PREVIEWFIG, APPLYTOFIG, RESTOREFIG, PRINT. - -% Copyright 2000-2009 The MathWorks, Inc. - -if (nargin < 2) - error('Too few input arguments'); -end - -% exportfig(H, filename, [options,] ...) -H = varargin{1}; -if ~LocalIsHG(H,'figure') - error('First argument must be a handle to a figure.'); -end -filename = varargin{2}; -if ~ischar(filename) - error('Second argument must be a string.'); -end -paramPairs = {varargin{3:end}}; -if nargin > 2 - if isstruct(paramPairs{1}) - pcell = LocalToCell(paramPairs{1}); - paramPairs = {pcell{:}, paramPairs{2:end}}; - end -end -verstr = version; -majorver = str2num(verstr(1)); -defaults = []; -if majorver > 5 - if ispref('exportfig','defaults') - defaults = getpref('exportfig','defaults'); - end -elseif exist('getappdata') - defaults = getappdata(0,'exportfigdefaults'); -end -if ~isempty(defaults) - dcell = LocalToCell(defaults); - paramPairs = {dcell{:}, paramPairs{:}}; -end - -% Do some validity checking on param-value pairs -if (rem(length(paramPairs),2) ~= 0) - error(['Invalid input syntax. Optional parameters and values' ... - ' must be in pairs.']); -end - -auto.format = 'eps'; -auto.preview = 'none'; -auto.width = -1; -auto.height = -1; -auto.color = 'bw'; -auto.defaultfontsize=10; -auto.fontsize = -1; -auto.fontmode='scaled'; -auto.fontmin = 8; -auto.fontmax = 60; -auto.defaultlinewidth = 1.0; -auto.linewidth = -1; -auto.linemode=[]; -auto.linemin = 0.5; -auto.linemax = 100; -auto.fontencoding = 'latin1'; -auto.renderer = []; -auto.resolution = []; -auto.stylemap = []; -auto.applystyle = 0; -auto.refobj = -1; -auto.bounds = 'tight'; -auto.boundscode = 'internal'; -explicitbounds = 0; -auto.lockaxes = 1; -auto.separatetext = 0; -opts = auto; - -% Process param-value pairs -args = {}; -for k = 1:2:length(paramPairs) - param = lower(paramPairs{k}); - if ~ischar(param) - error('Optional parameter names must be strings'); - end - value = paramPairs{k+1}; - - switch (param) - case 'format' - opts.format = LocalCheckAuto(lower(value),auto.format); - if strcmp(opts.format,'preview') - error(['Format ''preview'' no longer supported. Use PREVIEWFIG' ... - ' instead.']); - end - case 'preview' - opts.preview = LocalCheckAuto(lower(value),auto.preview); - if ~strcmp(opts.preview,{'none','tiff'}) - error('Preview must be ''none'' or ''tiff''.'); - end - case 'width' - opts.width = LocalToNum(value, auto.width); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.width) - error('Width must be a numeric scalar > 0'); - end - end - case 'height' - opts.height = LocalToNum(value, auto.height); - if ~ischar(value) | ~strcmp(value,'auto') - if(~LocalIsPositiveScalar(opts.height)) - error('Height must be a numeric scalar > 0'); - end - end - case 'color' - opts.color = LocalCheckAuto(lower(value),auto.color); - if ~strcmp(opts.color,{'bw','gray','rgb','cmyk'}) - error('Color must be ''bw'', ''gray'',''rgb'' or ''cmyk''.'); - end - case 'fontmode' - opts.fontmode = LocalCheckAuto(lower(value),auto.fontmode); - if ~strcmp(opts.fontmode,{'scaled','fixed'}) - error('FontMode must be ''scaled'' or ''fixed''.'); - end - case 'fontsize' - opts.fontsize = LocalToNum(value,auto.fontsize); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.fontsize) - error('FontSize must be a numeric scalar > 0'); - end - end - case 'defaultfixedfontsize' - opts.defaultfontsize = LocalToNum(value,auto.defaultfontsize); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.defaultfontsize) - error('DefaultFixedFontSize must be a numeric scalar > 0'); - end - end - case 'fontsizemin' - opts.fontmin = LocalToNum(value,auto.fontmin); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.fontmin) - error('FontSizeMin must be a numeric scalar > 0'); - end - end - case 'fontsizemax' - opts.fontmax = LocalToNum(value,auto.fontmax); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.fontmax) - error('FontSizeMax must be a numeric scalar > 0'); - end - end - case 'fontencoding' - opts.fontencoding = LocalCheckAuto(lower(value),auto.fontencoding); - if ~strcmp(opts.fontencoding,{'latin1','adobe'}) - error('FontEncoding must be ''latin1'' or ''adobe''.'); - end - case 'linemode' - opts.linemode = LocalCheckAuto(lower(value),auto.linemode); - if ~strcmp(opts.linemode,{'scaled','fixed'}) - error('LineMode must be ''scaled'' or ''fixed''.'); - end - case 'linewidth' - opts.linewidth = LocalToNum(value,auto.linewidth); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.linewidth) - error('LineWidth must be a numeric scalar > 0'); - end - end - case 'defaultfixedlinewidth' - opts.defaultlinewidth = LocalToNum(value,auto.defaultlinewidth); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.defaultlinewidth) - error(['DefaultFixedLineWidth must be a numeric scalar >' ... - ' 0']); - end - end - case 'linewidthmin' - opts.linemin = LocalToNum(value,auto.linemin); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.linemin) - error('LineWidthMin must be a numeric scalar > 0'); - end - end - case 'linewidthmax' - opts.linemax = LocalToNum(value,auto.linemax); - if ~ischar(value) | ~strcmp(value,'auto') - if ~LocalIsPositiveScalar(opts.linemax) - error('LineWidthMax must be a numeric scalar > 0'); - end - end - case 'linestylemap' - opts.stylemap = LocalCheckAuto(value,auto.stylemap); - case 'renderer' - opts.renderer = LocalCheckAuto(lower(value),auto.renderer); - if ~ischar(value) | ~strcmp(value,'auto') - if ~strcmp(opts.renderer,{'painters','zbuffer','opengl'}) - error(['Renderer must be ''painters'', ''zbuffer'' or' ... - ' ''opengl''.']); - end - end - case 'resolution' - opts.resolution = LocalToNum(value,auto.resolution); - if ~ischar(value) | ~strcmp(value,'auto') - if ~(isnumeric(value) & (prod(size(value)) == 1) & (value >= 0)); - error('Resolution must be a numeric scalar >= 0'); - end - end - case 'applystyle' % means to apply the options and not export - opts.applystyle = 1; - case 'reference' - if ischar(value) - if strcmp(value,'auto') - opts.refobj = auto.refobj; - else - opts.refobj = eval(value); - end - else - opts.refobj = value; - end - if ~LocalIsHG(opts.refobj,'axes') - error('Reference object must evaluate to an axes handle.'); - end - case 'bounds' - opts.bounds = LocalCheckAuto(lower(value),auto.bounds); - explicitbounds = 1; - if ~strcmp(opts.bounds,{'tight','loose'}) - error('Bounds must be ''tight'' or ''loose''.'); - end - case 'boundscode' - opts.boundscode = LocalCheckAuto(lower(value),auto.boundscode); - if ~strcmp(opts.boundscode,{'internal','mcode'}) - error('BoundsCode must be ''internal'' or ''mcode''.'); - end - case 'lockaxes' - opts.lockaxes = LocalToNum(value,auto.lockaxes); - case 'separatetext' - opts.separatetext = LocalToNum(value,auto.separatetext); - otherwise - error(['Unrecognized option ' param '.']); - end -end - -% make sure figure is up-to-date -drawnow; - -allLines = findall(H, 'type', 'line'); -allText = findall(H, 'type', 'text'); -allAxes = findall(H, 'type', 'axes'); -allImages = findall(H, 'type', 'image'); -allLights = findall(H, 'type', 'light'); -allPatch = findall(H, 'type', 'patch'); -allSurf = findall(H, 'type', 'surface'); -allRect = findall(H, 'type', 'rectangle'); -allFont = [allText; allAxes]; -allColor = [allLines; allText; allAxes; allLights]; -allMarker = [allLines; allPatch; allSurf]; -allEdge = [allPatch; allSurf]; -allCData = [allImages; allSurf]; - -% need to explicitly get ticks for some invisible figures -get(allAxes,'yticklabel'); - -old.objs = {}; -old.prop = {}; -old.values = {}; - -% Process format -if strncmp(opts.format,'eps',3) & ~strcmp(opts.preview,'none') - args = {args{:}, ['-' opts.preview]}; -end - -hadError = 0; -oldwarn = warning; -try - - % lock axes limits, ticks and labels if requested - if opts.lockaxes - old = LocalManualAxesMode(old, allAxes, 'TickMode'); - old = LocalManualAxesMode(old, allAxes, 'TickLabelMode'); - old = LocalManualAxesMode(old, allAxes, 'LimMode'); - end - - % Process size parameters - figurePaperUnits = get(H, 'PaperUnits'); - oldFigureUnits = get(H, 'Units'); - oldFigPos = get(H,'Position'); - set(H, 'Units', figurePaperUnits); - figPos = get(H,'Position'); - refsize = figPos(3:4); - if opts.refobj ~= -1 - oldUnits = get(opts.refobj, 'Units'); - set(opts.refobj, 'Units', figurePaperUnits); - r = get(opts.refobj, 'Position'); - refsize = r(3:4); - set(opts.refobj, 'Units', oldUnits); - end - aspectRatio = refsize(1)/refsize(2); - if (opts.width == -1) & (opts.height == -1) - opts.width = refsize(1); - opts.height = refsize(2); - elseif (opts.width == -1) - opts.width = opts.height * aspectRatio; - elseif (opts.height == -1) - opts.height = opts.width / aspectRatio; - end - wscale = opts.width/refsize(1); - hscale = opts.height/refsize(2); - sizescale = min(wscale,hscale); - old = LocalPushOldData(old,H,'PaperPositionMode', ... - get(H,'PaperPositionMode')); - set(H, 'PaperPositionMode', 'auto'); - newPos = [figPos(1) figPos(2)+figPos(4)*(1-hscale) ... - wscale*figPos(3) hscale*figPos(4)]; - set(H, 'Position', newPos); - set(H, 'Units', oldFigureUnits); - - old = LocalPushOldData(old,H,'Color', ... - get(H,'Color')); - set(H,'Color','w'); - - % process line-style map - if ~isempty(opts.stylemap) & ~isempty(allLines) - oldlstyle = LocalGetAsCell(allLines,'LineStyle'); - old = LocalPushOldData(old, allLines, {'LineStyle'}, ... - oldlstyle); - newlstyle = oldlstyle; - if ischar(opts.stylemap) & strcmpi(opts.stylemap,'bw') - newlstyle = LocalMapColorToStyle(allLines); - else - try - newlstyle = feval(opts.stylemap,allLines); - catch - warning(['Skipping stylemap. ' lasterr]); - end - end - set(allLines,{'LineStyle'},newlstyle); - end - - % Process rendering parameters - switch (opts.color) - case {'bw', 'gray'} - if ~strcmp(opts.color,'bw') & strncmp(opts.format,'eps',3) - opts.format = [opts.format 'c']; - end - args = {args{:}, ['-d' opts.format]}; - - %compute and set gray colormap - oldcmap = get(H,'Colormap'); - newgrays = 0.30*oldcmap(:,1) + 0.59*oldcmap(:,2) + 0.11*oldcmap(:,3); - newcmap = [newgrays newgrays newgrays]; - old = LocalPushOldData(old, H, 'Colormap', oldcmap); - set(H, 'Colormap', newcmap); - - %compute and set ColorSpec and CData properties - old = LocalUpdateColors(allColor, 'color', old); - old = LocalUpdateColors(allAxes, 'xcolor', old); - old = LocalUpdateColors(allAxes, 'ycolor', old); - old = LocalUpdateColors(allAxes, 'zcolor', old); - old = LocalUpdateColors(allMarker, 'MarkerEdgeColor', old); - old = LocalUpdateColors(allMarker, 'MarkerFaceColor', old); - old = LocalUpdateColors(allEdge, 'EdgeColor', old); - old = LocalUpdateColors(allEdge, 'FaceColor', old); - old = LocalUpdateColors(allCData, 'CData', old); - - if strcmp(opts.color,'bw') - lcolor = LocalGetAsCell(allLines,'color'); - n = length(lcolor); - for k=1:n - if (lcolor{k}(1) < 1-eps) - set(allLines(k),'color','k'); - end - end - end - case {'rgb','cmyk'} - if strncmp(opts.format,'eps',3) - opts.format = [opts.format 'c']; - args = {args{:}, ['-d' opts.format]}; - if strcmp(opts.color,'cmyk') - args = {args{:}, '-cmyk'}; - end - else - args = {args{:}, ['-d' opts.format]}; - end - otherwise - error('Invalid Color parameter'); - end - if (~isempty(opts.renderer)) - args = {args{:}, ['-' opts.renderer]}; - end - if (~isempty(opts.resolution)) | ~strncmp(opts.format,'eps',3) - if isempty(opts.resolution) - opts.resolution = 0; - end - args = {args{:}, ['-r' int2str(opts.resolution)]}; - end - - % Process font parameters - if ~isempty(opts.fontmode) - oldfonts = LocalGetAsCell(allFont,'FontSize'); - oldfontunits = LocalGetAsCell(allFont,'FontUnits'); - set(allFont,'FontUnits','points'); - switch (opts.fontmode) - case 'fixed' - if (opts.fontsize == -1) - set(allFont,'FontSize',opts.defaultfontsize); - else - set(allFont,'FontSize',opts.fontsize); - end - case 'scaled' - if (opts.fontsize == -1) - scale = sizescale; - else - scale = opts.fontsize; - end - newfonts = LocalScale(oldfonts,scale,opts.fontmin,opts.fontmax); - set(allFont,{'FontSize'},newfonts); - otherwise - error('Invalid FontMode parameter'); - end - old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts); - old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits); - end - if strcmp(opts.fontencoding,'adobe') & strncmp(opts.format,'eps',3) - args = {args{:}, '-adobecset'}; - end - - % Process line parameters - if ~isempty(opts.linemode) - oldlines = LocalGetAsCell(allMarker,'LineWidth'); - old = LocalPushOldData(old, allMarker, {'LineWidth'}, oldlines); - switch (opts.linemode) - case 'fixed' - if (opts.linewidth == -1) - set(allMarker,'LineWidth',opts.defaultlinewidth); - else - set(allMarker,'LineWidth',opts.linewidth); - end - case 'scaled' - if (opts.linewidth == -1) - scale = sizescale; - else - scale = opts.linewidth; - end - newlines = LocalScale(oldlines, scale, opts.linemin, opts.linemax); - set(allMarker,{'LineWidth'},newlines); - end - end - - % adjust figure bounds to surround axes - if strcmp(opts.bounds,'tight') - if (~strncmp(opts.format,'eps',3) & LocalHas3DPlot(allAxes)) | ... - (strncmp(opts.format,'eps',3) & opts.separatetext) - if (explicitbounds == 1) - warning(['Cannot compute ''tight'' bounds. Using ''loose''' ... - ' bounds.']); - end - opts.bounds = 'loose'; - end - end - warning('off'); - if ~isempty(allAxes) - usemcode = 1; - if strncmp(opts.format,'eps',3) - if strcmp(opts.boundscode,'internal') - if strcmp(opts.bounds,'loose') - args = {args{:}, '-loose'}; - end - old = LocalPushOldData(old,H,'Position', oldFigPos); - usemcode = 0; - else - usemcode = 1; - end - end - if strcmp(opts.bounds,'tight') & usemcode - oldaunits = LocalGetAsCell(allAxes,'Units'); - oldapos = LocalGetAsCell(allAxes,'Position'); - oldtunits = LocalGetAsCell(allText,'units'); - oldtpos = LocalGetAsCell(allText,'Position'); - set(allAxes,'units','points'); - apos = LocalGetAsCell(allAxes,'Position'); - oldunits = get(H,'Units'); - set(H,'units','points'); - origfr = get(H,'position'); - fr = []; - for k=1:length(allAxes) - if ~strcmpi(get(allAxes(k),'Tag'),'legend') - axesR = apos{k}; - r = LocalAxesTightBoundingBox(axesR, allAxes(k)); - if ~isempty(r) - r(1:2) = r(1:2) + axesR(1:2); - fr = LocalUnionRect(fr,r); - end - end - end - if isempty(fr) - fr = [0 0 origfr(3:4)]; - end - for k=1:length(allAxes) - ax = allAxes(k); - r = apos{k}; - r(1:2) = r(1:2) - fr(1:2); - set(ax,'Position',r); - end - old = LocalPushOldData(old, allAxes, {'Position'}, oldapos); - old = LocalPushOldData(old, allText, {'Position'}, oldtpos); - old = LocalPushOldData(old, allText, {'Units'}, oldtunits); - old = LocalPushOldData(old, allAxes, {'Units'}, oldaunits); - old = LocalPushOldData(old, H, 'Position', oldFigPos); - old = LocalPushOldData(old, H, 'Units', oldFigureUnits); - r = [origfr(1) origfr(2)+origfr(4)-fr(4) fr(3:4)]; - set(H,'Position',r); - end - end - - % Process text in a separate file if needed - if opts.separatetext & ~opts.applystyle - % First hide all text and export - oldtvis = LocalGetAsCell(allText,'visible'); - set(allText,'visible','off'); - oldax = LocalGetAsCell(allAxes,'XTickLabel',1); - olday = LocalGetAsCell(allAxes,'YTickLabel',1); - oldaz = LocalGetAsCell(allAxes,'ZTickLabel',1); - null = cell(length(oldax),1); - [null{:}] = deal([]); - set(allAxes,{'XTickLabel'},null); - set(allAxes,{'YTickLabel'},null); - set(allAxes,{'ZTickLabel'},null); - print(H, filename, args{:}); - set(allText,{'Visible'},oldtvis); - set(allAxes,{'XTickLabel'},oldax); - set(allAxes,{'YTickLabel'},olday); - set(allAxes,{'ZTickLabel'},oldaz); - % Now hide all non-text and export as eps in painters - [path, name, ext] = fileparts(filename); - tfile = fullfile(path,[name '_t.eps']); - tfile2 = fullfile(path,[name '_t2.eps']); - foundRenderer = 0; - for k=1:length(args) - if strncmp('-d',args{k},2) - args{k} = '-deps'; - elseif strncmp('-zbuffer',args{k},8) | ... - strncmp('-opengl', args{k},6) - args{k} = '-painters'; - foundRenderer = 1; - end - end - if ~foundRenderer - args = {args{:}, '-painters'}; - end - allNonText = [allLines; allLights; allPatch; ... - allImages; allSurf; allRect]; - oldvis = LocalGetAsCell(allNonText,'visible'); - oldc = LocalGetAsCell(allAxes,'color'); - oldaxg = LocalGetAsCell(allAxes,'XGrid'); - oldayg = LocalGetAsCell(allAxes,'YGrid'); - oldazg = LocalGetAsCell(allAxes,'ZGrid'); - [null{:}] = deal('off'); - set(allAxes,{'XGrid'},null); - set(allAxes,{'YGrid'},null); - set(allAxes,{'ZGrid'},null); - set(allNonText,'Visible','off'); - set(allAxes,'Color','none'); - print(H, tfile2, args{:}); - set(allNonText,{'Visible'},oldvis); - set(allAxes,{'Color'},oldc); - set(allAxes,{'XGrid'},oldaxg); - set(allAxes,{'YGrid'},oldayg); - set(allAxes,{'ZGrid'},oldazg); - %hack up the postscript file - fid1 = fopen(tfile,'w'); - fid2 = fopen(tfile2,'r'); - line = fgetl(fid2); - while ischar(line) - if strncmp(line,'%%Title',7) - fprintf(fid1,'%s\n',['%%Title: ', tfile]); - elseif (length(line) < 3) - fprintf(fid1,'%s\n',line); - elseif ~strcmp(line(end-2:end),' PR') & ... - ~strcmp(line(end-1:end),' L') - fprintf(fid1,'%s\n',line); - end - line = fgetl(fid2); - end - fclose(fid1); - fclose(fid2); - delete(tfile2); - - elseif ~opts.applystyle - drawnow; - print(H, filename, args{:}); - end - warning(oldwarn); - -catch - warning(oldwarn); - hadError = 1; -end - -% Restore figure settings -if opts.applystyle - varargout{1} = old; -else - for n=1:length(old.objs) - if ~iscell(old.values{n}) & iscell(old.prop{n}) - old.values{n} = {old.values{n}}; - end - set(old.objs{n}, old.prop{n}, old.values{n}); - end -end - -if hadError - error(deblank(lasterr)); -end - -% -% Local Functions -% - -function outData = LocalPushOldData(inData, objs, prop, values) -outData.objs = {objs, inData.objs{:}}; -outData.prop = {prop, inData.prop{:}}; -outData.values = {values, inData.values{:}}; - -function cellArray = LocalGetAsCell(fig,prop,allowemptycell); -cellArray = get(fig,prop); -if nargin < 3 - allowemptycell = 0; -end -if ~iscell(cellArray) & (allowemptycell | ~isempty(cellArray)) - cellArray = {cellArray}; -end - -function newArray = LocalScale(inArray, scale, minv, maxv) -n = length(inArray); -newArray = cell(n,1); -for k=1:n - newArray{k} = min(maxv,max(minv,scale*inArray{k}(1))); -end - -function gray = LocalMapToGray1(color) -gray = color; -if ischar(color) - switch color(1) - case 'y' - color = [1 1 0]; - case 'm' - color = [1 0 1]; - case 'c' - color = [0 1 1]; - case 'r' - color = [1 0 0]; - case 'g' - color = [0 1 0]; - case 'b' - color = [0 0 1]; - case 'w' - color = [1 1 1]; - case 'k' - color = [0 0 0]; - end -end -if ~ischar(color) - gray = 0.30*color(1) + 0.59*color(2) + 0.11*color(3); -end - -function newArray = LocalMapToGray(inArray); -n = length(inArray); -newArray = cell(n,1); -for k=1:n - color = inArray{k}; - if ~isempty(color) - color = LocalMapToGray1(color); - end - if isempty(color) | ischar(color) - newArray{k} = color; - else - newArray{k} = [color color color]; - end -end - -function newArray = LocalMapColorToStyle(inArray); -inArray = LocalGetAsCell(inArray,'Color'); -n = length(inArray); -newArray = cell(n,1); -styles = {'-','--',':','-.'}; -uniques = []; -nstyles = length(styles); -for k=1:n - gray = LocalMapToGray1(inArray{k}); - if isempty(gray) | ischar(gray) | gray < .05 - newArray{k} = '-'; - else - if ~isempty(uniques) & any(gray == uniques) - ind = find(gray==uniques); - else - uniques = [uniques gray]; - ind = length(uniques); - end - newArray{k} = styles{mod(ind-1,nstyles)+1}; - end -end - -function newArray = LocalMapCData(inArray); -n = length(inArray); -newArray = cell(n,1); -for k=1:n - color = inArray{k}; - if (ndims(color) == 3) & isa(color,'double') - gray = 0.30*color(:,:,1) + 0.59*color(:,:,2) + 0.11*color(:,:,3); - color(:,:,1) = gray; - color(:,:,2) = gray; - color(:,:,3) = gray; - end - newArray{k} = color; -end - -function outData = LocalUpdateColors(inArray, prop, inData) -value = LocalGetAsCell(inArray,prop); -outData.objs = {inData.objs{:}, inArray}; -outData.prop = {inData.prop{:}, {prop}}; -outData.values = {inData.values{:}, value}; -if (~isempty(value)) - if strcmp(prop,'CData') - value = LocalMapCData(value); - else - value = LocalMapToGray(value); - end - set(inArray,{prop},value); -end - -function bool = LocalIsPositiveScalar(value) -bool = isnumeric(value) & ... - prod(size(value)) == 1 & ... - value > 0; - -function value = LocalToNum(value,auto) -if ischar(value) - if strcmp(value,'auto') - value = auto; - else - value = str2num(value); - end -end - -%convert a struct to {field1,val1,field2,val2,...} -function c = LocalToCell(s) -f = fieldnames(s); -v = struct2cell(s); -opts = cell(2,length(f)); -opts(1,:) = f; -opts(2,:) = v; -c = {opts{:}}; - -function c = LocalIsHG(obj,hgtype) -c = 0; -if (length(obj) == 1) & ishandle(obj) - c = strcmp(get(obj,'type'),hgtype); -end - -function c = LocalHas3DPlot(a) -zticks = LocalGetAsCell(a,'ZTickLabel'); -c = 0; -for k=1:length(zticks) - if ~isempty(zticks{k}) - c = 1; - return; - end -end - -function r = LocalUnionRect(r1,r2) -if isempty(r1) - r = r2; -elseif isempty(r2) - r = r1; -elseif max(r2(3:4)) > 0 - left = min(r1(1),r2(1)); - bot = min(r1(2),r2(2)); - right = max(r1(1)+r1(3),r2(1)+r2(3)); - top = max(r1(2)+r1(4),r2(2)+r2(4)); - r = [left bot right-left top-bot]; -else - r = r1; -end - -function r = LocalAxesTightBoundingBox(axesR, a) -if strcmp(get(a,'handlevisibility'),'on') - r = get(a,'position'); -else - r = []; -end -atext = findall(a,'type','text','visible','on'); -if ~isempty(atext) - set(atext,'units','points'); - res=LocalGetAsCell(atext,'extent'); - for n=1:length(atext) - r = LocalUnionRect(r,res{n}); - end -end -if strcmp(get(a,'visible'),'on') - r = LocalUnionRect(r,[0 0 axesR(3:4)]); - oldunits = get(a,'fontunits'); - set(a,'fontunits','points'); - label = text(0,0,'','parent',a,... - 'units','points',... - 'fontsize',get(a,'fontsize'),... - 'fontname',get(a,'fontname'),... - 'fontweight',get(a,'fontweight'),... - 'fontangle',get(a,'fontangle'),... - 'visible','off'); - fs = get(a,'fontsize'); - - % handle y axis tick labels - ry = [0 0 0 axesR(4)]; - ylabs = get(a,'yticklabels'); - yticks = get(a,'ytick'); - maxw = 0; - if ~isempty(ylabs) - for n=1:size(ylabs,1) - if strcmp(get(a,'yscale'),'log') - set(label,'string',['10^' ylabs(n,:)]); - else - set(label,'string',ylabs(n,:)); - end - ext = get(label,'extent'); - maxw = max(maxw,ext(3)); - end - if strcmp(get(a,'yaxislocation'),'left') - ry(1) = -(maxw+5); - else - ry(1) = axesR(3); - end - ry(2) = -ext(4)/2; - ry(3) = maxw+5; - ry(4) = ry(4) + ext(4); - if ~isempty(yticks) - if ~strcmp(get(a,'yscale'),'log') & ... - ((str2num(ylabs(1,:)) ~= yticks(1)) | ... - (str2num(ylabs(end,:)) ~= yticks(end))) - set(label,'string','10^1'); - ext = get(label,'extent'); - if strcmp(get(a,'xaxislocation'),'bottom') - ry(4) = ry(4) + ext(4); - end - end - end - r = LocalUnionRect(r,ry); - end - - % handle x axis tick labels - rx = [0 0 axesR(3) 0]; - xlabs = get(a,'xticklabels'); - xticks = get(a,'xtick'); - if ~isempty(xlabs) - if strcmp(get(a,'xscale'),'log') - set(label,'string',['10^' xlabs(1,:)]); - else - set(label,'string',xlabs(1,:)); - end - ext1 = get(label,'extent'); - rx(1) = -ext1(3)/2; - if strcmp(get(a,'xscale'),'log') - set(label,'string',['10^' xlabs(size(xlabs,1),:)]); - else - set(label,'string',xlabs(size(xlabs,1),:)); - end - ext2 = get(label,'extent'); - rx(4) = max(ext1(4),ext2(4)); - rx(2) = -rx(4); - rx(3) = axesR(3) + (ext2(3) + ext1(3))/2; - if strcmp(get(a,'xaxislocation'),'top') - rx(2) = rx(2) + axesR(4); - end - if ~isempty(xticks) - if ~strcmp(get(a,'xscale'),'log') & ... - ((str2num(xlabs(1,:)) ~= xticks(1)) | ... - (str2num(xlabs(end,:)) ~= xticks(end))) - set(label,'string','10^1'); - ext = get(label,'extent'); - if strcmp(get(a,'xaxislocation'),'bottom') - rx(4) = rx(4) + ext(4)+5; - rx(2) = rx(2) - ext(4)-5; - else - rx(4) = rx(4) + ext(4)+5 + axesR(4); - rx(2) = rx(2) - ext(4)-5 - axesR(4); - end - end - end - r = LocalUnionRect(r,rx); - end - set(a,'fontunits',oldunits); - delete(label); -end - -function c = LocalManualAxesMode(old, allAxes, base) -xs = ['X' base]; -ys = ['Y' base]; -zs = ['Z' base]; -xlog = LocalGetAsCell(allAxes,'xscale'); -ylog = LocalGetAsCell(allAxes,'yscale'); -zlog = LocalGetAsCell(allAxes,'zscale'); -nonlogxflag = logical([]); -nonlogyflag = logical([]); -nonlogzflag = logical([]); -for k = 1:length(xlog) - nonlogxflag(k) = logical(~strcmp(xlog{k},'log')); - nonlogyflag(k) = logical(~strcmp(ylog{k},'log')); - nonlogzflag(k) = logical(~strcmp(zlog{k},'log')); -end -allNonLogXAxes = allAxes(nonlogxflag); -allNonLogYAxes = allAxes(nonlogyflag); -allNonLogZAxes = allAxes(nonlogzflag); -oldXMode = LocalGetAsCell(allNonLogXAxes,xs); -oldYMode = LocalGetAsCell(allNonLogYAxes,ys); -oldZMode = LocalGetAsCell(allNonLogZAxes,zs); -old = LocalPushOldData(old, allNonLogXAxes, {xs}, oldXMode); -old = LocalPushOldData(old, allNonLogYAxes, {ys}, oldYMode); -old = LocalPushOldData(old, allNonLogZAxes, {zs}, oldZMode); -set(allNonLogXAxes,xs,'manual'); -set(allNonLogYAxes,ys,'manual'); -set(allNonLogZAxes,zs,'manual'); -c = old; - -function val = LocalCheckAuto(val, auto) -if ischar(val) & strcmp(val,'auto') - val = auto; -end diff --git a/parameters.h b/parameters.h index 8cfde7d..61bd776 100644 --- a/parameters.h +++ b/parameters.h @@ -114,7 +114,8 @@ /*****************************************************************************************************/ /********************************** noise parameters **********************************/ /*****************************************************************************************************/ -#define phi_st 0.5E-3 +#define mphi_t 0E-3 +#define dphi_t 0.5E-3 #define phi_inp 0.0 /*****************************************************************************************************/ /********************************** end **********************************/ diff --git a/poster_plot.m b/poster_plot.m deleted file mode 100644 index 369aa28..0000000 --- a/poster_plot.m +++ /dev/null @@ -1,19 +0,0 @@ -Con = [10; % N_tr - 10; % N_rt - 50]; % N_rr -T = 10; % duration of the simulation -onset = 5; - -[Vt,~] = Thalamus(Con, T, onset); - -L = max(size(Vt)); -timeaxis = linspace(0,T,L); - -figure(1) -plot(timeaxis, Vt); -title('Membrane voltage of the thalamic relay population') -xlabel('time in s') -ylabel('V in mV') - -tightfig(); -savefig('Spindle', 'png') \ No newline at end of file diff --git a/randoms.h b/randoms.h deleted file mode 100644 index 821660d..0000000 --- a/randoms.h +++ /dev/null @@ -1,94 +0,0 @@ -/*****************************************************************************************************/ -/****************** functions for generation of noise and its evaluation in RK4 **********************/ -/*****************************************************************************************************/ -#pragma once -#include -#include -#include "MersenneTwister.h" -using std::vector; - -/*****************************************************************************************************/ -/***************************** parameters for integration of SRK4 **************************/ -/*****************************************************************************************************/ -// vectors needed for stochastic RK -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ - - -/*****************************************************************************************************/ -/********************************** random input generation **********************************/ -/*****************************************************************************************************/ -// function that creates a vector of random numbers with length N and SD sigma -vector rand_var (MTRand mtrand, int length, double mean, double sd) { - // length length of the simulation - // mean mean noise without stimulation - // sd its sd - - mtrand.seed(); - vector dW(length); - for (auto i=0; i rand_inp (MTRand mtrand, int length, double mean, double sd, int onset, double* var_stim) { - // res steps per secon - // length total length of the simulation - // onset time until data is saved - // var_stim[0] strength of stimulation - // var_stim[0] time between stimuli - // var_stim[0] length of stimuli - // mean mean noise without stimulation - // sd its sd - - extern const double res; - int mode = 0; - int count = 0; - - double strength = mean + var_stim[0]; - - // initializing the rando number generator - mtrand.seed(); - // generating random noise for the whole time length - vector dW(length); - for (auto i=0; i': Set compression for non-indexed bitmaps in PDFs - -% 0: lossless; 0.1: high quality; 0.5: medium; 1: high compression. -% '-r': Set resolution. -% '-crop': Removes points and line segments outside the viewing area -- permanently. -% Only use this on figures where many points and/or line segments are outside -% the area zoomed in to. This option will result in smaller vector files (has no -% effect on pixel files). -% '-dbg': Displays gs command line(s). -% -% EXAMPLE: -% savefig('nicefig', 'pdf', 'jpeg', '-cmyk', '-c0.1', '-r250'); -% Saves the current figure to nicefig.pdf and nicefig.png, both in cmyk and at 250 dpi, -% with high quality lossy compression. -% -% REQUIREMENT: Ghostscript. Version 8.57 works, probably older versions too, but '-dEPSCrop' -% must be supported. I think version 7.32 or newer is ok. -% -% HISTORY: -% Version 1.0, 2006-04-20. -% Version 1.1, 2006-04-27: -% - No 'epstopdf' stuff anymore! Using '-dEPSCrop' option in gs instead! -% Version 1.2, 2006-05-02: -% - Added a '-dbg' option (see options, above). -% - Now looks for a global variable 'savefig_defaults' (see options, above). -% - More detailed Ghostscript options (user will not really notice). -% - Warns when there is no device for a file-type/color-model combination. -% Version 1.3, 2006-06-06: -% - Added a check to see if there actually is a figure handle. -% - Now works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too). -% - Now compatible with Ghostscript 8.54, released 2006-06-01. -% Version 1.4, 2006-07-20: -% - Added an option '-soft' that enables anti-aliasing on pixel graphics (on by default). -% - Added an option '-hard' that don't do anti-aliasing on pixel graphics. -% Version 1.5, 2006-07-27: -% - Fixed a bug when calling with a figure handle argument. -% Version 1.6, 2006-07-28: -% - Added a crop option, see above. -% Version 1.7, 2007-03-31: -% - Fixed bug: calling print with invalid renderer value '-none'. -% - Removed GhostScript argument '-dUseCIEColor' as it sometimes discoloured things. -% Version 1.8, 2008-01-03: -% - Added MacIntel: 'MACI'. -% - Added 64bit PC (I think, can't test it myself). -% - Added option '-nointerpolate' (use it to prevent blurring of pixelated). -% - Removed '-hard' and '-soft'. Use '-nointerpolate' for '-hard', default for '-soft'. -% - Fixed the gs 8.57 warning on UseCIEColor (it's now set). -% - Added '-gray' for pdf, but gs 8.56 or newer is needed. -% - Added '-gray' and '-cmyk' for eps, but you a fairly recent gs might be needed. -% Version 1.9, 2008-07-27: -% - Added lossless compression, see option '-lossless', above. Works on most formats. -% - Added lossy compression, see options '-c...', above. Works on 'pdf'. -% Thanks to Olly Woodford for idea and implementation! -% - Removed option '-nointerpolate' -- now savefig never interpolates. -% - Fixed a few small bugs and removed some mlint comments. -% Version 2.0, 2008-11-07: -% - Added the possibility to include fonts into eps or pdf. -% -% TO DO: (Need Ghostscript support for these, so don't expect anything soon...) -% - svg output. -% - '-cmyk' for 'jpeg' and 'png'. -% - Preview in 'eps'. -% - Embedded vector fonts, not bitmap, in 'eps'. -% -% Copyright (C) Peder Axensten (peder at axensten dot se), 2006. - -% KEYWORDS: eps, pdf, jpg, jpeg, png, tiff, eps2pdf, epstopdf, ghostscript -% -% INSPIRATION: eps2pdf (5782), eps2xxx (6858) -% -% REQUIREMENTS: Works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too). - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - op_dbg= false; % Default value. - - % Compression - compr= [' -dUseFlateCompression=true -dLZWEncodePages=true -dCompatibilityLevel=1.6' ... - ' -dAutoFilterColorImages=false -dAutoFilterGrayImages=false ' ... - ' -dColorImageFilter=%s -dGrayImageFilter=%s']; % Compression. - lossless= sprintf (compr, '/FlateEncode', '/FlateEncode'); - lossy= sprintf (compr, '/DCTEncode', '/DCTEncode' ); - lossy= [lossy ' -c ".setpdfwrite << /ColorImageDict << /QFactor %g ' ... - '/Blend 1 /HSample [%s] /VSample [%s] >> >> setdistillerparams"']; - - % Create gs command. - cmdEnd= ' -sDEVICE=%s -sOutputFile="%s"'; % Essential. - epsCmd= ''; - epsCmd= [epsCmd ' -dSubsetFonts=true -dNOPLATFONTS']; % Future support? - epsCmd= [epsCmd ' -dUseCIEColor=true -dColorConversionStrategy=/UseDeviceIndependentColor']; - epsCmd= [epsCmd ' -dProcessColorModel=/%s']; % Color conversion. - pdfCmd= [epsCmd ' -dAntiAliasColorImages=false' cmdEnd]; - epsCmd= [epsCmd cmdEnd]; - - % Get file name. - if((nargin < 1) || isempty(fname) || ~ischar(fname)) % Check file name. - error('No file name specified.'); - end - [pathstr, namestr] = fileparts(fname); - if(isempty(pathstr)), fname= fullfile(cd, namestr); end - - % Get handle. - fighdl= get(0, 'CurrentFigure'); % See gcf. % Get figure handle. - if((nargin >= 2) && (numel(varargin{1}) == 1) && isnumeric(varargin{1})) - fighdl= varargin{1}; - varargin= {varargin{2:end}}; - end - if(isempty(fighdl)), error('There is no figure to save!?'); end - set(fighdl, 'Units', 'centimeters') % Set paper stuff. - sz= get(fighdl, 'Position'); - sz(1:2)= 0; - set(fighdl, 'PaperUnits', 'centimeters', 'PaperSize', sz(3:4), 'PaperPosition', sz); - - % Set up the various devices. - % Those commented out are not yet supported by gs (nor by savefig). - % pdf-cmyk works due to the Matlab '-cmyk' export being carried over from eps to pdf. - device.eps.rgb= sprintf(epsCmd, 'DeviceRGB', 'epswrite', [fname '.eps']); - device.jpeg.rgb= sprintf(cmdEnd, 'jpeg', [fname '.jpeg']); -% device.jpeg.cmyk= sprintf(cmdEnd, 'jpegcmyk', [fname '.jpeg']); - device.jpeg.gray= sprintf(cmdEnd, 'jpeggray', [fname '.jpeg']); - device.pdf.rgb= sprintf(pdfCmd, 'DeviceRGB', 'pdfwrite', [fname '.pdf']); - device.pdf.cmyk= sprintf(pdfCmd, 'DeviceCMYK', 'pdfwrite', [fname '.pdf']); - device.pdf.gray= sprintf(pdfCmd, 'DeviceGray', 'pdfwrite', [fname '.pdf']); - device.png.rgb= sprintf(cmdEnd, 'png16m', [fname '.png']); -% device.png.cmyk= sprintf(cmdEnd, 'png???', [fname '.png']); - device.png.gray= sprintf(cmdEnd, 'pnggray', [fname '.png']); - device.tiff.rgb= sprintf(cmdEnd, 'tiff24nc', [fname '.tiff']); - device.tiff.cmyk= sprintf(cmdEnd, 'tiff32nc', [fname '.tiff']); - device.tiff.gray= sprintf(cmdEnd, 'tiffgray', [fname '.tiff']); - - % Get options. - global savefig_defaults; % Add global defaults. - if( iscellstr(savefig_defaults)), varargin= {savefig_defaults{:}, varargin{:}}; - elseif(ischar(savefig_defaults)), varargin= {savefig_defaults, varargin{:}}; - end - varargin= {'-r300', '-lossless', '-rgb', varargin{:}}; % Add defaults. - res= ''; - types= {}; - fonts= 'false'; - crop= false; - for n= 1:length(varargin) % Read options. - if(ischar(varargin{n})) - switch(lower(varargin{n})) - case {'eps','jpeg','pdf','png','tiff'}, types{end+1}= lower(varargin{n}); - case '-rgb', color= 'rgb'; deps= {'-depsc2'}; - case '-cmyk', color= 'cmyk'; deps= {'-depsc2', '-cmyk'}; - case '-gray', color= 'gray'; deps= {'-deps2'}; - case '-fonts', fonts= 'true'; - case '-lossless', comp= 0; - case '-crop', crop= true; - case '-dbg', op_dbg= true; - otherwise - if(regexp(varargin{n}, '^\-r[0-9]+$')), res= varargin{n}; - elseif(regexp(varargin{n}, '^\-c[0-9.]+$')), comp= str2double(varargin{n}(3:end)); - else warning('pax:savefig:inputError', 'Unknown option in argument: ''%s''.', varargin{n}); - end - end - else - warning('pax:savefig:inputError', 'Wrong type of argument: ''%s''.', class(varargin{n})); - end - end - types= unique(types); - if(isempty(types)), error('No output format given.'); end - - if (comp == 0) % Lossless compression - gsCompr= lossless; - elseif (comp <= 0.1) % High quality lossy - gsCompr= sprintf(lossy, comp, '1 1 1 1', '1 1 1 1'); - else % Normal lossy - gsCompr= sprintf(lossy, comp, '2 1 1 2', '2 1 1 2'); - end - - % Generate the gs command. - switch(computer) % Get gs command. - case {'MAC','MACI'}, gs= '/usr/local/bin/gs'; - case {'PCWIN','PCWIN64'}, gs= 'C:\"Program Files"\gs\gs8.64\bin\gswin32c.exe'; - otherwise, gs= 'gs'; - end - gs= [gs ' -q -dNOPAUSE -dBATCH -dEPSCrop']; % Essential. - gs= [gs ' -dPDFSETTINGS=/prepress -dEmbedAllFonts=' fonts]; % Must be first? - gs= [gs ' -dUseFlateCompression=true']; % Useful stuff. - gs= [gs ' -dAutoRotatePages=/None']; % Probably good. - gs= [gs ' -dHaveTrueTypes']; % Probably good. - gs= [gs ' ' res]; % Add resolution to cmd. - - if(crop && ismember(types, {'eps', 'pdf'})) % Crop the figure. - fighdl= do_crop(fighdl); - end - - % Output eps from Matlab. - renderer= ['-' lower(get(fighdl, 'Renderer'))]; % Use same as in figure. - if(strcmpi(renderer, '-none')), renderer= '-painters'; end % We need a valid renderer. - print(fighdl, deps{:}, '-noui', renderer, res, [fname '-temp']); % Output the eps. - - % Convert to other formats. - for n= 1:length(types) % Output them. - if(isfield(device.(types{n}), color)) - cmd= device.(types{n}).(color); % Colour model exists. - else - cmd= device.(types{n}).rgb; % Use alternative. - if(~strcmp(types{n}, 'eps')) % It works anyways for eps (VERY SHAKY!). - warning('pax:savefig:deviceError', ... - 'No device for %s using %s. Using rgb instead.', types{n}, color); - end - end - cmp= lossless; - if (strcmp(types{n}, 'pdf')), cmp= gsCompr; end % Lossy compr only for pdf. - if (strcmp(types{n}, 'eps')), cmp= ''; end % eps can't use lossless. - cmd= sprintf('%s %s %s -f "%s-temp.eps"', gs, cmd, cmp, fname);% Add up. - status= system(cmd); % Run Ghostscript. - if (op_dbg || status), display (cmd), end - end - delete([fname '-temp.eps']); % Clean up. -end - - -function fig= do_crop(fig) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Remove line segments that are outside the view. -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - haxes= findobj(fig, 'Type', 'axes', '-and', 'Tag', ''); - for n=1:length(haxes) - xl= get(haxes(n), 'XLim'); - yl= get(haxes(n), 'YLim'); - lines= findobj(haxes(n), 'Type', 'line'); - for m=1:length(lines) - x= get(lines(m), 'XData'); - y= get(lines(m), 'YData'); - - inx= (xl(1) <= x) & (x <= xl(2)); % Within the x borders. - iny= (yl(1) <= y) & (y <= yl(2)); % Within the y borders. - keep= inx & iny; % Within the box. - - if(~strcmp(get(lines(m), 'LineStyle'), 'none')) - crossx= ((x(1:end-1) < xl(1)) & (xl(1) < x(2:end))) ... % Crossing border x1. - | ((x(1:end-1) < xl(2)) & (xl(2) < x(2:end))) ... % Crossing border x2. - | ((x(1:end-1) > xl(1)) & (xl(1) > x(2:end))) ... % Crossing border x1. - | ((x(1:end-1) > xl(2)) & (xl(2) > x(2:end))); % Crossing border x2. - crossy= ((y(1:end-1) < yl(1)) & (yl(1) < y(2:end))) ... % Crossing border y1. - | ((y(1:end-1) < yl(2)) & (yl(2) < y(2:end))) ... % Crossing border y2. - | ((y(1:end-1) > yl(1)) & (yl(1) > y(2:end))) ... % Crossing border y1. - | ((y(1:end-1) > yl(2)) & (yl(2) > y(2:end))); % Crossing border y2. - crossp= [( (crossx & iny(1:end-1) & iny(2:end)) ... % Crossing a x border within y limits. - | (crossy & inx(1:end-1) & inx(2:end)) ... % Crossing a y border within x limits. - | crossx & crossy ... % Crossing a x and a y border (corner). - ), false ... - ]; - crossp(2:end)= crossp(2:end) | crossp(1:end-1); % Add line segment's secont end point. - - keep= keep | crossp; - end - set(lines(m), 'XData', x(keep)) - set(lines(m), 'YData', y(keep)) - end - end -end diff --git a/saves.h b/saves.h index 576cac3..7b0f479 100644 --- a/saves.h +++ b/saves.h @@ -7,30 +7,17 @@ using std::vector; // saving file for the mex compilation - -// function to copy a std::vector into a matlab matrix -mxArray * getMexArray(const std::vector& v){ - mxArray * mx = mxCreateDoubleMatrix(1,v.size(), mxREAL); - std::copy(v.begin(), v.end(), mxGetPr(mx)); - return mx; -} - - // saving the fluctuations of the populations -inline void get_data(int counter, Thalamic_Column& Col, vector& Vt, vector& Vr, - vector& Cat, vector& Car, vector& Phi_tr, vector& Phi_rt, - vector& I_T_t, vector& I_T_r, vector& I_KCa, vector& I_CAN, - vector& I_h) { +inline void get_data(int counter, Thalamic_Column& Col, double* Vt) { Vt [counter] = Col.Vt [0]; - Vr [counter] = Col.Vr [0]; - Cat [counter] = Col.Cat [0]; - Car [counter] = Col.Car [0]; - Phi_tr [counter] = Col.Phi_tr [0]; - Phi_rt [counter] = Col.Phi_rt [0]; - I_T_t [counter] = Col.I_T_t (0); - I_T_r [counter] = Col.I_T_r (0); - I_KCa [counter] = Col.I_KCa (0); - I_CAN [counter] = Col.I_CAN (0); - I_h [counter] = Col.I_h (0); +} + +// function to create a 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; } diff --git a/tightfig.m b/tightfig.m deleted file mode 100644 index 67f3055..0000000 --- a/tightfig.m +++ /dev/null @@ -1,124 +0,0 @@ -function hfig = tightfig(hfig) -% tightfig: Alters a figure so that it has the minimum size necessary to -% enclose all axes in the figure without excess space around them. -% -% Note that tightfig will expand the figure to completely encompass all -% axes if necessary. If any 3D axes are present which have been zoomed, -% tightfig will produce an error, as these cannot easily be dealt with. -% -% hfig - handle to figure, if not supplied, the current figure will be used -% instead. - - if nargin == 0 - hfig = gcf; - end - - % There can be an issue with tightfig when the user has been modifying - % the contnts manually, the code below is an attempt to resolve this, - % but it has not yet been satisfactorily fixed -% origwindowstyle = get(hfig, 'WindowStyle'); - set(hfig, 'WindowStyle', 'normal'); - - % 1 point is 0.3528 mm for future use - - % get all the axes handles note this will also fetch legends and - % colorbars as well - hax = findall(hfig, 'type', 'axes'); - - % get the original axes units, so we can change and reset these again - % later - origaxunits = get(hax, 'Units'); - - % change the axes units to cm - set(hax, 'Units', 'centimeters'); - - % get various position parameters of the axes - if numel(hax) > 1 -% fsize = cell2mat(get(hax, 'FontSize')); - ti = cell2mat(get(hax,'TightInset')); - pos = cell2mat(get(hax, 'Position')); - else -% fsize = get(hax, 'FontSize'); - ti = get(hax,'TightInset'); - pos = get(hax, 'Position'); - end - - % ensure very tiny border so outer box always appears - ti(ti < 0.1) = 0.15; - - % we will check if any 3d axes are zoomed, to do this we will check if - % they are not being viewed in any of the 2d directions - views2d = [0,90; 0,0; 90,0]; - - for i = 1:numel(hax) - - set(hax(i), 'LooseInset', ti(i,:)); -% set(hax(i), 'LooseInset', [0,0,0,0]); - - % get the current viewing angle of the axes - [az,el] = view(hax(i)); - - % determine if the axes are zoomed - iszoomed = strcmp(get(hax(i), 'CameraViewAngleMode'), 'manual'); - - % test if we are viewing in 2d mode or a 3d view - is2d = all(bsxfun(@eq, [az,el], views2d), 2); - - if iszoomed && ~any(is2d) - error('TIGHTFIG:haszoomed3d', 'Cannot make figures containing zoomed 3D axes tight.') - end - - end - - % we will move all the axes down and to the left by the amount - % necessary to just show the bottom and leftmost axes and labels etc. - moveleft = min(pos(:,1) - ti(:,1)); - - movedown = min(pos(:,2) - ti(:,2)); - - % we will also alter the height and width of the figure to just - % encompass the topmost and rightmost axes and lables - figwidth = max(pos(:,1) + pos(:,3) + ti(:,3) - moveleft); - - figheight = max(pos(:,2) + pos(:,4) + ti(:,4) - movedown); - - % move all the axes - for i = 1:numel(hax) - - set(hax(i), 'Position', [pos(i,1:2) - [moveleft,movedown], pos(i,3:4)]); - - end - - origfigunits = get(hfig, 'Units'); - - set(hfig, 'Units', 'centimeters'); - - % change the size of the figure - figpos = get(hfig, 'Position'); - - - set(hfig, 'Position', [figpos(1), figpos(2), figwidth, figheight]); - - % change the size of the paper - set(hfig, 'PaperUnits','centimeters'); - set(hfig, 'PaperSize', [figwidth, figheight]); - set(hfig, 'PaperPositionMode', 'manual'); - set(hfig, 'PaperPosition',[0 0 figwidth figheight]); - - % reset to original units for axes and figure - if ~iscell(origaxunits) - origaxunits = {origaxunits}; - end - - for i = 1:numel(hax) - set(hax(i), 'Units', origaxunits{i}); - end - - set(hfig, 'Units', origfigunits); - - - - -% set(hfig, 'WindowStyle', origwindowstyle); - -end \ No newline at end of file