From 30534a99d6db408485a7b78db08addb0d6d7774a Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Sat, 1 Mar 2014 11:57:56 +0100 Subject: [PATCH] Moved RNGs into the population to reduce memory consumption. Removed the unneeded plotting functions. --- Main.cpp | 12 +- MersenneTwister.h | 462 ---------------------- ODE.h | 22 +- Spindle.png | Bin 26283 -> 0 bytes Thalamic_Column.cpp | 22 +- Thalamic_Column.h | 30 +- Thalamus.cpp | 40 +- Thalamus.mexa64 | Bin 41427 -> 0 bytes exportfig.m | 1057 --------------------------------------------------- parameters.h | 3 +- poster_plot.m | 19 - randoms.h | 94 ----- savefig.m | 288 -------------- saves.h | 33 +- tightfig.m | 124 ------ 15 files changed, 74 insertions(+), 2132 deletions(-) delete mode 100644 MersenneTwister.h delete mode 100644 Spindle.png delete mode 100755 Thalamus.mexa64 delete mode 100644 exportfig.m delete mode 100644 poster_plot.m delete mode 100644 randoms.h delete mode 100644 savefig.m delete mode 100644 tightfig.m 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 51971cc7bb9e0b7f6e13e0776acb612816b589e2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 26283 zcmdSB1z1(x);2shjkF*oNJ)2>0!k?$AO*Csbt-aQqV~#n-J?=5a+Pu)#R3$k@e+mFVa`lRm z4gh#{0N_RvV8Jt5Y#)9CKqz6OsHkmWX?6GULuZpm_sj%r%*_FCONx!WryUbO?ccfW zJxD~fP~@!V5zS+XnybE=pN^%z%%@IK$YO@W8Wbp|9Fc}bW$2v0tyD-(MOB!lf1ZF# zm{$3w_S$Dz%9)YXW^bR1>m9N)Nr~7B6u3Y^8Kt9x>m^U)sd>ek7QcoiDn8L}W^$5n z^w%L5GXUGSy=&L+aCHkOUI2P{XfI(`b|U`0652m0n}HXVZA|f z%s!<^c*6ekDP`F+NgM^Phu~a;;;9uD{I$FIRF?=Bo;}mX-t5BR?;_24eYWHCnKIJS zaeXI8x3Vvxbr91TtN@8O8b~a7yD8TUOMmE3uikzbStsuz6Xm zM9rQX$=@63e4n0@K9Z)XE&r{-Wfx5swwt~1J8u==@A>`}tH`T_zJwcm zSnAJyHK|#|7Wi0^6i%}@u!pceJ^zUPKD)pbf__@VYo%xmA60{np@3rrgNyW%pO_2lV3gQhyikefHfo*EP4(`if1N&6G{CO~`|9;P#0tFt!ajbxnU&oZBLS905Qvp(E663bqf*rIMpoJ|hPi&)+`BOcXbo#8t#yje8#_`atYQL!OO+s7=cZRc^q%PoJ=)yJU)C@L*MC;y(4Z582CrTGh~ICoNm8 zntS;)zke+~#ys$PKQPeyX?SF8aoQWR|2EuT%fB$Dosr-q4!uaO8=C5p$ocTjxb9q;X;f_~sLt18$MR|;_$Gqxird%JR{?Z`z6%TIA|$dcDcE z1{dV73*1i8m2=0;E9m0pq?b($*M;KS-=^37!IaMMRpxx&THf!sx8INGRvS^8h#EXF z8NAM^Tbg6H>aZ~$$yTB|F#l^_{kJ^hHL1OVx~$i^*Q*QHXvg)o3Tdugx`r+9=zHHL z&PKrI$m_yp*^1xjizDifadEfHe-ud+^XJ-bS?qkNl$Sjk>iIA@(p9&_C5@X?q*;By zJ`sLv+9X&~s;c)3SdsmW|#4blxYlbrsbJ~53A7wAPwl8my=F`g9zdgJ( zdhU*_%mMRC-RHs;;qlM20ixkiOx5m5hi-FUwj2-KR&R9Y#CA@lh#!5}^bF{_>2|l? z(QKvsc%3u4>Wi+g#ZlpfTRTsW(jOT(R#-M!Paov1FAT@e8*5gp`F;}9-N^Zw`Z!CfK$oF*##eJ?he} z4Zq%X8l-yGx_Qx=T-Dd##TQ7AsM4V^IfO)t(BX={6958i$RElrOWp;ZWOP>6bG~73 z00Elm70E@5d%mS)P;P|()VRZ$oyS;T__gOxg84}0VT5E3JyK;8qb1CcicB}E+M zot2&yc_Z&8N5MCApg#d}s3fL@yfLSjDHKX9?D&x#t$VYZg{&l5CM!&xpewZvfOq3erB+X;W5FsEcZ-zl)4q&yG%!KrNcnvnt@rFvIhdju(~Cgh`1l2Xpr(!*)HK zFr}^Tcdkr-esnvz*3a#M{Ncj_UEe#yIa$BtSIq;sot{V&y{%{#iU0BN?8LKJ-+OHn zs*5?FaCl`(?1==X4pQFkg}EKZzthTDIGeM;Gm!Lqza}T_@adakVH28^+N<)5ZP@nJ zc(r2_&tnFnumx*b^GkJSI)5e_O7%9MH5G<}-IE_q-t`KREg8`Yt~lQJmah9P(sxF1 z{Px7zeJwJn+KQxKR2MYDOKb<@y^cu7q&@StB!gI07Y}m9?wPAChEllfqs(jANM+ZZ zzoy)$&r5HAdz5pl)(&%3+BZ6PkM-foJ2xGh2c5Oa_O%ObI;rHpZG0`lZ51O&ou4f@ zrCoGNMD5L`|84hHl2NPHt4QfFF5@eX%3%g?`&3?KNYvZO zn})hJOjghzER26DmVb5TLfyHYx9iQ-V?J;8G~_Q9^^QuP5}0b^E^RUJxpRk|WMnxj z?!ufPjVXB(3s%SyrnxMK$$js~Y7f>nK0eXzMuJJ-=40j3p3f4K7uWaDdefB0#qGFc zJ%-_m6huKLg!o~a1!qlkSGVrs*BJYFP}g!q!99phb=iO7-d{(*P9w4=o4(q^?fLG7 zWPS_JKcFEEiXf@dDlF5SC+#j40VRSuo+H;F!F6ZK3lW&CYsa9 z^e&WYBYyYOIL0g0cF@@O_Z*Yck0Tcr=~ucI%FO~(8=E|$~fJC&ErW`|X*=4SZQ%g#RI_PQH$ha+ps zKgVuSDdJ zzlVF2OJU4XxZdfgU@jMlJ3!bz3m@Y5d9E9Amo?d7n6ASYH;_AE&GQhN#(nF7C-xb>f%a1 z{5jKN#zF6+XKr&#(sV?<;M3?r^n-`%p6(y`cZ;`-bMdU(q?FJpdp=7^+c%c8j@*nl zg<`9zXA&x^41yL4=Un#9oUZ2#9zPg7w=oj;+She^?{LSZdYH4tq@7@y)YKL+mfW4q zfq2(x{nAu!BMU<5T)iElTdO&7A1goXPT!wAwO>-V?-*C=HMz(zb2=`+W5{=|2hU`D zSgb~WBt)yq4Xq`nguWhg$RRY>{KE>bUp`ba<9BI0Ug6qHRsIq^HsiVSA#VSk#J!02 z9W%YT^Ae4e`ZfXBi#CBEc~__}p}P8&wXR0-TL;OgdT=tg61#kP~pIG&U2d8$A<;OiZ_vwLw)u@D9ci2{qo=P>%-7Li|7xb=(eHOQy zC{`att&k4MCpsx|H6~8YKPt0qvg}pE4BFzqp~T&~Jxnxm4s*?@qwTirDkLWtWEo~# zy`hJmJ*J$o&Mr&*GG=u^-97fEz08nS|J6RncVlJtb-UfG>LaDzNsH1Y?8=%)lM6U# zHMsP`A^FCXW4>1({Wiv?e$#W6MNWEswRcH)nv?WVQpN2#gWUZ_GS`oYQ`M)= zEYEj8bQyY3OQ1_1|I(=!uRjV;#Z5GExDn;?^IGfE`xgCicp1R=&>34Q4CD!Z$ zo>ujazH@x^z^k4(Y@t49DXc0MjSXjKiP<~9U6!;Xx;=Rub}U`A@MC&8QLTb0yKTnp zp)9XW0czJZ<*i4r|G=WHv#;vSZ->g2ZM+s)J&)?3z82YtFF*7;-D-4YrJ#Q_gy!Hq zY6HeO;+!$gDy?Apb;T(H!$Hf=V)u!eTH7!ObeOF|$9Aqss`AHb;RU<=nNwQjo-RyE zKKs^rt-%g&E;_`QGfkMVE4x@l@=Mp!HH*m==}JAvLgK`vjKs^axg2A`pS|<-SD$#F zO)1%WJgIKGIL|rP6tDa0mr&!N{xYb(_4^I3r|ifP^!PW!BTanEfhOWo-<{=OC2zz~ zn1t_tarA92nd)i63GbBbW07?^dh2*vwl?B3=dlK3_yVr)18@fklHlmUV2KWsm#e@jg;t*C7+1%|n4TT(K2I4E zP4SH6_OHKpbn>@Z^PF44lE0HlIn?JBi~0^n>2Mym{i2aBY7U!2*NUSzX#)rfvXjUH z6}B+h8POAMF($o|RX&Ll&>rRHl)K~yvWW?*Vy_go zT|DXYHqdISoKA(f;RZsvyuSAhKW6-xz9!s@cXPOuKUMSDjdtsMAEPU32cGeli(i=h z#Z**;zrFaagDD1E{Y@b<$4vOcxBLt<-Z*?JkoxTSS_+2aYTv5n-P^o~ANG+QEEthd z0KK6vub%_+M_;kh*!$pRAL_Sl0gyuOC>z#D#?GFI!9{{TC@9M@%&uY#wGLzzsyWp7 z1y6fITmyObgC))RL#uuEY#!*l<-Y_#HPjU_x;)+YaJ1>k3n@!fSBPUEz}y`>!!E-)c!yZtg6t<4?AOUw$1FpHMq4cJ3R>IGYaf zA9HI3npY7fmnQ6#NYp|YC{5O?uQ?41WAGH!BzdPtkVpRh`aJfdb8+S5$XTuqsfMVL z?CLT+D*;t8WVkdMlGT`IW?GAyd|aV>(&5K;DhM3{?aZVCh!mV~C34f}pQi z(OesHRJ`5ii9K!glNmW^6(#)}3FWsKsi8y8;aLbcx!1*>^2Al5(L!d@Q`v9GWGU%k zXMPueTo7E2pBnf6i_x^p)CJ_&%7Ve{26t9ipQ1B&lY4pQvrzqI4Eyp!JGbsfOs!a%5t4wu4)- z^g#&s#VCdR$&6_j1N8h_BN4qh;L_Lz5}mKbO*!GiH-zGM%XH2xJ4&xUr}rGLJ_|2F z^TL_5x8WzC7wiA?N`Elt%NAmIq##N^!JXI>Q5&YyrXTh?47tMDvD@?A&7ylX)r2+v zwMi^~+*xTn^nS?v3hML|`hz9a{`Fesc+U>Di<7X`kz3I8ZI%MeALqTjf9~QVClnm7 z#_`N9j@4m%?))U4%@2N*hZs)w58lm$RwJFaY@A!ju@LEh$DsbdpB4Pw>`rX`{|9Nz z{}F)=uFS>hEQdn#P=qNP`uCOx7aJB3V*5+NW)v?y4J+VCMknx{kpwO#9s=BG6`{A+PFV8FciN_MI(b4x~V(;gx=EiwVqz9FQ?%# zdP}Tf7P@tCeh)&p+y$l73>&zt(s9j!TnQ%?cg~kcT=u(azb#2SY8+xpgYzw zuI1}@(wn(EvU}P3%C|e7^Yu6ZC`ZcTtV#vX0S+^U z35y0*!B3X{EJdlH_iA4$^o6p^2?tMsLXzkGPA(kYg@&MZSzHs2GY9 zeMCT9jSS|jz=G8i03e&rx!_-AN_F=Frl93Aaw;t)!ODJ8uuq=o5ef^!~DM?IH@$29>wrTBua-+#M|W#iy!(gMPl~ ztKdvX34p&DgBe}cG^4LUnDYrh3(}ii%~DG>BLJBeNDRynia!Zt|B1>9F-w-ATD2q= za8pHs%--Pq>ktBf{p~Ns7iX&l7C@}|XAJt^;;K+F|1~=Ce}C4oAf4aUxG*Gl{X-QB z%G`9T_(13)1T(b3AIPS0?Z2s*n~^7=l1Ba)WQ6XGRKZ@VeWB8YW*Y*%t@vZ#=TwpO zX}Hne&?ObvbeI9qQ2P&)qUXxBh>U_tFM#44SSxapbxEZ*)frg>UC;mTqSn8!+Zm|; z86sM;p`;2oJTF3TLhJum!x2sFi2~7A;f)v>i%S5!x%3~@@}ILbOw7q9f^wvG#70(DeQ&6# zz(@%J53=!3NpwOxFA05 zz+a2caL90!+9|-SWqyDND{PUAHPn-~iF6!O^gsq*hQ;(4Whw@tLONa?2hgTI#ReEQ zlEPCcgoUCltpghqj1YQ=i8}JI$FOE6G$!dw2MqfFw_jRr0PvX}ytZj#p%zkBfQtrb zbFqdf9H<)rewf3cGZ`*Ed^HI<#QqvCQLPr^7?_jLcrxCaScAwPSKfZyT@EV?-Xp=| zN*@!Dg^U#hOEgS8@|fr`C%%J>YkF#01RHcc<%X*|p~JXMKw#wru@mZ*&zeyrp0Q7@ zKbTvE$_~HEjT`(J7m4^_K=kjJcS9r+fj_jBSOd!ii8?LHxN8=8@I)FBiKCKMxlq=| zjv(Gi9l0yY^IgcqY6nkgX9g+m$zte5NgbI@_S4WO$m-d|a(ny6PNn?LWAb#5Kf!h) zkm_V#n(hjmBpm;&R6RTW%h&?|w!Z@Df1!$vNz>seNzw_%0%I=`8Wp>^sjY$o#)SSO zSipt-@8Ozgxjzt3;T$qaQ*6{2DcC@;1e5j$)FJo-T2H6IAb9yd>u8w8X^9^Ir2fNe zpuGRf0*7#AjTpjSLim3b+KANtgdRNv6*2e!SQr&jU5BrK-&N2{G!(@KNr+bdydIsJ zW#PpHFwYVAI&;F&efw@I13MYoeaAbggW;y>0#4YxtWiGnB56t-n9vQ{`Ls~ZLII9h<;VF^1jre6dNVryd%`&}VeeT# zKDr+R_-v|ktA&IvdB3%u+g#&G%|SH>o^5(6b29ZnP&`z;_t z>NwI$1B(O->ihYOlKBvFc_DY^V-;0e`o7#VSu2N(1!P3x-I<$CxV#WEi6cjNWOte} z^2I?9{nsKngHuv4-mp2jeL5|AyT0S` zadA-K&A;V>=r3h2U|JFSZ{ZN5GvtH^(xm<+-WO2y)nfzn#9y8ABApyDnCyXsRZ1CP z0lgH&OBfoyM(Y%jq0r9|pA-3q&%wg=e^5w}V>=`j$f~=ZzN{?yn%Y7MB*eo7Dri-l2AG|H zqFOnT@AO~{GG2_5P8385G%36&;q#vaQ}{0%y@GB21cfI0&q{{zw38zOVCWy4gC?Xs z{WCI|?4%F^mw|*CbiEZ*ndJT!HZ0TqTaB(@Zzm;2X5fLO%zs2WIMBaHC?}FX!kY9S z&jS2=`v1>|*@!t|^E*gk8@(0-t&shni;0laQ{TdMEBi|pF<7-nRG>HRTYEm6NI5+- zz^I8e6yiv)9hs<|1X!58Ux-Lr4XF)thg$~F(IJX_b8XCZkaHKmX+FNjO&3bw;5bbM zL6n~^Kq|dYiHjZ!Uf31vFXh*K*FW{q_+;&?ZiL3ds_Ds)P??=Pvgf;!Ea6F@lEJbS zZ^Clz-U-1(62k;^Z88cds3H|`#_dZ_Q~mK@!Ws+&>C6s^hKq?d7OPVe`^N((HB^n9 zpHfP~!m4U-Cy%~V`K7ZVF4FaFXSzM5dq3VH(oJGUA=-(~v#es3mb zKbZ!zeFxVpG6-(!(z2wEWEPirjod6ICO?v<442r3hcAdzIaVJ+8-O(hBnXrl7nZm6 zg}m6+n&LjmGcbH#O~_VX3R3{9aV-NVU;bCghC4nQ4^2*MkrKwpApx{3EQrK$Tb#q# zLNEP~>UrB6WyPTvCsiGxs2QdmSn>wu(Fw4V)DW$ zD>f(7j8yZgbJGdG$1EK6cw*a?4)SM-c{0B8-bbbxi%yv&7yLeaE_GqoriBA008k4W z0gi!lsaj+|OKjs*;-qQzAMm&9cVtEJ!g^BVw|0D6u=8$Csv%WpAS6W~d}nKS`a*L~ zE~^~>y{&2-K<<2zOQL_xAW&C+cm`@q`_2>qqMB z8Wbk$=EJ}TIiFSQ-hA#0Xe0}yDQY{1NMXwQ>+HniZ{yKg?5jC5h88R^IQkLoR&Q1q z635}7ql3g7TC+cMDyK{LAPrVno>jI5AE zPY-P`0km^jP6Z;u7Wk{co`s_#Y!K9*g+RXcPEu*}q^g*Gd)4FQ5)9Tq!7|#geQ^{< z`<-ONRQw5pXv@2;=qw8ZENz_p1~!m?D@E$q8OY6^aZ)BT)$zR*Ujv^wAoFbxR1v`f z7HV_T8IZ?q_Ew)uA3y1-2vr7=Tvs@8QrzDAfmkQLmiKyTlYH5z|3sweTiYIik|)_J z;Lxr91DHE0`lPpNE-e9RqI0z5#+|a zI`}M4+S$ohP|>F#IAC6$ABbC6EM1@mdj9_;*!>Uj1m6G8o_Io!X7+j#0d&bf(fU6J z&VOVPEhdX+GGoFH$4YE1OA$>1sk9}jhR=w~KD5nP15>UuXVYO~>k!p>12UqtA4b!C zXluKZ#k9TI+a$OD4ec8;+PppQM1)kKMRguoVgs4n5yb3nf0LS7%#Ef%&7K78gl8cH z3W9$J#S^={^4@>;HV-mJf&cGG@@QTcCmVeAJ7d6(REG57jmNN2 z02QI{8yk1FxV|l9HYf-=A%wMUi2&-~6UDdU@%^mq&guDt{6}HywW8uYI(o{I17^6espGPE*2qdf`}z0uSKpN+dGEke%8K)R-(CD{wBc zwR(?P*)Rt}lJIb2BkA!gA83+yzmTKK)aV%DGKb)BM(H7>(84^Rfbhkv>`$N|tcLhi z_IX7dG?!li08{Hmj|`$tGlmV22^#M0LGHP+IWsXjdX)HM zXnDK!?1=~XIf*NJo#wy-RX0>l0Sr!%WuR95$!rN%cZpZ8CJ-PqQGGPCBuF=y|-8019$JHOCcc;)W1t$ zkvk-&kztOX7TiHM@oa+}Ff(6=nkKs_+8~|&nblLc`86^feFSTl zaRBk^J>%&tRY z42elN1t%k^yKB&EP=Pn&_x%&1qdr8tcn|gLANp{q3#`GC`lUWn1@2DIpzTa-ZNXfsg17vX<<2*s;>W`1ARFAqB40dr%i=Z^falHE(S9YCyOCrlKz)}F zdXHzLNJCb-BC=pT)hEP`4I3+F;-i6;P~*g^bHd=#!ScjteD*KiM|lhm$W`=}FA-ZI z5BjYPf^v6^r&}s+uAgyd403%3Ue89uD{dE9`|Eub?XEnnM!kE2zHb@~rbT*bu`nr~ z^iUop|DP>#Oc(Ch>rI_ZJVs78!M@KdcwIR`K!^{n4!*Qs%&^35tXX^xe>>W_SI5c( zi@mOcE~59b!E_qw>br2qkrLa0qFUq&_!4g+Cn_6VmM%Bma}D}a+u!KmU%7}#HpIIBl0-zXym>&b%5uT|yG7h@` ziZqJrUn}{4k^!K9Etu$g)@~1y?yU}(0V$<_xjlf{@6d_xVlF9Q{_C3LR(=?u2THDe z;Rl~jWGUX#AZrYSF8SAyJamEx6uh8-dVw}Dra3%?dxds)VYaRtK2lm#I#I=~g;WZgr)sL2@@4gBlxY6=RuAX7(usBcX%!fh znmI7NgnWtqNMEotB;eFY+%wtN%=3k+()Gkmu9oO!g z$Zp#JWm-iM6x=40-6B|SN56ubZSSSO7tIFuoA2_8V!9HSa|kuz$2QDl+ne!*=uLF^ zMjOaLVkjm;Dd|Gs46s&!=7?{_22&v8g(i*D`QZUipv}ISpJo06o5<8NjdlP97e`?6 zR{Ha)eF6msI~8DY61>Qt>rpauKoP1$JfR)_JPFO5i=+69M6@6w0SYv#K?J|J4fV>* zRfn*5a=0Ay^;iRyisQ=tPH1@e@OI9(lsG_BGy>2DBkqMxQwjx_pbVyXKhmMbHsQy> z%db(O^${h`$Q&I3;_KzalkiG$DA=q5ZSzaDY86Nh=>I_u&|M|h3Zb`w%jU z7dx+Va>8RqcKo&?dcsbwz(lkhGPvjn`(i|H!b0 zX?JOJGJs(Z2e7(1&+NfR;Z2=I8lq}&13Q@x+fezwSm4qEQ&ieIq!s>F5dL252(G(+ zfTuadX7wMBfTXd+DA4elv@pcGcw2`jA#tzzz+ghL#P&ug5~UPGf0Gkk!3hkP??1k) znYpv7DDaRO9P=2SeuE)5wc|md#xD8meVfY^Ouk%mTz;k<;0SLoc+{#cRYOLr;qP|K zDv1;2cjn<>1^2Tzxah|6LzIwYry}3F69WAU){n|zvAA7BB8z9g9GA_BaW7pI*tU7gkT%0N^t~RW5%j}XvCD&dvftAZFX90-(LQdq_iO=MF`Zd;d z7g&W=qBb4CXm@0;B9Qjvur*i;>L`Rq-jsgn zEsSG+uf==xocY%q)f60A;3W=WrS`7S5wrsc_;-H8VF7JCaUtV(n1exDZQcE;?aoA1 zQNB*(mALfF5l`Nf;==&JO(=x(I53`yaBa<`f8lut#lu&g7dDD(!9BH`+xzL!FefYu zDa?;%p=KL1yfpo4ESU(7#!H+2)gmzTyd0ew+rHnb$7cHTHj3G&FRZ>s)_Hm2fLq)T z1L8z-HTLMu*yx@Tuf9u)ZIM@}f{dj=1UfF-hPM0qiyOhb7J<)g2TREIjgPkDbT99$ z=^D!^QLrj_HQfW+L6tm;4(TG*%A>7DKMs9-wGL|zt)>H#cKpW9$%-vZE8extCzw0~ zf+|OkZ%o-y=|-I0meF|Hl5>w)?s1^yl#-{s@64W{pJgoj!uJ|hBWs`_#}=Tb@f7{# z4YV7&mFjr%>3hE6qg88!tGnl3QS6e^fiW?bmsiP!P;_SQ#+Ag;g;5W#y?j-C=;c;B z1rx(;3w}{VnJ`%U8Qs_TR$}`?(b=2HMDFhMXV5X$2}Pnt1(OMjuC}T96Htx4=;rMz zGg>WLqXB4vkl)hRPSA?|g$;VQh5Si9ae`EH;(F*QSdZ=yZpvn+V^#P%_v%BXILLDQ zaRyM$#CgOJNB%a7>*;BLlOk+9NtaDFx)}+E-7VyE3K|H>=QA6$x^mefFyvm?16Ws~ zp4v^p+)wcyo=r|6j;zygf=_Kg2cVPr8o}~G8J7tjASg~FVuCv=CW}RUDOI%Cdwns! z8(t=Tk|t=!3$p2^#U{QYcD!CuLQo8`2?$ns1SVJn&L`X*j~juS0BIK?qQO^>PQ*<= zVhXu8-L0E+j@5kMrSf#7P(;(r zJHb$s1Z`ATWX{m(R8dN|2o@0*3g{w242dG-iNua~macXgZ_h7TVprIc0iojSw~xg} zI#|L8d7XkDg%L-7nPYh19K;D2>N^A z$Ip*$oC^!Ml)AFcB^@G(|738nWWZGoJo61cJTfc{}v&m&o$U@o1frvwPxZE)zbuSMrSM*lQ^R^0?b#uAM%Su~B`}{xFG*qodZ8uYT!+FZzYt zJ&PX3D%0B|_u3~eq#ugWu^ty{im#m+_Dt$FYaU1szX$jBlL+rp{DJf&vv#TFu0*0VTuI+ZYG+M!~G&=X5TVjB9&qQ{`ev3A>8 z73}EY=#;0#DkxL4T~|YsUBL%2EQu|iYP;|If8iuNW{rbXBSE#*apa}zwH>P%eQNC& zgRJdSy@rexl;yS*KVC$CT-bex-oEUAMZ3o&wHohbH|fJ-hy>^Ys-bF6X{s;{Y|^zN zi$Lh>0$@x=+1SxK_zdl*>`jlgfzkpEPjS;1om~tgFg4IncXN&OiG!gPj75C6OSX)b zRa{EAEfp==C{MLo>cAnQ#H^mweK+6tVh3b3%vk1^$L>bBTa@|VqFied8pPoeL%lzAfjVElU(MFx7-`pM2Ov8#)G0AW-~AlGBpUw zxf*A-PznjMdBXl+E9FLqMw2+HL-GrT%{(WSYwJ;?j->|?PL<^?bx$TN0*CCZmVbpj zJ|CZK!}sYU;IqrCBMjDaq4Ic+9otzbb2t0dbIP=HbgWFhm2b}qCm69BhfWzDB_>khB2hP)0A`jYoqlDG=j0LHs=uD^#y=#;S$ zlFs${3KifOS znD3E_?_HQ5+?p>_@Nc=GQmfmoY&AO=HsDlh{a4=`g_MJg>bM7V!x2@E;PKq?5a`o zqzN(_GZ684idr-Ltl}NH@LzV+)t|1vo%=Kl8CpZdtqUeL~3-i=t&T!2q?@a{_; z5J`DkoA)!?cUICho28VI(y`9WF{bOg*p$Db2HMZstU2UWLpEP}KCin!5%_%gu;A<# zx!3xA*~0x`ra}peKq12F^Snaw*m`k>7-dpZM=YY28IQzg#EvXYr+}cw8?>Uv6Z9pz zb%z+kbITTiy*eubXIA`aYUM-AZ^3pw3jYcp@6peMCrNHuIdTyMw`U1017n;wWV)AO z7UC?dJw@TQjwhm>q5DfJ3Bu|vjfkI`1w_?S&fatxRZbePOxEq6?BC!3^Xl|6kqiMj z?9`njnbi*zeon+dGtdQ_ac}+ui-gZ0tDRx-dc^B?R~e|W-;<%7H&x{u-aDGTmQmuw z7QCFCDkMBCD4B!0`NQXC&zMQSmXOxL9kz?v8&aigqgO2h8OTMiGDLcw)4(LrrQucc z67%>PuRpS3s`zH&({J(aosRl1u^VA5cJ3ldAHFdbJHMf1N+(N~x*-Tyg*+`pH0Njn zNG~6O%`QiuT`%$cvMc@Ii_m-sGv*zDb{t~O9PDR5g zpck5N#nvLv%*5OBh_*D}bm+yse10;w7>huKyVV|g8P%-BX{Z$E7v$QnF29yAvzuz< zFCSH3YHe*}ou=^I(KovN_*rx4j8x`+F1JNsrqX=8F6Hhyx*Pp0L<*|#*$6T6`zuB^ zd(p>*(ga)*@^pPl~46;WTteTabq`%{0{FhZv?AC6ApGv0&Sa-OV; zM!m7TGg-a!zkJz-QkxcOxL5~ZA}tiO{iP)GyAFM4roU7rEHrsgj~+_+WG9(~c?FCq z`nT9KnHX&blFE*lyLrqxzS~;cJZA?({;IfJuKncdffHW$!V_;x-(6OxHpS^Tb3_8e z*azDMrFPhKtojw_{Xci7VttOB#>Pmhlgtg?^HOIV`5u!@e>T5~MbW6irSh<>7i((_ zC23k&5#`0ZXE*nQG0v#O+tA+dxU4Xo*{_rNQ>wucahqWy`2e$)ISKN0qoW1HL1g z^>)%#LSuO7@6KAMT4Lli=VgOlxe@mtj_(@-?W-LujArQ@j-RnEnB{K1SN1RaP&B&L zo53sJobZ^y_9g7hi2wzMrt5cbE-ZMf_F;!EKC>Sg^QqlP$D^D0b~|ha@3%Y6s-CsB z>{9QN=Tj@R_}uXe(-`Qp^}!lW3yogciJ85IXXMnvudDs+frRG);UL|`PqplME{eZ9 zQ21z_>=)+-q$COGDpt~0qmRaB-mU03p{h+_?L)@ld1?2mlu4Cv3oHbs5oqlhM#pD4 z3wwWp+_vmp(>2Slq%@ihZfcEU&@jiUuvdyUIUmegQlncBqR*U7n)K}t57<xN%x|OW>t?Ks&kuS&$X!iU-0^3aU z?pFrmTp!L2zX*`lUT68t#7JN4Ny4qxH&Lr4rx0PI7Z3-RE#7}$hzO)&%2@Xhp z2a9uuXoBgb;Pl1hPi>7;RHx7JOIV1H{@_-fm)e3ox;fhqrFxz*iXO?h%NXSf3TE4- z*HnNBrq||7q%~A$4Kmxj==*NwM;|(dDOUo_nLd-92g8nHf}!ck6TIJ5%-HA8Lryz@ z{9k}bZg6m>IsGFsy6)fR_Cv-Wa7evROP_!vjD5#&?0F;i66GwmX?iTZ<|8{+ zZS)-DAy(z8y6)ES^>|w}p_`9^Dy%A`_z=}I)@Mt3O@Fw&@3bv)K)6SCMu8`UXyM1%aZzHqfo1M z8Po4dx_py+`N%|2{*X@GhZvvTf4%%$=B8%A45#cQa{&HB>|8p>`t_v=YT;}1-lmC%N&N09eB=%F8g(r`et8*p(nLnY!#MpqMRbUCxT(3SkP;+eqfw zybx7$k8kJP!C@Yg*HLVaAPhdq3vlUJA=#Q3;;MLL6mv#32 zqx+mmJLHd~+%|p@4!*sWb}NN?~G5a$uaqF$1xP(YW#UR?VR^R_TS{e%y%5E!8@Oq4@v4HA{EUwe%NIgz?A#pqF%I`l-c?{R^D*9pRkD@3IYIMHxumgT z6q;bmb&92?@+JB@d5=?+C?y@MQPxf(P7XdSF*Q=5^*ZmxaKRdY1|76)eCpsX~Scgyk+4y;kxh7Zi1-6uqgnN} zB2Hkvfxp1%<|@yE2jv2<^IR@}D5-60b|-?*1Hk_ggZW9^HK@+CUex2WNQ51^hXTDB zq!eQXOTTQ4gfHdd>fhnYfAR|U1N8ifDaMHB!xW_lsC7&Mm>^nf7rvm^N~}N-I$8DZ zILZ1n+J1-)(8kCpp=(3tm>2@o81&GVjPHF*RnDQ&u(9uD)j|ql%#2W&S$LG!kvGkX zPo5_6*d6|U5DsqAT5}&}Mfggh1*CmaNJR~GWZq9BAUU!sP9_4pCOZhGdd$0?!)}9` zNG7(-fa|4w)*k)iJ0K`@DsjE+1Llx|MmtOr&z4{^_((FXF=saE$u4}OBzqRt( z=&fHAalMv6^!QE?(e3H_94Tm^-U&Z}f^E$NRvKcikJ+;Y;Sc#DzQk|bgbEzX z3|FL_JcEAa8loL-$RClmm*9a^77N74WG;NYj&;&dQi+nc6|KXzN-uyv4MCc0-mu#S zn-~SJKV>TQ=S+@Dg3Nf?akL-5?U+47A>rz zXe=lf zE$&1hF|ULL))?Vq4n+pXr&P(t7)HA&k!LF6tkh~rZL5onMYTjnGBg>4@#)Yo2L01K$!uhp`-rN$0>m5J0}sF$ zE0i!o`DYyi&-bRFH`E&so!^sTP91HGz?srNawJ=3X}QH?gL&9A0L0mjDOup7Fg^o( zZ~@-u@Ir*)5nl?i=Is{gI)Rh2F$$cb6=ut&c({40_3Ke(w6YI;q85=497WM%V^rZL zfV5`?n<>U+#9+A!{`f`j&QE4+dc&E6h5# zmW`>DfFY2&#t1b!6S{NxK5PskwKra!FdWPBY+!1K^wbL9;#k<7i#r7-4szn>00U@|-iW>YTh~^Ph+)8Ds!Ki)SK=nbpbk|7q{c!>R1n2EMQ< zLu_`5%#@If`z;c7L>b;~j+CKoNC_bk8#{@eWF8`gQVE$dBx7%8LQy*zHkqd!L)q9k z>(Tpt-#ORyUf=oW{Biy`&*i!v&tpC7dDgwwy6?5t{kxYMPj@JoN5|kN{^>~%b=I5j zPp~*LB+9g13VXaMQyP=2OyBnoaR41s&~r5!C&O5miiq@=rx{GvnQ^Av`B+uF%mvUf zurBU@)(DJQ#2MtPz|c5&4Wx!r1lXWkPXiH#?ODS&8RBoL8*(@lNfnV)>#InA!oJ-G9R z@H1lrS1Sj$E0MI2z-A^ZB`4Q$p^6%iq~2i@1f%nwLmBhk>AJU5QtG{cZyN(i4#9-z z2Q0*>6|WlLaGpaS$~o`5dJo#W>MG8Dcm5fWjrgle>TG*R4Q|K`ZWTM|O>*pE865kitqC3Y$bb7@S?HjE7Z=Fbr#bWx-o&L65nE*(F|k zB?w;7@FGB`D^hvXytTCGXDHCT0{xBv6pijUl&hoQFjsJ0@ayG!rk8^))!ac>9IPOF zyu-`r@wCKq*T({>W)d904OTaaDgJ7=kGQvWe#?UkAcFx?7NoxSA1xWZTh#~VgW0bL zVc^CMK4ehHmiew;NuYdQ6;0>hsJUVEgwnsOdI1zF!7>zPC(CM=a7X<30YE0y#z2{N z|IzVebuHZL@X(KFHbTN23iu%cv`XfI!5PhPk(-qAGAZ`qm-o#^$#=`XuC1}Ch%o>b zsOd;zLe}Hb&6CYmlTD! z0(F<_r1-#Eagu6GT2AAk#AwjfKGj=uAH1Go1ERukh%fk(jVMW)4_-JRNL;VGrG^jVJGO}S1DBWNtsWzi$TJ2zH-l(_TJ?E&1hA4IWKgME6s`RSNf%)>LmwpN7s8W5ba7}16^Y?--Jt2Y(K@Bt~)Law} zWO`>>@O(&Wf(6u>DAg|#c*PiTxbV0WMWU8|+zN2Mc{kswhHvSG9ExNrS@}XE*#y(& zcu)Wa(@die@%+YE7t&M702I=7X3#T5G5jC~M*J(J=qT5->9?WefQ-Y%tv#T_azQ@>q{Hk*r_te#V%vVm_suW_Gj(B1_Tm|s z8Y-y0O7#RbXQI3a+}~u=XgBQhP~*i-i0OLPQddq4m!C8xZKrdXMl-FZa6v5v6Sx=& z!4M$Enx^ooy<+;1-vVG#CyXz)9wDqqDtW({N9YQGm)o=ra{fy&>3@lv{|D9i@5BFR zf1oV>M3nw3uz%*z{v%r$eX5$F4z63MZ>*aZz0JhqWsByCaBYTOg{h;;TOOz;PBC**HhCgET3`0XBx`1=Ja z7JtF~lc&SJq>67hqegM3UnsJfdZ*TC~)u6 z2&E=bP9NGXq|=A-8Tbz`jXaOr&dB;J%Jj$RUNPZ7XF^c!jeQbGCJEM!{?lmJ3ZB90 zbO!nKGPghM;BByP0iQS>6&5HGymgJ>S4(}mjY}(E5}`?5%q7?z@pBRLUQ4#i;Vufze>KEv>s>qL`>Xd<%fh7q{LqaRfRyZ0O<2}FhWaP1Ww1#s?Um5GR zTAaj(@ilA&b+YbY-mWQZSxXk9zr1FTA@cg=|KTfK$HQIpHdwYI$bwbRTV z7wpx#H`6*4S@Y3tB1f>H)AefKkvSFYfb6mpianDq`1Z~IaC=FmZ10heUwutoQw}D0 z=cGuU5a=s&vUSgI9V!`FE)hb~dOAkgFYS1Ws`ZxvLE%c+GWjF&+wGIHdhqi-CQqY-NKP$o{6O`Aw6Opvl}H{ z(j2XU{>JjkkE4^ba=lbU9<+b)zn02zOdmxqFAsQ{>%+)cKjfubNDWFageR2eMruwn zL^_|x3E!;Vs`##vH&^TL!Rb>WbE|F2x?59Rg(8H?l(4xe=gXP@nTC(Fw|P{3K7I3O zTL?5k7LJ2uH&xE(blHw$8)bbp-t^N-K8RzntH-%K@%{{;`GVd+M%L-2)f>CsoIaS* z#+;pVeO^X8?bN3IT2QTKaBsYs!`c_AaRcBqWyz5GVN7dp`1tzbO*uLI);VW{bUdvX zlbX1L6X&b-^tUOmQ(|~VTR8a4Q5X#|d2Jp3)1}74B%$%PXdZ#(K8b@nL)n7Q+qmZl zFJ%Dx@%L5Q3j;4%kv?+?EoJ#!ArFteLv1?KiyU}$X#`?r^t>gaBVN`;USK6!I;>$M zPEuqzg;T}rhi4t?A6)31B^K#xKn%#JII_j1o_7A4nZHMCGJhOTThMAWRRXRxUQhE8PV09BuB%Yag? zNuJvX%DZ(t$u_?=THDcbOv+Y?U1*k^#S`yo5j8WG@@oBAPn36VGV?|dUglj(g`4XN z_l;1x6`m`j)WdLfR}ORdv*$wOTrbI#moKEl?Jtd@t`#cBa&b)8XOx^c7a;sCH0Tlh zp8(nKkN$0`WmeUvu&W*FxQcJ*6SZM!A~>&0F5dLPy14p~mCe#QjMJRr&@qagu+!`wInZ{6h%>YRpnIFUBpE!H?;FU*kvnGY|iJAMFU!z$t^6o)fADW_n zr#CDY1TMLa$a5093$iz9-FMaEUMVJOh(MNR-W&q~X92rOm{uI76MsGPh!g_06vx+(| zqWN}>KH zzQPRdtAWVkD#@PxPx3nq`+Mrt`!;|heug2CXnYK@aL{>{LZC*;eq#ZX2!oANIg~J7 zT~9VzJh^R(89+vKI`?py+oAdw@VRMdPXWDnA` zhSoP8arXSI?RC7{8bEF8i-g%3Y9s4U!=qWwPY$}K^l7d^55D&o4t=-sDXNrj84{lS9StqTNyHN=0 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 d52e85493a1439d27febd6a1404805ee706e3437..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 41427 zcmeHw3wTu3x%M8Cfp7^kXt2Q>Q;QZAV*&v~m6DJQ?7<0*5OhL;Nk|3~o69trK)ki# zgecQ7ibvG+w6-~FYpSiL)ml{4gj<4I4Tu;k*3_cT7>W_mrrzd%zjfJ{Nd_{&cuKU_+@5%hw>0UU>Vo6ZMm7siAu^^qD%zd_x^NK!>{Fo1lF#{c3QSQga55|TEEi2Ah1IDZtRNUu@oH3~gR9VCED(j=Z$ z5<%&I4>=xzPUEm5Daq1R@K&IV(rpziJB+%8JxQhB70{zNzw9OBTq^vSy77`e2`B#!CSm?q9tWrDAIZ)+as2w@$S;i} z|D8CvJx)9`;@}H`pM(FgS7{vkuRw65a&3%*uaAS@6bHX04o(d`lHa3o@YCYNe}5eL zH{#&Gi-RY`vHyoS^4G_~e-As8l=q`xo33%p{siUl85~ZLh4*^bRn=5`>jEXU0k2o_ zx(f5XWxiV9x{A7hueLBhx3Z?%S6H&9(#LrtbG)SuCEm3a)g_e`*Ftdddhb$UA-A%m zuFhAdcvqBF7OYp+*7$wZN@-FHoB@&yWNl`3CV zRn2-BTwcAQqPh&x)YZco3fs`&_4{h;YADpeHQx2m@l}=juK{UTSYKTlsHmytW|ozf zRF>4j^^#?QtcsdC#8ch0wmT6-lJbQAp2s~|qPep%qjZ;*fK;U8cb9<3Y@Sm6AB$C8q$ zyn#GjFg5I~ZTe+c)z3wivv# zgo!ErG5BC3mk-C_?P*+oFa}TG$mNG(@E&2mN7$3`$b6fS=Pk``(`UN@H|}dX4fr@+ zTU_Ze_eWsKV^; z<&|Z?$%nkM4Y;g3;;9B49a6gr47jme9s@p6lU9^s18y_m}R;HXWzl6G>5=l&7u~X=C|ul&7H)*~0RpC{JA~(#Y~3BeZhr zN)bQHAEP{Vp-3^yze{=QI*|gFf0Obwq$1fY{};;BP>Q6p{C>*Q5Q;ch{#nXX7mL_g zeh1~LYef{6e}eMVr6Pl$1332&l&7u~>1X)|C{JA|(#!JqQ2rFkce4C{P@cL>q=V)E zo$}OGB5f@H6UtAf{1)WZ7Hca2wPhaR%hmIT)Zklc)6m2XqgL%iLzUI?hB~$?+E48P z`4TR_!XV!@l%NK8s9o>QS1rA2aF4q0^LnN~c?Q$pAAtUwq~AF2eCR3l-&ULEeFHqW z489eZpoZpsj9mBzgcL4Eu~+S}DiE}MuyZY4X|}1G5HfX@d5twt3rt0D_y3Rwcj|J) ze>Gs(iuSIGTm=J>rTAF8h0Nz%m2+jz@|?nD=zLkO;HPTvgPh>OW41}oc6{Wwslk7S z(Fv6N(xa~K#A)Ql)txtuO&gS8<2b+hAGS^Bp&zN+V3@9{uDQF{Tf-M1INK=NZTX?hY5A6;$82`zwzamj&g;fevMv12&!SNh=y10* zrK0OT{zgaAPoE%NwdKi%Q^|z=teGB0lZhgd!)tBRk2F2G8C=EIe1K$weQL;`boH-dse8JV(7Q1!536p-Z6Ku zN6i>igH5E<-9+W_V>VEOKWjuX^auf{(jZH~G@-7L#nsf6MaZ&1T zW8k^*${BZSfQ^N+=RFBTPC=?+%rVt=LFd3^QrZ_mK(4V}$6U4xx{ldAXyzho*OgS^ zlwN6lilrwxd(h1P-&S_hgH+iTHS`EovKkt#Htk7LgFJ^xq3h^))>KHnJ!a$ERL5H0 zw()I6+4G*vA3X&rt%})>jH1J)(+tit0s`+n4vC&iCU#>f6(c(XhM|Lu8+= zDHVeQsRP%J!{^JkO-XF9*qxo>{*Q?>ofMO&EvLp>oe%9R_Cj#)b*E7qCUY3Vw`;<7 zCuzWEk+TE7Fp1>YCXIGA?M}*{`I6vodE@&gvZR zs-3YDmB$tQ3=>sw`L?{E^HFE;MH=8+@R8q2h1;BG>9RfiY%qU2&0ddT>2;U-w!6BW z7&A0|ImPVd6qDCs+oq|2jFZ!rfv&Y8A!oY7_j9ueQvQeuY^eCH<%USO^mi^%H2UlK|2#_d5$NY9r$LNYHu?ScPlLO~ z7~MrJbaMK`O7(wt`g8Eh_;`JZjn@OVO^Mog%{Zoy*Lr4c*mh+# zhcR9czb}Pj$7>tguH-yTx6rPpU6#DLzU{VKS7HqvNQ<|Bh)e{Vuhq8ZDcHY6E(hQa zhOvq9#BO46pB4N>4Q)tLFT6zU`aFScAf6$2;faXK72Kr;yTh{vaq}6PrG}OzsplV4 zyFN}(oA={JuCPR+NS0KyAeq>v~_oTA6eU{c16xs zGoDo|_WOckDpkvl;C?l@b$17Gn>sD(Ji584f8lW=d}p-HqU5kti$AKCex!!Bf5@au z2lGOkUZx;imR?uzU%1CO+nw>ES`kfNJ4I1vaWS{0I|#o#Lk_UNnpXqkKRQ(nHBX1O z8f+$eaDzZMeoT&Cw%q+q5lel`<0RQq0Kfay(4#MNzqkSar`q)%it)Kx@sTg&S*6Yr zyz>tRKT?A{184Cak@fqx!6WX!n0vK_ksNkwnL`}Z_)|S^AihumHu?MbEu$Jrr;R++ z#=fBLJF3p=V0VmfTHL|CE?ZuY8azg|8?@6MBi-<;(XF-SEQ1S`@o(?3^ zhvo>|jBxYX?K8x<1q+0gb9<$?m3;uZ;l$83+KypCdXz*iswi*4%*~2R+XPeRIY`N@AXGPThZAIb3A+*)X24HZj zl}0`5;q&NDxc+d6b)%rO=Lzac=kEyis;z3vrrk^~=y0gLjY#K97pSG5roRvK2S1 z*drzypuXYInK+A$kt!Yz-8UQx4~J&p95=?6e>n8$aA;^4gxR}|Eh){dAkOBf?UsJb z85OTlr3UZ*6pNR*VQwBmnm2DVHO7)qCVE=4oh4`VKS>k6JNU;iq|i@Vh`NKn6WL~> z!N1Ce`_l}ZcI~t@ePXSDK1N82JTVvd<#T^U0oZQMhHb%xYq#p>hPe1-J##vQFkuO>+CQ=79a-|NP8Hk z0)g}q1X8d+&;$xb5P0W&gMfbof$J}j0&2qDTOeQxpq=d7CM7hndc#rTE{k*SJGNUd zW9R8l6(Iibv5(NxLwE2!$|!{HC-jgJU62HL3&FchupMv*!4H|>RKUFi-(`XwfCmX~ zG{MsW+YxxU!UU%SPA7P&3C;psK=7p|I2*8^;AtjU1-ymeWD{JVh8V0iNAX;z4u`e? zNSy;~I6!VVa0Lg_H6WVRJI-Ub64BAg)|rDhYTLr{Pj1nZZ6X;tcLSk}?W(kB@G6J&&Ai}IiPZo7Z!iE;{rX;iReQ2(~oy&R*tF z!jmAawoLxRQU{JEW{(Giw+eoK1vuYeocvzEN3w8kgQ*m}v;5%EWEKgKv}7dq}xyS2{bw z!Xf#5S{cXDW4`cJV6FOKC%q5B9{aI(6E+@c#l0dOowdx~Gm|2kw*w@+7{dHOtC%gX zILz9IgZ9W6+)1h}&VIawVB7QtA%DV$8Y=9@D;y)AV-@ZhVn+%55bMB+4zVu#K;;Km zCh_kCAD;<11sSCuqyK7u|C4e9G-Fr9S>6h8V(0h zY){t0a$o8oipzhAgQ(gi4)nIssyK)WzM-E*ovOCfCGocZ3f=Ml^%mTx?u)3}>nUtD z6DNvJ^l+Th;dA6h99V`)(LHSM*Y+wJ5--V3kEF0y2HvqvRl*Mvh|2Cxpq-yq(UAl$ z5)OeM8TbrB8+Xx!!sbkVhEDH4;iwj8^Uu&x1mCGYhnZImJxTB7a^MLLbZEeDIMAvA zG^ss?^_YdRYz04Y;WWM*lNL(Yy@gxJ50jFz$Z^fRhJqqKCTP|ZKYvEUr8L-k!i1{x83pGD*Ynfqp2Y+BE{qb%(`0ozS zgCww|fm7LIm^Z?UG5gT^bsrdN2vLR*We8D*>WB97mr}%=bJTNO!Q*iPeg@045?(@9 zXf7_FPM)725X;#mGbpg`uvd6jyK_@nu-&;1K}k~KM#e%)i|+xPsZe~0%yMq3Fg;zG zK2w+`K3qCrI$fB|5;U7AH~r8-n4ZK|)9X^hOJF;)9rdmTbJK<4Y-u<_7=92cFq{R$ zsxVj}s7E)v++a98{2_*4WCb^@3d0_0_}vS6Wh~`}3t+fd7%UgmPn1Xe>C1Y=YWU}5 zc!(`?xy8b;UmD&g4F3x%h`1bv8@R!yMxwaX?xxK=CCuM#Fkc*=N9MP1^9|hmgCs9K zJS;ri#XLON$SrOW7GbE4NOD)t9zczyur_{@jEC2ew{~`@zzhg)q$RywDZD)mU3c(7vJJbPAao(w z!&8*(6-n-GrNJ9T7yj2L@Rr8CbpjCHNK1N~F1!_RZ)6*G`#|VIvY#bbGAxqZ+tHWw zvUP$J4w1J)|I+#4AX-bhP&dtP`8aBpNAc4=0p3rV}4Of`8s#o&#MhabipcactPxKTDTAiR;5 z^!5wkZ8Y~rw&8*n6}pg|t|!w)l9%ni{d!H3@o)}#TgbhU0pX3bq&L6t)(TzJB-w@w zTGNELY(1%pB=_c@9%KfDH`0>c76@+_aBpNA-U@_WBt3evSR}c( zo_%`RD7vtZye;M4$bj%hTGCsJ@U|VgC>z;^i*gWAHou;15J~QBg~1yc5C3xr-ume7 zB*+X1Z=@x?y^RBS>XS>jH?j?HjlwRHoAu-tk>uWv^y+1!=)ytrwwik*1Hv0=NpIVP zw|?lNY-AfQT0unF+Vo_*NOEsKGI%57;WNnFVWdztG9bK>mh|=$;jNZ?Bir!S0V2F@ z)03Sd$-UVP-pF|PcmIU9tC4~?G9bK>mh`q#cuVBo$Tqz7fCz8Bda_R>xwm^>)N7KA zhv$zPWnk3us)(;}Qh4tj1NOEtP0S&^F(^MgC|nM_frB#4Zj4XP0|~G zMxM5!^jb0B^EZ8-5IDc-^Tfb+nm!*D7^{ZH;K80Uki`?$dDH1pR#7&p(W4BUBhYUP z6txKyU(R4k$YiJD!)oVjL02;xUiQ3v;mDKx;G1e8Ngm>v<`6Y!T9osd2~RSo zrfBmh+I1A|DWZL-ni^QFR|~}Y;<>zv48f9!jC6sNyI;r|tCGp_I8T=zWL+{hT9mgj zZBulMDLOwzH(Bbcp_^#%u>k9&rm`uS65#wRgONxZ1iISomi>>{);E8BD-Kj2^si zYwU?99J8N%i^?FQ(Z=jy0rPN%kJ%Iv&-1888BY?$GhW0KVg^Qx=LB_AWr{}D%~hsg zhSg0gXIqJeL*2YS#xrn4|0k@EfpbQzj{^dZS0B@uyVX{r;V~kHo0!tD7;t)Er@3?% z^=cRO>WO^hq9G2uY&+q(PQzHkHdw z<-&i`wkBOw&{{*=Dy1mZ4xLXtP`%j`byIg&*4$Tdo;` zO*DKLW80TgJ54?Z2qpkp>5@l78R;Fe%_a)1hk%>Q^~G(rf2v_%#I)ICf@dh&iQ4Q| zUK7@Lpceki&e*awFt$C}&q~7Qx=LW-hnc9M&CX+fd4I%{8|s@@ij$Y^3-(SH5j0ja z!`^A&Y+g)Ol;L})eS~%Jkbrq{hOdLugeRhUc~0mlII^eOFYxpS zQFHfza!+!;Eya9m+jNV9R|hfYe$%*?V3$ODY1rFm_xh1FNQ?^kfvdn z!v+Ud2%SIUi;6I&Z<1qmv>9M<;5;o5xgcq_DH;I|hWFyu4mID2|JJv7C9nmF`t-`= z*EoVM$?id9nSHFxw}~>N^BL9|xS{(yz0f2kmzVLP1Q=|5n=ebiRZOK`DU}YM9;et+ z$x!TX>#-Qhtp}u)dzDsh_I+bRncrbC7RfTxJXX!hd#_$zx}~$Dq;wZdB}I~Mj*(1? z(8icoe)5?-@&YREK2*9vg>D+@Yp@<&w}?7}CusmXezKqZ*D1_m0g5+Fo3L1Is)nQi zF?ms@0Jln7h{$JCUdD;Mnu?FT_zz1()<3aeJ(J&Kv}lLnzH|5;07Xr=`c&cpc?Xa@ zgGX#C8<&!0lXn2g79naZh!o`ok#_*_!6TRcMC3(?yeQelqqcHtx(gUN>egXVuTN#= zFh$L!hDCk35H&_EMdL&*{fVeO*85q>S@ha{!chyp{DQiiydQsN1b>_1jC2Qzy&C<+ znuMepLar!r5ZUv!?1VwlaM(8$VijhHonRIsnj3D12B%4*1m-o9Q-;&KC zaf`^bvN+UaO>)y;WH-edGZVeMJ-_R^(hM|4T3`0Yi7%k=&GxN+%XXcsP z(9AP8rkQ7MRWr}ryk?%coy|ORgPVD#mNzV3?GgIgy=?x|@O=$cSxz6ECLuO}v=qH}PVc-^7b) zeiJXI`Axi-<~Q+Tn%~5WX@2lLG;Lt?KY+PFG(Qlr>5xsXBeJ;kK=n*A7bgjmtro>x4f6fUqY-u#&6&{Mcx{ryobJVHpwZ$pD=Ke zvrP#7S}1jheD4Qbse_axc<+Zap<4sf9SD!Fv6va`%VF8GybVVsli#1jCXIZkdFV(H zrwA))n}$E1f^YAPafhx;)dOvf3G}{^0=?*6UK$z0pM{bP{EXwo7cq2D3{@hAwzx5* z>oJ6xxgu+?od@1Z_Ycpp%JB|gdk0;AA*-a7&D3{?>ED5u z0-;ZT=4Lf(1{?HgqoM2HDfFq257W=U7fXiur(KDze^BVtSQw_im+Kq+!sD)#)Rbe3!Z$3RqyCuCq@AH_G|22U}LW7XoS=mejl|}{r zNDq7I*NN{8!T$1g2ux05j*IT#LiQG+nB2;Xwu*W~AIghAW3&>3gWbc3!GU7%{n<*~ zbQR%m8zHn53XQgC2#v~6#LYNc!&m*DMP&V+sI%|!Nb@NQDlSVTyEtkXf@&zY7;O4; z95gCf@qslEw%LTS*K`POG{a(l%qCU?2Z45CD&m2%Jv(&S9$wVbUN%}*vc&>=*1K6% z*^@af%11=9QTSZ(2yL_Z-#+6{e+`wMf^Wl+bJXWbwwEW1mkW!W7UiQn28*NJp+^Qa ziy`&^neGxqcuR~Hvw3@4#|n#_7UiQVOcoO!OhV(Jpg*Gr(C8Jbq`4^G@zny%66T&@ z&0$aGw1{mrOctH^y8kXJw&AsReulwru4yyo!Jn#{PPgTBe{R=qh8f}jlI&wnRI zL}bT`o(i09$UXa^`xj8!&Dje8d+_J8JTKvY0d!3|>0d)f{Te#OUqgp_>sKs~ywm%t zYNEG)RkcpE;g=86e%6+8wDp?&s>1=bqSNu8{?HlC5XesgOKyR(sTn)@jLEwIJ-MON zASRvdKr>>vYiiWv5y)3n6W#i&s)=3JS5>3c>}v!`i`dst6L-{KHAq^izJ@n(XZ}^a z(R%bXyzv#}D>pYSF<;f&Dtws*f14Apm*e%5o#}zTpt~ zxOG(q=W#oOQoi%k-@Ly+e@(W?_jB5~Fnj;OH?Cc3{p=Yoe@-@^*I(7<_gU=PDdexn zW^#@B*?)eXu2h~^7K{9uJfB@C@Lu6h=2JM{Se(twK9Vl<(s_Q?%0Ii-ESoeJC_PNP3RYid|LL@=7B4=o&EBBC0Mhs%bFzYtrT+0MR~?@yz*+n|E|zi zIrU#7@-9xZmvh>=RN$v&U%fqWdurli^B-FDw`)_%M^RG)G&Q@OD z&a1dyk+5pp5cB1^U3=a+kP>h zw`+l@XFJp1^vZ^untLtVt)~P;d|mUeXYp+kaqJcHvi+1sMO>CM5-xRPLoprs^n0RF z+TAtY7mfCTp8uO@lx~%d;`TKQH17}5=yaSqd+<7n2Q&$<*ldBE1C+KvoA*Ye#U%f1 zG`bmd6khS{0eu>D2=sbtpfmB^png#HUGzl!X|=v+)B&1~ugYbE zZU-#}b^itWpwEN0ffl?5ebDKLpbwh(I=2htT&(!It1DXy1F0w zpu0djLDTTP_kPe;Q2Ozwn{hkn0ByuwW;W5o&1a2J;aS0VZB_&??TXmmIHj41y|60;31E6)+ohaQeb zuOvX)&%HUJDRI;iDM;l`p+16F=6qCcDd%iTNcg@*im`HE2L9@vk;_f>LYL**0ynGK zF+ZzMaVTfV`la~E*G%YBIT9~Lp{1|FrqMTBTM}mL z%Z{;9oYl`J$fuyZ^xHm(bq4tYL!Cw^F`a3Uoq#?+)ESLFjI#4iNDcz-l{tBL;;jUNEN7yO@^_**spR`3_U6pjAQ#BbO5kAmL^{@o^i ztHwV7{+0MH)bCCF-)j8h;D^DF7k>)Ib|J=oy!bQ0-v$1UP4=(V>=%MR4Rb}4iT@*w z9{_(d_;ePd)$^#AH2zlbvoXgUGVx#0_>Y4BDfsd94}kwD=A?M~$HDhwj*1t53c5-! z_?;&IPip=%!S`cMTWjK%X#7I(A3Fg*0RE@oUv1L&Y5H5iUx>Lf9{*ABw}AhFN&l#( ze*pY4%(3zK$HC9SoEtUif2`@JU~WytT>M89zeD3^g8vZs@%m>W_?KbscAE6FnLdr< z0QhUbKNkmpTKgNvq&=C*k~)p2`yfwWqyH)5*5nhmYLj6!Y8kiu6#QD!uPl&6uGw<2 z^*XpYNWZZJ=tYb5;ZaKaDC_0~Wn+T%odl&n!TMsN(vxU?!m2!Cwc^{1|2f)v=NRSI zG1g5<%4cJ&$C8wRBy0az<<+s)7se^wWW|H#IsXUUh8GV?n-W+G`oua%y&ic$08ADW^(F!6)L z55A$iIVDq3cARCsdy4YYv#mEzQ3lVpev+zmPqn_As{F!XeKl42&|&RKRi6K*^-rnF zz2{heovM6%j`faIW#6}~n^KjbZ^6&W>tqd_tbvm?aIywY*1*XcI9UTHYv5!JoUDQW zoi!kTM^*lAs{EZ)87&>T(6vq^WS*Tt0hf9D?N+*Eo_BWI(}j2G*(HD9mF~yrq9Z4|KKwjdLv&F1({4vZLL9RxsW@|> z3ycQcRZ+Q=~uq zEgKepx={Q~a2-HOsRDDp={LH{N4vYcG%f9})1(ZP(%-_QLb zPlKXf?UzAflSf{qp)`7}Y#6LgNCJ))el9QcXo?>MF}UsGQlsK3~aALg#Dab%`l zk~Zt2jCuxUoX_((Uaq^QDp0ZpG*HWFxlC2p1bk`hs_WC%)K^rNT~tw~uv~dbUAdB0 zc1<-jIStfu&U*Z`HvND$6B8M)(2`k**7_<-NMEGsHK0za&+q~T|?@rhp&uK=&6 ztRzsPr1{FdYimoYeBSah=*k>?;PB!?1Le$yx1_eVz5_@%zAzZ+uo<$NuvoZm^G+9_85CSX*jay{t~^N^%+oq`Cv zs0_yI-;got_X>SU<@zdhqt1w}}A*dw(iG)eNf^*sgrDux9E7t$RU>fxY+2J(SKS6iN@gR)I_LDXr zkI{eB&pDDZ<5F128nViB$lz^UHU+G>f*71+>5>%RThM$nhWHq@kta@-f&!1pGCOkY8TDI(Jm;s*iR&!NmAsxVPqv>S1#89Gd@ncZ!zQJ6*+I4@d=8YkIi_pcE4lB zCo0u^o)VS!*`q_FEN#mzc$*?Ne!9OL zg}cz$cA@N?SPskG(Ix3b@4(ryL(^89ux|X}u=<@_rjxcJXy|c)M%^ccgOd z25!O1(_g~%r*p#ABf=eP%=kVbUmSze`a|+1%k0PCV+7t7gC`5TKL)4yl+};VY=4ThAA?Vo_G9pK zrTrNELTR5m04^DShO{4pUn1?t;Pa&Y7(5p^^&fdoL;D}Pyi9(IGJT_ivM3|Js)ZfL zI-Yk3JRtDw82nn`BlS13Hxhmevtw5}{KBvxKWF&xaYyWt?6gCk$~D8y4X}3vklZ8i zmAM?hl@le*?58TdjShw+2|RHWw{uG-2T}w+P2l$ke7nFs0v}w-0ro5($@K#7UBK~> zkZ%=upUQzwfo~J|UcA#o7kl>!^d*K5pXXl{@{Ku6T)CSQq;{og0#SI_1+y%R;((*kdG@x18P%Aq*&9|=4wo6FPuNY|-}Be(A?ftM}ja$@_h zTp{rN*&LVkSq+@(Q=U7?`urioQ(^mTzh`#lCL>>U&&KNI-%i#UFW6J@HE*Mnms2iUuLNX{1c^jREd z@5&*$1UTNx8&<9rLOym}+#+y~=;vZPuRJ91vlnpt>|HJ-pA&d&|2Yht{AO)*P(r&R zLVlAduTyZ6NAq&sw~zy}KVK~H&tv+}5`mvNpUcZ}S|#vyfy;U0W`VyD(@!20cx=7x z0#5xTeWQapko~-$;S&{kJ}L3P3%pm13p5|Q-T{ucmxk588eqI;{cgMj!aqu6+!JFdXcf`TlzgZ(T)Qz*p-HRC!A)YpQ*9$_S9xTUO&;S6Q>Bq|#d!sHv^< zmee;WI8pXj`U1YPv^ke#WWn+XLS8!U<|ox&oIclHqwv%2vihp3YoKDtdFhbcBwAQr zQdv?}Uzb+u_bX}MjG1%kb5=V0%wV50+2<_wnaMt9v(Gu~^Ah%%#XgyCCezJix|vKj zlj&wM-AtyN$#gTBZYI;soU3>jxxcp{$L;mj)vxioyr#&!E0zV^Szd2N4L|1hW~4I{ z>6a*6uA0HnGbJjDiWaQ3tLT3(@sse_tMK`HLZ-LIz%a^4i_7zKpg*&IDtjxRga%b5W z#2biLqDulL_1-FPIpb&1HzQF(|D627rASuLRaOS$>?z8KRg&rQ`r_6XYBM)y32ZOO zh!Y@pU+9e+;4BvpPV|WHEnK+FS?DdyS>SfEDp?kYsYY49aoUy~m5*4Lxf#xaIQ58f zxMhgL+5qhqPmd}`lq_y{$l%>!#Bkm5yCO_ikm&_8(WwjbbF~g$R#)RKN8{lk)`$&; zvRs~9!rD6{W1TOM?`z1R$6V-Su~kV1L|4eH^ZCknq{Av;PMt46oe*7RWS4WPmCLn( z#k<-1l3Ll#xZW_kd8}%@c%J613zXCb$oQ}r(Q02bYKB+sC43z7UO&PtK$R;x^)QD# zDqnTHvR2nr)cL#(U#dRh_WC*KaF@#-H$471XgBW_@wFK{?x^I(1$%g{^Z7|{d3r?a ztgLj-DXUn!)>o^P1+ZwaJ$1 z)OliDy9)EY`py6&nf3*R>>;$5cyfba<}jz2mDxgJcI&U%U1*3RtgT#qMbfy 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