def comp_conv(f,a,b,exact_value):
    nmax = 100
    
    x = np.arange(1.0, nmax, 1.0)

    t = np.arange(1.0, nmax, 1.0)
    u = np.arange(1.0, nmax, 1.0)
    v = np.arange(1.0, nmax, 1.0)
    w = np.arange(1.0, nmax, 1.0)

    for i in np.arange(0.0,t.size):
        t[i] = rect_meth(a,b,x[i],f)

    for i in np.arange(0.0,u.size):
        u[i] = simpson_meth(a,b,x[i],f)

    for i in np.arange(0.0,v.size):
        v[i] = midpoint_meth(a,b,x[i],f)

    for i in np.arange(0.0,v.size):
        w[i] = exact_value

    mp.clf()
    mp.semilogx()
    mp.plot(x, t, linewidth=1.0, label='Rectangle')
    mp.plot(x, u, linewidth=1.0, label='Simpson')
    mp.plot(x, v, linewidth=1.0, label='Midpoint')
    mp.plot(x, w, linewidth=1.0, label='Value of the integral')
    
    
    mp.title("Illustration of the convergence of the three integration methods");
    mp.xlabel('Number of subdivision points')
    mp.ylabel('Integral of f between a and b')
    mp.legend(loc='upper right')
    mp.show()

def plot_validation_curve(model, X, y, scorer, param_name, param_range=np.linspace(0.1, 1, 5), cv=None, n_jobs=5,
    ylim=None, title="Xval. validation curve"):
    ''' Plot learning curve for model on data '''

    df = pd.DataFrame()
    df['param_range'] = param_range
    train_scores, test_scores = validation_curve(model, X, y, param_name=param_name, param_range=param_range,
        cv=cv, scoring=scorer, n_jobs=n_jobs)
    df['train_mean'] = 1 - np.mean(train_scores, axis=1)
    df['train_std'] = np.std(train_scores, axis=1)
    df['test_mean'] = 1 - np.mean(test_scores, axis=1)
    df['test_std'] = np.std(test_scores, axis=1)

    plt.figure()
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Parameter value")
    plt.ylabel("Error (1-score)")
    plt.grid()
    plt.semilogx(param_range, df.train_mean, color="r", label="Training")
    plt.fill_between(param_range, df.train_mean - df.train_std, df.train_mean + df.train_std, alpha=0.1, color="r")
    plt.semilogx(param_range, df.test_mean, color="g", label="Test")
    plt.fill_between(param_range, df.test_mean - df.test_std, df.test_mean + df.test_std, alpha=0.1, color="g")
    plt.legend(loc="best")
    plt.show()
    return df, plt

def demon(n, d, g, kgiven=0):
    ks = np.arange(n + 1)
    c1, c2 = [], []
    for k in ks:
        a = hyp2f1(1, n - k + 1, n + 2, 1 - d)
        b = 0 if g == 0 else hyp2f1(1, n - k + 1, n + 2, 1 - d * g)
        c1.append(a * d - b * d * g)
        c2.append(a)
    plt.plot(ks / n, c1, label='c1')
    plt.plot(ks / n, c2, label='c2')
    plt.xlabel('k / n')
    plt.ylabel('Effective n / n')
    plt.legend()

    plt.figure()
    S = [100 * generalised_entropy_gain(n, k, d, g) for k in ks]
    plt.plot(ks / n, S, label='S')
    plt.xlabel('k / n')
    plt.ylabel('information gain')

    plt.figure()
    dif = np.logspace(-1, 1, 100)
    S = [100 * generalised_entropy_gain(n, kgiven, d, g) for d in dif]
    plt.semilogx(dif, S, label='S')
    plt.xlabel('difficulty')
    plt.ylabel('information gain')

def main():
    edges = [30, 60, 120, 240]
    corners = zip(edges[:-1], edges[1:])
    centres = [(a + b) / 2 for a, b in corners]

#c = [get_linkwitz_riley_coeffs(1, b, a, edges[-1] * 2) for b, a in corners]
    sr = 2000
    c = [get_peak_coeffs(-24, i, sr, 1) for i in centres]
    c.append([[1, 0, 0], [1, 0, 0]])

    bm = [BiquadMemory(0, 0) for _ in c]
    bc = [BiquadCoefficients(b0, b1, b2, a1, a2)
          for [b0, b1, b2], [a0, a1, a2] in c]

    c.append(series_coeffs(c))

#    c.append(impedance_filter(c[-1]))

    wh = [signal.freqz(b, a) for b, a in c]

    plt.subplot(111)
    plt.title("Frequency response - reflection filter")
    for w, h in wh:
        plt.semilogx(w, 20 * np.log10(np.abs(h)))
    plt.ylabel('Amplitude Response (dB)')
    plt.xlabel('Frequency (rad/sample)')
    plt.grid()

    plt.show()

    def plotPredictions(self):

        plt.semilogx(self.freqs, self.sPLs, label = 'ISO target')
        plt.semilogx(self.freqs, self.predictions, color = 'r', 
                linestyle = '--', label = 'Predicted')
        plt.legend()
        plt.show()

def plot_containment_radii(fraction):
    """Plotting script for 68% and 95% containment radii."""

    psf_gc = FermiGalacticCenter.psf()
    gtpsf_table_gc = get_psf_table(psf_gc,  10000, 300000, 15)

    psf_vela = FermiVelaRegion.psf()
    gtpsf_table_vela = get_psf_table(psf_vela, 10000, 300000, 15)
    
    if fraction == 68:
        true_table_rep = load_lat_psf_performance('P7REP_SOURCE_V15_68')
        true_table = load_lat_psf_performance('P7SOURCEV6_68')
        rad = 'CONT_68'
    elif fraction == 95:
        true_table_rep = load_lat_psf_performance('P7REP_SOURCE_V15_95')
        true_table = load_lat_psf_performance('P7SOURCEV6_95')
        rad = 'CONT_95'
    
    plt.plot(gtpsf_table_gc['ENERGY'], gtpsf_table_gc[rad],
             color='red',label='Fermi Tools PSF @ Galactic Center')
    plt.plot(gtpsf_table_vela['ENERGY'], gtpsf_table_vela[rad],
             color='blue', label='Fermi Tools PSF @ Vela Region')
    plt.plot(true_table_rep['energy'], true_table_rep['containment_angle'],
             color='green', linestyle='--', label='P7REP_SOURCE_V15')
    plt.plot(true_table['energy'], true_table['containment_angle'],
             color='black', linestyle='--', label='P7SOURCEV6')

    plt.xlim([10000, 300000])

    plt.legend()
    plt.semilogx()
    plt.xlabel('Energy/MeV')
    plt.ylabel('PSF Containment Radius/deg')

    return plt

def plotGroupXS(group_xs, title='', filename=''):
    
    global subdirectory

    directory = pinspec.get_output_directory() + subdirectory

    # Make directory if it does not exist
    if not os.path.exists(directory):
            os.makedirs(directory)

    # Plot Resonance Integrals
    fig = plt.figure()
    bins = group_xs.bin_edges
    plt.semilogx(bins[0:-1], group_xs.groupXS[:,:], drawstyle='steps-post')
    plt.xlabel('Energy [eV]')
    plt.ylabel('Group XS')
    plt.grid()

    if title is '':
        plt.title('Group XS')
    else:
        plt.title(title.title() + ' Group XS')

    if filename is '':
        filename = directory + '/group-xs.png'
    else:
        filename = directory + filename.replace(' ', '-').lower() + \
                                                        '-group-xs.png'

    plt.savefig(filename)
    plt.close(fig)

def plotValidationCurve(c, training_avgs, validation_avgs, model_name,
                        x_label = 'C'):
    """                                                                                                                                                       
    This function plots the training and validation averages                                                                                                  
    Against the parameter we are trying to optimize                                                                                                           
    """
    plt.clf()
    
    p1, = plt.plot(c, training_avgs, 'ro-', label = 'training')
    p2, = plt.plot(c, validation_avgs, 'go-', label = 'validation')

    plt.xlabel(x_label, fontsize = 16)
    plt.ylabel('Score', fontsize = 16)

    plt.title(model_name + ' Validation Curve',
              fontdict = {'fontsize': 16})
    plt.legend(loc = 0)
    plt.semilogx()

    # now save the figure                                                                                   
    model_name = re.sub(r'\W', '_', model_name)
    plt.savefig(model_name + '_validation.png',
                format = 'png')
    
    return None

 def plotDimensionsUpdateFrequencyEstimation(self, returnAxisValuesOnly=True):
     '''
     numberOfTimeUnits=10*24*12
     Experts stream 12
     Houston stream 2
     '''
     dataDistribution = defaultdict(list)
     for line in FileIO.iterateJsonFromFile(self.dimensionsUpdateFrequencyFile):
         for k, v in line[ParameterEstimation.dimensionsUpdateFrequencyId].iteritems():
             k = int(k) / self.timeUnitInSeconds.seconds
             if k not in dataDistribution: dataDistribution[k] = [0., 0.]
             dataDistribution[k][0] += v; dataDistribution[k][1] += 1
     x, y = [], []; [(x.append(k), y.append((dataDistribution[k][0] / dataDistribution[k][1]))) for k in sorted(dataDistribution)]
     x1, y1 = [], []; [(x1.append(k), y1.append((dataDistribution[k][0] / dataDistribution[k][1]) / k)) for k in sorted(dataDistribution)]
     x = x[:numberOfTimeUnits]; y = y[:numberOfTimeUnits]; x1 = x1[:numberOfTimeUnits]; y1 = y1[:numberOfTimeUnits]
     def subPlot(id):
         plt.subplot(id)
         inactivityCorordinates = max(zip(x1, y1), key=itemgetter(1))
         plt.semilogx(x1, y1, '-', color=self.stream_settings['plot_color'], label=getLatexForString(self.stream_settings['plot_label'] + ' (Update frequency=%d TU)' % inactivityCorordinates[0]), lw=2)
         plt.subplot(id).yaxis.set_major_formatter(FuncFormatter(lambda x, i: '%0.1f' % (x / 10. ** 3)))
         plt.semilogx([inactivityCorordinates[0]], [inactivityCorordinates[1]], 'o', alpha=0.7, color='r')
         plt.subplot(id).yaxis.set_major_formatter(FuncFormatter(lambda x, i: '%0.1f' % (x / 10. ** 3)))
         plt.yticks((min(y1), max(y1)))
         print self.stream_settings['plot_label'], inactivityCorordinates[0]
     plt.subplot(311)
     plt.title(getLatexForString('Dimensions update frequency estimation'))
     plt.semilogx(x, y, '-', color=self.stream_settings['plot_color'], label=getLatexForString(self.stream_settings['plot_label']), lw=2)
     plt.subplot(311).yaxis.set_major_formatter(FuncFormatter(lambda x, i: '%0.1f' % (x / 10. ** 5)))
     plt.text(0.0, 1.01, getLatexForString('10^5'), transform=plt.gca().transAxes)
     plt.ylabel(getLatexForString('\# of decayed dimensions'))
     if self.stream_settings['stream_id'] == 'experts_twitter_stream': subPlot(312)
     else: subPlot(313); plt.xlabel(getLatexForString(xlabelTimeUnits))
     plt.ylabel(getLatexForString('Rate of DD (10^3)'))
     plt.legend(loc=3)
     if returnAxisValuesOnly: plt.show()

def Pl(Outfile,line='-.',label=''):
    dd=np.loadtxt(Outfile+'PS_DD')
    dk=np.loadtxt(Outfile+'PS_DK')
    kk=np.loadtxt(Outfile+'PS_KK')
    halo=np.loadtxt('/home/mtx/data/tide/outdata/test/PS_haloDD')
    #plt.figure('Power Spectrum')
    #n=np.ones_like(kk[:,1])
    #plt.loglog(dd[:,0],dd[:,1],'b.-',label='P_DD')
    #plt.loglog(halo[:,0],halo[:,1],'m.-',label='P_halo')
    #plt.loglog(dk[:,0],dk[:,1],'g.-',label='P_DK')
    #plt.loglog(kk[:,0],kk[:,1],'r.-',label='P_KK')
    #plt.loglog(kk[:,0],n*noise3,'k-.',alpha=0.8,label='noise')
    #plt.title('Power Spectrum')
    #plt.xlabel('k $(h/Mpc)$')
    #plt.ylabel('P(k) $(Mpc^{3}/h^{3})$')
    #plt.legend()
    #plt.ylim(10**1,6*10**4)
    #plt.xlim([9*10**-3,1])
    ##plt.savefig(name+'PS.eps')
    
    plt.figure('correlation coefficient')
    plt.title('correlation coefficient')
    plt.semilogx(kk[:,0],dk[:,1]/np.sqrt(dd[:,1]*kk[:,1]),line,label=label)
    plt.xlabel('k $(h/Mpc)$')
    plt.ylabel('r')

 def plot_fft(self, plotname=None, window='hann', normalise=True, **kwargs):
     """Make a plot (in the frequency domain) of all channels"""
     
     ymin = kwargs.get('ymin', -160) #dB
     
     freq, mag = self.fft(window=window, normalise=normalise)
     
     fig_id = 1
     plt.figure(fig_id)
     
     #plt.semilogx(freq, mag, **kwargs)   # plots all channel directly
     plt.hold(True)
     for ch in range(self.ch):
         plt.semilogx(freq, mag[:,ch], label='ch%2i' %(ch+1))
     plt.hold(False)
     
     plt.xlim(xmin=1)    # we're not interested in freqs. below 1 Hz
     plt.ylim(ymin=ymin)
     
     plt.xlabel('Frequency [Hz]')
     plt.ylabel('Magnitude [dB]')
     
     plt.legend(loc='best')
     plt.grid(True)
     
     if plotname is None:
         plt.show()
     else:
         plt.savefig(plotname)
         plt.close(fig_id)

def make_plot():
    def u(t):
        return sym.exp(-a*t)

    a, t, dt, p = sym.symbols('a t dt p')
    dudt = sym.diff(u(t), t)

    from numpy import logspace, exp
    from matplotlib.pyplot import (
        semilogx, legend, show, loglog, savefig)

    # Map operator function name to logical names
    operator2name = dict(
        D_f='forward', D_b='backward', D_c='central')
    legends = []
    for operator in D_f, D_b, D_c:
        E = operator(u, dt, t)/dudt
        # Expand, set p=a*dt, simplify
        E = sym.expand(E)
        E = E.subs(a*dt, p)
        E = sym.simplify(E)
        print '%s E:' % operator2name[operator.__name__], E
        print 'Taylor series:', E.series(p, 0, 3)
        latex_expr = sym.latex(E)

        E = sym.lambdify([p], E, modules='numpy')
        p_values = logspace(-6, -0.5, 101)
        y = E(p_values)
        semilogx(p_values, y)
        legends.append(operator2name[operator.__name__] +
                       ': $' + latex_expr + '$')
    legend(legends, loc='lower left')
    savefig('tmp.png'); savefig('tmp.pdf')
    show()

 def plotActiveFilters(self, colour='r'):
     if self.initialized:
         plt.semilogx(self.freqs, 10 *
                      np.log10(self.wActive + 1e-10), colour)
         plt.xlabel('Frequency, Hz')
         plt.ylabel('Response, dB')
         plt.show()

    def test_frequency(self):
        import scipy.signal as ss
        nsteps = 10000
        dt = 0.1

        t,w = self.regular_integrator.run(self.regular_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.regular_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_regular.png"))

        # ----------------------------------------------------------------------
        t,w = self.chaotic_integrator.run(self.chaotic_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.chaotic_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_chaotic.png"))

def cmc(cmc_scores, logx = True, **kwargs):
  """Plots the (cumulative) match characteristics curve and returns the maximum rank.

  This function plots a CMC curve using the given CMC scores, which can be read from the our score files using the :py:func:`bob.measure.load.cmc_four_column` or :py:func:`bob.measure.load.cmc_five_column` methods.
  The structure of the ``cmc_scores`` parameter is relatively complex.
  It contains a list of pairs of lists.
  For each probe object, a pair of list negative and positive scores is required.

  **Parameters:**

  ``cmc_scores`` : [(array_like(1D, float), array_like(1D, float))]
    See :py:func:`bob.measure.cmc`

  ``logx`` : bool
    Plot the rank axis in logarithmic scale using :py:func:`matplotlib.pyplot.semilogx` or in linear scale using :py:func:`matplotlib.pyplot.plot`? (Default: ``True``)

  ``kwargs`` : keyword arguments
    Extra plotting parameters, which are passed directly to :py:func:`matplotlib.pyplot.plot` or :py:func:`matplotlib.pyplot.semilogx`.

  **Returns:**

  The number of classes (clients) in the given scores.
  """
  from matplotlib import pyplot
  from . import cmc as calc

  out = calc(cmc_scores)

  if logx:
    pyplot.semilogx(range(1, len(out)+1), out * 100, **kwargs)
  else:
    pyplot.plot(range(1, len(out)+1), out * 100, **kwargs)

  return len(out)

def bode(G, w1, w2, label='Figure', margin=False):
    """ 
    Shows the bode plot for a plant model
    
    Parameters
    ----------
    G : tf
        plant transfer function
    w1 : real
        start frequency
    w2 : real
        end frequency
    label : string
        title for the figure (optional)
    margin : boolean
        show the cross over frequencies on the plot (optional)        
          
    Returns
    -------
    GM : array containing a real number      
        gain margin
    PM : array containing a real number           
        phase margin         
    """

    GM, PM, wc, w_180 = margins(G)

    # plotting of Bode plot and with corresponding frequencies for PM and GM
#    if ((w2 < numpy.log(w_180)) and margin):
#        w2 = numpy.log(w_180)  
    w = numpy.logspace(w1, w2, 1000)
    s = 1j*w

    plt.figure(label)
    # Magnitude of G(jw)
    plt.subplot(211)
    gains = numpy.abs(G(s))
    plt.loglog(w, gains)
    if margin:
        plt.axvline(w_180, color='black')
        plt.text(w_180, numpy.average([numpy.max(gains), numpy.min(gains)]), r'$\angle$G(jw) = -180$\degree$')
    plt.axhline(1., color='red')
    plt.grid()
    plt.ylabel('Magnitude')

    # Phase of G(jw)
    plt.subplot(212)
    phaseangle = phase(G(s), deg=True)
    plt.semilogx(w, phaseangle)
    if margin:
        plt.axvline(wc, color='black')
        plt.text(wc, numpy.average([numpy.max(phaseangle), numpy.min(phaseangle)]), '|G(jw)| = 1')
    plt.axhline(-180., color='red')
    plt.grid()
    plt.ylabel('Phase')
    plt.xlabel('Frequency [rad/unit time]')
    
    plt.show()

    return GM, PM

def plot_validation_curve(estimator, title, X, y, param_name, param_range,
							cv=10, scoring='accuracy', n_jobs=2):
	from sklearn.learning_curve import validation_curve
	train_scores, test_scores = validation_curve(
	    estimator, X, y, param_name, param_range,
	    cv=cv, scoring=scoring, n_jobs=n_jobs)
	train_scores_mean = np.mean(train_scores, axis=1)
	train_scores_std = np.std(train_scores, axis=1)
	test_scores_mean = np.mean(test_scores, axis=1)
	test_scores_std = np.std(test_scores, axis=1)

	plt.figure()
	plt.title(title)
	plt.xlabel(param_name)
	plt.ylabel("Score")
	plt.ylim(0.0, 1.1)
	plt.semilogx(param_range, train_scores_mean, label="Training score", color="r")
	plt.fill_between(param_range, train_scores_mean - train_scores_std,
	                 train_scores_mean + train_scores_std, alpha=0.2, color="r")
	plt.semilogx(param_range, test_scores_mean, label="Cross-validation score",
	             color="g")
	plt.fill_between(param_range, test_scores_mean - test_scores_std,
	                 test_scores_mean + test_scores_std, alpha=0.2, color="g")
	plt.legend(loc="best")
	plt.show()

def main():
	# define paths
	dir			= 'media/results'
	sizefile	= os.path.join(dir,'sizes.csv')
	timefile	= os.path.join(dir,'times.csv')

	# load data
	sizedata = np.loadtxt(open(sizefile,"rb"),delimiter=",")
	timedata = np.loadtxt(open(timefile,"rb"),delimiter=",")

	timedata = 1000*timedata
	
	# calculate mean and standard deviation
	mean_time = timedata.mean(axis=0);
	stdev_time = timedata.std(axis=0);
	
	size	= [elem[0]*elem[1] for elem in sizedata]
	
	t_max	= mean_time+stdev_time
	t_min	= mean_time-stdev_time
	# show data
	plt.figure(facecolor='white')
	#plt.fill_between(size, t_min, t_max, facecolor='blue', interpolate=True,alpha = 0.5)
	plt.fill_between(size, t_min, t_max, facecolor='#8080FF',color='blue',linewidth=0.5, interpolate=True)
	plt.semilogx(size,mean_time, color='blue',linewidth=1)
	plt.title('Computational time vs number of pixels')
	plt.ylim([0,1])
	plt.xlim([min(size),max(size)])
	plt.xlabel("Size [pixels]")
	plt.ylabel("Time [ms]")
	#plt.show()

	# save plot
	plt.savefig('media/plots/size-time.eps',format='eps', facecolor="white")

def plotRI(RI, title='', filename=''):
    
    global subdirectory

    directory = pinspec.get_output_directory() + subdirectory

    # Make directory if it does not exist
    if not os.path.exists(directory):
            os.makedirs(directory)

    # Plot Resonance Integrals
    fig = plt.figure()
    bins = RI.bin_edges
    plt.semilogx(bins[0:-1], RI.RIs, drawstyle='steps-post')
    plt.xlabel('Energy [eV]')
    plt.ylabel('RI')
    plt.grid()

    if title is '':
        plt.title('Resonance Integrals')
    else:
        plt.title(title.title() + ' Resonance Integrals')

    if filename is '':
        filename = directory + 'RI.png'
    else:
        filename = directory + filename.replace(' ', '-').lower() +'-RI.png'

    plt.savefig(filename)

def hypers_demo():
    steps, mc = 500, 20
    ssm = UNGM()  # initialize UNGM model
    x, z = ssm.simulate(steps, mc_sims=mc)  # generate some data
    lscale = [1e-3, 3e-3, 1e-2, 3e-2, 1e-1, 3e-1, 1, 3, 1e1, 3e1]  # , 1e2, 3e2]
    sigmas_ut = Unscented.unit_sigma_points(ssm.xD, kappa=0.0)
    mean_f, cov_f = np.zeros((ssm.xD, steps, mc, len(lscale))), np.zeros((ssm.xD, ssm.xD, steps, mc, len(lscale)))
    for iel, el in enumerate(lscale):
        # initialize BHKF with current lenghtscale
        f = GPQuadKalman(ssm, usp_dyn=sigmas_ut, usp_meas=sigmas_ut,
                         hyp_dyn={'sig_var': 1.0, 'lengthscale': el * np.ones(ssm.xD, ), 'noise_var': 1e-8},
                         hyp_meas={'sig_var': 1.0, 'lengthscale': el * np.ones(ssm.xD, ), 'noise_var': 1e-8})
        # filtering
        for s in range(mc):
            mean_f[..., s, iel], cov_f[..., s, iel] = f.forward_pass(z[..., s])

    # compute average (over MC sims) RMSE and NCI
    rmseVsEl = rmse(x, mean_f).mean(axis=1)
    nciVsEl = nci(x, mean_f, cov_f).mean(axis=1)
    # plot influence of changing lengthscale on the RMSE and NCI filter performance
    plt.figure()
    plt.semilogx(lscale, rmseVsEl.squeeze(), color='k', ls='-', lw=2, marker='o', label='RMSE')
    plt.semilogx(lscale, nciVsEl.squeeze(), color='k', ls='--', lw=2, marker='o', label='NCI')
    plt.grid(True)
    plt.legend()
    plt.show()
    return lscale, rmseVsEl, nciVsEl