MOOSE - Multiscale Object Oriented Simulation Environment
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
testBiophysics.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 ** This program is part of 'MOOSE', the
3 ** Messaging Object Oriented Simulation Environment.
4 ** Copyright (C) 2003-2013 Upinder S. Bhalla. and NCBS
5 ** It is made available under the terms of the
6 ** GNU Lesser General Public License version 2.1
7 ** See the file COPYING.LIB for the full notice.
8 **********************************************************************/
9 
10 
11 #include "../basecode/header.h"
12 #include "../basecode/global.h"
13 #include "../shell/Shell.h"
14 #include "../utility/testing_macros.hpp"
15 
16 #include "CompartmentBase.h"
17 #include "Compartment.h"
18 
19 extern void testCompartment(); // Defined in Compartment.cpp
20 extern void testCompartmentProcess(); // Defined in Compartment.cpp
21 extern void testMarkovRateTable(); //Defined in MarkovRateTable.cpp
22 extern void testVectorTable(); //Defined in VectorTable.cpp
23 
24 
25 #ifdef DO_UNIT_TESTS
26 
27 // Use a larger value of runsteps when benchmarking
28 void testIntFireNetwork( unsigned int runsteps = 5 )
29 {
30  static const double thresh = 0.8;
31  static const double Vmax = 1.0;
32  static const double refractoryPeriod = 0.4;
33  static const double weightMax = 0.02;
34  static const double timestep = 0.2;
35  static const double delayMax = 4;
36  static const double delayMin = 0;
37  static const double connectionProbability = 0.1;
38  static const unsigned int NUM_TOT_SYN = 104831;
39  unsigned int size = 1024;
40  string arg;
41  Eref sheller( Id().eref() );
42  Shell* shell = reinterpret_cast< Shell* >( sheller.data() );
43 
44  Id fire = shell->doCreate( "IntFire", Id(), "network", size );
45  assert( fire.element()->getName() == "network" );
46 
47  Id i2 = shell->doCreate( "SimpleSynHandler", fire, "syns", size );
48  assert( i2.element()->getName() == "syns" );
49 
50  Id synId( i2.value() + 1 );
51  Element* syn = synId.element();
52  assert( syn->getName() == "synapse" );
53 
54  DataId di( 1 ); // DataId( data, field )
55  Eref syne( syn, di );
56 
57  ObjId mid = shell->doAddMsg( "Sparse", fire, "spikeOut",
58  ObjId( synId, 0 ), "addSpike" );
59 
60  SetGet2< double, long >::set( mid, "setRandomConnectivity",
61  connectionProbability, 5489UL );
62 
63  mid = shell->doAddMsg( "OneToOne", i2, "activationOut",
64  fire, "activation" );
65  assert( !mid.bad() );
66 
67  unsigned int nd = syn->totNumLocalField();
68  if ( Shell::numNodes() == 1 )
69  {
70  EXPECT_EQ( nd, NUM_TOT_SYN, "" );
71  }
72  else if ( Shell::numNodes() == 2 )
73  assert( nd == 52446 );
74  else if ( Shell::numNodes() == 3 )
75  //assert( nd == 34969 );
76  assert( nd == 35087 );
77  else if ( Shell::numNodes() == 4 )
78  assert( nd == 26381 );
79 
81  // Checking access to message info through SparseMsg on many nodes.
83  vector< ObjId > tgts;
84  vector< string > funcs;
85  ObjId oi( fire, 123 );
86  tgts = LookupField< string, vector< ObjId > >::
87  get( oi, "msgDests", "spikeOut" );
88  funcs = LookupField< string, vector< string > >::
89  get( oi, "msgDestFunctions", "spikeOut" );
90  assert( tgts.size() == funcs.size() );
91  /*
92  assert( tgts.size() == 116 );
93  assert( tgts[0] == ObjId( synId, 20, 11 ) );
94  assert( tgts[1] == ObjId( synId, 27, 15 ) );
95  assert( tgts[2] == ObjId( synId, 57, 14 ) );
96  assert( tgts[90] == ObjId( synId, 788, 15 ) );
97  assert( tgts[91] == ObjId( synId, 792, 12 ) );
98  assert( tgts[92] == ObjId( synId, 801, 17 ) );
99  */
100  for ( unsigned int i = 0; i < funcs.size(); ++i )
101  assert( funcs[i] == "addSpike" );
102 
104  // Here we have an interesting problem. The mtRand might be called
105  // by multiple threads if the above Set call is not complete.
106 
107  vector< double > origVm( size, 0.0 );
108 
109  // NOTE: From
110  moose::mtseed( 5489UL );
111  for ( unsigned int i = 0; i < size; ++i )
112  origVm[i] = moose::mtrand() * Vmax;
113 
114  double origVm100 = origVm[100];
115  double origVm900 = origVm[900];
116 
117  vector< double > temp;
118  temp.clear();
119  temp.resize( size, thresh );
120  bool ret = Field< double >::setVec( fire, "thresh", temp );
121  assert( ret );
122  temp.clear();
123  temp.resize( size, refractoryPeriod );
124  ret = Field< double >::setVec( fire, "refractoryPeriod", temp );
125  assert( ret );
126 
127  // cout << Shell::myNode() << ": fieldSize = " << fieldSize << endl;
128  vector< unsigned int > numSynVec;
129  Field< unsigned int >::getVec( i2, "numSynapses", numSynVec );
130  assert ( numSynVec.size() == size );
131  unsigned int numTotSyn = 0;
132  for ( unsigned int i = 0; i < size; ++i )
133  numTotSyn += numSynVec[i];
134  assert( numTotSyn == NUM_TOT_SYN );
135 
136  vector< vector< double > > weight( size );
137  for ( unsigned int i = 0; i < size; ++i )
138  {
139  weight[i].resize( numSynVec[i], 0.0 );
140  vector< double > delay( numSynVec[i], 0.0 );
141  for ( unsigned int j = 0; j < numSynVec[i]; ++j )
142  {
143  weight[i][ j ] = moose::mtrand() * weightMax;
144  delay[ j ] = delayMin + moose::mtrand() * ( delayMax - delayMin );
145  }
146  ret = Field< double >::setVec( ObjId( synId, i ), "weight", weight[i] );
147  assert( ret );
148  ret = Field< double >::setVec( ObjId( synId, i ), "delay", delay );
149  assert( ret );
150  }
151 
152  for ( unsigned int i = 0; i < size; ++i )
153  {
154  vector< double > retVec(0);
155  Field< double >::getVec( ObjId( synId, i ), "weight", retVec );
156  assert( retVec.size() == numSynVec[i] );
157  for ( unsigned int j = 0; j < numSynVec[i]; ++j )
158  {
159  ASSERT_DOUBLE_EQ("", retVec[j], weight[i][j] );
160  }
161  }
162 
163  // We have to have the SynHandlers called before the network of
164  // IntFires since the 'activation' message must be delivered within
165  // the same timestep.
166  shell->doUseClock("/network/syns", "process", 0 );
167  shell->doUseClock("/network", "process", 1 );
168  shell->doSetClock( 0, timestep );
169  shell->doSetClock( 1, timestep );
170  shell->doSetClock( 9, timestep );
171  shell->doReinit();
172  ret = Field< double >::setVec( fire, "Vm", origVm );
173  assert( ret );
174 
175  double retVm100 = Field< double >::get( ObjId( fire, 100 ), "Vm" );
176  double retVm900 = Field< double >::get( ObjId( fire, 900 ), "Vm" );
177  assert( fabs( retVm100 - origVm100 ) < 1e-6 );
178  assert( fabs( retVm900 - origVm900 ) < 1e-6 );
179 
180  shell->doStart( static_cast< double >( timestep * runsteps) + 0.0 );
181  if ( runsteps == 5 ) // default for unit tests, others are benchmarks
182  {
183  retVm100 = Field< double >::get( ObjId( fire, 100 ), "Vm" );
184  double retVm101 = Field< double >::get( ObjId( fire, 101 ), "Vm" );
185  double retVm102 = Field< double >::get( ObjId( fire, 102 ), "Vm" );
186  double retVm99 = Field< double >::get( ObjId( fire, 99 ), "Vm" );
187  retVm900 = Field< double >::get( ObjId( fire, 900 ), "Vm" );
188  double retVm901 = Field< double >::get( ObjId( fire, 901 ), "Vm" );
189  double retVm902 = Field< double >::get( ObjId( fire, 902 ), "Vm" );
190 
191  /*
192  ASSERT_DOUBLE_EQ("", retVm100, 0.00734036 );
193  ASSERT_DOUBLE_EQ("", retVm101, 0.246818 );
194  ASSERT_DOUBLE_EQ("", retVm102, 0.200087 );
195  ASSERT_DOUBLE_EQ("", retVm99, 0.0095779083 );
196  ASSERT_DOUBLE_EQ("", retVm900, 0.1150573482 );
197  ASSERT_DOUBLE_EQ("", retVm901, 0.289321534 );
198  ASSERT_DOUBLE_EQ("", retVm902, 0.01011172486 );
199  ASSERT_DOUBLE_EQ("", retVm100, 0.008593194687366486 );
200  ASSERT_DOUBLE_EQ("", retVm101, 0.24931678857743744 );
201  ASSERT_DOUBLE_EQ("", retVm102, 0.19668269662484533 );
202  ASSERT_DOUBLE_EQ("", retVm99, 0.00701607616202429 );
203  ASSERT_DOUBLE_EQ("", retVm900, 0.12097053045094018 );
204  ASSERT_DOUBLE_EQ("", retVm901, 0.2902593120492995 );
205  ASSERT_DOUBLE_EQ("", retVm902, 0.00237157280699805 );
206  ASSERT_DOUBLE_EQ("", retVm100, 0.015766608829826119 );
207  ASSERT_DOUBLE_EQ("", retVm101, 0.24405557875013356 );
208  ASSERT_DOUBLE_EQ("", retVm102, 0.20878261213859917 );
209  ASSERT_DOUBLE_EQ("", retVm99, 0.0081746848675747306 );
210  ASSERT_DOUBLE_EQ("", retVm900, 0.12525297735741736 );
211  ASSERT_DOUBLE_EQ("", retVm901, 0.28303358631241327 );
212  ASSERT_DOUBLE_EQ("", retVm902, 0.0096374021108587178 );
213  */
214 
215 #if 0
216  cout << endl;
217  cout << std::setprecision(12) << retVm100 << endl;
218  cout << std::setprecision(12) << retVm101 << endl;
219  cout << std::setprecision(12) << retVm102 << endl;
220  cout << std::setprecision(11) << retVm99 << endl;
221  cout << std::setprecision(12) << retVm900 << endl;
222  cout << std::setprecision(12) << retVm901 << endl;
223  cout << std::setprecision(12) << retVm902 << endl;
224 #endif
225  ASSERT_DOUBLE_EQ("", retVm100, 0.0752853031478);
226  ASSERT_DOUBLE_EQ("", retVm101, 0.226731547886 );
227  ASSERT_DOUBLE_EQ("", retVm102, 0.204294350789 );
228  ASSERT_DOUBLE_EQ("", retVm99, 0.2814616871 );
229  ASSERT_DOUBLE_EQ("", retVm900, 0.194820080944 );
230  ASSERT_DOUBLE_EQ("", retVm901, 0.0490452677121);
231  ASSERT_DOUBLE_EQ("", retVm902, 0.214295483534 );
232 
233  }
234  /*
235  cout << "testIntFireNetwork: Vm100 = " << retVm100 << ", " <<
236  retVm101 << ", " << retVm102 << ", " << retVm99 <<
237  ", " << Vm100 << endl;
238  cout << "Vm900 = " << retVm900 << ", "<< retVm901 << ", " <<
239  retVm902 << ", " << Vm900 << endl;
240  */
241 
242  cout << "." << flush;
243  shell->doDelete( fire );
244 }
245 
246 
247 static const double EREST = -0.07;
248 
249 
250 
251 #if 0
252 void testHHGateCreation()
253 {
254  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
255  vector< int > dims( 1, 1 );
256  Id nid = shell->doCreate( "Neutral", Id(), "n", dims );
257  Id comptId = shell->doCreate( "Compartment", nid, "compt", dims );
258  Id chanId = shell->doCreate( "HHChannel", nid, "Na", dims );
259  MsgId mid = shell->doAddMsg( "Single", ObjId( comptId ), "channel",
260  ObjId( chanId ), "channel" );
261  assert( mid != Msg::bad );
262 
263  // Check gate construction and removal when powers are assigned
264  vector< Id > kids;
265  HHChannel* chan = reinterpret_cast< HHChannel* >(chanId.eref().data());
266  assert( chan->xGate_ == 0 );
267  assert( chan->yGate_ == 0 );
268  assert( chan->zGate_ == 0 );
269  kids = Field< vector< Id > >::get( chanId, "children" );
270  assert( kids.size() == 3 );
271  assert( kids[0] == Id( chanId.value() + 1 ) );
272  assert( kids[1] == Id( chanId.value() + 2 ) );
273  assert( kids[2] == Id( chanId.value() + 3 ) );
274  assert( kids[0]()->dataHandler()->localEntries() == 0 );
275  assert( kids[1]()->dataHandler()->localEntries() == 0 );
276  assert( kids[2]()->dataHandler()->localEntries() == 0 );
277 
278  Field< double >::set( chanId, "Xpower", 1 );
279  assert( chan->xGate_ != 0 );
280  assert( chan->yGate_ == 0 );
281  assert( chan->zGate_ == 0 );
282  kids = Field< vector< Id > >::get( chanId, "children" );
283  assert( kids.size() == 3 );
284  assert( kids[0]()->dataHandler()->localEntries() == 1 );
285  assert( kids[1]()->dataHandler()->localEntries() == 0 );
286  assert( kids[2]()->dataHandler()->localEntries() == 0 );
287  // Read the size of the gateIds.
288 
289  Field< double >::set( chanId, "Xpower", 2 );
290  assert( chan->xGate_ != 0 );
291  assert( chan->yGate_ == 0 );
292  assert( chan->zGate_ == 0 );
293  assert( kids[0]()->dataHandler()->localEntries() == 1 );
294  assert( kids[1]()->dataHandler()->localEntries() == 0 );
295  assert( kids[2]()->dataHandler()->localEntries() == 0 );
296 
297  Field< double >::set( chanId, "Xpower", 0 );
298  assert( chan->xGate_ == 0 );
299  assert( chan->yGate_ == 0 );
300  assert( chan->zGate_ == 0 );
301  kids = Field< vector< Id > >::get( chanId, "children" );
302  assert( kids.size() == 3 );
303  // Even though gate was deleted, its Id is intact.
304  assert( kids[0] == Id( chanId.value() + 1 ) );
305  assert( kids[0]()->dataHandler()->localEntries() == 0 );
306  assert( kids[1]()->dataHandler()->localEntries() == 0 );
307  assert( kids[2]()->dataHandler()->localEntries() == 0 );
308 
309  Field< double >::set( chanId, "Xpower", 2 );
310  assert( chan->xGate_ != 0 );
311  assert( chan->yGate_ == 0 );
312  assert( chan->zGate_ == 0 );
313  kids = Field< vector< Id > >::get( chanId, "children" );
314  assert( kids.size() == 3 );
315  // New gate was created but has orig Id
316  assert( kids[0] == Id( chanId.value() + 1 ) );
317  assert( kids[0]()->dataHandler()->localEntries() == 1 );
318  assert( kids[1]()->dataHandler()->localEntries() == 0 );
319  assert( kids[2]()->dataHandler()->localEntries() == 0 );
320 
321  Field< double >::set( chanId, "Ypower", 3 );
322  assert( chan->xGate_ != 0 );
323  assert( chan->yGate_ != 0 );
324  assert( chan->zGate_ == 0 );
325  kids = Field< vector< Id > >::get( chanId, "children" );
326  assert( kids.size() == 3 );
327  assert( kids[0]()->dataHandler()->localEntries() == 1 );
328  assert( kids[1]()->dataHandler()->localEntries() == 1 );
329  assert( kids[2]()->dataHandler()->localEntries() == 0 );
330 
331  Field< double >::set( chanId, "Zpower", 1 );
332  assert( chan->xGate_ != 0 );
333  assert( chan->yGate_ != 0 );
334  assert( chan->zGate_ != 0 );
335  kids = Field< vector< Id > >::get( chanId, "children" );
336  assert( kids.size() == 3 );
337  assert( kids[0] == Id( chanId.value() + 1 ) );
338  assert( kids[1] == Id( chanId.value() + 2 ) );
339  assert( kids[2] == Id( chanId.value() + 3 ) );
340  assert( kids[0]()->dataHandler()->localEntries() == 1 );
341  assert( kids[1]()->dataHandler()->localEntries() == 1 );
342  assert( kids[2]()->dataHandler()->localEntries() == 1 );
343 
344  shell->doDelete( nid );
345  cout << "." << flush;
346 }
347 
348 
349 // AP measured in millivolts wrt EREST, at a sample interval of
350 // 100 usec
351 static double actionPotl[] = { 0,
352  1.2315, 2.39872, 3.51917, 4.61015, 5.69088, 6.78432, 7.91934, 9.13413,
353  10.482, 12.0417, 13.9374, 16.3785, 19.7462, 24.7909, 33.0981, 47.967,
354  73.3833, 98.7377, 105.652, 104.663, 101.815, 97.9996, 93.5261, 88.6248,
355  83.4831, 78.2458, 73.0175, 67.8684, 62.8405, 57.9549, 53.217, 48.6206,
356  44.1488, 39.772, 35.4416, 31.0812, 26.5764, 21.7708, 16.4853, 10.6048,
357  4.30989, -1.60635, -5.965, -8.34834, -9.3682, -9.72711,
358  -9.81085, -9.78371, -9.71023, -9.61556, -9.50965, -9.39655,
359  -9.27795, -9.15458, -9.02674, -8.89458, -8.75814, -8.61746,
360  -8.47254, -8.32341, -8.17008, -8.01258, -7.85093, -7.68517,
361  -7.51537, -7.34157, -7.16384, -6.98227, -6.79695, -6.60796,
362  -6.41542, -6.21944, -6.02016, -5.81769, -5.61219, -5.40381,
363  -5.19269, -4.979, -4.76291, -4.54459, -4.32422, -4.10197,
364  -3.87804, -3.65259, -3.42582, -3.19791, -2.96904, -2.7394,
365  -2.50915, -2.27848, -2.04755, -1.81652, -1.58556, -1.3548,
366  -1.12439, -0.894457, -0.665128, -0.436511, -0.208708, 0.0181928,
367  };
368 
369 // y(x) = (A + B * x) / (C + exp((x + D) / F))
370 // So: A = 0.1e6*(EREST+0.025), B = -0.1e6, C = -1.0, D = -(EREST+0.025)
371 // F = -0.01
372 
373 double Na_m_A( double v )
374 {
375  if ( fabs( EREST + 0.025 - v ) < 1e-6 )
376  v += 1e-6;
377  return 0.1e6 * ( EREST + 0.025 - v ) / ( exp ( ( EREST + 0.025 - v )/ 0.01 ) - 1.0 );
378 }
379 
380 // A = 4.0e3, B = 0, C = 0.0, D = -EREST, F = 0.018
381 double Na_m_B( double v )
382 {
383  return 4.0e3 * exp ( ( EREST - v ) / 0.018 );
384 }
385 
386 // A = 70, B = 0, C = 0, D = -EREST, F = 0.02
387 double Na_h_A( double v )
388 {
389  return 70.0 * exp ( ( EREST - v ) / 0.020 );
390 }
391 
392 // A = 1e3, B = 0, C = 1.0, D = -(EREST + 0.03 ), F = -0.01
393 double Na_h_B( double v )
394 {
395  return 1.0e3 / ( exp ( ( 0.030 + EREST - v )/ 0.01 ) + 1.0 );
396 }
397 
398 // A = 1e4 * (0.01 + EREST), B = -1e4, C = -1.0, D = -(EREST + 0.01 ), F = -0.01
399 double K_n_A( double v )
400 {
401  if ( fabs( EREST + 0.025 - v ) < 1e-6 )
402  v += 1e-6;
403 
404  return ( 1.0e4 * ( 0.01 + EREST - v ) ) / ( exp ( ( 0.01 + EREST - v ) / 0.01 ) - 1.0 );
405 }
406 
407 // A = 0.125e3, B = 0, C = 0, D = -EREST ), F = 0.08
408 double K_n_B( double v )
409 {
410  return 0.125e3 * exp ( (EREST - v ) / 0.08 );
411 }
412 
413 void testHHGateLookup()
414 {
415  Id shellId = Id();
416  HHGate gate( shellId, Id( 1 ) );
417  Eref er = Id(1).eref();
418  Qinfo q;
419  gate.setMin( er, &q, -2.0 );
420  gate.setMax( er, &q, 2.0 );
421  gate.setDivs( er, &q, 1 );
422  assert( gate.A_.size() == 2 );
423  assert( gate.B_.size() == 2 );
424  assert( gate.getDivs( er, &q ) == 1 );
425  ASSERT_DOUBLE_EQ("", gate.invDx_, 0.25 );
426  gate.A_[0] = 0;
427  gate.A_[1] = 4;
428  gate.lookupByInterpolation_ = 0;
429  ASSERT_DOUBLE_EQ("", gate.lookupA( -3 ), 0 );
430  ASSERT_DOUBLE_EQ("", gate.lookupA( -2 ), 0 );
431  ASSERT_DOUBLE_EQ("", gate.lookupA( -1.5 ), 0 );
432  ASSERT_DOUBLE_EQ("", gate.lookupA( -1 ), 0 );
433  ASSERT_DOUBLE_EQ("", gate.lookupA( -0.5 ), 0 );
434  ASSERT_DOUBLE_EQ("", gate.lookupA( 0 ), 0 );
435  ASSERT_DOUBLE_EQ("", gate.lookupA( 1 ), 0 );
436  ASSERT_DOUBLE_EQ("", gate.lookupA( 2 ), 4 );
437  ASSERT_DOUBLE_EQ("", gate.lookupA( 3 ), 4 );
438  gate.lookupByInterpolation_ = 1;
439  ASSERT_DOUBLE_EQ("", gate.lookupA( -3 ), 0 );
440  ASSERT_DOUBLE_EQ("", gate.lookupA( -2 ), 0 );
441  ASSERT_DOUBLE_EQ("", gate.lookupA( -1.5 ), 0.5 );
442  ASSERT_DOUBLE_EQ("", gate.lookupA( -1 ), 1 );
443  ASSERT_DOUBLE_EQ("", gate.lookupA( -0.5 ), 1.5 );
444  ASSERT_DOUBLE_EQ("", gate.lookupA( 0 ), 2 );
445  ASSERT_DOUBLE_EQ("", gate.lookupA( 1 ), 3 );
446  ASSERT_DOUBLE_EQ("", gate.lookupA( 2 ), 4 );
447  ASSERT_DOUBLE_EQ("", gate.lookupA( 3 ), 4 );
448 
449  gate.B_[0] = -1;
450  gate.B_[1] = 1;
451  double x = 0;
452  double y = 0;
453  gate.lookupBoth( -3 , &x, &y );
454  ASSERT_DOUBLE_EQ("", x, 0 );
455  ASSERT_DOUBLE_EQ("", y, -1 );
456  gate.lookupBoth( -2 , &x, &y );
457  ASSERT_DOUBLE_EQ("", x, 0 );
458  ASSERT_DOUBLE_EQ("", y, -1 );
459  gate.lookupBoth( -0.5, &x, &y );
460  ASSERT_DOUBLE_EQ("", x, 1.5 );
461  ASSERT_DOUBLE_EQ("", y, -0.25 );
462  gate.lookupBoth( 0, &x, &y );
463  ASSERT_DOUBLE_EQ("", x, 2 );
464  ASSERT_DOUBLE_EQ("", y, 0 );
465  gate.lookupBoth( 1.5, &x, &y );
466  ASSERT_DOUBLE_EQ("", x, 3.5 );
467  ASSERT_DOUBLE_EQ("", y, 0.75 );
468  gate.lookupBoth( 100000, &x, &y );
469  ASSERT_DOUBLE_EQ("", x, 4 );
470  ASSERT_DOUBLE_EQ("", y, 1 );
471 
472  cout << "." << flush;
473 }
474 
475 void testHHGateSetup()
476 {
477  Id shellId = Id();
478  HHGate gate( shellId, Id( 1 ) );
479  Eref er = Id(1).eref();
480  Qinfo q;
481 
482  vector< double > parms;
483  // Try out m-gate of NA.
484 // For the alpha:
485 // A = 0.1e6*(EREST*0.025), B = -0.1e6, C= -1, D= -(EREST+0.025), F = -0.01
486 // For the beta: A = 4.0e3, B = 0, C = 0.0, D = -EREST, F = 0.018
487 // xdivs = 100, xmin = -0.1, xmax = 0.05
488  unsigned int xdivs = 100;
489  double xmin = -0.1;
490  double xmax = 0.05;
491  parms.push_back( 0.1e6 * ( EREST + 0.025 ) ); // A alpha
492  parms.push_back( -0.1e6 ); // B alpha
493  parms.push_back( -1 ); // C alpha
494  parms.push_back( -(EREST + 0.025 ) ); // D alpha
495  parms.push_back( -0.01 ); // F alpha
496 
497  parms.push_back( 4e3 ); // A beta
498  parms.push_back( 0 ); // B beta
499  parms.push_back( 0 ); // C beta
500  parms.push_back( -EREST ); // D beta
501  parms.push_back( 0.018 ); // F beta
502 
503  parms.push_back( xdivs );
504  parms.push_back( xmin );
505  parms.push_back( xmax );
506 
507  gate.setupAlpha( er, &q, parms );
508  assert( gate.A_.size() == xdivs + 1 );
509  assert( gate.B_.size() == xdivs + 1 );
510 
511  double x = xmin;
512  double dx = (xmax - xmin) / xdivs;
513  for ( unsigned int i = 0; i <= xdivs; ++i )
514  {
515  double ma = Na_m_A( x );
516  double mb = Na_m_B( x );
517  ASSERT_DOUBLE_EQ("", gate.A_[i], ma );
518  ASSERT_DOUBLE_EQ("", gate.B_[i], ma + mb );
519  x += dx;
520  }
521 
522  cout << "." << flush;
523 }
524 
526 // Check construction and result of HH squid simulation
528 
530 Id makeSquid()
531 {
532  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
533  vector< int > dims( 1, 1 );
534  Id nid = shell->doCreate( "Neutral", Id(), "n", dims );
535  Id comptId = shell->doCreate( "Compartment", nid, "compt", dims );
536  Id naId = shell->doCreate( "HHChannel", comptId, "Na", dims );
537  MsgId mid = shell->doAddMsg( "Single", ObjId( comptId ), "channel",
538  ObjId( naId ), "channel" );
539  assert( mid != Msg::bad );
540  Id kId = shell->doCreate( "HHChannel", comptId, "K", dims );
541  mid = shell->doAddMsg( "Single", ObjId( comptId ), "channel",
542  ObjId( kId ), "channel" );
543  assert( mid != Msg::bad );
544 
546  // set up compartment properties
548 
549  Field< double >::set( comptId, "Cm", 0.007854e-6 );
550  Field< double >::set( comptId, "Ra", 7639.44e3 ); // does it matter?
551  Field< double >::set( comptId, "Rm", 424.4e3 );
552  Field< double >::set( comptId, "Em", EREST + 0.010613 );
553  Field< double >::set( comptId, "inject", 0.1e-6 );
554  Field< double >::set( comptId, "initVm", EREST );
555 
557  // set up Na channel properties
559  Field< double >::set( naId, "Gbar", 0.94248e-3 );
560  Field< double >::set( naId, "Ek", EREST + 0.115 );
561  Field< double >::set( naId, "Xpower", 3.0 );
562  Field< double >::set( naId, "Ypower", 1.0 );
563 
565  // set up K channel properties
567  Field< double >::set( kId, "Gbar", 0.282743e-3 );
568  Field< double >::set( kId, "Ek", EREST - 0.012 );
569  Field< double >::set( kId, "Xpower", 4.0 );
570 
572  // set up m-gate of Na.
574  vector< Id > kids; // These are the HHGates.
575  kids = Field< vector< Id > >::get( naId, "children" );
576  assert( kids.size() == 3 );
577  vector< double > parms;
578 // For the alpha:
579 // A = 0.1e6*(EREST*0.025), B = -0.1e6, C= -1, D= -(EREST+0.025), F = -0.01
580 // For the beta: A = 4.0e3, B = 0, C = 0.0, D = -EREST, F = 0.018
581 // xdivs = 100, xmin = -0.1, xmax = 0.05
582  unsigned int xdivs = 150;
583  double xmin = -0.1;
584  double xmax = 0.05;
585  parms.push_back( 0.1e6 * ( EREST + 0.025 ) ); // A alpha
586  parms.push_back( -0.1e6 ); // B alpha
587  parms.push_back( -1 ); // C alpha
588  parms.push_back( -(EREST + 0.025 ) ); // D alpha
589  parms.push_back( -0.01 ); // F alpha
590 
591  parms.push_back( 4e3 ); // A beta
592  parms.push_back( 0 ); // B beta
593  parms.push_back( 0 ); // C beta
594  parms.push_back( -EREST ); // D beta
595  parms.push_back( 0.018 ); // F beta
596 
597  parms.push_back( xdivs );
598  parms.push_back( xmin );
599  parms.push_back( xmax );
600 
601  SetGet1< vector< double > >::set( kids[0], "setupAlpha", parms );
602  Field< bool >::set( kids[0], "useInterpolation", 1 );
603 
605  // set up h-gate of NA.
607  // Alpha rates: A = 70, B = 0, C = 0, D = -EREST, F = 0.02
608  // Beta rates: A = 1e3, B = 0, C = 1.0, D = -(EREST + 0.03 ), F = -0.01
609  parms.resize( 0 );
610  parms.push_back( 70 );
611  parms.push_back( 0 );
612  parms.push_back( 0 );
613  parms.push_back( -EREST );
614  parms.push_back( 0.02 );
615 
616  parms.push_back( 1e3 ); // A beta
617  parms.push_back( 0 ); // B beta
618  parms.push_back( 1 ); // C beta
619  parms.push_back( -( EREST + 0.03 ) ); // D beta
620  parms.push_back( -0.01 ); // F beta
621 
622  parms.push_back( xdivs );
623  parms.push_back( xmin );
624  parms.push_back( xmax );
625 
626  SetGet1< vector< double > >::set( kids[1], "setupAlpha", parms );
627  Field< bool >::set( kids[1], "useInterpolation", 1 );
628 
629  // Check parameters
630  vector< double > A = Field< vector< double > >::get( kids[1], "tableA");
631  vector< double > B = Field< vector< double > >::get( kids[1], "tableB");
632  double x = xmin;
633  double dx = (xmax - xmin) / xdivs;
634  for ( unsigned int i = 0; i <= xdivs; ++i )
635  {
636  double ha = Na_h_A( x );
637  double hb = Na_h_B( x );
638  ASSERT_DOUBLE_EQ("", A[i], ha );
639  ASSERT_DOUBLE_EQ("", B[i], ha + hb );
640  x += dx;
641  }
642 
644  // set up n-gate of K.
646  // Alpha rates: A = 1e4 * (0.01 + EREST), B = -1e4, C = -1.0,
647  // D = -(EREST + 0.01 ), F = -0.01
648  // Beta rates: A = 0.125e3, B = 0, C = 0, D = -EREST ), F = 0.08
649  kids = Field< vector< Id > >::get( kId, "children" );
650  parms.resize( 0 );
651  parms.push_back( 1e4 * (0.01 + EREST) );
652  parms.push_back( -1e4 );
653  parms.push_back( -1.0 );
654  parms.push_back( -( EREST + 0.01 ) );
655  parms.push_back( -0.01 );
656 
657  parms.push_back( 0.125e3 ); // A beta
658  parms.push_back( 0 ); // B beta
659  parms.push_back( 0 ); // C beta
660  parms.push_back( -EREST ); // D beta
661  parms.push_back( 0.08 ); // F beta
662 
663  parms.push_back( xdivs );
664  parms.push_back( xmin );
665  parms.push_back( xmax );
666 
667  SetGet1< vector< double > >::set( kids[0], "setupAlpha", parms );
668  Field< bool >::set( kids[0], "useInterpolation", 1 );
669 
670  // Check parameters
671  A = Field< vector< double > >::get( kids[0], "tableA" );
672  B = Field< vector< double > >::get( kids[0], "tableB" );
673  x = xmin;
674  for ( unsigned int i = 0; i <= xdivs; ++i )
675  {
676  double na = K_n_A( x );
677  double nb = K_n_B( x );
678  if ( i != 40 && i != 55 )
679  {
680  // Annoying near-misses due to different ways of handling
681  // singularity. I think the method in the HHGate is better.
682  ASSERT_DOUBLE_EQ("", A[i], na );
683  ASSERT_DOUBLE_EQ("", B[i], na + nb );
684  }
685  x += dx;
686  }
687  return nid;
688 }
689 
690 void testHHChannel()
691 {
692  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
693  vector< int > dims( 1, 1 );
694 
695  Id nid = makeSquid();
696  Id comptId( "/n/compt" );
697 
698  Id tabId = shell->doCreate( "Table", nid, "tab", dims );
699  MsgId mid = shell->doAddMsg( "Single", ObjId( tabId, 0 ),
700  "requestOut", ObjId( comptId, 0 ), "get_Vm" );
701  assert( mid != Msg::bad );
703  // Schedule, Reset, and run.
705 
706  shell->doSetClock( 0, 1.0e-5 );
707  shell->doSetClock( 1, 1.0e-5 );
708  shell->doSetClock( 2, 1.0e-5 );
709  shell->doSetClock( 3, 1.0e-4 );
710 
711  shell->doUseClock( "/n/compt", "init", 0 );
712  shell->doUseClock( "/n/compt", "process", 1 );
713  // shell->doUseClock( "/n/compt/##", "process", 2 );
714  shell->doUseClock( "/n/compt/Na,/n/compt/K", "process", 2 );
715  shell->doUseClock( "/n/tab", "process", 3 );
716 
717  shell->doReinit();
718  shell->doReinit();
719  shell->doStart( 0.01 );
720 
722  // Check output
724  vector< double > vec = Field< vector< double > >::get( tabId, "vector" );
725  // assert( vec.size() == 101 );
726  double delta = 0;
727  for ( unsigned int i = 0; i < 100; ++i )
728  {
729  double ref = EREST + actionPotl[i] * 0.001;
730  delta += (vec[i] - ref) * (vec[i] - ref);
731  }
732  assert( sqrt( delta/100 ) < 0.001 );
733 
735  // Clear it all up
737  shell->doDelete( nid );
738  cout << "." << flush;
739 }
740 
741 #endif //#if 0
742 // Markov Channel unit tests.
745 
746 //Sample current obtained from channel in Chapter 20, Sakmann & Neher, Pg. 603.
747 //The current is sampled at intervals of 10 usec.
748 static double sampleCurrent[] = {0.0, // This is to handle the value sent by ChanBase on reinit
749  0.0000000e+00, 3.0005743e-26, 1.2004594e-25, 2.7015505e-25, 4.8036751e-25, 7.5071776e-25,
750  1.0812402e-24, 1.4719693e-24, 1.9229394e-24, 2.4341850e-24, 3.0057404e-24, 3.6376401e-24,
751  4.3299183e-24, 5.0826095e-24, 5.8957481e-24, 6.7693684e-24, 7.7035046e-24, 8.6981913e-24,
752  9.7534627e-24, 1.0869353e-23, 1.2045897e-23, 1.3283128e-23, 1.4581082e-23, 1.5939791e-23,
753  1.7359292e-23, 1.8839616e-23, 2.0380801e-23, 2.1982878e-23, 2.3645883e-23, 2.5369850e-23,
754  2.7154813e-23, 2.9000806e-23, 3.0907863e-23, 3.2876020e-23, 3.4905309e-23, 3.6995766e-23,
755  3.9147423e-23, 4.1360317e-23, 4.3634480e-23, 4.5969946e-23, 4.8366751e-23, 5.0824928e-23,
756  5.3344511e-23, 5.5925535e-23, 5.8568033e-23, 6.1272040e-23, 6.4037589e-23, 6.6864716e-23,
757  6.9753453e-23, 7.2703835e-23, 7.5715897e-23, 7.8789672e-23, 8.1925194e-23, 8.5122497e-23,
758  8.8381616e-23, 9.1702584e-23, 9.5085435e-23, 9.8530204e-23, 1.0203692e-22, 1.0560563e-22,
759  1.0923636e-22, 1.1292913e-22, 1.1668400e-22, 1.2050099e-22, 1.2438013e-22, 1.2832146e-22,
760  1.3232502e-22, 1.3639083e-22, 1.4051894e-22, 1.4470937e-22, 1.4896215e-22, 1.5327733e-22,
761  1.5765494e-22, 1.6209501e-22, 1.6659757e-22, 1.7116267e-22, 1.7579032e-22, 1.8048057e-22,
762  1.8523345e-22, 1.9004900e-22, 1.9492724e-22, 1.9986821e-22, 2.0487195e-22, 2.0993849e-22,
763  2.1506786e-22, 2.2026010e-22, 2.2551524e-22, 2.3083331e-22, 2.3621436e-22, 2.4165840e-22,
764  2.4716548e-22, 2.5273563e-22, 2.5836888e-22, 2.6406527e-22, 2.6982483e-22, 2.7564760e-22,
765  2.8153360e-22, 2.8748287e-22, 2.9349545e-22, 2.9957137e-22, 3.0571067e-22
766  };
767 
768 void testMarkovGslSolver()
769 {
770  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
771  unsigned size = 1;
772 
773  Id nid = shell->doCreate( "Neutral", Id(), "n", size );
774  Id comptId = shell->doCreate( "Compartment", nid, "compt", size );
775  Id rateTabId = shell->doCreate( "MarkovRateTable", comptId, "rateTab", size );
776  Id mChanId = shell->doCreate( "MarkovChannel", comptId, "mChan", size );
777  Id gslSolverId = shell->doCreate( "MarkovGslSolver", comptId, "gslSolver", size );
778 
779  Id tabId = shell->doCreate( "Table", nid, "tab", size );
780 
781  ObjId mid = shell->doAddMsg( "Single", ObjId( comptId ), "channel",
782  ObjId( mChanId ), "channel" );
783  assert( !mid.bad() );
784 
785  mid = shell->doAddMsg( "Single", ObjId( comptId ), "channel",
786  ObjId( rateTabId ), "channel" );
787  assert( !mid.bad() );
788  mid = shell->doAddMsg( "Single", ObjId( gslSolverId ), "stateOut",
789  ObjId( mChanId ), "handleState" );
790  assert( !mid.bad() );
791 
792  mid = shell->doAddMsg("Single", ObjId( rateTabId ), "instratesOut",
793  ObjId( gslSolverId ), "handleQ" );
794 
795  mid = shell->doAddMsg( "Single", ObjId( tabId, 0 ), "requestOut",
796  ObjId( mChanId, 0 ), "getIk" );
797  assert( !mid.bad() );
798 
800  // set up compartment properties
802 
803  Field< double >::set( comptId, "Cm", 0.007854e-6 );
804  Field< double >::set( comptId, "Ra", 7639.44e3 ); // does it matter?
805  Field< double >::set( comptId, "Rm", 424.4e3 );
806  Field< double >::set( comptId, "Em", -0.1 );
807  Field< double >::set( comptId, "inject", 0 );
808  Field< double >::set( comptId, "initVm", -0.1 );
809 
811  //
812  //Setup of Markov Channel.
813  //This is a simple 5-state channel model taken from Chapter 20, "Single-Channel
814  //Recording", Sakmann & Neher.
815  //All the transition rates are constant.
816  //
818 
819  //Setting number of states, number of open states.
820  Field< unsigned int >::set( mChanId, "numStates", 5 );
821  Field< unsigned int >::set( mChanId, "numOpenStates", 2 );
822 
823  //Setting initial state of system.
824  vector< double > initState;
825 
826  initState.push_back( 0.0 );
827  initState.push_back( 0.0 );
828  initState.push_back( 0.0 );
829  initState.push_back( 0.0 );
830  initState.push_back( 1.0 );
831 
832  Field< vector< double > >::set( mChanId, "initialState", initState );
833 
834  vector< string > stateLabels;
835 
836  stateLabels.push_back( "O1" );
837  stateLabels.push_back( "O2" );
838  stateLabels.push_back( "C1" );
839  stateLabels.push_back( "C2" );
840  stateLabels.push_back( "C3" );
841 
842  Field< vector< string > >::set( mChanId, "labels", stateLabels );
843 
844  vector< double > gBars;
845 
846  gBars.push_back( 40e-12 );
847  gBars.push_back( 50e-12 );
848 
849  Field< vector< double > >::set( mChanId, "gbar", gBars );
850 
851  //Setting up rate tables.
852  SetGet1< unsigned int >::set( rateTabId, "init", 5 );
853 
854  //Filling in values into one parameter rate table. Note that all rates here
855  //are constant.
856  SetGet3< unsigned int, unsigned int, double >::
857  set( rateTabId, "setconst", 1, 2, 0.05 );
858  SetGet3< unsigned int, unsigned int, double >::
859  set( rateTabId, "setconst", 1, 4, 3 );
860  SetGet3< unsigned int, unsigned int, double >::
861  set( rateTabId, "setconst", 2, 1, 0.00066667 );
862  SetGet3< unsigned int, unsigned int, double >::
863  set( rateTabId, "setconst", 2, 3, 0.5 );
864  SetGet3< unsigned int, unsigned int, double >::
865  set( rateTabId, "setconst", 3, 2, 15 );
866  SetGet3< unsigned int, unsigned int, double >::
867  set( rateTabId, "setconst", 3, 4, 4 );
868  SetGet3< unsigned int, unsigned int, double >::
869  set( rateTabId, "setconst", 4, 1, 0.015 );
870  SetGet3< unsigned int, unsigned int, double >::
871  set( rateTabId, "setconst", 4, 3, 0.05 );
872  SetGet3< unsigned int, unsigned int, double >::
873  set( rateTabId, "setconst", 4, 5, 2.0 );
874  SetGet3< unsigned int, unsigned int, double >::
875  set( rateTabId, "setconst", 5, 4, 0.01 );
876 
877  //Setting initial state of the solver. Once this is set, the solver object
878  //will send out messages containing the updated state to the channel object.
879  SetGet1< vector< double > >::set( gslSolverId, "init", initState );
880 
881  shell->doSetClock( 0, 1.0e-5 );
882  shell->doSetClock( 1, 1.0e-5 );
883  shell->doSetClock( 2, 1.0e-5 );
884  shell->doSetClock( 3, 1.0e-5 );
885  shell->doSetClock( 4, 1.0e-5 );
886 
887  //Voltage is clamped to -100 mV in the example. Hence, we skip running the
888  //process function.
889  shell->doUseClock( "/n/compt", "init", 0 );
890  shell->doUseClock( "/n/compt", "process", 1 );
891  shell->doUseClock( "/n/compt/rateTab", "process", 1 );
892  shell->doUseClock( "/n/compt/gslSolver", "process", 1 );
893  shell->doUseClock( "/n/compt/mChan,/n/tab", "process", 2 );
894 
895  shell->doReinit( );
896  shell->doReinit( );
897  shell->doStart( 1.0e-3 );
898 
899  vector< double > vec = Field< vector< double > >::get( tabId, "vector" );
900 
901  for ( unsigned i = 0; i < 101; ++i )
902  {
903  if (!doubleEq( sampleCurrent[i] * 1e25, vec[i] * 1e25 ))
904  {
905  cout << "testMarkovGslSolver: sample=" << sampleCurrent[i]*1e25 << " calculated=" << vec[i]*1e25 << endl;
906  }
907  ASSERT_DOUBLE_EQ("", sampleCurrent[i] * 1e25, vec[i] * 1e25 );
908  }
909  //Currents involved here are incredibly small. Scaling them up is necessary
910  //for the doubleEq function to do its job.
911 
912  shell->doDelete( nid );
913  cout << "." << flush;
914 }
915 
917 //The testMarkovGslSolver() function includes the MarkovChannel object, but
918 //is a rather trivial case, in that the rates are all constant.
919 //This test simultaneously tests the MarkovChannel, MarkovGslSolver,
920 //MarkovSolverBase and MarkovSolver classes.
921 //This test involves simulating the 4-state NMDA channel model specified
922 //in the following paper :
923 //"Voltage Dependence of NMDA-Activated Macroscopic Conductances Predicted
924 //by Single-Channel Kinetics", Craig E. Jahr and Charles F. Stevens, The Journal
925 //of Neuroscience, 1990, 10(9), pp. 3178-3182.
926 //It is expected that the MarkovGslSolver and the MarkovSolver objects will
927 //give the same answer.
928 //
929 //Note that this is different from the NMDAChan test which involves synapses.
931 void testMarkovChannel()
932 {
933  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
934  unsigned size = 1;
935 
936  Id nid = shell->doCreate( "Neutral", Id(), "n", size );
937 
938  Id gslComptId = shell->doCreate( "Compartment", nid, "gslCompt", size );
939  Id exptlComptId = shell->doCreate( "Compartment", nid,
940  "exptlCompt", size );
941 
942  Id gslRateTableId = shell->doCreate( "MarkovRateTable", gslComptId,
943  "gslRateTable", size );
944  Id exptlRateTableId = shell->doCreate( "MarkovRateTable", exptlComptId,
945  "exptlRateTable", size );
946 
947  Id mChanGslId = shell->doCreate( "MarkovChannel", gslComptId,
948  "mChanGsl", size );
949  Id mChanExptlId = shell->doCreate( "MarkovChannel", exptlComptId,
950  "mChanExptl", size );
951 
952  Id gslSolverId = shell->doCreate( "MarkovGslSolver", gslComptId,
953  "gslSolver", size );
954  Id exptlSolverId = shell->doCreate( "MarkovSolver", exptlComptId,
955  "exptlSolver", size );
956 
957  Id gslTableId = shell->doCreate( "Table", nid, "gslTable", size );
958  Id exptlTableId = shell->doCreate( "Table", nid, "exptlTable", size );
959 
960  Id int2dTableId = shell->doCreate( "Interpol2D", nid, "int2dTable", size );
961  Id vecTableId = shell->doCreate( "VectorTable", nid, "vecTable", size );
962 
963  vector< double > table1d;
964  vector< vector< double > > table2d;
965 
967  //Setting up the messaging.
969 
970  //Connecting Compartment and MarkovChannel objects.
971  //Compartment sends Vm to MarkovChannel object. The MarkovChannel,
972  //via its ChanBase base class, sends back the conductance and current through
973  //it.
974  ObjId mid = shell->doAddMsg( "Single", ObjId( gslComptId ), "channel",
975  ObjId( mChanGslId ), "channel" );
976  assert( !mid.bad() );
977 
978  mid = shell->doAddMsg( "Single", ObjId( exptlComptId ), "channel",
979  ObjId( mChanExptlId ), "channel" );
980  assert( !mid.bad() );
981 
983  //Connecting up the MarkovGslSolver.
985 
986  //Connecting Compartment and MarkovRateTable.
987  //The MarkovRateTable's job is to send out the instantaneous rate matrix,
988  //Q, to the solver object(s).
989  //In order to do so, the MarkovRateTable object needs information on
990  //Vm and ligand concentration to look up the rate from the table provided
991  //by the user. Hence, the need of the connection to the Compartment object.
992  //However, unlike a channel object, the MarkovRateTable object does not
993  //return anything to the Compartment directly, and communicates only with the
994  //solvers.
995  mid = shell->doAddMsg( "Single", ObjId( gslComptId ), "channel",
996  ObjId( gslRateTableId ), "channel" );
997  assert( !mid.bad() );
998 
999  //Connecting the MarkovRateTable with the MarkovGslSolver object.
1000  //As mentioned earlier, the MarkovRateTable object sends out information
1001  //about Q to the MarkovGslSolver. The MarkovGslSolver then churns out
1002  //the state of the system for the next time step.
1003  mid = shell->doAddMsg("Single", ObjId( gslRateTableId ), "instratesOut",
1004  ObjId( gslSolverId ), "handleQ" );
1005 
1006  //Connecting MarkovGslSolver with MarkovChannel.
1007  //The MarkovGslSolver object, upon computing the state of the channel,
1008  //sends this information to the MarkovChannel object. The MarkovChannel
1009  //object will compute the expected conductance of the channel and send
1010  //this information to the compartment.
1011  mid = shell->doAddMsg( "Single", ObjId( gslSolverId ), "stateOut",
1012  ObjId( mChanGslId ), "handleState" );
1013  assert( !mid.bad() );
1014 
1016  //Connecting up the MarkovSolver class.
1018 
1019  //Connecting the MarkovSolver and Compartment.
1020  //The MarkovSolver derives from the MarkovSolverBase class.
1021  //The base class need Vm and ligand concentration information to
1022  //perform lookup and interpolation on the matrix exponential lookup
1023  //tables.
1024  mid = shell->doAddMsg( "Single", ObjId( exptlComptId ), "channel",
1025  ObjId( exptlSolverId ), "channel" );
1026  assert( !mid.bad() );
1027 
1028  mid = shell->doAddMsg( "Single", ObjId( exptlSolverId ), "stateOut",
1029  ObjId( mChanExptlId ), "handleState" );
1030  assert( !mid.bad() );
1031 
1033  //Connecting up the Table objects to cross-check values.
1035 
1036  //Get the current values from the GSL solver based channel.
1037  mid = shell->doAddMsg( "Single", ObjId( gslTableId ), "requestOut",
1038  ObjId( mChanGslId ), "getIk" );
1039  assert( !mid.bad() );
1040 
1041  //Get the current values from the matrix exponential solver based channel.
1042  mid = shell->doAddMsg( "Single", ObjId( exptlTableId ), "requestOut",
1043  ObjId( mChanExptlId ), "getIk" );
1044  assert( !mid.bad() );
1045 
1047  //Compartment properties. Identical to ones used in testHHChannel()
1048  //barring a few modifications.
1050 
1051  Field< double >::set( gslComptId, "Cm", 0.007854e-6 );
1052  Field< double >::set( gslComptId, "Ra", 7639.44e3 ); // does it matter?
1053  Field< double >::set( gslComptId, "Rm", 424.4e3 );
1054  Field< double >::set( gslComptId, "Em", EREST + 0.1 );
1055  Field< double >::set( gslComptId, "inject", 0 );
1056  Field< double >::set( gslComptId, "initVm", EREST );
1057 
1058  Field< double >::set( exptlComptId, "Cm", 0.007854e-6 );
1059  Field< double >::set( exptlComptId, "Ra", 7639.44e3 ); // does it matter?
1060  Field< double >::set( exptlComptId, "Rm", 424.4e3 );
1061  Field< double >::set( exptlComptId, "Em", EREST + 0.1 );
1062  Field< double >::set( exptlComptId, "inject", 0 );
1063  Field< double >::set( exptlComptId, "initVm", EREST );
1064 
1066  //Setup of rate tables.
1067  //Refer paper mentioned at the header of the unit test for more
1068  //details.
1070 
1071  //Number of states and open states.
1072  Field< unsigned int >::set( mChanGslId, "numStates", 4 );
1073  Field< unsigned int >::set( mChanExptlId, "numStates", 4 );
1074 
1075  Field< unsigned int >::set( mChanGslId, "numOpenStates", 1 );
1076  Field< unsigned int >::set( mChanExptlId, "numOpenStates", 1 );
1077 
1078  vector< string > stateLabels;
1079 
1080  //In the MarkovChannel class, the opening states are listed first.
1081  //This is in line with the convention followed in Chapter 20, Sakmann &
1082  //Neher.
1083  stateLabels.push_back( "O" ); //State 1.
1084  stateLabels.push_back( "B1" ); //State 2.
1085  stateLabels.push_back( "B2" ); //State 3.
1086  stateLabels.push_back( "C" ); //State 4.
1087 
1088  Field< vector< string > >::set( mChanGslId, "labels", stateLabels );
1089  Field< vector< string > >::set( mChanExptlId, "labels", stateLabels );
1090 
1091  //Setting up conductance value for single open state. Value chosen
1092  //is quite arbitrary.
1093  vector< double > gBar;
1094 
1095  gBar.push_back( 5.431553e-9 );
1096 
1097  Field< vector< double > >::set( mChanGslId, "gbar", gBar );
1098  Field< vector< double > >::set( mChanExptlId, "gbar", gBar );
1099 
1100  //Initial state of the system. This is really an arbitrary choice.
1101  vector< double > initState;
1102 
1103  initState.push_back( 0.00 );
1104  initState.push_back( 0.20 );
1105  initState.push_back( 0.80 );
1106  initState.push_back( 0.00 );
1107 
1108  Field< vector< double > >::set( mChanGslId, "initialState", initState );
1109  Field< vector< double > >::set( mChanExptlId, "initialState", initState );
1110 
1111  //This initializes the GSL solver object.
1112  SetGet1< vector< double > >::set( gslSolverId, "init", initState );
1113 
1114  //Initializing MarkovRateTable object.
1115  double v;
1116  double conc;
1117 
1118  SetGet1< unsigned int >::set( gslRateTableId, "init", 4 );
1119  SetGet1< unsigned int >::set( exptlRateTableId, "init", 4 );
1120 
1121  //Setting up lookup tables for the different rates.
1122  //Please note that the rates should be in sec^(-1).
1123 
1124  //Transition from "O" to "B1" i.e. r12 or a1.
1125  Field< double >::set( vecTableId, "xmin", -0.10 );
1126  Field< double >::set( vecTableId, "xmax", 0.10 );
1127  Field< unsigned int >::set( vecTableId, "xdivs", 200 );
1128 
1129  v = Field< double >::get( vecTableId, "xmin" );
1130  for ( unsigned int i = 0; i < 201; ++i )
1131  {
1132  table1d.push_back( 1e3 * exp( -16 * v - 2.91 ) );
1133  v += 0.001;
1134  }
1135 
1136  Field< vector< double > >::set( vecTableId, "table", table1d );
1137 
1138  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1139  gslRateTableId, "set1d", 1, 2, vecTableId, 0 );
1140  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1141  exptlRateTableId, "set1d", 1, 2, vecTableId, 0 );
1142 
1143  table1d.erase( table1d.begin(), table1d.end() );
1144 
1145  //Transition from "B1" back to O i.e. r21 or b1
1146  v = Field< double >::get( vecTableId, "xmin" );
1147  for ( unsigned int i = 0; i < 201; ++i )
1148  {
1149  table1d.push_back( 1e3 * exp( 9 * v + 1.22 ) );
1150  v += 0.001;
1151  }
1152 
1153  Field< vector< double > >::set( vecTableId, "table", table1d );
1154  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1155  gslRateTableId, "set1d", 2, 1, vecTableId, 0 );
1156  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1157  exptlRateTableId, "set1d", 2, 1, vecTableId, 0 );
1158 
1159  table1d.erase( table1d.begin(), table1d.end() );
1160 
1161  //Transition from "O" to "B2" i.e. r13 or a2
1162  //This is actually a 2D rate. But, there is no change in Mg2+ concentration
1163  //that occurs. Hence, I create a 2D lookup table anyway but I manually
1164  //set the concentration on the rate table object anyway.
1165 
1166  Field< double >::set( gslRateTableId, "ligandConc", 24e-6 );
1167  Field< double >::set( exptlRateTableId, "ligandConc", 24e-6 );
1168 
1169  Field< double >::set( int2dTableId, "xmin", -0.10 );
1170  Field< double >::set( int2dTableId, "xmax", 0.10 );
1171  Field< double >::set( int2dTableId, "ymin", 0 );
1172  Field< double >::set( int2dTableId, "ymax", 30e-6 );
1173  Field< unsigned int >::set( int2dTableId, "xdivs", 200 );
1174  Field< unsigned int >::set( int2dTableId, "ydivs", 30 );
1175 
1176  table2d.resize( 201 );
1177  v = Field< double >::get( int2dTableId, "xmin" );
1178  for ( unsigned int i = 0; i < 201; ++i )
1179  {
1180  conc = Field< double >::get( int2dTableId, "ymin" );
1181  for ( unsigned int j = 0; j < 31; ++j )
1182  {
1183  table2d[i].push_back( 1e3 * conc * exp( -45 * v - 6.97 ) );
1184  conc += 1e-6;
1185  }
1186  v += 1e-3;
1187  }
1188 
1189  Field< vector< vector< double > > >::set( int2dTableId, "tableVector2D",
1190  table2d );
1191 
1192  SetGet3< unsigned int, unsigned int, Id >::set( gslRateTableId,
1193  "set2d", 1, 3, int2dTableId );
1194  SetGet3< unsigned int, unsigned int, Id >::set( exptlRateTableId,
1195  "set2d", 1, 3, int2dTableId );
1196 
1197  //There is only one 2D rate, so no point manually erasing the elements.
1198 
1199  //Transition from "B2" to "O" i.e. r31 or b2
1200  v = -0.10;
1201  for ( unsigned int i = 0; i < 201; ++i )
1202  {
1203  table1d.push_back( 1e3 * exp( 17 * v + 0.96 ) );
1204  v += 0.001;
1205  }
1206 
1207  Field< vector< double > >::set( vecTableId, "table", table1d );
1208  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1209  gslRateTableId, "set1d", 3, 1, vecTableId, 0 );
1210  SetGet4< unsigned int, unsigned int, Id, unsigned int >::set(
1211  exptlRateTableId, "set1d", 3, 1, vecTableId, 0 );
1212 
1213  table1d.erase( table1d.begin(), table1d.end() );
1214 
1215  //Transition from "O" to "C" i.e. r14
1216  SetGet3< unsigned int, unsigned int, double >::set( gslRateTableId,
1217  "setconst", 1, 4, 1e3 * exp( -2.847 ) );
1218  SetGet3< unsigned int, unsigned int, double >::set( exptlRateTableId,
1219  "setconst", 1, 4, 1e3 * exp( -2.847 ) );
1220 
1221  //Transition from "B1" to "C" i.e. r24
1222  SetGet3< unsigned int, unsigned int, double >::set( gslRateTableId,
1223  "setconst", 2, 4, 1e3 * exp( -0.693 ) );
1224  SetGet3< unsigned int, unsigned int, double >::set( exptlRateTableId,
1225  "setconst", 2, 4, 1e3 * exp( -0.693 ) );
1226 
1227  //Transition from "B2" to "C" i.e. r34
1228  SetGet3< unsigned int, unsigned int, double >::set( gslRateTableId,
1229  "setconst", 3, 4, 1e3* exp( -3.101 ) );
1230  SetGet3< unsigned int, unsigned int, double >::set( exptlRateTableId,
1231  "setconst", 3, 4, 1e3* exp( -3.101 ) );
1232 
1233  //Once the rate tables have been set up, we can initialize the
1234  //tables in the MarkovSolver class.
1235  SetGet2< Id, double >::set( exptlSolverId, "init", exptlRateTableId, 1.0e-3 );
1236  SetGet1< double >::set( exptlSolverId, "ligandConc", 24e-6 );
1237  Field< vector< double > >::set( exptlSolverId, "initialState",
1238  initState );
1239 
1240  Field< double >::set( gslSolverId, "relativeAccuracy", 1e-24 );
1241  Field< double >::set( gslSolverId, "absoluteAccuracy", 1e-24 );
1242  Field< double >::set( gslSolverId, "internalDt", 1e-24 );
1243 
1244  shell->doSetClock( 0, 1.0e-3 );
1245  shell->doSetClock( 1, 1.0e-3 );
1246  shell->doSetClock( 2, 1.0e-3 );
1247  shell->doSetClock( 3, 1.0e-3 );
1248  shell->doSetClock( 4, 1.0e-3 );
1249  shell->doSetClock( 5, 1.0e-3 );
1250 
1251  shell->doUseClock( "/n/gslCompt,/n/exptlCompt", "init", 0 );
1252  shell->doUseClock( "/n/gslCompt,/n/exptlCompt", "process", 1 );
1253  shell->doUseClock( "/n/gslCompt/gslRateTable,/n/exptlCompt/exptlRateTable",
1254  "process", 2 );
1255  shell->doUseClock( "/n/gslCompt/gslSolver,/n/exptlCompt/exptlSolver",
1256  "process", 3 );
1257  shell->doUseClock( "/n/gslCompt/mChanGsl,/n/gslTable","process", 4 );
1258  shell->doUseClock( "/n/exptlCompt/mChanExptl,/n/exptlTable", "process", 5 );
1259 
1260  shell->doReinit();
1261  // shell->doReinit();// why twice? - subha
1262  shell->doStart( 1.0 );
1263 
1264  vector< double > gslVec = Field< vector< double > >::get( gslTableId, "vector" );
1265  vector< double > exptlVec = Field< vector< double > >::get( exptlTableId, "vector");
1266 
1267  assert( gslVec.size() == 1001 );
1268  assert( exptlVec.size() == 1001 );
1269 
1270  for ( unsigned int i = 0; i < 1001; ++i )
1271  assert( doubleApprox( gslVec[i], exptlVec[i] ) );
1272 
1273  shell->doDelete( nid );
1274  cout << "." << flush;
1275 }
1276 
1277 
1279 // Unit tests for SynChan
1281 
1282 // #include "SpikeGen.h"
1283 
1291 void testSynChan()
1292 {
1293  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
1294 
1295  Id nid = shell->doCreate( "Neutral", Id(), "n", 1 );
1296 
1297  Id synChanId = shell->doCreate( "SynChan", nid, "synChan", 1 );
1298  Id synHandlerId = shell->doCreate( "SimpleSynHandler", synChanId, "syns", 1 );
1299  Id synId( synHandlerId.value() + 1 );
1300  Id sgId1 = shell->doCreate( "SpikeGen", nid, "sg1", 1 );
1301  Id sgId2 = shell->doCreate( "SpikeGen", nid, "sg2", 1 );
1302  ProcInfo p;
1303  p.dt = 1.0e-4;
1304  p.currTime = 0;
1305  bool ret;
1306  assert( synId.element()->getName() == "synapse" );
1307  ret = Field< double >::set( synChanId, "tau1", 1.0e-3 );
1308  assert( ret );
1309  ret = Field< double >::set( synChanId, "tau2", 1.0e-3 );
1310  assert( ret );
1311  ret = Field< double >::set( synChanId, "Gbar", 1.0 );
1312  assert( ret );
1313 
1314  // This is a hack, should really inspect msgs to automatically figure
1315  // out how many synapses are needed.
1316  ret = Field< unsigned int >::set( synHandlerId, "numSynapse", 2 );
1317  assert( ret );
1318 
1319  Element* syne = synId.element();
1320  assert( syne->totNumLocalField() == 2 );
1321 
1322  ObjId mid = shell->doAddMsg( "single",
1323  ObjId( sgId1, 0 ), "spikeOut", ObjId( synId, 0, 0 ), "addSpike" );
1324  assert( mid != Id() );
1325  mid = shell->doAddMsg( "single",
1326  ObjId( sgId2, 0 ), "spikeOut", ObjId( synId, 0, 1 ), "addSpike" );
1327  assert( mid != Id() );
1328  mid = shell->doAddMsg( "single",
1329  synHandlerId, "activationOut", synChanId, "activation" );
1330  assert( mid != Id() );
1331 
1332  ret = Field< double >::set( sgId1, "threshold", 0.0 );
1333  ret = Field< double >::set( sgId1, "refractT", 1.0 );
1334  ret = Field< bool >::set( sgId1, "edgeTriggered", 0 );
1335  ret = Field< double >::set( sgId2, "threshold", 0.0 );
1336  ret = Field< double >::set( sgId2, "refractT", 1.0 );
1337  ret = Field< bool >::set( sgId2, "edgeTriggered", 0 );
1338 
1339  ret = Field< double >::set( ObjId( synId, 0, 0 ), "weight", 1.0 );
1340  assert( ret);
1341  ret = Field< double >::set( ObjId( synId, 0, 0 ), "delay", 0.001 );
1342  assert( ret);
1343  ret = Field< double >::set( ObjId( synId, 0, 1 ), "weight", 1.0 );
1344  assert( ret);
1345  ret = Field< double >::set( ObjId( synId, 0, 1 ), "delay", 0.01 );
1346  assert( ret);
1347 
1348  double dret;
1349  dret = Field< double >::get( ObjId( synId, 0, 0 ), "weight" );
1350  ASSERT_DOUBLE_EQ("", dret, 1.0 );
1351  dret = Field< double >::get( ObjId( synId, 0, 0 ), "delay" );
1352  ASSERT_DOUBLE_EQ("", dret, 0.001 );
1353  dret = Field< double >::get( ObjId( synId, 0, 1 ), "weight" );
1354  ASSERT_DOUBLE_EQ("", dret, 1.0 );
1355  dret = Field< double >::get( ObjId( synId, 0, 1 ), "delay" );
1356  ASSERT_DOUBLE_EQ("", dret, 0.01 );
1357 
1358  dret = SetGet1< double >::set( sgId1, "Vm", 2.0 );
1359  dret = SetGet1< double >::set( sgId2, "Vm", 2.0 );
1360  dret = Field< double >::get( synChanId, "Gk" );
1361  ASSERT_DOUBLE_EQ("", dret, 0.0 );
1362 
1364 
1365  shell->doSetClock( 0, 1e-4 );
1366  shell->doSetClock( 1, 1e-4 );
1367  shell->doSetClock( 2, 1e-4 );
1368  // shell->doUseClock( "/n/##", "process", 0 );
1369  // shell->doUseClock( "/n/synChan/syns,/n/sg1,/n/sg2", "process", 0 );
1370  //It turns out that the order of setting of the spikes (sg1, sg2)
1371  //does not affect the outcome. The one thing that is critical is that
1372  //the 'process' call for the 'syns' should be before that of the
1373  //synChan. This is because the 'activation' message from the syns to
1374  //the synChan should proceed within a given timestep otherwise the
1375  //apparent arrival time of the event is delayed.
1376  shell->doUseClock( "/n/sg1,/n/sg2", "process", 0 );
1377  shell->doUseClock( "/n/synChan/syns", "process", 1 );
1378  shell->doUseClock( "/n/synChan", "process", 2 );
1379  // shell->doStart( 0.001 );
1380  shell->doReinit();
1381  shell->doReinit();
1382 
1383  shell->doStart( 0.001 );
1384  dret = Field< double >::get( synChanId, "Gk" );
1385  assert( doubleApprox( dret, 0.0 ) );
1386 
1387  shell->doStart( 0.0005 );
1388  dret = Field< double >::get( synChanId, "Gk" );
1389  assert( doubleApprox( dret, 0.825 ) );
1390 
1391  shell->doStart( 0.0005 );
1392  dret = Field< double >::get( synChanId, "Gk" );
1393  assert( doubleApprox( dret, 1.0 ) );
1394 
1395  shell->doStart( 0.001 );
1396  dret = Field< double >::get( synChanId, "Gk" );
1397  assert( doubleApprox( dret, 0.736 ) );
1398 
1399  shell->doStart( 0.001 );
1400  dret = Field< double >::get( synChanId, "Gk" );
1401  assert( doubleApprox( dret, 0.406 ) );
1402 
1403  shell->doStart( 0.007 );
1404  dret = Field< double >::get( synChanId, "Gk" );
1405  // assert( doubleApprox( dret, 0.997 ) );
1406  assert( doubleApprox( dret, 1.002 ) );
1407 
1408  shell->doDelete( nid );
1409  cout << "." << flush;
1410 }
1411 
1412 #if 0
1413 
1414 void testNMDAChan()
1415 {
1416  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
1417 
1418  vector< int > dims( 1, 1 );
1419  Id nid = shell->doCreate( "Neutral", Id(), "n", dims );
1420 
1421  Id synChanId = shell->doCreate( "NMDAChan", nid, "nmdaChan", dims );
1422  Id synId( synChanId.value() + 1 );
1423  Id sgId1 = shell->doCreate( "SpikeGen", nid, "sg1", dims );
1424  ProcInfo p;
1425  p.dt = 1.0e-4;
1426  p.currTime = 0;
1427  bool ret;
1428  assert( synId()->getName() == "synapse" );
1429  ret = Field< double >::set( synChanId, "tau1", 130.5e-3 );
1430  assert( ret );
1431  ret = Field< double >::set( synChanId, "tau2", 5.0e-3 );
1432  assert( ret );
1433  ret = Field< double >::set( synChanId, "Gbar", 1.0 );
1434  assert( ret );
1435 
1436  // This is a hack, should really inspect msgs to automatically figure
1437  // out how many synapses are needed.
1438  ret = Field< unsigned int >::set( synChanId, "num_synapse", 1 );
1439  assert( ret );
1440 
1441  Element* syne = synId();
1442  assert( syne->dataHandler()->localEntries() == 1 );
1443  dynamic_cast< FieldDataHandlerBase* >( syne->dataHandler() )->setNumField( synChanId.eref().data(), 1 );
1444 
1445  assert( syne->dataHandler()->totalEntries() == 1 );
1446  assert( syne->dataHandler()->numDimensions() == 1 );
1447  assert( syne->dataHandler()->sizeOfDim( 0 ) == 1 );
1448 
1449  MsgId mid = shell->doAddMsg( "single",
1450  ObjId( sgId1, DataId( 0 ) ), "spikeOut",
1451  ObjId( synId, DataId( 0 ) ), "addSpike" );
1452  assert( mid != Id() );
1453 
1454  ret = Field< double >::set( sgId1, "threshold", 0.0 );
1455  ret = Field< double >::set( sgId1, "refractT", 1.0 );
1456  ret = Field< bool >::set( sgId1, "edgeTriggered", 0 );
1457 
1458  ret = Field< double >::set( ObjId( synId, DataId( 0 ) ),
1459  "weight", 1.0 );
1460  assert( ret);
1461  ret = Field< double >::set( ObjId( synId, DataId( 0 ) ),
1462  "delay", 0.001 );
1463  assert( ret);
1464 
1465  double dret;
1466  dret = Field< double >::get( ObjId( synId, DataId( 0 ) ), "weight" );
1467  ASSERT_DOUBLE_EQ("", dret, 1.0 );
1468  dret = Field< double >::get( ObjId( synId, DataId( 0 ) ), "delay" );
1469  ASSERT_DOUBLE_EQ("", dret, 0.001 );
1470 
1471  dret = SetGet1< double >::set( sgId1, "Vm", 2.0 );
1472  dret = Field< double >::get( synChanId, "Gk" );
1473  ASSERT_DOUBLE_EQ("", dret, 0.0 );
1474 
1476 
1477  shell->doSetClock( 0, 1e-4 );
1478  // shell->doUseClock( "/n/##", "process", 0 );
1479  shell->doUseClock( "/n/synChan,/n/sg1", "process", 0 );
1480  // shell->doStart( 0.001 );
1481  shell->doReinit();
1482  shell->doReinit();
1483 
1484  shell->doStart( 0.001 );
1485  dret = Field< double >::get( synChanId, "Gk" );
1486  assert( doubleApprox( dret, 0.0 ) );
1487 
1488  shell->doStart( 0.0005 );
1489  dret = Field< double >::get( synChanId, "Gk" );
1490  cout << "Gk:" << dret << endl;
1491  assert( doubleApprox( dret, 1.0614275017053588e-07 ) );
1492 
1493  // shell->doStart( 0.0005 );
1494  // dret = Field< double >::get( synChanId, "Gk" );
1495  // cout << "Gk:" << dret << endl;
1496  // assert( doubleApprox( dret, 1.0 ) );
1497 
1498  // shell->doStart( 0.001 );
1499  // dret = Field< double >::get( synChanId, "Gk" );
1500  // cout << "Gk:" << dret << endl;
1501  // assert( doubleApprox( dret, 0.736 ) );
1502 
1503  // shell->doStart( 0.001 );
1504  // dret = Field< double >::get( synChanId, "Gk" );
1505  // cout << "Gk:" << dret << endl;
1506  // assert( doubleApprox( dret, 0.406 ) );
1507 
1508  // shell->doStart( 0.007 );
1509  // dret = Field< double >::get( synChanId, "Gk" );
1510  // cout << "Gk:" << dret << endl;
1511  // assert( doubleApprox( dret, 0.997 ) );
1512 
1513  shell->doDelete( nid );
1514  cout << "." << flush;
1515 
1516 }
1517 #endif
1518 
1519 static Id addCompartment( const string& name,
1520  Id neuron, Id parent,
1521  double dx, double dy, double dz, double dia )
1522 {
1523  static Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1524  double x0 = 0.0;
1525  double y0 = 0.0;
1526  double z0 = 0.0;
1527  Id compt = shell->doCreate( "Compartment", neuron, name, 1 );
1528  if ( parent != Id() )
1529  {
1530  ObjId mid = shell->doAddMsg( "single",
1531  parent, "axial", compt, "raxial" );
1532  assert( mid != Id() );
1533  x0 = Field< double >::get( parent, "x" );
1534  y0 = Field< double >::get( parent, "y" );
1535  z0 = Field< double >::get( parent, "z" );
1536  }
1537  Field< double >::set( compt, "x0", x0 );
1538  Field< double >::set( compt, "y0", y0 );
1539  Field< double >::set( compt, "z0", z0 );
1540  Field< double >::set( compt, "x", x0 + dx );
1541  Field< double >::set( compt, "y", y0 + dy );
1542  Field< double >::set( compt, "z", z0 + dz );
1543  double length = sqrt( dx*dx + dy*dy + dz*dz );
1544  Field< double >::set( compt, "length", length );
1545  Field< double >::set( compt, "diameter", dia );
1546  Field< double >::set( compt, "Rm", 1.0 / ( PI * length * dia ) );
1547  Field< double >::set( compt, "Cm", 0.01 * ( PI * length * dia ) );
1548  Field< double >::set( compt, "Ra", 1 * length / ( PI*0.25*dia*dia ) );
1549  return compt;
1550 }
1551 
1552 #include "../utility/Vec.h"
1553 #include "SwcSegment.h"
1554 #include "Spine.h"
1555 #include "Neuron.h"
1556 #include "../shell/Wildcard.h"
1557 static void testNeuronBuildTree()
1558 {
1559  Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
1560 
1561  Id nid = shell->doCreate( "Neuron", Id(), "n", 1 );
1562  double somaDia = 5e-6;
1563  double dendDia = 2e-6;
1564  double branchDia = 1e-6;
1565  static double len[] = { 10e-6, 100e-6, 200e-6, 500e-6 };
1566  static double dia[] = { somaDia, dendDia, branchDia, branchDia };
1567  Id soma = addCompartment ( "soma", nid, Id(), 10e-6, 0, 0, somaDia );
1568  Id dend1 = addCompartment ( "dend1", nid, soma, 100e-6, 0, 0, dendDia);
1569  Id branch1 = addCompartment ( "branch1", nid, dend1, 0, 200e-6, 0, branchDia );
1570  Id branch2 = addCompartment ( "branch2", nid, dend1, 0, -500e-6, 0, branchDia );
1571  static double x[] = { 10e-6, 110e-6, 110e-6, 110e-6 };
1572  static double y[] = {0, 0, 200e-6, -500e-6};
1573  static double z[] = {0, 0, 0, 0};
1574 
1575  SetGet0::set( nid, "buildSegmentTree" );
1576 
1577  vector< double > e = Field< vector< double > >::get(
1578  nid, "electrotonicDistanceFromSoma" );
1579  vector< double > g = Field< vector< double > >::get(
1580  nid, "geometricalDistanceFromSoma" );
1581  vector< double > p = Field< vector< double > >::get(
1582  nid, "pathDistanceFromSoma" );
1583  assert( e.size() == 4 );
1584  ASSERT_DOUBLE_EQ("", e[0], 0.0 );
1585  double dL = 100e-6 / sqrt( dendDia /4.0 );
1586  ASSERT_DOUBLE_EQ("", e[1], dL );
1587  double bL1 = dL + 200e-6 / sqrt( branchDia /4.0 );
1588  ASSERT_DOUBLE_EQ("", e[2], bL1 );
1589  double bL2 = dL + 500e-6 / sqrt( branchDia/4.0 );
1590  ASSERT_DOUBLE_EQ("", e[3], bL2 );
1591 
1592  ASSERT_DOUBLE_EQ("", p[0], 0.0 );
1593  ASSERT_DOUBLE_EQ("", p[1], 100.0e-6 );
1594  ASSERT_DOUBLE_EQ("", p[2], 300.0e-6 ); // 100 + 200 microns
1595  ASSERT_DOUBLE_EQ("", p[3], 600.0e-6 ); // 100 + 500 microns
1596 
1597  ASSERT_DOUBLE_EQ("", g[0], 0.0 );
1598  ASSERT_DOUBLE_EQ("", g[1], 100.0e-6 );
1599  ASSERT_DOUBLE_EQ("", g[2], sqrt(5.0) * 100.0e-6 ); // 100 + 200 microns
1600  ASSERT_DOUBLE_EQ("", g[3], sqrt(26.0) * 100.0e-6 ); // 100 + 500 microns
1601 
1603  // Here we test Neuron::evalExprForElist, which uses the muParser
1604  // Note that the wildcard list starts with the spine which is not
1605  // a compartment. So the indexing of the arrays e, p and g needs care.
1606  unsigned int nuParserNumVal = 13;
1607  vector< ObjId > elist;
1608  wildcardFind( "/n/#[ISA=Compartment]", elist );
1609  Neuron* n = reinterpret_cast< Neuron* >( nid.eref().data() );
1610  vector< double > val;
1611  n->evalExprForElist( elist, "p + g + L + len + dia + H(1-L)", val );
1612  assert( val.size() == nuParserNumVal * elist.size() );
1613  double maxP = 0.0;
1614  double maxG = 0.0;
1615  double maxL = 0.0;
1616  for ( unsigned int i = 0; i < elist.size(); ++i )
1617  {
1618  if ( maxP < p[i] ) maxP = p[i];
1619  if ( maxG < g[i] ) maxG = g[i];
1620  if ( maxL < e[i] ) maxL = e[i];
1621  }
1622  unsigned int j = 0;
1623  for ( unsigned int i = 0; i < elist.size(); ++i )
1624  {
1625  if ( !elist[i].element()->cinfo()->isA( "CompartmentBase" ) )
1626  continue;
1627  assert( val[i * nuParserNumVal] ==
1628  p[j] + g[j] + e[j] + len[j] + dia[j] + ( 1.0 - e[j] > 0 ) );
1629  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 1], p[j] );
1630  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 2], g[j] );
1631  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 3], e[j] );
1632  assert( doubleEq( val[i * nuParserNumVal + 4], len[j] ));
1633  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 5], dia[j] );
1634  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 6], maxP );
1635  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 7], maxG );
1636  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 8], maxL );
1637  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 9], x[j] );
1638  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 10], y[j] );
1639  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 11], z[j] );
1640  ASSERT_DOUBLE_EQ("", val[i * nuParserNumVal + 12], 0.0 );
1641  j++;
1642  }
1644  // Here we test Neuron::makeSpacingDistrib, which uses the muParser
1645  // n->evalExprForElist( elist, "H(p-50e-6)*5e-6", val );
1646  n->evalExprForElist( elist, "H(p-100e-6)*5e-6", val );
1647  vector< unsigned int > seglistIndex;
1648  vector< unsigned int > elistIndex;
1649  vector< double > pos;
1650  vector< string > line; // empty, just use the default spacingDistrib=0
1651  n->makeSpacingDistrib( elist, val, seglistIndex, elistIndex, pos, line );
1652  // Can't do this now, it is not determinisitic.
1653  // assert( pos.size() == ((800 - 100)/5) );
1654  // ASSERT_DOUBLE_EQ("", pos[0], 2.5e-6 );
1655  // ASSERT_DOUBLE_EQ("", pos.back(), 500e-6 - 7.5e-6 );
1656  assert( seglistIndex[0] == 2 );
1657  assert( seglistIndex.back() == 3 );
1658  assert( elistIndex[0] == 2 );
1659  assert( elistIndex.back() == 3 );
1660 
1661  shell->doDelete( nid );
1662 }
1663 
1664 
1665 // This tests stuff without using the messaging.
1666 void testBiophysics()
1667 {
1668  testCompartment();
1669  testVectorTable();
1670  testNeuronBuildTree();
1671 #if 0
1672  testMarkovSolverBase();
1673  testMarkovSolver();
1674  testHHGateCreation();
1675  testHHGateLookup();
1676  testHHGateSetup();
1677  testSpikeGen();
1678  testCaConc();
1679  testNernst();
1680 #endif
1681 }
1682 
1683 // This is applicable to tests that use the messaging and scheduling.
1684 void testBiophysicsProcess()
1685 {
1686  // testSynChan();
1687  testIntFireNetwork();
1688  testCompartmentProcess();
1689  // testMarkovGslSolver();
1690  // testMarkovChannel();
1691 #if 0
1692  testHHChannel();
1693 #endif
1694 }
1695 
1696 #else // ifdef DO_UNIT_TESTS
1697 void testBiophysics()
1698 {
1699  ;
1700 }
1701 void testBiophysicsProcess()
1702 {
1703  ;
1704 }
1705 void testIntFireNetwork( unsigned int unsteps = 5 )
1706 {
1707  ;
1708 }
1709 #endif
wildcardFind
int wildcardFind(const string &path, vector< ObjId > &ret)
Definition: Wildcard.cpp:169
testIntFireNetwork
void testIntFireNetwork(unsigned int unsteps=5)
Definition: testBiophysics.cpp:1705
moose::mtseed
void mtseed(unsigned int x)
Set the global seed or all rngs.
Definition: global.cpp:89
Shell::doStart
void doStart(double runtime, bool notify=false)
Definition: Shell.cpp:332
HHChannel::zGate_
HHGate * zGate_
HHGate data structure for the yGate.
Definition: HHChannel.h:206
Eref::data
char * data() const
Definition: Eref.cpp:41
HHGate
Definition: HHGate.h:31
Neuron
Definition: Neuron.h:18
ObjId::data
char * data() const
Definition: ObjId.cpp:113
Id::element
Element * element() const
Synonym for Id::operator()()
Definition: Id.cpp:113
Neuron::makeSpacingDistrib
void makeSpacingDistrib(const vector< ObjId > &elist, const vector< double > &val, vector< unsigned int > &seglistIndex, vector< unsigned int > &elistIndex, vector< double > &pos, const vector< string > &line) const
Definition: Neuron.cpp:1841
Shell::doSetClock
void doSetClock(unsigned int tickNum, double dt)
Definition: Shell.cpp:377
ProcInfo::currTime
double currTime
Definition: ProcInfo.h:19
ObjId::bad
bool bad() const
Definition: ObjId.cpp:18
Id::value
unsigned int value() const
Definition: Id.cpp:197
Field
Definition: SetGet.h:236
ObjId
Definition: ObjId.h:20
Neuron::evalExprForElist
void evalExprForElist(const vector< ObjId > &elist, const string &expn, vector< double > &val) const
Definition: Neuron.cpp:1478
Id::eref
Eref eref() const
Definition: Id.cpp:125
testCompartment
void testCompartment()
Field::set
static bool set(const ObjId &dest, const string &field, A arg)
Definition: SetGet.h:245
testCompartmentProcess
void testCompartmentProcess()
Shell::doCreate
Id doCreate(string type, ObjId parent, string name, unsigned int numData, NodePolicy nodePolicy=MooseBlockBalance, unsigned int preferredNode=1)
Definition: Shell.cpp:181
testVectorTable
void testVectorTable()
testMarkovRateTable
void testMarkovRateTable()
doubleApprox
bool doubleApprox(double x, double y)
Definition: doubleEq.cpp:24
doubleEq
bool doubleEq(double x, double y)
Definition: doubleEq.cpp:16
ProcInfo::dt
double dt
Definition: ProcInfo.h:18
SetGet4::set
static bool set(const ObjId &dest, const string &field, A1 arg1, A2 arg2, A3 arg3, A4 arg4)
Definition: SetGet.h:687
HHChannel::yGate_
HHGate * yGate_
HHGate data structure for the yGate.
Definition: HHChannel.h:203
Eref
Definition: Eref.h:26
EXPECT_EQ
#define EXPECT_EQ(a, b, token)
Definition: testing_macros.hpp:61
HHChannel::xGate_
HHGate * xGate_
Definition: HHChannel.h:200
Field::setVec
static bool setVec(ObjId destId, const string &field, const vector< A > &arg)
Definition: SetGet.h:252
SetGet0::set
static bool set(const ObjId &dest, const string &field)
Definition: SetGet.h:103
Element::totNumLocalField
virtual unsigned int totNumLocalField() const =0
testBiophysicsProcess
void testBiophysicsProcess()
Definition: testBiophysics.cpp:1701
Shell::doReinit
void doReinit()
Definition: Shell.cpp:362
name
static char name[]
Definition: mfield.cpp:401
Shell::doAddMsg
ObjId doAddMsg(const string &msgType, ObjId src, const string &srcField, ObjId dest, const string &destField)
Definition: Shell.cpp:269
SetGet1::set
static bool set(const ObjId &dest, const string &field, A arg)
Definition: SetGet.h:153
SetGet3::set
static bool set(const ObjId &dest, const string &field, A1 arg1, A2 arg2, A3 arg3)
Definition: SetGet.h:626
PI
const double PI
Definition: consts.cpp:12
Shell::doDelete
bool doDelete(ObjId oid)
Definition: Shell.cpp:259
moose::mtrand
double mtrand(void)
Generate a random double between 0 and 1.
Definition: global.cpp:97
Id
Definition: Id.h:17
testBiophysics
void testBiophysics()
Definition: testBiophysics.cpp:1697
Shell::doUseClock
void doUseClock(string path, string field, unsigned int tick)
Definition: Shell.cpp:382
Field::get
static A get(const ObjId &dest, const string &field)
Definition: SetGet.h:284
Shell::numNodes
static unsigned int numNodes()
Definition: ShellThreads.cpp:76
Field::getVec
static void getVec(ObjId dest, const string &field, vector< A > &vec)
Definition: SetGet.h:317
ASSERT_DOUBLE_EQ
#define ASSERT_DOUBLE_EQ(token, a, b)
Definition: testing_macros.hpp:149
Element::getName
const string & getName() const
Definition: Element.cpp:56
DataId
unsigned int DataId
Definition: header.h:47
Shell
Definition: Shell.h:43
SetGet2::set
static bool set(const ObjId &dest, const string &field, A1 arg1, A2 arg2)
Definition: SetGet.h:365