def get_response_content(fs):
    # make the laplacian matrix for the graph
    weight = 1 / fs.edge_length
    n = fs.nvertices
    L = np.zeros((n,n))
    # set the diagonal
    for i in range(n):
        L[i,i] = 2 * weight
    L[0,0] = weight
    L[-1,-1] = weight
    # complete the tridiagonal
    for i in range(n-1):
        L[i+1,i] = -weight
        L[i,i+1] = -weight
    # define other matrices
    L_pinv = np.linalg.pinv(L)
    HDH = -2*L_pinv
    v = np.diag(HDH)
    e = np.ones(n)
    D = HDH - (np.outer(v, e) + np.outer(e, v))/2
    # show some matrices
    out = StringIO()
    np.set_printoptions(linewidth=300)
    print >> out, 'Laplacian matrix:'
    print >> out, L
    print >> out
    print >> out, 'HDH:'
    print >> out, HDH
    print >> out
    print >> out, 'EDM:'
    print >> out, D
    print >> out
    return out.getvalue()

	def updateA(self, featureVector, decay=None, current_time=None):
		if decay:
			assert decay <= 1 and decay >=0
			self.DD = self.decayAverage(decay, self.DD, np.outer(featureVector, featureVector), current_time)
			self.A = self.DD + self.identityMatrix
		else:
			self.A += np.outer(featureVector, featureVector)

def plot_surface(image, center=None, size=15, output=False, ds9_indexing=False, 
                 **kwargs):
    """
    Create a surface plot from image.
    
    By default, the whole image is plotted. The 'center' and 'size' attributes 
    allow to crop the image.
        
    Parameters
    ----------
    image : numpy.array
        The image as a numpy.array.
    center : tuple of 2 int (optional, default=None)
        If None, the whole image will be plotted. Otherwise, it grabs a square
        subimage at the 'center' (Y,X) from the image.
    size : int (optional, default=15)
        It corresponds to the size of a square in the image.
    output : {False, True}, bool optional
        Whether to output the grids and intensities or not.
    ds9_indexing : {False, True}, bool optional 
        If True the coordinates are in X,Y convention and in 1-indexed format.
    kwargs:
        Additional attributes are passed to the matplotlib figure() and 
        plot_surface() method.        
    
    Returns
    -------
    out : tuple of 3 numpy.array
        x and y for the grid, and the intensity
        
    """        
    if not center:
        size = image.shape[0]
        x = np.outer(np.arange(0,size,1), np.ones(size))
        y = x.copy().T 
        z = image
    else: 
        if ds9_indexing:
            center = (center[0]-1,center[1]-1) 
            cx, cy = center
        else:
            cy, cx = center
        if size % 2:            # if size is odd             
            x = np.outer(np.arange(0,size,1), np.ones(size))
        else:                   # otherwise, size is even
            x = np.outer(np.arange(0,size+1,1), np.ones(size+1))
        y = x.copy().T            
        z = image[cy-size//2:cy+size//2+1,cx-size//2:cx+size//2+1]           
    
    figure(figsize=kwargs.pop('figsize',(5,5)))
    ax = axes(projection='3d')
    ax.plot_surface(x, y, z, rstride=1, cstride=1, linewidth=0, **kwargs) 
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$I(x,y)$')
    ax.set_title('Data')
    show()
    
    if output:
        return (x,y,z)

def unscented_obs_model (f, sigmaPoints, meanWeight, covWeight, state, stateAngleMask, obsModelAngleMask, **kwargs):
    n = sigmaPoints.shape[1]
    stateLen = len (state)

    y = f (sigmaPoints[0:stateLen,0], **kwargs) + sigmaPoints[stateLen:,0]
    m = len (y)

    y = zeros ((m, n))

    for i in range (n):
        y[:,i] = f (sigmaPoints[0:stateLen,i], **kwargs) + sigmaPoints[stateLen:,i]

    muPrime = sum (y * meanWeight, axis=1)
    for i, mask in enumerate (obsModelAngleMask):
        if mask:
            muPrime[i] = circularMean (y[i,:], weights=meanWeight)
            muPrime[i] = minimizedAngle (muPrime[i])

    SigmaPrime = zeros ((m, m))
    for i in range (n):
        SigmaPrime += covWeight[i] * outer (minimizedAngle (y[:,i] - muPrime, obsModelAngleMask), 
                                            minimizedAngle (y[:,i] - muPrime, obsModelAngleMask))

    crossCov = zeros ((stateLen, m))
    for i in range (n):
        diffState = minimizedAngle (sigmaPoints[0:stateLen,i] - state, stateAngleMask)
        diffMeas = minimizedAngle (y[:,i] - muPrime, obsModelAngleMask)
        crossCov += covWeight[i] * outer (diffState, diffMeas)

    return muPrime, SigmaPrime, crossCov

def Cramer(var1, var2):
	"""
	Compute Cramer's V statistic for two Pandas series

	Parameters:
	----------
	var1, var2: Pandas series

	Returns:
	--------
	v : float
		The Cramer's V statistic of two categorical-variable series

	Status:
	-------	
	Cramer's V Implementation
	Author: Jesse Lund, [email protected]
	Date: 9/12/2015

	##Round 1##
	Comments: Thomas Roderick, [email protected]
	Date: 9/13/2015

	"""

	table = crosstab(var1,var2) #For Pandas: must have an index, can't just feed in two lists. This could be a sticking point. Might be better to do a check or roll our own crosstab implementation
	l,w = table.shape #save on a (small) function call here--reads in both outputs 
	df = min(l-1, w-1)
	colsum, rowsum = table.sum(0), table.sum(1) 
	n = float(l*w)
	expectmat = outer(rowsum,colsum)/n
	outmat = outer(table.sum(0),table.sum(1))/n #this works if same size
	return  sqrt((((table - expectmat)**2)/(expectmat*n*df)).sum().sum())

def generate_ODGD_spec_chirped(F1, F2, Fs, lengthOdgd=2048, Nfft=2048, \
                               Ot=0.5, t0=0.0, \
                               analysisWindowType='sinebell'):
    """
    generateODGDspecChirped:
    
    generates a waveform ODGD and the corresponding spectrum,
    using as analysis window the -optional- window given as
    argument.
    """
    
    # converting input to double:
    F1 = np.double(F1)
    F2 = np.double(F2)
    F0 = np.double(F1 + F2) / 2.0
    Fs = np.double(Fs)
    Ot = np.double(Ot)
    t0 = np.double(t0)
    
    # compute analysis window of given type:
    if analysisWindowType == 'sinebell':
        analysisWindow = sinebell(lengthOdgd)
    else:
        if analysisWindowType == 'hanning' or \
               analysisWindowType == 'hann':
            analysisWindow = hann(lengthOdgd)
    
    # maximum number of partials in the spectral comb:
    partialMax = np.floor((Fs / 2) / np.max(F1, F2))
    
    # Frequency numbers of the partials:
    frequency_numbers = np.arange(1,partialMax + 1)
    
    # intermediate value
    temp_array = 1j * 2.0 * np.pi * frequency_numbers * Ot
    
    # compute the amplitudes for each of the frequency peaks:
    amplitudes = F0 * 27 / 4 * \
                 (np.exp(-temp_array) \
                  + (2 * (1 + 2 * np.exp(-temp_array)) / temp_array) \
                  - (6 * (1 - np.exp(-temp_array)) \
                     / (temp_array ** 2))) \
                  / temp_array
    
    # Time stamps for the time domain ODGD
    timeStamps = np.arange(lengthOdgd) / Fs + t0 / F0
    
    # Time domain odgd:
    odgd = np.exp(2.0 * 1j * np.pi \
                  * (np.outer(F1 * frequency_numbers,timeStamps) \
                     + np.outer((F2 - F1) \
                                * frequency_numbers,timeStamps ** 2) \
                     / (2 * lengthOdgd / Fs))) \
                     * np.outer(amplitudes,np.ones(lengthOdgd))
    odgd = np.sum(odgd,axis=0)
    
    # spectrum:
    odgdSpectrum = np.fft.fft(real(odgd * analysisWindow), n=Nfft)
    
    return odgd, odgdSpectrum

    def _set_expected_stats(self,smoothed_mus,smoothed_sigmas,E_xtp1_xtT):
        assert not np.isnan(E_xtp1_xtT).any()
        assert not np.isnan(smoothed_mus).any()
        assert not np.isnan(smoothed_sigmas).any()

        data = self.data
        EyyT = data.T.dot(data)
        EyxT = data.T.dot(smoothed_mus)
        ExxT = smoothed_sigmas.sum(0) + smoothed_mus.T.dot(smoothed_mus)

        E_xt_xtT = \
            ExxT - (smoothed_sigmas[-1]
                    + np.outer(smoothed_mus[-1],smoothed_mus[-1]))
        E_xtp1_xtp1T = \
            ExxT - (smoothed_sigmas[0]
                    + np.outer(smoothed_mus[0], smoothed_mus[0]))

        E_xtp1_xtT = E_xtp1_xtT.sum(0)

        def is_symmetric(A):
            return np.allclose(A,A.T)

        assert is_symmetric(ExxT)
        assert is_symmetric(E_xt_xtT)
        assert is_symmetric(E_xtp1_xtp1T)

        self.E_emission_stats = np.array([EyyT, EyxT, ExxT, self.T])
        self.E_dynamics_stats = np.array([E_xtp1_xtp1T, E_xtp1_xtT, E_xt_xtT, self.T-1])

 def calculate(self, w, ionization_data, beta_rad, t_electrons, t_rad,
     beta_electron, delta_input, chi_0):
     if delta_input is None:
         if self.departure_coefficient is None:
             departure_coefficient = 1. / w
         else:
             departure_coefficient = self.departure_coefficient
         radiation_field_correction = -np.ones((len(ionization_data), len(
             beta_rad)))
         less_than_chi_0 = (
             ionization_data.ionization_energy < chi_0).values
         factor_a = (t_electrons / (departure_coefficient * w * t_rad))
         radiation_field_correction[~less_than_chi_0] = factor_a * \
             np.exp(np.outer(ionization_data.ionization_energy.values[
             ~less_than_chi_0], beta_rad - beta_electron))
         radiation_field_correction[less_than_chi_0] = 1 - np.exp(np.outer(
             ionization_data.ionization_energy.values[less_than_chi_0],
             beta_rad) - beta_rad * chi_0)
         radiation_field_correction[less_than_chi_0] += factor_a * np.exp(
             np.outer(ionization_data.ionization_energy.values[
             less_than_chi_0],beta_rad) - chi_0 * beta_electron)
     else:
         radiation_field_correction = np.ones((len(ionization_data),
             len(beta_rad))) * delta_input
     delta = pd.DataFrame(radiation_field_correction,
         columns=np.arange(len(t_rad)), index=ionization_data.index)
     return delta

def test_strucdamping(use_GPU):
    n_inputs = 3
    sig_len = 5

    inputs = np.outer(np.linspace(0.1, 0.9, n_inputs),
                      np.ones(sig_len))[:, :, None]
    targets = np.outer(np.linspace(0.1, 0.9, n_inputs),
                       np.linspace(0, 1, sig_len))[:, :, None]
    inputs = inputs.astype(np.float32)
    targets = targets.astype(np.float32)

    optimizer = HessianFree(CG_iter=100)

    rnn = hf.RNNet(
        shape=[1, 5, 1],
        loss_type=[hf.loss_funcs.SquaredError(),
                   hf.loss_funcs.StructuralDamping(0.1, optimizer=optimizer)],
        debug=True, use_GPU=use_GPU)

    rnn.run_epochs(inputs, targets, optimizer=optimizer,
                   max_epochs=30, print_period=None)

    outputs = rnn.forward(inputs, rnn.W)

    assert rnn.loss.batch_loss(outputs, targets) < 1e-4

 def trainNN(self, imagesTrainSet, labelsTrainSet, etha):
     self.reset_weights()
     trainingSetSize = labelsTrainSet.shape[0];
     j = 0
     while j < 30:
         i = 0
         # print("Round: " + str(j + 1))
         while i < trainingSetSize :
             x = imagesTrainSet[i].ravel()  # Convert 28x28 pixel image into a (784,) vector
             x = np.array([ 0 if val == 0 else 1 for val in x ])
             x_a = np.insert(x, 0, values=1, axis=0)  # Augmented Feature vector
             net_hidd = np.dot(self.w1, x_a)
             y = self.signum(net_hidd)
             y_a = np.insert(y, 0, values=1, axis=0)  # Augmented Feature vector
             
             net_out = np.dot(self.w2, y_a)
             z = self.signum(net_out)
             lab = np.array([ 1 if k == self.labels[i] else 0 for k in range(10) ])
             
             J = z - lab;
             J = np.sum(0.5 * J * J);
             if J < 1 and self.enableWeightDecay:
                 break;
             out_sensitivity = (lab - z) * self.signum_prime(net_out)
             net_hidd_prime = self.signum_prime(net_hidd) 
             hid_sensitivity = np.dot(self.w2.T, out_sensitivity) * np.insert(net_hidd_prime, 0, 1)
             
             grad_hidd_out = etha * np.outer(out_sensitivity, y_a.T)
             grad_in_hidd = etha * np.outer(hid_sensitivity[1:] , x_a.T) 
             
             self.update_weights_bias(grad_in_hidd, grad_hidd_out)
             i += 1
         j += 1
         
     return self.w1, self.w2

    def _create_displacement_matrix(self,
                                    disp_pairs,
                                    site_symmetry,
                                    rot_atom_map):
        rot_disp1s = []
        rot_disp2s = []
        rot_pair12 = []
        rot_pair21 = []
        rot_pair11 = []
        rot_pair22 = []

        for disp_pairs_u1 in disp_pairs:
            for rot_atom_num, ssym in zip(rot_atom_map, site_symmetry):
                ssym_c = similarity_transformation(self._lattice, ssym)
                for (u1, u2) in disp_pairs_u1[rot_atom_num]:
                    Su1 = np.dot(ssym_c, u1)
                    Su2 = np.dot(ssym_c, u2)
                    rot_disp1s.append(Su1)
                    rot_disp2s.append(Su2)
                    rot_pair12.append(np.outer(Su1, Su2).flatten() / 2)
                    rot_pair21.append(np.outer(Su2, Su1).flatten() / 2)
                    rot_pair11.append(np.outer(Su1, Su1).flatten() / 2)
                    rot_pair22.append(np.outer(Su2, Su2).flatten() / 2)
    
        ones = np.ones(len(rot_disp1s)).reshape((-1, 1))

        return np.hstack((ones, rot_disp1s, rot_disp2s,
                          rot_pair12, rot_pair21, rot_pair11, rot_pair22))

    def test_extract_signals_pca(self):
        # Prepare some testing data
        x = np.linspace(0, 2*np.pi, num=int(2*np.pi*50))
        signal0 = np.sin(2*x)
        signal1 = np.sin(3*x)
        in_weights = np.array([[0, 0.25, 0.5, 0.75, 1],
                               [1, 0.75, 0.5, 0.25, 0]])
        features = np.outer(in_weights[0], signal0)
        features += np.outer(in_weights[1], signal1)
        self.assertEqual(features.shape, (5, 314))
        # Extract the signals
        comps, weights = extract_signals_pca(spectra=features,
                                          n_components=2)
        weights = np.swapaxes(weights, 0, 1)
        # Check the results
        new_features = np.outer(weights[0], comps[0])
        new_features += np.outer(weights[1], comps[1])
        self.assertEqual(comps.shape, (2, len(x)))
        self.assertEqual(weights.shape, in_weights.shape)

        plt.plot(x, comps[0], label='c0')
        plt.plot(x, comps[1], label='c1')
        plt.plot(x, features[0], label='f0')
        plt.plot(x, features[1], label='f1')
        plt.plot(x, new_features[0], label='nf0')
        plt.plot(x, new_features[1], label='nf1')
        plt.legend()
        plt.show()

 def TwoSampleTest(self,sample1,sample2,numShuffles=1000,method='vanilla',blockSize=20):
     """
     Compute the p-value associated to the MMD between two samples
     method determines the null approximation procedure:
     ----'vanilla': standard permutation test
     ----'block': block permutation test
     ----'wild': wild bootstrap
     ----'wild-center': wild bootstrap with empirical degeneration
     """
     n1=shape(sample1)[0]
     n2=shape(sample2)[0]
     merged = concatenate( [sample1, sample2], axis=0 )
     merged_len=shape(merged)[0]
     numBlocks = merged_len/blockSize
     K=self.kernel(merged)
     mmd = mean(K[:n1,:n1])+mean(K[n1:,n1:])-2*mean(K[n1:,:n1])
     null_samples = zeros(numShuffles)
     
     if method=='vanilla':
         for i in range(numShuffles):
             pp = permutation(merged_len)
             Kpp = K[pp,:][:,pp]
             null_samples[i] = mean(Kpp[:n1,:n1])+mean(Kpp[n1:,n1:])-2*mean(Kpp[n1:,:n1])
             
     elif method=='block':
         blocks=reshape(arange(merged_len),(numBlocks,blockSize))
         for i in range(numShuffles):
             pb = permutation(numBlocks)
             pp = reshape(blocks[pb],(merged_len))
             Kpp = K[pp,:][:,pp]
             null_samples[i] = mean(Kpp[:n1,:n1])+mean(Kpp[n1:,n1:])-2*mean(Kpp[n1:,:n1])
             
     elif method=='wild' or method=='wild-center':
         if n1!=n2:
             raise ValueError("Wild bootstrap MMD available only on the same sample sizes")
         alpha = exp(-1/float(blockSize))
         coreK = K[:n1,:n1]+K[n1:,n1:]-K[n1:,:n1]-K[:n1,n1:]
         for i in range(numShuffles):
             """
             w is a draw from the Ornstein-Uhlenbeck process
             """
             w = HelperFunctions.generateOU(n=n1,alpha=alpha)
             if method=='wild-center':
                 """
                 empirical degeneration (V_{n,2} in Leucht & Neumann)
                 """
                 w = w - mean(w)
             null_samples[i]=mean(outer(w,w)*coreK)
     elif method=='wild2':
         
         alpha = exp(-1/float(blockSize))
         for i in range(numShuffles):
             wx=HelperFunctions.generateOU(n=n1,alpha=alpha)
             wx = wx - mean(wx)
             wy=HelperFunctions.generateOU(n=n2,alpha=alpha)
             wy = wy - mean(wy)
             null_samples[i]=mean(outer(wx,wx)*K[:n1,:n1])+mean(outer(wy,wy)*K[n1:,n1:])-2*mean(outer(wx,wy)*K[:n1,n1:])
     else:
         raise ValueError("Unknown null approximation method")
     return sum(mmd<null_samples)/float(numShuffles)

def hess(A):
    """Computes the upper Hessenberg form of A using Householder reflectors.
    input:  A, mxn array
    output: Q, orthogonal mxm array
            H, upper Hessenberg
            s.t. Q.dot(H).dot(Q.T) = A
    """
    # similar approach as the householder function.
    # again, not perfectly optimized, but good enough.
    Q = np.eye(A.shape[0]).T
    H = np.array(A, order="C")
    # initialize m and n for convenience
    m, n = H.shape
    # avoid reallocating v in the for loop
    v = np.empty(A.shape[1]-1)
    for k in xrange(n-2):
        # get a slice of the temporary array
        vk = v[k:]
        # fill it with corresponding values from R
        vk[:] = H[k+1:,k]
        # add in the term that makes the reflection work
        vk[0] += copysign(la.norm(vk), vk[0])
        # normalize it so it's an orthogonal transform
        vk /= la.norm(vk)
        # apply projection to H on the left
        H[k+1:,k:] -= 2 * np.outer(vk, vk.dot(H[k+1:,k:]))
        # apply projection to H on the right
        H[:,k+1:] -= 2 * np.outer(H[:,k+1:].dot(vk), vk)
        # Apply it to Q
        Q[k+1:] -= 2 * np.outer(vk, vk.dot(Q[k+1:]))
    return Q, H

    def evolveDraw(self, r):
        """ This can contain triggers for things to be drawn, e.g. the shark."""
        if self.mode == 3 or self.mode == 4 or self.mode == 7:
            # Draw shark
            self.updateShark(r)

            self.axis.scatter(
                self.predatorLocation[:, 0],
                self.predatorLocation[:, 1],
                self.predatorLocation[:, 2],
                color="r",
                s=4 * self.length,
            )

        if self.mode == 1 or self.mode == 2 or self.mode == 3 or self.mode == 6:

            if self.noSphere == True:
                u = np.linspace(0, 2 * np.pi, 100)
                v = np.linspace(0, np.pi, 100)

                self.sphere_x = self.habitatSize * np.outer(np.cos(u), np.sin(v))
                self.sphere_y = self.habitatSize * np.outer(np.sin(u), np.sin(v))
                self.sphere_z = self.habitatSize * np.outer(np.ones(np.size(u)), np.cos(v))

                self.noSphere = False
            self.axis.plot_wireframe(
                self.sphere_x, self.sphere_y, self.sphere_z, rstride=13, cstride=13, color="r", alpha=0.3
            )

    def _get_H(self, debug=False):
        """
        returns H_t as defined in algorithm 2
        
        Reference:
        https://en.wikipedia.org/wiki/Limited-memory_BFGS
        http://www.ccms.or.kr/data/pdfpaper/jcms21_1/21_1_117.pdf
        https://homes.cs.washington.edu/~galen/files/quasi-newton-notes.pdf
        """
        I = np.identity(len(self.w))
        
        if min(len(self.s), len(self.y)) == 0:
                print "Warning: No second order information used!"
                return I
            
        assert len(self.s) > 0, "s cannot be empty."
        assert len(self.s) == len(self.y), "s and y must have same length"
        assert self.s[0].shape == self.y[0].shape, \
            "s and y must have same shape"
        assert abs(self.y[-1]).sum() != 0, "latest y entry cannot be 0!"
        assert 1/np.inner(self.y[-1], self.s[-1]) != 0, "!"

        I = np.identity(len(self.s[0]))
        H = np.dot((np.inner(self.s[-1], self.y[-1]) / np.inner(self.y[-1],
                   self.y[-1])), I)

        for (s_j, y_j) in itertools.izip(self.s, self.y):
            rho = 1.0/np.inner(y_j, s_j)
            V = I - np.multiply(rho, np.outer(s_j, y_j))
            H = (V).dot(H).dot(V.T)
            H += np.multiply(rho, np.outer(s_j, s_j))

        return H

    def learn(self, docs, alpha=0.1, tau=5):
        index = numpy.arange(len(docs))
        numpy.random.shuffle(index)
        for i in index:
            doc = docs[i]
            pre_s = [numpy.zeros(self.K)]
            pre_w = [0] # <s>
            for w in doc:
                s = 1 / (numpy.exp(- numpy.dot(self.W, pre_s[-1]) - self.U[:, pre_w[-1]]) + 1)
                z = numpy.dot(self.V, s)
                y = numpy.exp(z - z.max())
                y = y / y.sum()

                # calculate errors
                y[w] -= 1  # -e0
                eh = [numpy.dot(y, self.V) * s * (s - 1)] # eh[t]
                for t in xrange(min(tau, len(pre_s)-1)):
                    st = pre_s[-1-t]
                    eh.append(numpy.dot(eh[-1], self.W) * st * (1 - st))

                # update parameters
                pre_w.append(w)
                pre_s.append(s)
                self.V -= numpy.outer(y, s * alpha)
                for t in xrange(len(eh)):
                    self.U[:, pre_w[-1-t]] += eh[t] * alpha
                    self.W += numpy.outer(pre_s[-2-t], eh[t]) * alpha

def Dgrid_zone(zone,Px):
    """ Distance point to zone

    A zone is a quadrilateral zone.

    Parameters
    ----------

    zone : dictionnary
           xmin xmax Nx
           ymin ymax Ny

    Px : np.array
         point

    Build the distance matrix between Tx and points in the zone

    Notes
    -----

    use broadcasting instead

    """

    rx = np.linspace(zone['xmin'],zone['xmax'],zone['Nx'])
    ry = np.linspace(zone['ymin'],zone['ymax'],zone['Ny'])

    R_x = np.outer(np.ones(len(ry)),rx)
    R_y = np.outer(ry,np.ones(len(rx)))

    Dx = R_x - Px[0]
    Dy = R_y - Px[1]
    D = np.sqrt(Dx*Dx+Dy*Dy)
    return (D)

    def backProp(self,node,error=None):
        # Clear nodes
        node.fprop = False
        errorCur = node.probs - make_onehot(node.label,len(self.bs))
        self.dWs += np.outer(errorCur, node.hActs1)
        self.dbs += errorCur

        errorCur = errorCur.dot(self.Ws)
        if error is not None:
            errorCur += error

        if node.isLeaf == True:
            self.dL[node.word] += errorCur
            return

        errorCur = errorCur*self.df(node.hActs1)
        LR = np.hstack([node.left.hActs1, node.right.hActs1])
        self.dW += np.outer(errorCur,LR)
        self.db += errorCur

        S = np.zeros(len(LR))
        for i in range(len(self.V)):
            self.dV[i] += errorCur[i]*np.outer(LR,LR)
            S += (self.V[i]+self.V[i].T).dot(LR)*errorCur[i]
        
        errorDown = errorCur.dot(self.W) + S        
        self.backProp(node.left,errorDown[:self.wvecDim])
        self.backProp(node.right,errorDown[self.wvecDim:])

    def bptt(self, x, y):
        # forward pass
        # save all hidden states and outpts at each time step
        # because need them during back propagation.
        # add one additional 0s as the initial hidden state
        s = np.zeros((len(x) + 1, self.hidden_dim))
        s[-1] = np.zeros(self.hidden_dim)
        o = np.zeros((len(x), self.word_dim))
        for t in range(len(x)):
            # Note that we are indxing U by x[t]. This is the same as
            # multiplying U with a one-hot vector.
            s[t] = np.tanh(self.U[:,x[t]] + self.W.dot(s[t-1]))
            o[t] = softmax(self.V.dot(s[t]))

        # backward pass
        dLdU = np.zeros(self.U.shape)
        dLdV = np.zeros(self.V.shape)
        dLdW = np.zeros(self.W.shape)
        # dLdy = o - y
        dLdy = o
        dLdy[np.arange(len(y)), y] -= 1.
        # For each output backwards...
        for t in np.arange(len(y))[::-1]:
            dLdV += np.outer(dLdy[t], s[t].T)
            # Initial delta calculation
            dLdz = self.V.T.dot(dLdy[t]) * (1 - (s[t] ** 2))
            # Backpropagation through time (for at most self.bptt_truncate steps)
            for bptt_step in np.arange(max(0, t-self.bptt_truncate), t+1)[::-1]:
                # print "Backpropagation step t=%d bptt step=%d " % (t, bptt_step)
                dLdW += np.outer(dLdz, s[bptt_step-1])
                dLdU[:,x[bptt_step]] += dLdz
                # Update delta for next step
                dLdz = self.W.T.dot(dLdz) * (1 - s[bptt_step-1] ** 2)
        return [dLdU, dLdV, dLdW]