def backward_pass(self, a1L, a1R, a2L, a2LR, a2R, a3, z1Lb, z1LRb, z1Rb, z2b, xLb, xRb, t):
	 	
		# Third Layer
		if self.k == 2 :
			r3=	-t*self.sigmoid(-t*a3)
		else :
			r3 = a3 - t

		grad3=sp.dot(r3,z2b.T)
		
		# Second Layer
		r3w3T = sp.dot(self.w3[:,:-1].T, r3)

		r2L=r3w3T*a2LR*self.sigmoid(a2R)*self.divsigmoid(a2L)
		r2R=r3w3T*a2LR*self.sigmoid(a2L)*self.divsigmoid(a2R)
		r2LR=r3w3T*self.sigmoid(a2L)*self.sigmoid(a2R)
		
		grad2L = sp.dot(r2L, z1Lb.T)
		grad2LR = sp.dot(r2LR, z1LRb.T)
		grad2R = sp.dot(r2R, z1Rb.T)
		
		# First Layer
		r1L = sp.power(1.0/sp.cosh(a1L),2)*(sp.dot(self.w2l[:,:-1].T, r2L)+sp.dot(self.w2lr[:,:self.H1].T, r2LR))
		r1R = sp.power(1.0/sp.cosh(a1R),2)*(sp.dot(self.w2r[:,:-1].T, r2R)+sp.dot(self.w2lr[:,self.H1:-1].T, r2LR))
		
		grad1L = sp.dot(r1L, xLb.T)
		grad1R = sp.dot(r1R, xRb.T)
		
		
		return grad3, grad2L, grad2LR, grad2R, grad1L, grad1R

 def Function(self, x, param):
     Sp, alpha, beta, Ta  = param
     S0 = self.CalcS0(Ta)
     x2 = (scipy.greater_equal(x, Ta) * (x - Ta))
     #y  = Sp * numpy.abs(scipy.power((scipy.e / (alpha*beta)), alpha)) * numpy.abs(scipy.power(x2, alpha)) * scipy.exp(-x2/beta) + S0
     y  = Sp * numpy.abs(scipy.power((scipy.e / (alpha*beta)), alpha) * scipy.power(x2, alpha) * scipy.exp(-x2/beta)) + S0
     return y

 def DiffFunction(self, x, param):
     Sp, alpha, beta, Ta  = param
     if x < Ta:
         return 0
     x2 = x - Ta
     y = Sp * scipy.power((scipy.e / (alpha*beta)), alpha) * (alpha * scipy.power(x2, alpha - 1) * scipy.exp(-x2/beta) - scipy.power(x2, alpha) * scipy.exp(-x2/beta) / beta)
     return y

def setTransitionsFork(dimension):
    """ setting transitions in the state space """

    space_size = np.power(2, dimension)
    transition = np.ndarray(shape=(space_size, space_size), dtype=bool)
    transition.fill(False)

    state1 = [0 for i in range(dimension)]
    state2 = state1[:]
    for i in range(dimension-1):
        state2[1+i] = 1
        state1[0] = 1
        state2[0] = 1
        transition[st2Ind(state1)][st2Ind(state2)] = True  # forward transition
        transition[st2Ind(state2)][st2Ind(state1)] = True  # backward transition
        state1 = state2[:]

    state1 = [0 for i in range(dimension)]
    state2 = state1[:]
    for i in range(dimension-1):
        state2[dimension-i-1] = 1
        state1[0] = 1
        state2[0] = 1
        transition[st2Ind(state1)][st2Ind(state2)] = True  # forward transition
        transition[st2Ind(state2)][st2Ind(state1)] = True  # backward transition
        state1 = state2[:]

    transition[0][np.power(2, (dimension-1))] = True
    transition[np.power(2, (dimension-1))][0] = True

    print 'Fork transitions'
    printTransitions(transition, dimension)

    return transition

 def calc_volume(self):
     """calculate the volume over which the compound can move. We have
        Cavg = mass/volume
     """
     return sp.sum((sp.power(self.grid_edge[1:],2) - 
                         sp.power(self.grid_edge[:-1],2)) 
                     ) * sp.pi 

def r_ion_neutral(s,t,Ni,Nn,Ti,Tn):
    """ This will calculate resonant ion - neutral reactions collision frequencies. See
    table 4.5 in Schunk and Nagy.
    Inputs
    s - Ion name string
    t - neutral name string
    Ni - Ion density cm^-3
    Nn - Neutral density cm^-3
    Ti - Ion tempreture K
    Tn - Neutral tempreture K
    Outputs
    nu_ineu - collision frequency s^-1
    """
    Tr = (Ti+Tn)*0.5
    sp1 = (s,t)
    # from Schunk and Nagy table 4.5
    nudict={('H+','H'):[2.65e-10,0.083],('He+','He'):[8.73e-11,0.093], ('N+','N'):[3.84e-11,0.063],
            ('O+','O'):[3.67e-11,0.064], ('N2+','N'):[5.14e-11,0.069], ('O2+','O2'):[2.59e-11,0.073],
            ('H+','O'):[6.61e-11,0.047],('O+','H'):[4.63e-12,0.],('CO+','CO'):[3.42e-11,0.085],
            ('CO2+','CO'):[2.85e-11,0.083]}
    A = nudict[sp1][0]
    B = nudict[sp1][1]
    if sp1==('O+','H'):

        nu_ineu = A*Nn*sp.power(Ti/16.+Tn,.5)
    elif sp1==('H+','O'):
        nu_ineu = A*Nn*sp.power(Ti,.5)*(1-B*sp.log10(Ti))**2
    else:
        nu_ineu = A*Nn*sp.power(Tr,.5)*(1-B*sp.log10(Tr))**2
    return nu_ineu

 def __init__(self, *args, **kwargs):
     MultiModalFunction.__init__(self, *args, **kwargs)
     self._opts = (rand((self.numPeaks, self.xdim)) - 0.5) * 9.8
     self._opts[0] = (rand(self.xdim) - 0.5) * 8
     alphas = [power(self.maxCond, 2 * i / float(self.numPeaks - 2)) for i in range(self.numPeaks - 1)]
     shuffle(alphas)
     self._covs = [generateDiags(alpha, self.xdim, shuffled=True) / power(alpha, 0.25) for alpha in [self.optCond] + alphas]
     self._R = orth(rand(self.xdim, self.xdim))
     self._ws = [10] + [1.1 + 8 * i / float(self.numPeaks - 2) for i in range(self.numPeaks - 1)]

def gvf(x, Sp, alpha, beta, Ta, S0):
    global scipy, numpy
    y = (
        Sp
        * numpy.abs(scipy.power((scipy.e / (alpha * beta)), alpha))
        * numpy.abs(scipy.power((x - Ta), alpha))
        * scipy.exp(-(x - Ta) / beta)
        + S0
    )
    return y

def freqs_b_by_freqs(numFilters, lowFreq, highFreq, c, d, order = 1):
    # Find the center frequencies of the filters from the begin and end frequencies
    EarQ = c
    minBW = d
    vec = scipy.arange(numFilters, 0, -1)
    freqs = -(EarQ*minBW) + scipy.exp(vec*(-scipy.log(highFreq + EarQ*minBW) + scipy.log(lowFreq + EarQ*minBW))/numFilters) * (highFreq + EarQ*minBW);

    ERB = scipy.power((scipy.power((freqs/EarQ), order) + scipy.power(minBW, order)), (1.0 / order))
    B = 1.019 * 2.0 * scipy.pi * ERB

    return (freqs, B)

def plotdata(ionofile_in,ionofile_fit,madfile,time1):


    fig1,axmat =plt.subplots(2,2,facecolor='w',figsize=(10,10))
    axvec = axmat.flatten()
    paramlist = ['ne','te','ti','vo']
    paramlisti = ['Ne','Te','Ti','Vi']
    paramlistiname = ['$N_e$','$T_e$','$T_i$','$V_i$']
    paramunit = ['$m^{-3}$','$^\circ$ K','$^\circ$ K','m/s']
    boundlist = [[0.,7e11],[500.,3200.],[500.,2500.],[-500.,500.]]
    IonoF = IonoContainer.readh5(ionofile_fit)
    IonoI = IonoContainer.readh5(ionofile_in)
    gfit = GeoData(readIono,[IonoF,'spherical'])
    ginp = GeoData(readIono,[IonoI,'spherical'])
    data1 = GeoData(readMad_hdf5,[madfile,['nel','te','ti','vo','dnel','dte','dti','dvo']])
    data1.data['ne']=sp.power(10.,data1.data['nel'])
    data1.data['dne']=sp.power(10.,data1.data['dnel'])
    
    t1,t2 = data1.timelisting()[340]
    handlist = []
    for inum,iax in enumerate(axvec):
        ploth = rangevsparam(data1,data1.dataloc[0,1:],time1,gkey=paramlist[inum],fig=fig1,ax=iax,it=False)
        handlist.append(ploth[0])
        ploth = rangevsparam(ginp,ginp.dataloc[0,1:],0,gkey=paramlisti[inum],fig=fig1,ax=iax,it=False)
        handlist.append(ploth[0])
        ploth = rangevsparam(gfit,gfit.dataloc[0,1:],0,gkey=paramlisti[inum],fig=fig1,ax=iax,it=False)
        handlist.append(ploth[0])
        iax.set_xlim(boundlist[inum])
        iax.set_ylabel('Altitude in km')
        iax.set_xlabel(paramlistiname[inum]+' in '+paramunit[inum])
    # with error bars
    plt.tight_layout()
    fig1.suptitle('Comparison Without Error Bars\nPFISR Data Times: {0} to {1}'.format(t1,t2))
    plt.subplots_adjust(top=0.9)
    plt.figlegend( handlist[:3], ['PFISR', 'SimISR Input','SimISR Fit'], loc = 'lower center', ncol=5, labelspacing=0. )
    fig2,axmat2 =plt.subplots(2,2,facecolor='w',figsize=(10,10))
    axvec2 = axmat2.flatten()
    handlist2 = []
    for inum,iax in enumerate(axvec2):
        ploth = rangevsparam(data1,data1.dataloc[0,1:],time1,gkey=paramlist[inum],gkeyerr='d'+paramlist[inum],fig=fig2,ax=iax,it=False)
        handlist2.append(ploth[0])
        ploth = rangevsparam(ginp,ginp.dataloc[0,1:],0,gkey=paramlisti[inum],fig=fig2,ax=iax,it=False)
        handlist2.append(ploth[0])
        ploth = rangevsparam(gfit,gfit.dataloc[0,1:],0,gkey=paramlisti[inum],gkeyerr='n'+paramlisti[inum],fig=fig2,ax=iax,it=False)
        handlist2.append(ploth[0])
        iax.set_xlim(boundlist[inum])
        iax.set_ylabel('Altitude in km')
        iax.set_xlabel(paramlistiname[inum]+' in '+paramunit[inum])
    plt.tight_layout()
    fig2.suptitle('Comparison With Error Bars\nPFISR Data Times: {0} to {1}'.format(t1,t2))
    plt.subplots_adjust(top=0.9)
    plt.figlegend( handlist2[:3], ['PFISR', 'SimISR Input','SimISR Fit'], loc = 'lower center', ncol=5, labelspacing=0. )
    return (fig1,axvec,handlist,fig2,axvec2,handlist2)

 def calc_mass(self, conc_r, split=0):
     """calculate the mass of component present given value in cell center
     This is given by 2 \pi int_r1^r2 C(r)r dr
     
     conc_r: concentration in self.grid
     """
     if split == 1:
         grid = self.grid_edge_sp1
     elif split == 2:
         grid = self.grid_edge_sp2
     else:
         grid = self.grid_edge
     return sp.sum(conc_r * (sp.power(grid[1:], 2) - 
                             sp.power(grid[:-1], 2)) 
                     ) * sp.pi 

def f_L(q,a_1,a_2,a_3):
    '''integrand of the static geometrical tensor

    Parameters
    ----------
    'q' = free variable for the geometrical tensor
    'a_1,a_2,a_3' = three axes of the ellipsoid (in nm)

    Returns
    -------
    'L' = integrand of the static geometrical tensor'''

    L = sp.power(q+a_1**2,-1.5)*sp.power(q+a_2**2,-0.5)*sp.power(q+a_3**2,-0.5)

    return L

	def setNewScale(self, state):
		Names = ( 'Y_XScale', 'XLogScale', 'YLogScale', 'LogScale' )
		Types = {'c' : 0, 's' : 1, 'r' : 2} 
		senderName = self.sender().objectName()
		t, Type = senderName[0], Types[senderName[0]]
		data = self.getData(Type)
		Scale = data.Scale()
		ui_obj = self.findUi([t + i for i in Names])
		if senderName[1:] == Names[0]:
			#ui_obj = getattr(self.ui, t + "LogScale")
			if state:
				Scale[1] = 2
				data[:,1] = data[:,1] / data[:,0]
			else:
				Scale[1] = 0
				data[:,1] = data[:,1] * data[:,0]
			ui_obj[3].setEnabled(not ui_obj[0].isChecked())
		else:
			index = bool(senderName[1] != "X")
			#ui_obj = getattr(self.ui, t + Names[0])
			if Scale[index] != state:
				if state == 1:
					data[:,index] = sp.log10(data[:,index])
				else:
					data[:,index] = sp.power(10.,data[:,index])
				Scale[index] = int(state)
				ui_obj[0].setEnabled(not (ui_obj[1].isChecked() or ui_obj[2].isChecked()))
		self.updateData(array = Array(data, Type = Type, scale = Scale))

    def filter(self, mode='soft'):

        if self.level > self.max_dec_level():
            clevel = self.max_dec_level()
        else:
            clevel = self.level

        # decompose
        coeffs = pywt.wavedec(self.sig, pywt.Wavelet(self.wt), \
                              mode=self.mode, \
                              level=clevel)

        # threshold evaluation
        th = sqrt(2 * log(len(self.sig)) * power(self.sigma, 2))

        # thresholding
        for (i, cAD) in enumerate(coeffs):
            if mode == 'soft':
                coeffs[i] = pywt.thresholding.soft(cAD, th)
            elif mode == 'hard':
                coeffs[i] = pywt.thresholding.hard(cAD, th)

        # reconstruct
        rec_sig = pywt.waverec(coeffs, pywt.Wavelet(self.wt), mode=self.mode)
        if len(rec_sig) == (len(self.sig) + 1):
            self.sig = rec_sig[:-1]

  def generateGaborMotherWavelet(self):
    pitch = 440.0
    sigma = 6.
    NL = 48
    NU = 39
    print 'sampling rate:', self.fs, 'Hz'
    fs = float(self.fs)
    self.sample_duration = 10.
    #asigma = 0.3
    limit_t = 0.1
    #zurashi = 1.

    #NS = NL + NU + 1
    f = sp.array([2**(i/12.) for i in range(NL+NU+1)]) * pitch*2**(-NL/12.)
    f = f[:, sp.newaxis]
    sigmao = sigma*10**(-3)*sp.sqrt(fs/f)
    t = sp.arange(-limit_t, limit_t+1/fs, 1/fs)

    inv_sigmao = sp.power(sigmao, -1)
    inv_sigmao_t = inv_sigmao * t
    t_inv_sigmao2 = sp.multiply(inv_sigmao_t, inv_sigmao_t)
    omega_t = 2*sp.pi*f*t
    gabor = (1/sp.sqrt(2*sp.pi))
    gabor = sp.multiply(gabor, sp.diag(inv_sigmao))
    exps = -0.5*t_inv_sigmao2+sp.sqrt(-1)*omega_t
    self.gabor = gabor*sp.exp(exps)

def minowski(h1, h2, p = 2): # 46..45..14,11..43..44 / 45 us for p=int(-inf..-24..-1,1..24..inf) / float @array, +20 us @list \w 100 bins
    r"""
    Minowski distance.
    
    With :math:`p=2` equal to the Euclidean distance, with :math:`p=1` equal to the Manhattan distance,
    and the Chebyshev distance implementation represents the case of :math:`p=\pm inf`.
    
    The Minowksi distance between two histograms :math:`H` and :math:`H'` of size :math:`m` is
    defined as:
    
    .. math::
    
        d_p(H, H') = \left(\sum_{m=1}^M|H_m - H'_m|^p  
            \right)^{\frac{1}{p}}

    *Attributes:*
    
    - a real metric
    
    *Attributes for normalized histograms:*
    
    - :math:`d(H, H')\in[0, \sqrt[p]{2}]`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-normalized histograms:*
    
    - :math:`d(H, H')\in[0, \infty)`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-equal histograms:*
    
    - not applicable
    
    Parameters
    ----------
    h1 : sequence
        The first histogram.
    h2 : sequence
        The second histogram.
    p : float
        The :math:`p` value in the Minowksi distance formula.
    
    Returns
    -------
    minowski : float
        Minowski distance.
    
    Raises
    ------
    ValueError
        If ``p`` is zero.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    if 0 == p: raise ValueError('p can not be zero')
    elif int == type(p):
        if p > 0 and p < 25: return __minowski_low_positive_integer_p(h1, h2, p)
        elif p < 0 and p > -25: return __minowski_low_negative_integer_p(h1, h2, p)
    return math.pow(scipy.sum(scipy.power(scipy.absolute(h1 - h2), p)), 1./p)

 def addReward(self):
     """ A filtered mapping towards performAction of the underlying environment. """
     # by default, the cumulative reward is just the sum over the episode
     if self.discount:
         self.cumreward += power(self.discount, self.samples) * self.getReward()
     else:
         self.cumreward += self.getReward()

def gauss_legendre(n):
    k = sp.arange(1.0, n)
    a_band = sp.zeros((2, n))
    a_band[1, 0:n - 1] = k / sp.sqrt(4 * k * k - 1)
    x, V = sp.linalg.eig_banded(a_band, lower=True)
    w = 2 * sp.real(sp.power(V[0, :], 2))
    return x, w

def asymmetrify(x, beta=0.2):
    res = x.copy()
    dim = len(x)
    for i, xi in enumerate(x):
        if xi > 0:
            res[i] = power(xi, 1+beta*i/(dim-1.)*sqrt(xi))
    return res

def multivariateNormalPdf(z, x, sigma):
    """ The pdf of a multivariate normal distribution (not in scipy).
    The sample z and the mean x should be 1-dim-arrays, and sigma a square 2-dim-array. """
    assert len(z.shape) == 1 and len(x.shape) == 1 and len(x) == len(z) and sigma.shape == (len(x), len(z))
    tmp = -0.5 * dot(dot((z - x), inv(sigma)), (z - x))
    res = (1. / power(2.0 * pi, len(z) / 2.)) * (1. / sqrt(det(sigma))) * exp(tmp)
    return res