/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
**           Copyright (C) 2003-2014 Upinder S. Bhalla. and NCBS
** It is made available under the terms of the
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/
/*
#include <vector>
#include <algorithm>
#include <cassert>
#include <functional>
#include <iostream>
#include <iomanip>
//#include "/usr/include/gsl/gsl_linalg.h"
using namespace std;
*/
#ifdef USE_GSL
#include <gsl/gsl_linalg.h>
#endif
#include "header.h"
#include "../basecode/SparseMatrix.h"
#include "FastMatrixElim.h"
#include "../shell/Shell.h"



double checkAns( 
	const double* m, unsigned int numCompts, 
	const double* ans, const double* rhs )
{
	vector< double > check( numCompts, 0.0 );
	for ( unsigned int i = 0; i < numCompts; ++i ) {
		for ( unsigned int j = 0; j < numCompts; ++j )
			check[i] += m[i*numCompts + j] * ans[j];
	}
	double ret = 0.0;
	for ( unsigned int i = 0; i < numCompts; ++i )
		ret += (check[i] - rhs[i]) * (check[i] - rhs[i] );
	return ret;
}



void testFastMatrixElim()
{


/*
2    11
 1  4
   3    10
    9  5
     6
     7
     8

        1 2 3 4 5 6 7 8 9 10 11
1       x x x x
2       x x          
3       x   x x         x
4       x   x x              x
5               x x     x x
6               x x x   x
7                 x x x
8                   x x  
9           x   x x     x  
10              x         x   
11            x              x
	static double test[] = {
		1,  2,  3,  4,  0,  0,  0,  0,  0,  0,  0,
		5,  6,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		7,  0,  8,  9,  0,  0,  0,  0, 10,  0,  0,
		11, 0, 12, 13,  0,  0,  0,  0,  0,  0, 14,
		0,  0,  0,  0, 15, 16,  0,  0, 17, 18,  0,
		0,  0,  0,  0, 19, 20, 21,  0, 22,  0,  0,
		0,  0,  0,  0,  0, 23, 24, 25,  0,  0,  0,
		0,  0,  0,  0,  0,  0, 26, 27,  0,  0,  0,
		0,  0, 28,  0, 29, 30,  0,  0, 31,  0,  0,
		0,  0,  0,  0, 32,  0,  0,  0,  0, 33,  0,
		0,  0,  0, 34,  0,  0,  0,  0,  0,  0, 35,
	};
	const unsigned int numCompts = 11;
//	static unsigned int parents[] = { 3,1,9,3,6,7,8,-1,6,5,4 };
	static unsigned int parents[] = { 2,0,8,2,5,6,7,-1,5,4,3 };
*/

/*
1   3
 2 4
   7   5
    8 6
     9
     10
     11

        1 2 3 4 5 6 7 8 9 10 11
1       x x
2       x x         x
3           x x
4           x x     x
5               x x
6               x x     x
7         x   x     x x
8                   x x x
9                 x   x x x
10                      x x  x
11                        x  x
	static double test[] = {
		1,  2,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		3,  4,  0,  0,  0,  0,  5,  0,  0,  0,  0,
		0,  0,  6,  7,  0,  0,  0,  0,  0,  0,  0,
		0,  0,  8,  9,  0,  0, 10,  0,  0,  0,  0,
		0,  0,  0,  0, 11, 12,  0,  0,  0,  0,  0,
		0,  0,  0,  0, 13, 14,  0,  0, 15,  0,  0,
		0, 16,  0, 17,  0,  0, 18, 19,  0,  0,  0,
		0,  0,  0,  0,  0,  0, 20, 21, 22,  0,  0,
		0,  0,  0,  0,  0, 23,  0, 24, 25, 26,  0,
		0,  0,  0,  0,  0,  0,  0,  0, 27, 28, 29,
		0,  0,  0,  0,  0,  0,  0,  0,  0, 30, 31,
	};
	const unsigned int numCompts = 11;
	static unsigned int parents[] = { 1,6,3,6,5,8,7,8,9,10,-1};
*/

/*
Linear cable, 12 segments.
*/

	static double test[] = {
		1,  2,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		3,  4,  5,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		0,  6,  7,  8,  0,  0,  0,  0,  0,  0,  0,  0,
		0,  0,  9, 10, 11,  0,  0,  0,  0,  0,  0,  0,
		0,  0,  0, 12, 13, 14,  0,  0,  0,  0,  0,  0,
		0,  0,  0,  0, 15, 16, 17,  0,  0,  0,  0,  0,
		0,  0,  0,  0,  0, 18, 19, 20,  0,  0,  0,  0,
		0,  0,  0,  0,  0,  0, 21, 22, 23,  0,  0,  0,
		0,  0,  0,  0,  0,  0,  0, 24, 25, 26,  0,  0,
		0,  0,  0,  0,  0,  0,  0,  0, 27, 28, 29,  0,
		0,  0,  0,  0,  0,  0,  0,  0,  0, 30, 31, 32,
		0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 33, 34,
	};
	const unsigned int numCompts = 12;
	static unsigned int parents[] = { 1,2,3,4,5,6,7,8,9,10,11, unsigned(-1)};

		/*
	static double test[] = {
		1,  2,
		3,  4
	};
	const unsigned int numCompts = 2;
	static double test[] = {
		1, 2, 0, 0,
		3, 4, 5, 0,
		0, 6, 7, 8,
		0, 0, 9, 10
	};
	const unsigned int numCompts = 4;
	static double test[] = {
		1,  2,  0,  0,  0,  0,
		3,  4,  5,  0,  0,  0,
		0,  6,  7,  8,  0,  0,
		0,  0,  9, 10, 11,  0,
		0,  0,  0, 12, 13, 14,
		0,  0,  0,  0, 15, 16,
	};
	const unsigned int numCompts = 6;
	static double test[] = {
		1,  2,  0,  0,  0,  0,
		3,  4,  0,  0,  1,  0,
		0,  0,  7,  8,  0,  0,
		0,  0,  9, 10, 11,  0,
		0,  1,  0, 12, 13, 14,
		0,  0,  0,  0, 15, 16,
	};
	const unsigned int numCompts = 6;
	*/
	// testSorting(); // seems to work fine.
	FastMatrixElim fe;
	vector< Triplet< double > > fops;
	fe.makeTestMatrix( test, numCompts );
	// fe.print();
	vector< unsigned int > parentVoxel;
	vector< unsigned int > lookupOldRowsFromNew;
	parentVoxel.insert( parentVoxel.begin(), &parents[0], &parents[numCompts] );
	fe.hinesReorder( parentVoxel, lookupOldRowsFromNew );
	/*
	*/
	/*
	vector< unsigned int > shuf;
	for ( unsigned int i = 0; i < numCompts; ++i )
		shuf.push_back( i );
	shuf[0] = 1;
	shuf[1] = 0;
	fe.shuffleRows( shuf );
	*/
	// fe.print();
	FastMatrixElim foo = fe;

	vector< unsigned int > diag;
	vector< double > diagVal;
	fe.buildForwardElim( diag, fops );
	// fe.print();
	fe.buildBackwardSub( diag, fops, diagVal );
	vector< double > y( numCompts, 1.0 );
	vector< double > ones( numCompts, 1.0 );
	FastMatrixElim::advance( y, fops, diagVal );
	/*
	for ( unsigned int i = 0; i < numCompts; ++i )
		cout << "y" << i << "]=	" << y[i] << endl;
		*/

	// Here we verify the answer
	
	vector< double > alle;
	for( unsigned int i = 0; i < numCompts; ++i ) {
		for( unsigned int j = 0; j < numCompts; ++j ) {
			alle.push_back( foo.get( i, j ) );
		}
	}
	// cout << "myCode: " << checkAns( &alle[0], numCompts, &y[0], &ones[0] ) << endl;

	assert(	checkAns( &alle[0], numCompts, &y[0], &ones[0] ) < 1e-25 );

#if USE_GSL
	/////////////////////////////////////////////////////////////////////
	// Here we do the gsl test.
	vector< double > temp( &test[0], &test[numCompts*numCompts] );
	gsl_matrix_view m = gsl_matrix_view_array( &temp[0], numCompts, numCompts );

	vector< double > z( numCompts, 1.0 );
	gsl_vector_view b = gsl_vector_view_array( &z[0], numCompts );
	gsl_vector* x = gsl_vector_alloc( numCompts );
	int s;
	gsl_permutation* p = gsl_permutation_alloc( numCompts );
	gsl_linalg_LU_decomp( &m.matrix, p, &s );
	gsl_linalg_LU_solve( &m.matrix, p, &b.vector, x);
	vector< double > gslAns( numCompts );
	for ( unsigned int i = 0; i < numCompts; ++i ) {
		gslAns[i] = gsl_vector_get( x, i );
		// cout << "x[" << i << "]=	" << gslAns[i] << endl;
	}
	// cout << "GSL: " << checkAns( test, numCompts, &gslAns[0], &ones[0] ) << endl;
	assert( checkAns( test, numCompts, &gslAns[0], &ones[0] ) < 1e-25 );
	gsl_permutation_free( p );
	gsl_vector_free( x );
	cout << "." << flush;
#endif
}

void testSorting()
{
	static unsigned int k[] = {20,40,60,80,100,10,30,50,70,90};
	static double d[] = {1,2,3,4,5,6,7,8,9,10};
	vector< unsigned int > col;
	col.insert( col.begin(), k, k+10);
	vector< double > entry;
	entry.insert( entry.begin(), d, d+10);
	sortByColumn( col, entry );

	for ( unsigned int i = 0; i < col.size(); ++i )
		assert( col[i] == (i + 1) * 10 );

	assert( entry[0] == 6 );
	assert( entry[1] == 1 );
	assert( entry[2] == 7 );
	assert( entry[3] == 2 );
	assert( entry[4] == 8 );
	assert( entry[5] == 3 );
	assert( entry[6] == 9 );
	assert( entry[7] == 4 );
	assert( entry[8] == 10 );
	assert( entry[9] == 5 );

	/*
	cout << "testing sorting\n";
	for ( unsigned int i = 0; i < col.size(); ++i ) {
		cout << "d[" << i << "]=	" << k[i] << 
		   ", col[" << i << "]= " <<	col[i] << ", e=" << entry[i] << endl;
	}
	cout << endl;
	*/
	cout << "." << flush;
}

void testSetDiffusionAndTransport()
{
	static double test[] = {
		0,  2,  0,  0,  0,  0,
		1,  0,  2,  0,  0,  0,
		0,  1,  0,  2,  0,  0,
		0,  0,  1,  0,  2,  0,
		0,  0,  0,  1,  0,  2,
		0,  0,  0,  0,  1,  0,
	};
	const unsigned int numCompts = 6;
	FastMatrixElim fm;
	fm.makeTestMatrix( test, numCompts );
	vector< unsigned int > parentVoxel( numCompts );
	parentVoxel[0] = -1;
	parentVoxel[1] = 0;
	parentVoxel[2] = 1;
	parentVoxel[3] = 2;
	parentVoxel[4] = 3;
	parentVoxel[5] = 4;

	// cout << endl;
	// fm.print();
	// cout << endl;
	// fm.printInternal();
	fm.setDiffusionAndTransport( parentVoxel, 1, 10, 0.1 );
	// cout << endl;
	// fm.print();
	// cout << endl;
	// fm.printInternal();

	for( unsigned int i =0; i < numCompts; ++i ) {
		unsigned int start = 0;
		if ( i > 0 )
			start = i - 1;
		for( unsigned int j = start ; j < i+1 && j < numCompts ; ++j ) {
			if ( i == j + 1 )
				assert( doubleEq( fm.get( i, j ), 0.1 ) );
			else if ( i + 1 == j ) {
				assert( doubleEq( fm.get( i, j ), 2.2 ) );
			} else if ( i == j ) {
				if ( i == 0 )
					assert( doubleEq( fm.get( i, j ), 0.8 ) );
				else if ( i == numCompts - 1 )
					assert( doubleEq( fm.get( i, j ), -0.1 ) );
				else 
					assert( doubleEq( fm.get( i, j ), -0.3 ) );
			}
		}
	}
	cout << "." << flush;
}

void testCylDiffn()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 1e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
	Field< double >::set( pool, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", cyl );
	// Next: build by doing reinit
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	Field< string >::set( dsolve, "path", "/model/cyl/pool" );
	// Then find a way to test it.
	vector< double > poolVec;
	Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );
   	Field< double >::getVec( pool, "nInit", poolVec );
	assert( poolVec.size() == ndc );
	assert( doubleEq( poolVec[0], 1.0 ) );
	assert( doubleEq( poolVec[1], 0.0 ) );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );

	s->doReinit();
	s->doStart( runtime );

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
   	Field< double >::getVec( pool, "n", poolVec );
	assert( nvec.size() == poolVec.size() );
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		assert( doubleEq( nvec[i], poolVec[i] ) );
	/*
	cout << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << "	";
	cout << endl;
	*/

	double dx = diffLength;
	double err = 0.0;
	double analyticTot = 0.0;
	double myTot = 0.0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		double x = i * dx + dx * 0.5;
		// This part is the solution as a func of x,t.
		double y = dx *  // This part represents the init n of 1 in dx
			( 1.0 / sqrt( PI * diffConst * runtime ) ) * 
			exp( -x * x / ( 4 * diffConst * runtime ) ); 
		err += ( y - nvec[i] ) * ( y - nvec[i] );
		//cout << i << "	" << x << "	" << y << "	" << conc[i] << endl;
		analyticTot += y;
		myTot += nvec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );
	// cout << "analyticTot= " << analyticTot << ", myTot= " << myTot << endl;
	assert( err < 1.0e-5 );


	s->doDelete( model );
	cout << "." << flush;
}

void testTaperingCylDiffn()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 2e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
	double diffConst = 1.0e-12; 
	// Should set explicitly, currently during creation of DiffPoolVec
	//double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
	Field< double >::set( pool, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", cyl );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build by setting the path of the dsolve.
	Field< string >::set( dsolve, "path", "/model/cyl/pool" );
	// Then find a way to test it.
	assert( pool.element()->numData() == ndc );
	Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );

	s->doReinit();
	s->doStart( runtime );

	double myTot = 0.0;
	vector< double > poolVec;
   	Field< double >::getVec( pool, "n", poolVec );
	for ( unsigned int i = 0; i < poolVec.size(); ++i ) {
		myTot += poolVec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );

	s->doDelete( model );
	cout << "." << flush;
}

/**
 * Cell looks like:
 *                soma
 *                dend[0]
 *                dend[1]
 *             b1[0]    b2[0]
 *            b1[1]      b2[1]
 *         t1[0]   t2[0]
 *        t1[1]      t2[1]
 *
 *
 * The matrix should look like:
 *
 *		s	d0	d1	b10	b11	t10	t11	t20	t21	b20	b21
 * s	#	1	0	0	0	0	0	0	0	0	0
 * d0	1	#	1	0	0	0	0	0	0	0	0
 * d1	0	1	#	1	0	0	0	0	0	1	0
 * b10	0	0	1	#	1	0	0	0	0	c	0
 * b11	0	0	0	1	#	1	0	1	0	0	0
 * t10	0	0	0	0	1	#	1	c	0	0	0
 * t11	0	0	0	0	0	1	#	0	0	0	0
 * t20	0	0	0	0	1	c	0	#	1	0	0
 * t21	0	0	0	0	0	0	0	1	#	0	0
 * b20	0	0	1	c	0	0	0	0	0	#	1
 * b21	0	0	0	0	0	0	0	0	0	1	#
 *                
 */
void testSmallCellDiffn()
{
	Id makeCompt( Id parentCompt, Id parentObj,
		string name, double len, double dia, double theta );
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 20e-6;
	double dia = 10e-6;
	double diffLength = 10e-6;
	double dt = 1.0e-1;
	double runtime = 100.0;
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id soma = makeCompt( Id(), model, "soma", dia, dia, 90 );
	Id dend = makeCompt( soma, model, "dend", len, 3e-6, 0 );
	Id branch1 = makeCompt( dend, model, "branch1", len, 2e-6, 45.0 );
	Id branch2 = makeCompt( dend, model, "branch2", len, 2e-6, -45.0 );
	Id twig1 = makeCompt( branch1, model, "twig1", len, 1.5e-6, 90.0 );
	Id twig2 = makeCompt( branch1, model, "twig2", len, 1.5e-6, 0.0 );

	Id nm = s->doCreate( "NeuroMesh", model, "neuromesh", 1 );
	Field< double >::set( nm, "diffLength", diffLength );
	Field< string >::set( nm, "geometryPolicy", "cylinder" );
	Field< string >::set( nm, "subTreePath", "/model/#" );
	//SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );
	unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
	assert( ns == 6 );
	unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
	assert( ndc == 11  );
	Id pool1 = s->doCreate( "Pool", nm, "pool1", 1 );
	Field< double >::set( pool1, "diffConst", diffConst );
	Id pool2 = s->doCreate( "Pool", nm, "pool2", 1 );
	Field< double >::set( pool2, "diffConst", diffConst );
	Id pool3 = s->doCreate( "Pool", nm, "pool3", 1 );
	Field< double >::set( pool3, "diffConst", 0.0 );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", nm );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build diffusion by setting path
	Field< string >::set( dsolve, "path", "/model/neuromesh/pool#" );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );
	assert( pool1.element()->numData() == ndc );
	Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
	Field< double >::set( ObjId( pool2, ndc - 1 ), "nInit", 2.0 );
	Field< double >::set( ObjId( pool3, 0 ), "nInit", 3.0 );

	s->doReinit();
	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( doubleEq( nvec[0], 1.0 ) );
	assert( doubleEq( nvec[1], 0.0 ) );
	s->doStart( runtime );

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	vector< double > pool1Vec;
	Field< double >::getVec( pool1, "n", pool1Vec );
	assert( pool1Vec.size() == ndc );

	vector< double > pool2Vec;
	Field< double >::getVec( pool2, "n", pool2Vec );
	assert( pool2Vec.size() == ndc );

	vector< double > pool3Vec;
	Field< double >::getVec( pool3, "n", pool3Vec );
	assert( pool3Vec.size() == ndc );
	double myTot1 = 0;
	double myTot2 = 0;
	double myTot3 = 0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		assert( doubleEq( pool1Vec[i], nvec[i] ) );
		myTot1 += nvec[i];
		myTot2 += pool2Vec[i];
		myTot3 += pool3Vec[i];
	}
	assert( doubleEq( myTot1, 1.0 ) );
	assert( doubleEq( myTot2, 2.0 ) );
	assert( doubleEq( myTot3, 3.0 ) );
	assert( doubleEq( pool3Vec[0], 3.0 ) );

	/*
	cout << "Small cell: " << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << ", " << pool2Vec[i] << endl;
	cout << endl;
	*/


	s->doDelete( model );
	cout << "." << flush;
}

void testCellDiffn()
{
	Id makeCompt( Id parentCompt, Id parentObj,
		string name, double len, double dia, double theta );
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 40e-6;
	double dia = 10e-6;
	double diffLength = 1e-6;
	double dt = 1.0e-1;
	double runtime = 100.0;
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id soma = makeCompt( Id(), model, "soma", dia, dia, 90 );
	Id dend = makeCompt( soma, model, "dend", len, 3e-6, 0 );
	Id branch1 = makeCompt( dend, model, "branch1", len, 2e-6, 45.0 );
	Id branch2 = makeCompt( dend, model, "branch2", len, 2e-6, -45.0 );
	Id twig1 = makeCompt( branch1, model, "twig1", len, 1.5e-6, 90.0 );
	Id twig2 = makeCompt( branch1, model, "twig2", len, 1.5e-6, 0.0 );

	Id nm = s->doCreate( "NeuroMesh", model, "neuromesh", 1 );
	Field< double >::set( nm, "diffLength", diffLength );
	Field< string >::set( nm, "geometryPolicy", "cylinder" );
	Field< string >::set( nm, "subTreePath", "/model/#" );
	//SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );
	unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
	assert( ns == 6 );
	unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
	assert( ndc == 210  );
	Id pool1 = s->doCreate( "Pool", nm, "pool1", 1 );
	Field< double >::set( pool1, "diffConst", diffConst );
	Id pool2 = s->doCreate( "Pool", nm, "pool2", 1 );
	Field< double >::set( pool2, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", nm );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build by setting path
	Field< string >::set( dsolve, "path", "/model/neuromesh/pool#" );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );
	assert( pool1.element()->numData() == ndc );
	Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
	Field< double >::set( ObjId( pool2, ndc - 1 ), "nInit", 2.0 );

	s->doReinit();
	s->doStart( runtime );

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	vector< double > pool1Vec;
	Field< double >::getVec( pool1, "n", pool1Vec );
	assert( pool1Vec.size() == ndc );

	vector< double > pool2Vec;
	Field< double >::getVec( pool2, "n", pool2Vec );
	assert( pool2Vec.size() == ndc );
	double myTot1 = 0;
	double myTot2 = 0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		assert( doubleEq( pool1Vec[i], nvec[i] ) );
		myTot1 += nvec[i];
		myTot2 += pool2Vec[i];
	}
	assert( doubleEq( myTot1, 1.0 ) );
	assert( doubleEq( myTot2, 2.0 ) );

	/*
	cout << endl;
	cout << "Big cell: " << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << ", " << pool2Vec[i] << endl;
	cout << endl;
	*/


	s->doDelete( model );
	cout << "." << flush;
}

void testCylDiffnWithStoich()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 1e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt0 = 0.1; // Used for diffusion. 0.2 is the highest dt for which the error is in bounds
	double dt1 = 1; // Used for chem.
	double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool1 = s->doCreate( "Pool", cyl, "pool1", 1 );
	Id pool2 = s->doCreate( "Pool", cyl, "pool2", 1 );
	Field< double >::set( pool1, "diffConst", diffConst );
	Field< double >::set( pool2, "diffConst", diffConst/2 );

	Id stoich = s->doCreate( "Stoich", model, "stoich", 1 );
	Id ksolve = s->doCreate( "Ksolve", model, "ksolve", 1 );
	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( stoich, "compartment", cyl );
	Field< Id >::set( stoich, "ksolve", ksolve );
	Field< Id >::set( stoich, "dsolve", dsolve );
	Field< string >::set( stoich, "path", "/model/cyl/#" );
	assert( pool1.element()->numData() == ndc );

	// Then find a way to test it.
	vector< double > poolVec;
	Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
	Field< double >::set( ObjId( pool2, 0 ), "nInit", 1.0 );
   	Field< double >::getVec( pool1, "nInit", poolVec );
	assert( poolVec.size() == ndc );
	assert( doubleEq( poolVec[0], 1.0 ) );
	assert( doubleEq( poolVec[1], 0.0 ) );

	vector< double > nvec = 
		LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
	assert( nvec.size() == ndc );

	// Next: build by doing reinit
	s->doUseClock( "/model/dsolve", "process", 0 );
	s->doUseClock( "/model/ksolve", "process", 1 );
	s->doSetClock( 0, dt0 );
	s->doSetClock( 1, dt1 );
	s->doReinit();
	s->doStart( runtime );

	nvec = LookupField< unsigned int, vector< double > >::get( 
						dsolve, "nVec", 0);
   	Field< double >::getVec( pool1, "n", poolVec );
	assert( nvec.size() == poolVec.size() );
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		assert( doubleEq( nvec[i], poolVec[i] ) );
	/*
	cout << endl;
	for ( unsigned int i = 0; i < nvec.size(); ++i )
		cout << nvec[i] << "	";
	cout << endl;
	*/

	double dx = diffLength;
	double err = 0.0;
	double analyticTot = 0.0;
	double myTot = 0.0;
	for ( unsigned int i = 0; i < nvec.size(); ++i ) {
		double x = i * dx + dx * 0.5;
		// This part is the solution as a func of x,t.
		double y = dx *  // This part represents the init n of 1 in dx
			( 1.0 / sqrt( PI * diffConst * runtime ) ) * 
			exp( -x * x / ( 4 * diffConst * runtime ) ); 
		err += ( y - nvec[i] ) * ( y - nvec[i] );
		//cout << i << "	" << x << "	" << y << "	" << conc[i] << endl;
		analyticTot += y;
		myTot += nvec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );
	// cout << "analyticTot= " << analyticTot << ", myTot= " << myTot << endl;
	assert( err < 1.0e-5 );


	s->doDelete( model );
	cout << "." << flush;
}
#if 0
void testBuildTree()
{
	static double test[] = {
		1,  2,  3,  4,  0,  0,  0,  0,  0,  0,  0,
		5,  6,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		7,  0,  8,  9,  0,  0,  0,  0, 10,  0,  0,
		11, 0, 12, 13,  0,  0,  0,  0,  0,  0, 14,
		0,  0,  0,  0, 15, 16,  0,  0, 17, 18,  0,
		0,  0,  0,  0, 19, 20, 21,  0, 22,  0,  0,
		0,  0,  0,  0,  0, 23, 24, 25,  0,  0,  0,
		0,  0,  0,  0,  0,  0, 26, 27,  0,  0,  0,
		0,  0, 28,  0, 29, 30,  0,  0, 31,  0,  0,
		0,  0,  0,  0, 32,  0,  0,  0,  0, 33,  0,
		0,  0,  0, 34,  0,  0,  0,  0,  0,  0, 35,
	};
	const unsigned int numCompts = 11;
	FastMatrixElim fe;
	fe.makeTestMatrix( test, numCompts );
	fe.print();
	vector< unsigned int > parentVoxel;
	bool ret = fe.buildTree( 0, parentVoxel );
	assert( ret );
	assert( parentVoxel[0] == static_cast< unsigned int >( -1 ) );
	assert( parentVoxel[1] == 0 );
	assert( parentVoxel[2] == 0 );
	assert( parentVoxel[3] == 0 );
	assert( parentVoxel[8] == 2 );
	assert( parentVoxel[10] == 3 );
	assert( parentVoxel[5] == 4 );
	assert( parentVoxel[9] == 4 );
	assert( parentVoxel[6] == 5 );
	assert( parentVoxel[7] == 6 );

	assert( parentVoxel[10] == 2 );

	/*
	 * This is the sequence of traversal. x means empty, s means sibling,
	 * . means below diagonal. The numbers are the sequence.
	static double traverseIndex[] = {
	// col  1   2   3   4   5   6   7   8   9   10
		#,  1,  2,  3,  x,  x,  x,  x,  x,  x,  x,
		.,  #,  x,  x,  x,  x,  x,  x,  x,  x,  x,
		.,  x,  #,  s,  x,  x,  x,  x,  4,  x,  x,
		.,  x,  .,  #,  x,  x,  x,  x,  x,  x,  5,
		x,  x,  x,  x,  #,  6,  x,  x,  s,  7,  x,
		x,  x,  x,  x,  .,  #,  8,  x,  s,  x,  x,
		x,  x,  x,  x,  x,  .,  #,  9,  x,  x,  x,
		x,  x,  x,  x,  x,  x,  .,  #,  x,  x,  x,
		x,  x,  .,  x,  .,  .,  x,  x,  #,  x,  x,
		x,  x,  x,  x,  .,  x,  x,  x,  x,  #,  x,
		x,  x,  x,  .,  x,  x,  x,  x,  x,  x,  #,
	};
	*/
}
#endif

void testCalcJunction()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	// Make a neuron with same-size dend and spine. PSD is tiny.
	// Put a, b, c in dend, b, c, d in spine, c, d, f in psd. No reacs.
	// See settling of all concs by diffusion, pairwise.
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id dend = s->doCreate( "Compartment", model, "dend", 1 );
	Id neck = s->doCreate( "Compartment", model, "spine_neck", 1 );
	Id head = s->doCreate( "Compartment", model, "spine_head", 1 );
	Field< double >::set( dend, "x", 10e-6 );
	Field< double >::set( dend, "diameter", 2e-6 );
	Field< double >::set( dend, "length", 10e-6 );
	Field< double >::set( neck, "x0", 9e-6 );
	Field< double >::set( neck, "x", 9e-6 );
	Field< double >::set( neck, "y", 1e-6 );
	Field< double >::set( neck, "diameter", 0.5e-6 );
	Field< double >::set( neck, "length", 1.0e-6 );
	Field< double >::set( head, "x0", 9e-6 );
	Field< double >::set( head, "x", 9e-6 );
	Field< double >::set( head, "y0", 1e-6 );
	Field< double >::set( head, "y", 11e-6 );
	Field< double >::set( head, "diameter", 2e-6 );
	Field< double >::set( head, "length", 10e-6 );
	s->doAddMsg( "Single", ObjId( dend ), "raxial", ObjId( neck ), "axial");
	s->doAddMsg( "Single", ObjId( neck ), "raxial", ObjId( head ), "axial");

	Id nm = s->doCreate( "NeuroMesh", model, "nm", 1 );
	Field< double >::set( nm, "diffLength", 10e-6 );
	Field< bool >::set( nm, "separateSpines", true );
	Id sm = s->doCreate( "SpineMesh", model, "sm", 1 );
	Id pm = s->doCreate( "PsdMesh", model, "pm", 1 );
	ObjId mid = s->doAddMsg( "Single", ObjId( nm ), "spineListOut", ObjId( sm ), "spineList" );
	assert( !mid.bad() );
	mid = s->doAddMsg( "Single", ObjId( nm ), "psdListOut", ObjId( pm ), "psdList" );
	Field< string >::set( nm, "subTreePath", "/model/#" );
	//SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );

	vector< Id > pools( 9 );
	static string names[] = {"a", "b", "c", "b", "c", "d", "c", "d", "e" };
	static Id parents[] = {nm, nm, nm, sm, sm, sm, pm, pm, pm};
	for ( unsigned int i = 0; i < 9; ++i ) {
		pools[i] = s->doCreate( "Pool", parents[i], names[i], 1 );
		assert( pools[i] != Id() );
		Field< double >::set( pools[i], "concInit", 1.0 + 1.0 * i );
		Field< double >::set( pools[i], "diffConst", 1e-11 );
		if ( i < 6 ) {
			double vol = Field< double >::get( pools[i], "volume" );
			assert( doubleEq( vol, 10e-6 * 1e-12 * PI ) );
		}
	}
	Id dendsolve = s->doCreate( "Dsolve", model, "dendsolve", 1 );
	Id spinesolve = s->doCreate( "Dsolve", model, "spinesolve", 1 );
	Id psdsolve = s->doCreate( "Dsolve", model, "psdsolve", 1 );
	Field< Id >::set( dendsolve, "compartment", nm );
	Field< Id >::set( spinesolve, "compartment", sm );
	Field< Id >::set( psdsolve, "compartment", pm );
	Field< string >::set( dendsolve, "path", "/model/nm/#" );
	Field< string >::set( spinesolve, "path", "/model/sm/#" );
	Field< string >::set( psdsolve, "path", "/model/pm/#" );
	assert( Field< unsigned int >::get( dendsolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( spinesolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( psdsolve, "numAllVoxels" ) == 1 );
	assert( Field< unsigned int >::get( dendsolve, "numPools" ) == 3 );
	assert( Field< unsigned int >::get( spinesolve, "numPools" ) == 3 );
	assert( Field< unsigned int >::get( psdsolve, "numPools" ) == 3 );
	SetGet2< Id, Id >::set( dendsolve, "buildNeuroMeshJunctions", 
					spinesolve, psdsolve );
	s->doSetClock( 0, 0.01 );
	s->doUseClock( "/model/#solve", "process", 0 );
	s->doReinit();
	s->doStart( 100 );

	/*
	for ( unsigned int i = 0; i < 9; ++i ) {
		double c = Field< double >::get( pools[i], "conc" );
		double n = Field< double >::get( pools[i], "n" );
		double v = Field< double >::get( pools[i], "volume" );
		cout << pools[i].path() << ": " << c << ", " << n << ", " <<
				n / v << ", " <<
				v << endl;
	}
	*/
	s->doDelete( model );
	cout << "." << flush;
}

void testDiffusion()
{
	testSorting();
	testFastMatrixElim();
	testSetDiffusionAndTransport();
	testCylDiffn();
	testTaperingCylDiffn();
	testSmallCellDiffn();
	testCellDiffn();
	testCylDiffnWithStoich();
	testCalcJunction();
}