def testAll(tests, allalgos, tolerant=True):
    countgood = 0
    for i, algo in enumerate(sorted(allalgos)):
        print "%d, %s:" % (i + 1, algo.__name__)
        print " " * int(log10(i + 1) + 2),
        good = True
        messages = []
        for t in tests:
            try:
                res = t(algo)
            except Exception, e:
                if not tolerant:
                    raise e
                res = e

            if res is True:
                print ".",
            else:
                good = False
                messages.append(res)
                print "F",
        if good:
            countgood += 1
            print "--- OK."
        else:
            print "--- NOT OK."
            for m in messages:
                if m is not None:
                    print " " * int(log10(i + 1) + 2), "->", m

def testAll(tests, allalgos, tolerant=True):
    countgood = 0
    for i, algo in enumerate(sorted(allalgos)):
        print("%d, %s:" % (i + 1, algo.__name__))
        print(' ' * int(log10(i + 1) + 2),)
        good = True
        messages = []
        for t in tests:
            try:
                res = t(algo)
            except Exception, e:
                if not tolerant:
                    raise e
                res = e

            if res is True:
                print('.',)
            else:
                good = False
                messages.append(res)
                print('F',)
        if good:
            countgood += 1
            print('--- OK.')
        else:
            print('--- NOT OK.')
            for m in messages:
                if m is not None:
                    print(' ' * int(log10(i + 1) + 2), '->', m)

 def residual_lmfit(self, pars, x, y):
     a = P4Rm()
     self.strain_DW(pars)
     res = f_Refl_fit(a.AllDataDict["geometry"], self.Data4f_Refl)
     y_cal = convolve(abs(res) ** 2, a.ParamDict["resol"], mode="same")
     y_cal = y_cal / y_cal.max() + a.AllDataDict["background"]
     return log10(y) - log10(y_cal)

    def run(self, npts=25, inv_points=None, access_limited=True, **kwargs):
        r"""
        Parameters
        ----------
        npts : int (default = 25)
            The number of pressure points to apply.  The list of pressures
            is logarithmically spaced between the lowest and highest throat
            entry pressures in the network.

        inv_points : array_like, optional
            A list of specific pressure point(s) to apply.

        """
        if 'inlets' in kwargs.keys():
            logger.info('Inlets recieved, passing to set_inlets')
            self.set_inlets(pores=kwargs['inlets'])
        if 'outlets' in kwargs.keys():
            logger.info('Outlets recieved, passing to set_outlets')
            self.set_outlets(pores=kwargs['outlets'])
        self._AL = access_limited
        if inv_points is None:
            logger.info('Generating list of invasion pressures')
            min_p = sp.amin(self['throat.entry_pressure']) * 0.98  # nudge down
            max_p = sp.amax(self['throat.entry_pressure']) * 1.02  # bump up
            inv_points = sp.logspace(sp.log10(min_p),
                                     sp.log10(max_p),
                                     npts)

        self._npts = sp.size(inv_points)
        # Execute calculation
        self._do_outer_iteration_stage(inv_points)

def getLogBins(first_point, last_point, log_step):
    """
    get the bin in log scale and the center bin value
    
    Parameters:
    ----------------
    first_point, last_point : number
    First and last point of the x-axis
    
    log_step : number
    Required log-distance between x-points
    
    Returns:
    -----------
    xbins : array of the x values at the center (in log-scale) of the bin
    bins : array of the x values of the bins 
    """
    log_first_point = scipy.log10(first_point)
    log_last_point = scipy.log10(last_point)
    # Calculate the bins as required by the histogram function, i.e. the bins edges including the rightmost one
    N_log_steps = scipy.floor((log_last_point - log_first_point) / log_step) + 1.0
    llp = N_log_steps * log_step + log_first_point
    bins_in_log_scale = np.linspace(log_first_point, llp, N_log_steps + 1)
    bins = 10 ** bins_in_log_scale
    center_of_bins_log_scale = bins_in_log_scale[:-1] + log_step / 2.0
    xbins = 10 ** center_of_bins_log_scale
    return xbins, bins

def powerlaw_fit(xdata, ydata, yerr):
    # Power-law fitting is best done by first converting
    # to a linear equation and then fitting to a straight line.
    #  y = a * x^b
    #  log(y) = log(a) + b*log(x)
    from scipy import log10
    from scipy import optimize

    powerlaw = lambda x, amp, index: amp*np.power(x,index)

    logx = log10(xdata)
    logy = log10(ydata)
    logyerr = yerr / ydata

    # define our (line) fitting function
    fitfunc = lambda p, x: p[0] + p[1] * x
    errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err

    pinit = [1.0, -1.0]
    out = optimize.leastsq(errfunc, pinit, args=(logx, logy, logyerr), full_output=1)	
    pfinal = out[0]
    covar = out[1]

    #y = amp * x^exponent
    exponent = pfinal[1] #index in original
    amp = 10.0**pfinal[0]
    
    exponentErr = np.sqrt( covar[0][0] )
    ampErr = np.sqrt( covar[1][1] ) * amp

    chisq = np.sum((((ydata - powerlaw(xdata, amp, exponent))/yerr)**2),axis=0)

    return exponent, amp, chisq

def degree_distrib(net, deg_type="total", node_list=None, use_weights=True,
                   log=False, num_bins=30):
    '''
    Computing the degree distribution of a network.
    
    Parameters
    ----------
    net : :class:`~nngt.Graph` or subclass
        the network to analyze.
    deg_type : string, optional (default: "total")
        type of degree to consider ("in", "out", or "total").
    node_list : list or numpy.array of ints, optional (default: None)
        Restrict the distribution to a set of nodes (default: all nodes).
    use_weights : bool, optional (default: True)
        use weighted degrees (do not take the sign into account: all weights
        are positive).
    log : bool, optional (default: False)
        use log-spaced bins.
    
    Returns
    -------
    counts : :class:`numpy.array`
        number of nodes in each bin
    deg : :class:`numpy.array`
        bins
    '''
    ia_node_deg = net.get_degrees(node_list, deg_type, use_weights)
    ra_bins = sp.linspace(ia_node_deg.min(), ia_node_deg.max(), num_bins)
    if log:
        ra_bins = sp.logspace(sp.log10(sp.maximum(ia_node_deg.min(),1)),
                               sp.log10(ia_node_deg.max()), num_bins)
    counts,deg = sp.histogram(ia_node_deg, ra_bins)
    ia_indices = sp.argwhere(counts)
    return counts[ia_indices], deg[ia_indices]

def plotHeatmap(fwrap, aclass, algoparams, trials, maxsteps):
    """ Visualizing performance across trials and across time 
    (iterations in powers of 2) """
    psteps = int(log2(maxsteps)) + 1
    storesteps = [0] + [2 ** x  for x in range(psteps)]
    ls = lossTraces(fwrap, aclass, dim=trials, maxsteps=maxsteps,
                    storesteps=storesteps, algoparams=algoparams,
                    minLoss=1e-10)
            
    initv = mean(ls[0])
    maxgain = exp(fwrap.stochfun.maxLogGain(maxsteps) + 1)
    maxneggain = (sqrt(maxgain))
    
    M = zeros((psteps, trials))
    for sid in range(psteps):
        # skip the initial values
        winfactors = clip(initv / ls[sid+1], 1. / maxneggain, maxgain)
        winfactors[isnan(winfactors)] = 1. / maxneggain
        M[sid, :] = log10(sorted(winfactors))
        
    pylab.imshow(M.T, interpolation='nearest', cmap=cm.RdBu, #@UndefinedVariable
                 aspect=psteps / float(trials) / 1,  
                 vmin= -log10(maxgain), vmax=log10(maxgain),
                 )   
    pylab.xticks([])
    pylab.yticks([])
    return ls

 def fit(self, kk=None):
     """
     Fit Fourier spectrum with the function set at class instantination
     ==> NB: fitting is done in logarithmic coordinates
     and fills plotting arrays with data
     --------
     Options:
     --------
     kk
        (k1,k2) <None> spectral interval for function fitting
        by default interval [ kk[1], kk[imax__kk] ] will be fitted
        ==> i.e. k=0 is excluded
     """
     # fitting interval
     if kk:
         ik_min=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[0]).nonzero()[0][-1]
         ik_max=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[1]).nonzero()[0][-1]
     else:
         ik_min=1;
         ik_max=self.fft_data.imax__kk
     # do fitting
     self.__popt,self.__pcov = scipy.optimize.curve_fit(self.__func_fit,
                                                        scipy.log(self.fft_data.kk[ik_min:ik_max]),
                                                        scipy.log(self.fft_data.Ik[ik_min:ik_max]) )
     # boundaries of fitted interval
     self.kmin = self.fft_data.kk[ik_min]
     self.kmax = self.fft_data.kk[ik_max]
     # fill plot arrays <===============
     self.kk_plot=scipy.logspace( scipy.log10(self.kmin),
                                  scipy.log10(self.kmax),
                                  self.nk_plot )
     self.Ik_plot=self.fitting_function(self.kk_plot)

def create_grid(r_in, r_out, nshell, space = 'powerlaw1', end = True):
    # function to create grid
    if space == 'log10':
        from scipy import log10, logspace
        # get the exponent of the start- and
        # stop-radius in input units
        start = [log10(r_in), 0][r_in == 0]
        stop = log10(r_out)
        radii = logspace(start, stop, num=nshell, endpoint=end)
    elif space == "powerlaw1":
        from scipy import arange
        radii = r_in * (r_out/r_in)**(arange(nshell)/(nshell - 1.0))
    elif space == 'linear':
        from scipy import linspace
        # linearly spaced grid
        radii = linspace(r_in, r_out, num=nshell, endpoint=end)
    elif space == 'powerlaw2':
        from scipy import linspace
        # first check if coefficients to the power-law was given
        #~ if 'exp' in kwargs:
            #~ p_exp = kwargs['exp']
        #~ else: # if not, set it to 2, i.e. r^2
            #~ p_exp = 2
        radii = r_in + (r_out - r_in)*(linspace(r_in, r_out, num=nshell, endpoint=end)/(r_out))**2
        #pr_int('Not implemented yet.')
        #raise ParError(spaced)
    else:
        raise Exception(space)
    return radii

def testPlot():
    """ Get/generate the data to play with """
    TIME_INC = 1e-6
    NUM_POINTS = 40000
    t = timeScale(TIME_INC,NUM_POINTS)
    noisy_sig = genData(TIME_INC,True, t)
    clean_sig = genData(TIME_INC,False,t)
    """ Get FFT of signal and the sampling frequency from the time intervals used to generate the signals"""
    freq, s_fft  = getFFT(noisy_sig, TIME_INC)
    freq2,s_fft2 = getFFT(clean_sig, TIME_INC)


    """ Show in 2 subplots the signals and their spectrums"""
    plb.subplot(211,axisbg='#FFFFCC')
    p.plot(t,clean_sig,'b')
    p.hold(True)
    p.grid(True)
    p.plot(t,noisy_sig,'r')
    plb.subplot(212,axisbg='#FFFFCC')
    #p.hold(False)
    p.plot(freq2, 20*sp.log10(s_fft2),'x-b')
    p.hold(True)
    p.plot(freq,  20*sp.log10(s_fft), '+-r')
    p.xticks([-10e4,-5e4,-4e4,-3e4,-2e4,-1e4,0,1e4,2e4,3e4,4e4,5e4,10e4])
    p.xlim([-1e5,1e5])
    p.grid(True)
    #p.show()
    q = ScrollingToolQT(p.gcf())
    return q   # WARNING: it's important to return this object otherwise

def plot_median_errors(RefinementLevels):
        for i in RefinementLevels[0].cases:
            x =[];
            y =[];
            print "Analyzing median error on: ", i ;
            for r in RefinementLevels:                
                x.append(r.LUT.D_dim*r.LUT.P_dim)
                r.get_REL_ERR_SU2(i)
                y.append(r.SU2[i].median_ERR*100)
            
            x = sp.array(x)
            y = sp.array(y)            
            y = y[sp.argsort(x)]
            x = x[sp.argsort(x)]
                                    
            LHM = sp.ones((len(x),2))
            RHS = sp.ones((len(x),1))            
            LHM[:,1] = sp.log10(x)
            RHS[:,0] = sp.log10(y)

            sols = sp.linalg.lstsq(LHM,RHS)
            b = -sols[0][1]
            plt.loglog(x,y, label='%s, %s'%(i,r'$O(\frac{1}{N})^{%s}$'%str(sp.around(b,2))), basex=10, basey=10, \
                       subsy=sp.linspace(10**(-5), 10**(-2),20),\
                       subsx=sp.linspace(10**(2), 10**(5),50))
            
            #for r in RefinementLevels:                
               # x.append(r.LUT.D_dim*r.LUT.P_dim)
              #  r.get_REL_ERR_SciPy(i)
             #   y.append(r.SciPy[i].median_ERR*100)
            #plt.plot(x,y, label='SciPy: %s'%i)
        plt.grid(which='both')
        plt.xlabel('Grid Nodes (N)')
        plt.ylabel('Median relative error [%]')
        return;

def testAll(tests, allalgos, tolerant=True):
    countgood = 0
    for i, algo in enumerate(sorted(allalgos)):
        print(("%d, %s:" % (i + 1, algo.__name__)))
        print((" " * int(log10(i + 1) + 2),))
        good = True
        messages = []
        for t in tests:
            try:
                res = t(algo)
            except Exception as e:
                if not tolerant:
                    raise e
                res = e

            if res is True:
                print((".",))
            else:
                good = False
                messages.append(res)
                print(("F",))
        if good:
            countgood += 1
            print("--- OK.")
        else:
            print("--- NOT OK.")
            for m in messages:
                if m is not None:
                    print((" " * int(log10(i + 1) + 2), "->", m))
    print()
    print(("Summary:", countgood, "/", len(allalgos), "of test were passed."))

def entropyloss(act, pred):
    epsilon = 1e-15
    pred = sp.maximum(epsilon, pred)
    pred = sp.minimum(1-epsilon, pred)
    el = sum(act*sp.log10(pred) + sp.subtract(1,act)*sp.log10(sp.subtract(1,pred)))
    el = el * -1.0/len(act)
    return el

 def test_permutation(self):
     #test permutation function
     for dn in self.datasets:
         D = data.load(os.path.join(self.dir_name,dn))
         perm = SP.random.permutation(D['X'].shape[0])
         #1. set permuattion
         lmm = dlimix.CLMM()
         lmm.setK(D['K'])
         lmm.setSNPs(D['X'])
         lmm.setCovs(D['Cov'])
         lmm.setPheno(D['Y'])
         if 1:
             #pdb.set_trace()
             perm = SP.array(perm,dtype='int32')#Windows needs int32 as long -> fix interface to accept int64 types
         lmm.setPermutation(perm)
         lmm.process()
         pv_perm1 = lmm.getPv().ravel()
         #2. do by hand
         lmm = dlimix.CLMM()
         lmm.setK(D['K'])
         lmm.setSNPs(D['X'][perm])
         lmm.setCovs(D['Cov'])
         lmm.setPheno(D['Y'])
         lmm.process()
         pv_perm2 = lmm.getPv().ravel()
         D2 = (SP.log10(pv_perm1)-SP.log10(pv_perm2))**2
         RV = SP.sqrt(D2.mean())
         self.assertTrue(RV<1E-6)

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 addqqplotinfo(qnull,M,xl='-log10(P) observed',yl='-log10(P) expected',xlim=None,ylim=None,alphalevel=0.05,legendlist=None,fixaxes=False):    
    distr='log10'
    pl.plot([0,qnull.max()], [0,qnull.max()],'k')
    pl.ylabel(xl)
    pl.xlabel(yl)
    if xlim is not None:
        pl.xlim(xlim)
    if ylim is not None:
        pl.ylim(ylim)        
    if alphalevel is not None:
        if distr == 'log10':
            betaUp, betaDown, theoreticalPvals = _qqplot_bar(M=M,alphalevel=alphalevel,distr=distr)
            lower = -sp.log10(theoreticalPvals-betaDown)
            upper = -sp.log10(theoreticalPvals+betaUp)
            pl.fill_between(-sp.log10(theoreticalPvals),lower,upper,color="grey",alpha=0.5)
            #pl.plot(-sp.log10(theoreticalPvals),lower,'g-.')
            #pl.plot(-sp.log10(theoreticalPvals),upper,'g-.')
    if legendlist is not None:
        leg = pl.legend(legendlist, loc=4, numpoints=1)
        # set the markersize for the legend
        for lo in leg.legendHandles:
            lo.set_markersize(10)

    if fixaxes:
        fix_axes()        

 def getAxis(self,X,Y):
     """
     return the proper axis limits for the plots
     """
     out = []
     mM = [(min(X),max(X)),(min(Y),max(Y))]
     for i,j in mM:
         #YJC: checking if values are negative, if yes, return 0 and break
         if j <0 or i <0:
             return 0
         log_i = scipy.log10(i)
         d, I = scipy.modf(log_i)
         if log_i < 0:
             add = 0.5 *(scipy.absolute(d)<0.5)
         else:
             add = 0.5 *(scipy.absolute(d)>0.5)
         m = scipy.floor(log_i) + add
         out.append(10**m)
         log_j = scipy.log10(j)
         d, I = scipy.modf(log_j)
         if log_j < 0:
             add = - 0.5 *(scipy.absolute(d)>0.5)
         else:
             add = - 0.5 *(scipy.absolute(d)<0.5)
         m = scipy.ceil(log_j) + add
         out.append(10**m)
     return tuple(out)

    def PlotManhattan(self,
                      xname=str,
                      yname=str,
                      Log=Logger):

        LogString = '**** Generating Manhattan plot ...'
        print LogString
        Log.Write(LogString+'\n')

        XName = ''
        YName = ''
        for Key in self.DataContainers.iterkeys():
            if(re.search(xname,Key)):
                XName = Key
            if(re.search(yname,Key)):
                YName = Key

        X = []
        for Entry in self.DataContainers[XName].GetDataArray():
            X.append(int(Entry))
        X = scipy.array(X)
        Y = []
        for Entry in self.DataContainers[YName].GetDataArray():
            if(Entry=='-1'):
                Y.append(float(1.0))
            else:
                Y.append(float(Entry))
        Y = -scipy.log10(scipy.array(Y))

        PylabParameters,\
        Rectangle         = self.PylabUpdateParams()
        pylab.rcParams.update(PylabParameters)
        PylabFigure = pylab.figure()
        PylabFigure.clf()
        PylabAxis = PylabFigure.add_axes(Rectangle)
        PylabAxis.scatter(x=X,
                          y=Y)
        XSign = scipy.array(PylabAxis.get_xlim())
        YSign = -scipy.log10(scipy.array([5.0e-8,5.0e-8]))
        PylabAxis.plot(XSign,
                       YSign,
                       linestyle='--',
                       color='grey',
                       linewidth=1.0)
        XSugg = scipy.array(PylabAxis.get_xlim())
        YSugg = -scipy.log10(scipy.array([1.0e-5,1.0e-5]))
        PylabAxis.plot(XSugg,
                       YSugg,
                       linestyle=':',
                       color='grey',
                       linewidth=1.0)
        PylabAxis.set_ylim([0.0,PylabAxis.get_ylim()[1]])
        PylabFigure.savefig('Manhattan.png')
        PylabAxis.clear()
        pylab.close(PylabFigure)
        del PylabFigure
        del PylabAxis

        return

def effective_order(x, y):
    '''Find slope of log log plot to find our effective order of accuracy'''

    logx = log10(x)
    logy = log10(y)
    out = curve_fit(linear_fit, logx, logy)

    return out[0][1]