def reconstructionError(sp, inputVectors, activeColumnVectors, threshold=0.):
  """
  Computes a reconstruction error. The reconstuction $r(x)$ of an input vector $x$
  is given by the sum of the active column's connected synapses vector of 
  the SDR representation $sdr(x)$ of $x$ normalized by $1/numActiveColumns$. 
  The error is the normalized sum over the "hamming distance" (i.e. distance 
  induced by L1 norm) of $x$ and its reconstruction $r(x)$, i.e. (mathy stuff in LaTex)
  \[
      Reconstruction Error = (1/batchSize) * \sum_{x \in InputBatch} \| x - r(x) \|_1 .
  \]
  Note that $r(x)$ can be expressed as
  \[
      r(x) = (1/numActiveColumns) * C * sdr(x) ,
  \]
  where we view $sdr(x)$ as a binary column vector and $C$ is the 
  binary matrix whose jth column encodes the synaptic connectivity of 
  the pooler's columns and the input bits, i.e. 
  \[
        c_{i,j} = 1   :<=>  column j has a stable synaptic 
                            connection to input bit i.
  \]
  Note: Turns out that in our setting (x and syn(i) binary vectors) we have 
  \[
        Reconstruction Error = Witness Error.
  \]
  It can be shown that the error is optimized by the Hebbian-like update rule 
  of the spatial pooler. 
  
  @param sp (SpatialPooler) the spatial pooler instance
  @param inputVectors (array) 2D numpy array of input vectors
  @param activeColumnVectors (array) 2D numpy array of activation history
  @param threshold (float) if set > 0 it serves as threshold for a step function 
                           applied to the reconstruction vectors (values smaller than 
                           threshold are set to zero, and values bigger to one)
  @return error (float) the reconstruction error
  """
  batchSize        = inputVectors.shape[0]
  connectionMatrix = getConnectedSyns(sp)

  reconstructionVectors = np.dot(activeColumnVectors, connectionMatrix)
  # numActiveColumns      = np.sum(activeColumnVectors, 1)[0]
  numActiveColumns = int(sp._localAreaDensity * sp._numColumns) + 0.0
  reconstructionVectors = reconstructionVectors/numActiveColumns

  if threshold > 0.:
    reconstructionVectors =np.where( 
                            reconstructionVectors > threshold, 
                            np.ones( reconstructionVectors.shape), 
                            np.zeros(reconstructionVectors.shape))

  Err = np.sum(np.absolute(reconstructionVectors  - inputVectors))

  return Err/batchSize