def relative_bin_deviation(h1, h2): # 79 us @array, 104 us @list \w 100 bins
    r"""
    Calculate the bin-wise deviation between two histograms.
    
    The relative bin deviation between two histograms :math:`H` and :math:`H'` of size
    :math:`m` is defined as:
    
    .. math::
    
        d_{rbd}(H, H') = \sum_{m=1}^M
            \frac{
                \sqrt{(H_m - H'_m)^2}
              }{
                \frac{1}{2}
                \left(
                    \sqrt{H_m^2} +
                    \sqrt{{H'}_m^2}
                \right)
              }
    
    *Attributes:*

    - semimetric (triangle equation satisfied?)
    
    *Attributes for normalized histograms:*

    - :math:`d(H, H')\in[0, \infty)`
    - :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, same bins as ``h1``.
    
    Returns
    -------
    relative_bin_deviation : float
        Relative bin deviation between the two histograms.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    numerator = scipy.sqrt(scipy.square(h1 - h2))
    denominator = (scipy.sqrt(scipy.square(h1)) + scipy.sqrt(scipy.square(h2))) / 2.
    old_err_state = scipy.seterr(invalid='ignore') # divide through zero only occurs when the bin is zero in both histograms, in which case the division is 0/0 and leads to (and should lead to) 0
    result = numerator / denominator
    scipy.seterr(**old_err_state)
    result[scipy.isnan(result)] = 0 # faster than scipy.nan_to_num, which checks for +inf and -inf also
    return scipy.sum(result)

    def run(self, phase=None, throats=None):
        logger.warning('This algorithm can take some time...')
        conduit_lengths = sp.sum(misc.conduit_lengths(network=self._net,
                                 mode='centroid'), axis=1)
        graph = self._net.create_adjacency_matrix(data=conduit_lengths,
                                                  sprsfmt='csr')

        if phase is not None:
            self._phase = phase
            if 'throat.occupancy' in self._phase.props():
                temp = conduit_lengths*(self._phase['throat.occupancy'] == 1)
                graph = self._net.create_adjacency_matrix(data=temp,
                                                          sprsfmt='csr',
                                                          prop='temp')
        path = spgr.shortest_path(csgraph=graph, method='D', directed=False)

        Px = sp.array(self._net['pore.coords'][:, 0], ndmin=2)
        Py = sp.array(self._net['pore.coords'][:, 1], ndmin=2)
        Pz = sp.array(self._net['pore.coords'][:, 2], ndmin=2)

        Cx = sp.square(Px.T - Px)
        Cy = sp.square(Py.T - Py)
        Cz = sp.square(Pz.T - Pz)
        Ds = sp.sqrt(Cx + Cy + Cz)

        temp = path / Ds

        temp[sp.isnan(temp)] = 0
        temp[sp.isinf(temp)] = 0

        return temp

def csr3(complex_n):
    ang = sp.angle(complex_n) # sp.arctan(a.imag/a.real) why it does not work?!?!
    r = sp.sqrt(sp.square(complex_n.real)+sp.square(complex_n.imag))
    if (sp.sin(ang/2)>=0): #sin>0
        return sp.sqrt(r)*(complex(sp.cos(ang/2),sp.sin(ang/2)))
    else:
        return sp.sqrt(r)*(complex(sp.cos((ang/2)+sp.pi),sp.sin((ang/2)+sp.pi)))

    def run(self,phase=None):
        r'''
        '''
        logger.warning('This algorithm can take some time...')
        graph = self._net.create_adjacency_matrix(data=self._net['throat.length'],sprsfmt='csr')

        if phase is not None:
            self._phase = phase
            if 'throat.occupancy' in self._phase.props():
                temp = self._net['throat.length']*(self._phase['throat.occupancy']==1)
                graph = self._net.create_adjacency_matrix(data=temp,sprsfmt='csr',prop='temp')

        #self._net.tic()
        path = spgr.shortest_path(csgraph = graph, method='D', directed = False)
        #self._net.toc()

        Px = sp.array(self._net['pore.coords'][:,0],ndmin=2)
        Py = sp.array(self._net['pore.coords'][:,1],ndmin=2)
        Pz = sp.array(self._net['pore.coords'][:,2],ndmin=2)

        Cx = sp.square(Px.T - Px)
        Cy = sp.square(Py.T - Py)
        Cz = sp.square(Pz.T - Pz)
        Ds = sp.sqrt(Cx + Cy + Cz)

        temp = path/Ds
        #temp = path

        temp[sp.isnan(temp)] = 0
        temp[sp.isinf(temp)] = 0

        return temp

def relative_deviation(h1, h2): # 18 us @array, 42 us @list \w 100 bins
    r"""
    Calculate the deviation between two histograms.
    
    The relative deviation between two histograms :math:`H` and :math:`H'` of size :math:`m` is
    defined as:
    
    .. math::
    
        d_{rd}(H, H') =
            \frac{
                \sqrt{\sum_{m=1}^M(H_m - H'_m)^2}
              }{
                \frac{1}{2}
                \left(
                    \sqrt{\sum_{m=1}^M H_m^2} +
                    \sqrt{\sum_{m=1}^M {H'}_m^2}
                \right)
              }
    
    *Attributes:*

    - semimetric (triangle equation satisfied?)
    
    *Attributes for normalized histograms:*

    - :math:`d(H, H')\in[0, \sqrt{2}]`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-normalized histograms:*

    - :math:`d(H, H')\in[0, 2]`
    - :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, same bins as ``h1``.
    
    Returns
    -------
    relative_deviation : float
        Relative deviation between the two histograms.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    numerator = math.sqrt(scipy.sum(scipy.square(h1 - h2)))
    denominator = (math.sqrt(scipy.sum(scipy.square(h1))) + math.sqrt(scipy.sum(scipy.square(h2)))) / 2.
    return numerator / denominator

def tsc(parameters,positions,values):
    
    values_tsc = sp.zeros((parameters.Ng,parameters.Ng,parameters.Ng))
    cellsize = parameters.boxsize/parameters.Ng


    for position,pvalue in zip(positions,values):
        
        position = sp.array(position)
        
        position_cellunits = position/cellsize



        cell_indices = sp.floor(position_cellunits)
        leftcell_indices = cell_indices - 1
        rightcell_indices = cell_indices + 1


        cell_position = cell_indices + 0.5
        leftcell_position = leftcell_indices + 0.5
        rightcell_position = rightcell_indices + 0.5
        
        particle_cell_distances = sp.absolute(position_cellunits - cell_position)
        particle_leftcell_distances = sp.absolute(position_cellunits - leftcell_position)        
        particle_rightcell_distances = sp.absolute(position_cellunits - rightcell_position)

        weights_cell = 0.75 - sp.square(particle_cell_distances)
        weights_leftcell = 0.5*sp.square(1.5 - particle_leftcell_distances)
        weights_rightcell = 0.5*sp.square(1.5 - particle_rightcell_distances)



        if periodic_boundaries:
            cell_indices = sp.mod(cell_indices,parameters.Ng)
            leftcell_indices = sp.mod(leftcell_indices,parameters.Ng)
            rightcell_indices = sp.mod(rightcell_indices,parameters.Ng)


        indices_x, weights_x = [cell_indices[0],leftcell_indices[0],rightcell_indices[0]],\
                                [weights_cell[0],weights_leftcell[0],weights_rightcell[0]]
        indices_y, weights_y = [cell_indices[1],leftcell_indices[1],rightcell_indices[1]],\
                                [weights_cell[1],weights_leftcell[1],weights_rightcell[1]]
        indices_z, weights_z = [cell_indices[2],leftcell_indices[2],rightcell_indices[2]],\
                                [weights_cell[2],weights_leftcell[2],weights_rightcell[2]]

        for index_x,weight_x in zip(indices_x, weights_x):
            for index_y,weight_y in zip(indices_y, weights_y):
                for index_z,weight_z in zip(indices_z, weights_z):
                    values_tsc[index_x][index_y][index_z] += pvalue*weight_x*weight_y*weight_z/cellsize**3                                    


                
    return values_tsc 

def cosine(h1, h2): # 17 us @array, 42 us @list \w 100 bins
    r"""
    Cosine simmilarity.
    
    Compute the angle between the two histograms in vector space irrespective of their
    length. The cosine similarity between two histograms :math:`H` and :math:`H'` of size
    :math:`m` is defined as:
    
    .. math::
    
        d_{\cos}(H, H') = \cos\alpha = \frac{H * H'}{\|H\| \|H'\|} = \frac{\sum_{m=1}^M H_m*H'_m}{\sqrt{\sum_{m=1}^M H_m^2} * \sqrt{\sum_{m=1}^M {H'}_m^2}}
    
    
    *Attributes:*

    - not a metric, a similarity
    
    *Attributes for normalized histograms:*

    - :math:`d(H, H')\in[0, 1]`
    - :math:`d(H, H) = 1`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-normalized histograms:*

    - :math:`d(H, H')\in[-1, 1]`
    - :math:`d(H, H) = 1`
    - :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, same bins as ``h1``.
    
    Returns
    -------
    cosine : float
        Cosine simmilarity.
        
    Notes
    -----
    The resulting similarity ranges from -1 meaning exactly opposite, to 1 meaning
    exactly the same, with 0 usually indicating independence, and in-between values
    indicating intermediate similarity or dissimilarity.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    return scipy.sum(h1 * h2) / math.sqrt(scipy.sum(scipy.square(h1)) * scipy.sum(scipy.square(h2)))

def relative_bin_deviation(h1, h2):  # 79 us @array, 104 us @list \w 100 bins
    """
    Calculate the bin-wise deviation between two histograms.
    The relative bin deviation between two histograms \f$H\f$ and \f$H'\f$ of size
    \f$m\f$ is defined as
    \f[
        d_{rbd}(H, H') = \sum_{m=1}^M
            \frac{
                \sqrt{(H_m - H'_m)^2}
              }{
                \frac{1}{2}
                \left(
                    \sqrt{H_m^2} +
                    \sqrt{{H'}_m^2}
                \right)
              }
    \f]
    
    Attributes:
    - semimetric (triangle equation satisfied?)
    
    Attributes for normalized histograms:
    - \f$d(H, H')\in[0, \infty)\f$
    - \f$d(H, H) = 0\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-normalized histograms:
    - \f$d(H, H')\in[0, \infty)\f$
    - \f$d(H, H) = 0\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-equal histograms:
    - not applicable 
    
    @param h1 the first histogram
    @type h1 array-like sequence
    @param h2 the second histogram, same bins as h1
    @type h2 array-like sequence
    
    @return relative bin deviation
    @rtype float
    """
    h1, h2 = __prepare_histogram(h1, h2)
    numerator = scipy.sqrt(scipy.square(h1 - h2))
    denominator = (scipy.sqrt(scipy.square(h1)) + scipy.sqrt(scipy.square(h2))) / 2.0
    old_err_state = scipy.seterr(
        invalid="ignore"
    )  # divide through zero only occurs when the bin is zero in both histograms, in which case the division is 0/0 and leads to (and should lead to) 0
    result = numerator / denominator
    scipy.seterr(**old_err_state)
    result[scipy.isnan(result)] = 0  # faster than scipy.nan_to_num, which checks for +inf and -inf also
    return scipy.sum(result)

		def __fit(index1, index2):
			from scipy import stats, sqrt, square
			
			# do the fit
			(cijFitted,intercept,r,tt,stderr) = stats.linregress(strain[:,index2-1],stress[:,index1-1])

			if (S.__version__ < '0.7.0'):
				# correct for scipy weirdness - see http://www.scipy.org/scipy/scipy/ticket/8
				# This was fixed before 0.7.0 release. Maybe in some versions of 0.6.x too - 
				# will report huge errors if the check is wrong
				stderr = S.sqrt((numsteps * stderr**2)/(numsteps-2))
				error  = stderr/sqrt(sum(square(strain[:,index2-1])))
			else:
				# Work out the error ourselves as I cannot get it from
				# stderr and this has been checked with gnuplot's fitter
				fit_str = ((strain[:,index2-1] * cijFitted) + intercept)
				error = sqrt((sum(square(stress[:,index1-1] - fit_str)) / \
				             (numsteps-2))/(sum(square(strain[:,index2-1]))))
			
			# print info about the fit
			print '\n'
			print     'Cij (gradient)          :    ', cijFitted
			print     'Error in Cij            :    ', error
			print     'Intercept               :    ', intercept
			if abs(r) > 0.9:
				print 'Correlation coefficient :    ',r
			else:
				print 'Correlation coefficient :    ',r, '     <----- WARNING'
			
			# if using graphics, add a subplot
			if options.graphics:
					
				# position this plot in a 6x6 grid
				sp = P.subplot(6,6,6*(index1-1)+index2)
				sp.set_axis_on()
				
				# change the labels on the axes
				xlabels = sp.get_xticklabels()
				P.setp(xlabels,'rotation',90,fontsize=7)
				ylabels = sp.get_yticklabels()
				P.setp(ylabels,fontsize=7)
			
				# colour the plot depending on the strain pattern
				sp.set_axis_bgcolor(colourDict[patt])

				# plot the data
				P.plot([strain[0,index2-1],strain[numsteps-1,index2-1]],[cijFitted*strain[0,index2-1]+intercept,cijFitted*strain[numsteps-1,index2-1]+intercept])
				P.plot(strain[:,index2-1],stress[:,index1-1],'ro')
			
			return cijFitted, error

def conduit_lengths(network, throats=None, mode='pore'):
    r"""
    Return the respective lengths of the conduit components defined by the throat
    conns P1 - T - P2

    Notes
    -----
    mode = 'pore' - uses pore coordinates
    mode = 'centroid' uses pore and throat centroids

    """
    if throats is None:
        throats = network.throats()
    Ps = network['throat.conns']
    pdia = network['pore.diameter']

    if mode == 'centroid':
        try:
            pcentroids = network['pore.centroid']
            tcentroids = network['throat.centroid']
            if _sp.sum(_sp.isnan(pcentroids)) + _sp.sum(_sp.isnan(tcentroids)) > 0:
                mode = 'pore'
            else:
                plen1 = _sp.sqrt(_sp.sum(_sp.square(pcentroids[Ps[:, 0]] -
                                         tcentroids), 1))-network['throat.length']/2
                plen2 = _sp.sqrt(_sp.sum(_sp.square(pcentroids[Ps[:, 1]] -
                                         tcentroids), 1))-network['throat.length']/2
        except KeyError:
            mode = 'pore'
    if mode == 'pore':
        # Find half-lengths of each pore
        pcoords = network['pore.coords']
        # Find the pore-to-pore distance, minus the throat length
        lengths = _sp.sqrt(_sp.sum(_sp.square(pcoords[Ps[:, 0]] -
                                   pcoords[Ps[:, 1]]), 1)) - network['throat.length']
        lengths[lengths < 0.0] = 2e-9
        # Calculate the fraction of that distance from the first pore
        try:
            fractions = pdia[Ps[:, 0]]/(pdia[Ps[:, 0]] + pdia[Ps[:, 1]])
            # Don't allow zero lengths
            # fractions[fractions == 0.0] = 0.5
            # fractions[fractions == 1.0] = 0.5
        except:
            fractions = 0.5
        plen1 = lengths*fractions
        plen2 = lengths*(1-fractions)

    return _sp.vstack((plen1, network['throat.length'], plen2)).T[throats]

	def get_chisq(self, pops, writeback=False):
		"""Calculate the chi-square goodness-of-fit between the experimental population abundances and the fitted ones.

		Arguments
		---------
		pops : ndarray
			The normalized abundances of each lattice+ligand stoichiometries at each of the provided component concentrations.
		writeback : boolean
			Update the experiment's simulated heat attribute (dQ_fit) with the provided Qs?
			
		Returns
		-------
		float
			The goodness of the fit, as a reduced chi-square.

		Notes
		-----
			This method takes advantage of the fact that ITCSim doesn't inspect the data that is returned by the model, and instead lets the associated experiment handle the goodness-of-fit.
			The only caveat of course is that model must return the same number of lattice+ligand stoichiometries as are present in the experimental results.
			Variances (sigma**2) must have been precomputed as self.PopSigmas (perhaps should be self.PopVariances?)
		"""

		assert self.PopIntens.shape == pops.shape

		self.chisq = scipy.sum(scipy.square(self.PopIntens - pops) / self.PopSigmas) / self.PopIntens.size

		if writeback:
			self.PopFits = pops		

		return self.chisq

    def predict_gmm(self, testSamples, tau=0):
        """
            Function that predict the label for testSamples using the learned model
            Inputs:
                testSamples: the samples to be classified
                tau:         regularization parameter
            Outputs:
                predLabels: the class
                scores:     the decision value for each class
        """
        # Get information from the data
        nbTestSpl = testSamples.shape[0] # Number of testing samples

        # Initialization
        scores = sp.empty((nbTestSpl,self.C))

        # Start the prediction for each class
        for c in xrange(self.C):
            testSamples_c = testSamples - self.mean[c,:]

            regvp = self.vp[c,:] + tau

            logdet        = sp.sum(sp.log(regvp))
            cst           = logdet - 2*sp.log(self.prop[c]) # Pre compute the constant term

            # compute ||lambda^{-0.5}q^T(x-mu)||^2 + cst for all samples
            scores[:,c] = sp.sum( sp.square( sp.dot( (self.Q[c,:,:][:,:]/sp.sqrt(regvp)).T, testSamples_c.T ) ), axis=0 ) + cst

            del testSamples_c

        # Assign the label to the minimum value of scores
        predLabels = sp.argmin(scores,1)+1

        return predLabels,scores

def objF(x):
	global INDEX
	modelname = 'Model-'+str(INDEX)
	INDEX += 1
	K,e0,n,rth,rz= x
	parts = [7.85e-9,210000.0,.3,K,e0,n,rth,rz]		
	paramFile = open('result.txt','a+')
	paramFile.write('%s %10.6E %10.6E %10.6E %10.6E %10.6E '%(str(INDEX),K,e0,n,rth,rz))
	paramFile.flush()
	material = {'part':parts}
	shapes = {'outDimater':72.5,'thick':3.65,'length':145}
	mesh={'tube':100}
	timelist = [.20,.21,.22,.24]
	amps = [.5547,.576,.5824,.5982]
	midpairs = [(.05,.0329),(.1,.2246),(.15,.4013),(.18,.4957)]
	args={'timelist':timelist,'amp':amps,'midpairs':midpairs}
	load=args
	inits={}
	BCs={}
	positions = {}
	material={'part':parts}
	meshSize = {'pressDie':shapes['thick']*3/0.6,'tube':shapes['thick']/0.6 * 2}
	t = TBFEA(modelname)
	t.setParameter(shapes,material,positions,inits,\
			len(timelist),BCs,load,meshSize,args)
	t.setModels()
	coords = t.getResults()
	npFEA = scipy.array(coords)
	npExp = scipy.array([37.00,37.55,38.53,40.58])
	npres = npFEA - npExp
	mse=scipy.sum(scipy.square(npres))
	paramFile.write('%10.6E\n'%(mse))
	paramFile.close()
	return mse

 def propup(self, X, eps=1e-8):
     #~ F = self.W.dot(X.T) 
     F = X.dot(self.W.T).T
     Fs = sqrt(square(F) + eps)
     NFs, L2Fs = l2row(Fs)
     Fhat, L2Fn = l2row(NFs.T)
     return F, Fs, NFs, L2Fs, Fhat, L2Fn

def __fit(index1,index2):
       from scipy import stats, sqrt, square
       print strain
       print stress
       (cijFitted,intercept,r,tt,stderr) = stats.linregress(strain[:,index2-1],stress[:,index1-1])
       if (S.__version__ < '0.7.0'):
           stderr = S.sqrt((numsteps * stderr**2)/(numsteps-2))
           error  = stderr/sqrt(sum(square(strain[:,index2-1])))
       else:
           fit_str = ((strain[index2-1,:] * cijFitted) + intercept)                
           error = sqrt((sum(square(stress[:,index1-1] - fit_str)) / \
                       (numsteps-2))/(sum(square(strain[:,index2-1]))))
       print 'Cij   ', cijFitted
       print 'Error   ', error
       print 'intercept   ', intercept
       return cijFitted, error

    def learn(self, X, t, tol=0.01, amax=1e10):
        u"""学習"""

        N = X.shape[0]
        a = sp.ones(N+1) # hyperparameter
        b = 1.0
        phi = sp.ones((N, N+1)) # design matrix
        phi[:,1:] = [[self._kernel(xi, xj) for xj in X] for xi in X]

        diff = 1
        while diff >= tol:
            sigma = spla.inv(sp.diag(a) + b * sp.dot(phi.T, phi))
            m = b * sp.dot(sigma, sp.dot(phi.T, t))
            gamma = sp.ones(N+1) - a * sigma.diagonal()
            anew = gamma / (m * m)
            bnew = (N -  gamma.sum()) / sp.square(spla.norm(t - sp.dot(phi, m)))
            anew[anew >= amax] = amax
            adiff, bdiff = anew - a, bnew - b
            diff = (adiff * adiff).sum() + bdiff * bdiff
            a, b = anew, bnew
            print ".",

        self._a = a
        self._b = b
        self._X = X
        self._m = m
        self._sigma = sigma
        self._amax = amax

def euclidean(h1, h2):  # 9 us @array, 33 us @list \w 100 bins
    """
    Equal to Minowski distance with p=2.
    @see minowski()
    """
    h1, h2 = __prepare_histogram(h1, h2)
    return math.sqrt(scipy.sum(scipy.square(scipy.absolute(h1 - h2))))

    def errorApproximation(self, ratio, dim=20):

        self.buildMatrix()

        sumNonzeros = (self.vxm !=0).sum()
        numTest = int(ratio*sumNonzeros)

        elementList = []

        nonZeroTuple = sp.nonzero(self.vxm)

        for x in range(int(numTest)):
            rInt = sp.random.randint(0,nonZeroTuple[0].size)
            randrow = nonZeroTuple[0][rInt]
            randcolumn = nonZeroTuple[1][rInt]

            valElementIndex = [randrow,randcolumn]
            elementList.append(valElementIndex)

        self.modvxm = sp.copy(self.vxm)

        for x in elementList:
            self.modvxm[x[0],x[1]] = 0

        self.modvmx = self.fillAverages(vxm = self.modvxm)
        self.newmodvxm = self.predict(dim,vxm=self.modvxm)

        sqDiff = 0
        for x in elementList:
            sqDiff += sp.square(self.newmodvxm[x[0],x[1]] - self.vxm[x[0],x[1]])
        self.rmse = sp.sqrt(sqDiff/len(elementList))

def relative_deviation(h1, h2):  # 18 us @array, 42 us @list \w 100 bins
    """
    Calculate the deviation between two histograms.
    The relative deviation between two histograms \f$H\f$ and \f$H'\f$ of size \f$m\f$ is
    defined as
    \f[
        d_{rd}(H, H') =
            \frac{
                \sqrt{\sum_{m=1}^M(H_m - H'_m)^2}
              }{
                \frac{1}{2}
                \left(
                    \sqrt{\sum_{m=1}^M H_m^2} +
                    \sqrt{\sum_{m=1}^M {H'}_m^2}
                \right)
              }
    \f]
    
    Attributes:
    - semimetric (triangle equation satisfied?)
    
    Attributes for normalized histograms:
    - \f$d(H, H')\in[0, \sqrt{2}]\f$
    - \f$d(H, H) = 0\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-normalized histograms:
    - \f$d(H, H')\in[0, 2]\f$
    - \f$d(H, H) = 0\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-equal histograms:
    - not applicable    
    
    @param h1 the first histogram
    @type h1 array-like sequence
    @param h2 the second histogram, same bins as h1
    @type h2 array-like sequence
    
    @return relative deviation
    @rtype float
    """
    h1, h2 = __prepare_histogram(h1, h2)
    numerator = math.sqrt(scipy.sum(scipy.square(h1 - h2)))
    denominator = (math.sqrt(scipy.sum(scipy.square(h1))) + math.sqrt(scipy.sum(scipy.square(h2)))) / 2.0
    return numerator / denominator

def correlate(h1, h2):  # 31 us @array, 55 us @list \w 100 bins
    """
    Compute the correlation between two histograms.
    The histogram correlation between two histograms \f$H\f$ and \f$H'\f$ of size \f$m\f$
    is defined as
    \f[
        d_{corr}(H, H') = 
        \frac{
            \sum_{m=1}^M (H_m-\bar{H}) \cdot (H'_m-\bar{H'})
        }{
            \sqrt{\sum_{m=1}^M (H_m-\bar{H})^2 \cdot \sum_{m=1}^M (H'_m-\bar{H'})^2}
        }
    \f]
    with \f$\bar{H}\f$ and \f$\bar{H'}\f$ being the mean values of \f$H\f$ resp. \f$H'\f$
        
    Attributes:
    - not a metric, a similarity
    
    Attributes for normalized histograms:
    - \f$d(H, H')\in[-1, 1]\f$
    - \f$d(H, H) = 1\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-normalized histograms:
    - \f$d(H, H')\in[-1, 1]\f$
    - \f$d(H, H) = 1\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-equal histograms:
    - not applicable
    
    @note returns 0 if one of h1 or h2 contains only zeros.
    
    @param h1 the first histogram
    @type h1 array-like sequence
    @param h2 the second histogram, same bins as h1
    @type h2 array-like sequence    
    """
    h1, h2 = __prepare_histogram(h1, h2)
    h1m = h1 - scipy.sum(h1) / float(h1.size)
    h2m = h2 - scipy.sum(h2) / float(h2.size)
    a = scipy.sum(scipy.multiply(h1m, h2m))
    b = math.sqrt(scipy.sum(scipy.square(h1m)) * scipy.sum(scipy.square(h2m)))
    return 0 if 0 == b else a / b