VectorMatrixCU32Accessor::VectorMatrixCU32Accessor(VectorMatrix &mat) : mat(mat)
{
    mat.writeLock(1); // 1 = CU32Device!
    data_x = static_cast<CU32Array*>(mat.getArray(1, 0))->ptr();
    data_y = static_cast<CU32Array*>(mat.getArray(1, 1))->ptr();
    data_z = static_cast<CU32Array*>(mat.getArray(1, 2))->ptr();
}

void gradient_cpu(double delta_x, double delta_y, double delta_z, const double *phi, VectorMatrix &field)
{
	const int dim_x = field.dimX();
	const int dim_y = field.dimY();
	const int dim_z = field.dimZ();

	const int phi_sx = 1;
	const int phi_sy = (dim_x+1);
	const int phi_sz = (dim_x+1)*(dim_y+1);

	VectorMatrix::accessor field_acc(field);

	for (int z=0; z<dim_z; ++z)
	for (int y=0; y<dim_y; ++y)
	for (int x=0; x<dim_x; ++x) {
		const int i = phi_sx*x + phi_sy*y + phi_sz*z;

		const double dx = (+ phi[i+1*phi_sx+0*phi_sy+0*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+0*phi_sz] - phi[i+0*phi_sx+1*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+0*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+1*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+1*phi_sy+1*phi_sz]) / (4.0 * delta_x);

		const double dy = (+ phi[i+0*phi_sx+1*phi_sy+0*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+0*phi_sz] - phi[i+1*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+0*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+1*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+1*phi_sx+0*phi_sy+1*phi_sz]) / (4.0 * delta_y);

		const double dz = (+ phi[i+0*phi_sx+0*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+0*phi_sy+1*phi_sz] - phi[i+1*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+0*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+1*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+1*phi_sx+1*phi_sy+0*phi_sz]) / (4.0 * delta_z);

		field_acc.set(x, y, z, Vector3d(dx, dy, dz));
	}
}

VectorVectorConvolution_Simple::VectorVectorConvolution_Simple(const VectorMatrix &lhs, int dim_x, int dim_y, int dim_z)
	: lhs(lhs), dim_x(dim_x), dim_y(dim_y), dim_z(dim_z)
{
	assert(lhs.getShape().getDim(0) == 3);
	exp_x = lhs.getShape().getDim(1);
	exp_y = lhs.getShape().getDim(2);
	exp_z = lhs.getShape().getDim(3);
}

 // must pass in a lvalue
 void copy_to( VectorMatrix<T>& vector_matrix_obj ) {
   vector_matrix_obj.resize( this->nrow_, this->ncol_ );
   const size_t length = this->nrow_ * this->ncol_;
   std :: vector<T> temp_store;
   temp_store.resize(length);
   std :: copy_n( this->store_.begin(), length, temp_store.begin() );
   vector_matrix_obj.set_store() = std :: move( temp_store );
 }

double rk_scaled_error_norm_cpu(double h, double eps_abs, double eps_rel, const VectorMatrix &y, const VectorMatrix &y_error)
{
	double norm;
	
	if (false) { // Maximum norm
		double max = 0.0;

		VectorMatrix::const_accessor y_acc(y), y_error_acc(y_error);

		const size_t s = y.size() * 3;
		for (size_t i=0; i<s; ++i) {
			const Vector3d value =       y_acc.get(i);
			const Vector3d error = y_error_acc.get(i);
			{
				const double D0 = eps_rel * std::fabs(value.x) + eps_abs;
				const double r  = std::fabs(error.x) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			} {
				const double D0 = eps_rel * std::fabs(value.y) + eps_abs;
				const double r  = std::fabs(error.y) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			} {
				const double D0 = eps_rel * std::fabs(value.z) + eps_abs;
				const double r  = std::fabs(error.z) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			}
		}
		norm = max;

	} else { // Euclidian norm
		double sum = 0.0;

		VectorMatrix::const_accessor y_acc(y), y_error_acc(y_error);
		const size_t s = y.size();
		for (size_t i=0; i<s; ++i) {
			const Vector3d value =       y_acc.get(i);
			const Vector3d error = y_error_acc.get(i);

			{
				const double D0 = eps_rel * std::fabs(value.x) + eps_abs;
				const double r  = std::fabs(error.x) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			} {
				const double D0 = eps_rel * std::fabs(value.y) + eps_abs;
				const double r  = std::fabs(error.y) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			} {
				const double D0 = eps_rel * std::fabs(value.z) + eps_abs;
				const double r  = std::fabs(error.z) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			}
		}
		norm = std::sqrt(sum / (s*3));
	}

	return norm;
}

Vector3d findExtremum_cpu(VectorMatrix &M, int z_slice, int component)
{
	if (M.getShape().getRank() != 3) {
		throw std::runtime_error("findExtremum: Fixme: Need matrix of rank 3");
	}

	if (component < 0 || component > 2) {
		throw std::runtime_error("findExtremum: Invalid 'component' value, must be 0, 1 or 2.");
	}

	const int dim_x = M.getShape().getDim(0);
	const int dim_y = M.getShape().getDim(1);

	VectorMatrix::const_accessor M_acc(M);

	// Find cell with maximum absolute value
	double max_val = -1.0;
	int max_x = -1, max_y = -1;
	for (int y=1; y<dim_y-1; ++y)
	for (int x=1; x<dim_x-1; ++x) {
		const int val = std::fabs(M_acc.get(x, y, z_slice)[component]);
		if (val > max_val) {
			max_val = val;
			max_x = x;
			max_y = y;
		}
	}
	assert(max_x > 0);
	assert(max_y > 0);
	
	// Refine maximum by fitting to sub-cell precision
	const double xdir_vals[3] = {
		M_acc.get(max_x-1, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+1, max_y+0, z_slice)[component]
	};

	const double ydir_vals[3] = {
		M_acc.get(max_x+0, max_y-1, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+1, z_slice)[component]
	};

	return Vector3d(
		fit(max_x-1, max_x+0, max_x+1, xdir_vals[0], xdir_vals[1], xdir_vals[2]),
		fit(max_y-1, max_y+0, max_y+1, ydir_vals[0], ydir_vals[1], ydir_vals[2]),
		static_cast<double>(z_slice)
	);
}

void VectorMatrix::add(const VectorMatrix &op, double factor)
{
	if (this == &op) {
		scale(1.0 + factor);
	} else if (isUniform() && op.isUniform()) {
		fill(getUniformValue() + op.getUniformValue() * factor);
	} else {
		const int dev = computeStrategy2(op);
		writeLock(dev); op.readLock(dev);
		for (int c=0; c<num_arrays; ++c) {
			matty::getDevice(dev)->add(getArray(dev, c), op.getArray(dev, c), factor);
		}
		writeUnlock(dev); op.readUnlock(dev);
	}
}

void VectorMatrix::assign(const VectorMatrix &op)
{
	if (this == &op) {
		return;
	} else if (op.isUniform()) {
		fill(op.getUniformValue());
	} else {
		const int dev = computeStrategy2(op);
		writeLock(dev); op.readLock(dev);
		for (int c=0; c<num_arrays; ++c) {
			matty::getDevice(dev)->assign(getArray(dev, c), op.getArray(dev, c));
		}
		writeUnlock(dev); op.readUnlock(dev);
	}
}

void rk_combine_result_cpu(
	double h, ButcherTableau &tab,
	const VectorMatrix &k0, const VectorMatrix &k1, const VectorMatrix &k2, const VectorMatrix &k3,
	VectorMatrix &y, VectorMatrix &y_error)
{
	const size_t s = y.size();
	VectorMatrix::accessor y_acc(y), y_error_acc(y_error);

	// tab
	const double  c0 = tab. c[0],  c1 = tab. c[1],  c2 = tab. c[2],  c3 = tab. c[3];
	const double ec0 = tab.ec[0], ec1 = tab.ec[1], ec2 = tab.ec[2], ec3 = tab.ec[3];

	VectorMatrix::const_accessor k0_acc(k0), k1_acc(k1), k2_acc(k2), k3_acc(k3);
	for (size_t i=0; i<s; ++i) {
		const Vector3d y_i = y_acc.get(i);
		y_acc.set(i, 
			y_i + h * (  c0*k0_acc.get(i) 
			           + c1*k1_acc.get(i)
				   + c2*k2_acc.get(i)
				   + c3*k3_acc.get(i))
		);
		y_error_acc.set(i,
			h * (  ec0*k0_acc.get(i) 
			     + ec1*k1_acc.get(i)
			     + ec2*k2_acc.get(i)
			     + ec3*k3_acc.get(i))
		);
	}
}

 void copy_from( const VectorMatrix<T>& vector_matrix_obj ) {
   try {
     const size_t obj_size = vector_matrix_obj.size();
     if( ( obj_size < STACK_DOUBLE_LIMIT ) == false ) {
       throw obj_size;
     }
     this->nrow_ = vector_matrix_obj.nrow();
     this->ncol_ = vector_matrix_obj.ncol();
     std :: vector<T> temp_store = vector_matrix_obj.store();
     std :: copy_n( temp_store.begin(), obj_size, this->store_.begin() );
   }
   catch( size_t n ) {
     std :: cout << " array load exception: " << n << " >= " << STACK_DOUBLE_LIMIT << std :: endl;
     abort();
   }
 }

/*
  @method DistanceList
  @discussion Allocate storage for an n-unit unitList.
  @throws KII_bad_alloc
*/
DistanceList::DistanceList(const FLOAT xlocs[], const FLOAT ylocs[], 
			   const VectorMatrix& radii) 
  : unitList(NULL), numUnits(radii.Size()), numOverlap(0)
{
  // Since this is a symmetric array, there's no need to store
  // information for unit n (all that information is already in the
  // lists for units 1..n-1
  if ((unitList = new ListItem[numUnits-1]) == NULL)
    throw KII_bad_alloc("Failed allocating memory for DistanceList.");

  // Fill the lists. For unit i, we need only store information for
  // connections to units i+1..n, as earlier units already have the
  // information for their connections with i.
  FLOAT tempDist, tempDist2, tempDelta, deltaX, deltaY;
  for (int u1=0; u1<numUnits-1; u1++) {
    unitList[u1].radius = radii[u1];
    for (int u2=u1+1; u2<numUnits; u2++) {
      deltaX = xlocs[u1]-xlocs[u2];
      deltaY = ylocs[u1]-ylocs[u2];
      tempDist2 = deltaX*deltaX + deltaY*deltaY;
      tempDist = sqrt(tempDist2);
      tempDelta = tempDist - (radii[u1]+radii[u2]);
      if (tempDelta < 0.0) {
	unitList[u1].overlapping.push_back(SublistItem(u2, radii[u2], tempDist,
						       tempDist2, tempDelta));
	numOverlap++;
      } else {
	unitList[u1].nonOverlapping.push_back(SublistItem(u2, radii[u2],
							  tempDist, tempDist2, tempDelta));
      }
    }
  }
}

void rk_combine_result(
	double h, ButcherTableau &tab,
	const VectorMatrix &k0, const VectorMatrix &k1, const VectorMatrix &k2,
	const VectorMatrix &k3, const VectorMatrix &k4, const VectorMatrix &k5,
	VectorMatrix &y, VectorMatrix &y_error)
{
	const int s = y.size();
	if (   s != y_error.size()
	    || s != k1.size()
	    || s != k2.size()
	    || s != k3.size()
	    || s != k4.size()
	    || s != k5.size()) throw std::runtime_error("rk_combine_result: Input matrix size mismatch.");
	if (!tab.num_steps == 6) throw std::runtime_error("Need num_steps == 6 in rk_combine_result");

	if (isCudaEnabled()) {
#ifdef HAVE_CUDA
		rk_combine_result_cuda(h, tab, k0, k1, k2, k3, k4, k5, y, y_error, isCuda64Enabled());
#else
		assert(0);
#endif
	} else {
		rk_combine_result_cpu(h, tab, k0, k1, k2, k3, k4, k5, y, y_error);
	}
}

VectorMatrix linearInterpolate(const VectorMatrix &src, Shape dest_dim)
{
	Shape src_dim = src.getShape();

	if (src_dim.getRank() != dest_dim.getRank()) {
		throw std::runtime_error("linearInterpolate: Source and destination matrices need to have the same rank.");
	}

	if (src_dim.getRank() != 3) {
		throw std::runtime_error("linearInterpolate: Fixme: Need to have matrix of rank 3");
	}

	VectorMatrix dest(dest_dim);

	VectorMatrix::      accessor dest_acc(dest);
	VectorMatrix::const_accessor  src_acc(src);

	const bool sing_x = (src_dim.getDim(0) == 1);
	const bool sing_y = (src_dim.getDim(1) == 1);
	const bool sing_z = (src_dim.getDim(2) == 1);

	Vector3d scale(1.0, 1.0, 1.0);
	if (!sing_x) scale.x = double(dest_dim.getDim(0)-1) / double(src_dim.getDim(0)-1);
	if (!sing_y) scale.y = double(dest_dim.getDim(1)-1) / double(src_dim.getDim(1)-1);
	if (!sing_z) scale.z = double(dest_dim.getDim(2)-1) / double(src_dim.getDim(2)-1);

	for (int k=0; k<dest_dim.getDim(2); ++k)
	for (int j=0; j<dest_dim.getDim(1); ++j)
	for (int i=0; i<dest_dim.getDim(0); ++i) {
		// (x,y,z): coordinates of point with indices (i,j,k) in dst matrix
		const double x = i / scale.x;
		const double y = j / scale.y;
		const double z = k / scale.z;

		const double u = x - std::floor(x);
		const double v = y - std::floor(y);
		const double w = z - std::floor(z);

		const int I = std::floor(x);
		const int J = std::floor(y);
		const int K = std::floor(z);

		Vector3d tmp(0.0, 0.0, 0.0);
		if (true                         ) tmp = tmp + (1.0-u) * (1.0-v) * (1.0-w) * src_acc.get(I  , J  , K  );
		if (                      !sing_z) tmp = tmp + (1.0-u) * (1.0-v) *      w  * src_acc.get(I  , J  , K+1);
		if (           !sing_y           ) tmp = tmp + (1.0-u) *      v  * (1.0-w) * src_acc.get(I  , J+1, K  );
		if (           !sing_y && !sing_z) tmp = tmp + (1.0-u) *      v  *      w  * src_acc.get(I  , J+1, K+1);
		if (!sing_x                      ) tmp = tmp +      u  * (1.0-v) * (1.0-w) * src_acc.get(I+1, J  , K  );
		if (!sing_x            && !sing_z) tmp = tmp +      u  * (1.0-v) *      w  * src_acc.get(I+1, J  , K+1);
		if (!sing_x && !sing_y           ) tmp = tmp +      u  *      v  * (1.0-w) * src_acc.get(I+1, J+1, K  );
		if (!sing_x && !sing_y && !sing_z) tmp = tmp +      u  *      v  *      w  * src_acc.get(I+1, J+1, K+1);
		dest_acc.set(i, j, k, tmp);
	}

	return dest;
}

void Transposer_CUDA::copy_pad(const VectorMatrix &M, float *out_x, float *out_y, float *out_z)
{
	// Ifdef HAVE_CUDA_64, we directly support input matrices that
	// are stored with 64 bit precision on the GPU.
#ifdef HAVE_CUDA_64
	const bool M_is_cuda64_bit = M.isCached(2); // 0 = CPU device, 2 = CUDA_64 device
	if (M_is_cuda64_bit) {
		VectorMatrix::const_cu64_accessor M_acc(M);
		// xyz, M -> s1
		cuda_copy_pad_r2r(dim_x, dim_y, dim_z, exp_x, M_acc.ptr_x(), M_acc.ptr_y(), M_acc.ptr_z(), out_x, out_y, out_z);
	} 
	else
#endif
	{
		VectorMatrix::const_cu32_accessor M_acc(M);
		// xyz, M -> s1
		cuda_copy_pad_r2r(dim_x, dim_y, dim_z, exp_x, M_acc.ptr_x(), M_acc.ptr_y(), M_acc.ptr_z(), out_x, out_y, out_z);
	}
}

ConstVectorMatrixAccessor::ConstVectorMatrixAccessor(const VectorMatrix &mat) : mat(mat) 
{
	mat.readLock(0);
	data_x = static_cast<CPUArray*>(mat.getArray(0, 0))->ptr();
	data_y = static_cast<CPUArray*>(mat.getArray(0, 1))->ptr();
	data_z = static_cast<CPUArray*>(mat.getArray(0, 2))->ptr();

	// Precalculate strides
	const int rank = mat.getShape().getRank();
	strides[0] = 1;
	strides[1] = strides[0] * (rank > 0 ? mat.getShape().getDim(0) : 1);
	strides[2] = strides[1] * (rank > 1 ? mat.getShape().getDim(1) : 1);
	strides[3] = strides[2] * (rank > 2 ? mat.getShape().getDim(2) : 1);
}

double VectorMatrix::dotSum(const VectorMatrix &other) const
{
	if (isUniform() && other.isUniform()) {
		const double x = uval[0], y = uval[1], z = uval[2];
		const double dot = x*x + y*y + z*z;
		return size() * dot;
	} else {
		const int dev = computeStrategy2(other);
		readLock(dev); if (this != &other) other.readLock(dev);
		const double sum = matty::getDevice(dev)->sumdot3(
			this->getArray(dev, 0), this->getArray(dev, 1), this->getArray(dev, 2),
			other.getArray(dev, 0), other.getArray(dev, 1), other.getArray(dev, 2)
		);
		if (this != &other) other.readUnlock(dev); readUnlock(dev);
		return sum;
	}
}

void rk_prepare_step_cpu(
	int step,
	double h,
	ButcherTableau &tab,

	const VectorMatrix &k0,
	const VectorMatrix &k1,
	const VectorMatrix &k2,
	const VectorMatrix &k3,
	const VectorMatrix &k4,
	const VectorMatrix &k5,

	const VectorMatrix &y,
	VectorMatrix &ytmp)
{
	const size_t s = y.size();

	VectorMatrix::const_accessor y_acc(y); 
	VectorMatrix::const_accessor k0_acc(k0), k1_acc(k1), k2_acc(k2), k3_acc(k3), k4_acc(k4), k5_acc(k5);
	VectorMatrix::accessor ytmp_acc(ytmp);

	const double b10 = tab.b[1][0];
	const double b20 = tab.b[2][0], b21 = tab.b[2][1];
	const double b30 = tab.b[3][0], b31 = tab.b[3][1], b32 = tab.b[3][2];
	const double b40 = tab.b[4][0], b41 = tab.b[4][1], b42 = tab.b[4][2], b43 = tab.b[4][3];
	const double b50 = tab.b[5][0], b51 = tab.b[5][1], b52 = tab.b[5][2], b53 = tab.b[5][3], b54 = tab.b[5][4];

	switch (step) {
		case 1:
			for (size_t i=0; i<s; ++i) {
				ytmp_acc.set(i, y_acc.get(i) + h * (b10*k0_acc.get(i)));
			}
			break;

		case 2:
			for (size_t i=0; i<s; ++i) {
				ytmp_acc.set(i, 
					y_acc.get(i) + h * (  b20*k0_acc.get(i) 
				                            + b21*k1_acc.get(i))
				);
			}
			break;

		case 3: 
			for (size_t i=0; i<s; ++i) {
				ytmp_acc.set(i, 
					y_acc.get(i) + h * (  b30*k0_acc.get(i) 
				                            + b31*k1_acc.get(i) 
			                                    + b32*k2_acc.get(i)) 
				);
			}
			break;

		case 4:
			for (size_t i=0; i<s; ++i) {
				ytmp_acc.set(i, 
					y_acc.get(i) + h * (  b40*k0_acc.get(i) 
					                    + b41*k1_acc.get(i)
					                    + b42*k2_acc.get(i)
					                    + b43*k3_acc.get(i))
				);
			}
			break;

		case 5:
			for (size_t i=0; i<s; ++i) {
				ytmp_acc.set(i, 
					y_acc.get(i) + h * (  b50*k0_acc.get(i) 
					                    + b51*k1_acc.get(i)
					                    + b52*k2_acc.get(i)
					                    + b53*k3_acc.get(i)
					                    + b54*k4_acc.get(i))
				);
			}
			break;

		default:
			throw std::runtime_error("Cant handle runge-kutta methods with more than 6 steps (not implemented)");
	}
}

/*
 * Writes simulation results to an output destination.
 *
 * @param  neurons the Neuron list to search from.
 **/
void XmlRecorder::saveSimData(vector<Cluster *> &vtClr, vector<ClusterInfo *> &vtClrInfo)
{
    // create Neuron Types matrix
    VectorMatrix neuronTypes(MATRIX_TYPE, MATRIX_INIT, 1, m_sim_info->totalNeurons, EXC);
    for (int i = 0; i < m_sim_info->totalNeurons; i++) {
        neuronTypes[i] = m_model->getLayout()->neuron_type_map[i];
    }

    // create neuron threshold matrix
    VectorMatrix neuronThresh(MATRIX_TYPE, MATRIX_INIT, 1, m_sim_info->totalNeurons, 0);
    for (CLUSTER_INDEX_TYPE iCluster = 0; iCluster < vtClr.size(); iCluster++) {
        AllIFNeurons *neurons = dynamic_cast<AllIFNeurons*>(vtClr[iCluster]->m_neurons);

        int neuronLayoutIndex = vtClrInfo[iCluster]->clusterNeuronsBegin;
        int totalClusterNeurons = vtClrInfo[iCluster]->totalClusterNeurons;
        for (int iNeurons = 0; iNeurons < totalClusterNeurons; iNeurons++, neuronLayoutIndex++) {
            neuronThresh[neuronLayoutIndex] = neurons->Vthresh[iNeurons];
        }
    }

    // Write XML header information:
    stateOut << "<?xml version=\"1.0\" standalone=\"no\"?>\n" << "<!-- State output file for the DCT growth modeling-->\n";
    //stateOut << version; TODO: version

    // Write the core state information:
    VectorMatrix* xloc = new VectorMatrix(MATRIX_TYPE, MATRIX_INIT, 1, m_sim_info->totalNeurons);
    VectorMatrix* yloc = new VectorMatrix(MATRIX_TYPE, MATRIX_INIT, 1, m_sim_info->totalNeurons);
    for (int i = 0; i < m_sim_info->totalNeurons; i++) {
        (*xloc)[i] = m_model->getLayout()->xloc[i];
        (*yloc)[i] = m_model->getLayout()->yloc[i];
    }

    stateOut << "<SimState>\n";
    stateOut << "   " << burstinessHist.toXML("burstinessHist") << endl;
    stateOut << "   " << spikesHistory.toXML("spikesHistory") << endl;
    stateOut << "   " << xloc->toXML("xloc") << endl;
    stateOut << "   " << yloc->toXML("yloc") << endl;
    stateOut << "   " << neuronTypes.toXML("neuronTypes") << endl;

    delete xloc;
    delete yloc;

    // create starter nuerons matrix
    int num_starter_neurons = static_cast<int>(m_model->getLayout()->num_endogenously_active_neurons);
    if (num_starter_neurons > 0)
    {
        VectorMatrix starterNeurons(MATRIX_TYPE, MATRIX_INIT, 1, num_starter_neurons);
        getStarterNeuronMatrix(starterNeurons, m_model->getLayout()->starter_map, m_sim_info);
        stateOut << "   " << starterNeurons.toXML("starterNeurons") << endl;
    }

    // Write neuron thresold
    stateOut << "   " << neuronThresh.toXML("neuronThresh") << endl;

    // write time between growth cycles
    stateOut << "   <Matrix name=\"Tsim\" type=\"complete\" rows=\"1\" columns=\"1\" multiplier=\"1.0\">" << endl;
    stateOut << "   " << m_sim_info->epochDuration << endl;
    stateOut << "</Matrix>" << endl;

    // write simulation end time
    stateOut << "   <Matrix name=\"simulationEndTime\" type=\"complete\" rows=\"1\" columns=\"1\" multiplier=\"1.0\">" << endl;
    stateOut << "   " << g_simulationStep * m_sim_info->deltaT << endl;
    stateOut << "</Matrix>" << endl;
    stateOut << "</SimState>" << endl;
}

/*
  @method Update
  @discussion Update the DistanceList information based on the given
  new radii.
  @param radii unit connectivity radii must have same number of
  elements as numUnits
  @throws KII_invalid_argument
*/
void DistanceList::Update(const VectorMatrix& radii)
{
  if (radii.Size() != numUnits)
    throw KII_invalid_argument("Wrong number of elements in radii for distance update.");

#ifdef DEBUG2
  cerr << "Updating DistanceList with radii " << radii << endl;
#endif
  FLOAT tempDelta;
  SublistItem tempItem;
  list<SublistItem>::iterator u2iter;
  // Iterate through all units
  for (int u1=0; u1<numUnits-1; u1++) {
#ifdef DEBUG2
    cerr << u1 << ", ";
#endif
    // Update unit connection radius
    unitList[u1].radius = radii[u1];
    // Iterate over all overlapping and non-overlapping other units
    // First, overlapping:
    u2iter = unitList[u1].overlapping.begin();
    while (u2iter != unitList[u1].overlapping.end()) {
      // Compute and update other unit's information
      u2iter->radius = radii[u2iter->otherUnit];
      tempDelta = u2iter->dist - (unitList[u1].radius+u2iter->radius);
      u2iter->Delta = tempDelta;
      // Move it if non-overlapping, if necessary (Delta>=0). If the
      // move occurs, the iterator will point to the next item. If
      // not, then the iterator must be advanced explicitly.
      if (tempDelta >= 0) {
	tempItem = *u2iter;
	u2iter = unitList[u1].overlapping.erase(u2iter);
	unitList[u1].nonOverlapping.push_back(tempItem);
	numOverlap--;
      } else
	u2iter++;
    }
    // Then, non-overlapping:
    u2iter = unitList[u1].nonOverlapping.begin();
    while (u2iter != unitList[u1].nonOverlapping.end()) {
      // Compute and update other unit's information
      u2iter->radius = radii[u2iter->otherUnit];
      tempDelta = u2iter->dist - (unitList[u1].radius+u2iter->radius);
      u2iter->Delta = tempDelta;
      // Move it if overlapping, if necessary (Delta<0). If the
      // move occurs, the iterator will point to the next item. If
      // not, then the iterator must be advanced explicitly.
      if (tempDelta < 0) {
	tempItem = *u2iter;
	u2iter = unitList[u1].nonOverlapping.erase(u2iter);
	unitList[u1].overlapping.push_back(tempItem);
	numOverlap++;
#ifdef DLDEBUG
	cerr << "\tDL::Update(): overlap of " << u1 << " and " << tempItem.otherUnit
	     << " with radii " << unitList[u1].radius << " and " << tempItem.radius
	     << ", respectively" << endl;
#endif
      } else
	u2iter++;
    }
  }  // end for (int u1-0; ...)
#ifdef DEBUG2
  cerr << "done." << endl;;
#endif
}