def main():
    SAMPLE_NUM = 10
    degree = 9
    x, y = sin_wgn_sample(SAMPLE_NUM)
    fig = pylab.figure(1)
    pylab.grid(True)
    pylab.xlabel('x')
    pylab.ylabel('y')
    pylab.axis([-0.1,1.1,-1.5,1.5])

    # sin(x) + noise
    # markeredgewidth mew
    # markeredgecolor mec
    # markerfacecolor mfc

    # markersize      ms
    # linewidth       lw
    # linestyle       ls
    pylab.plot(x, y,'bo',mew=2,mec='b',mfc='none',ms=8)

    # sin(x)
    x2 = linspace(0, 1, 1000)
    pylab.plot(x2,sin(2*x2*pi),'#00FF00',lw=2,label='$y = \sin(x)$')

    # polynomial fit
    reg = exp(-18)
    w = curve_poly_fit(x, y, degree,reg) #w = polyfit(x, y, 3)
    po = poly1d(w)      
    xx = linspace(0, 1, 1000)
    pylab.plot(xx, po(xx),'-r',label='$M = 9, \ln\lambda = -18$',lw=2)
    
    pylab.legend()
    pylab.show()
    fig.savefig("poly_fit9_10_reg.pdf")

def simulationWithDrug():

    """

    Runs simulations and plots graphs for problem 4.
    Instantiates a patient, runs a simulation for 150 timesteps, adds
    guttagonol, and runs the simulation for an additional 150 timesteps.
    total virus population vs. time and guttagonol-resistant virus population
    vs. time are plotted
    """

    maxBirthProb = .1
    clearProb = .05
    resistances = {'guttagonal': False}
    mutProb = .005

    total = [100]
    g = [0]

    badVirus = ResistantVirus(maxBirthProb, clearProb, resistances, mutProb)
    viruses = [badVirus]*total[0]
    maxPop = 1000

    Bob = Patient(viruses, maxPop)

    for i in range(150):
        Bob.update()
        gVirus = 0

        for v in Bob.viruses:
            if v.isResistantTo('guttagonal'):
                gVirus += 1

        #print "g = ", gVirus
        #print "t = ", len(Bob.viruses)
        #print
        g += [gVirus]
        total += [len(Bob.viruses)]

    Bob.addPrescription('guttagonal')

    for i in range(150):
        Bob.update()
        gVirus = 0

        for v in Bob.viruses:
            if v.isResistantTo('guttagonal'):
                gVirus += 1

        g += [gVirus]
        total += [len(Bob.viruses)]

    pylab.title("Number of Viruses with Different Resistances to Guttagonal")
    pylab.xlabel("Number of Timesteps")
    pylab.ylabel("Number of Viruses")

    pylab.plot(g, '-r', label = 'Resistant')
    pylab.plot(total, '-b', label = 'Total')
    pylab.legend(loc = 'lower right')
    pylab.show()

def embed_two_dimensions(data, vectorizer, size=10, n_components=5, colormap='YlOrRd'):
    if hasattr(data, '__iter__'):
        iterable = data
    else:
        raise Exception('ERROR: Input must be iterable')
    import itertools
    iterable_1, iterable_2 = itertools.tee(iterable)
    # get labels
    labels = []
    for graph in iterable_2:
        label = graph.graph.get('id', None)
        if label:
            labels.append(label)

    # transform iterable into sparse vectors
    data_matrix = vectorizer.transform(iterable_1)
    # embed high dimensional sparse vectors in 2D
    from sklearn import metrics
    distance_matrix = metrics.pairwise.pairwise_distances(data_matrix)

    from sklearn.manifold import MDS
    feature_map = MDS(n_components=n_components, dissimilarity='precomputed')
    explicit_data_matrix = feature_map.fit_transform(distance_matrix)

    from sklearn.decomposition import TruncatedSVD
    pca = TruncatedSVD(n_components=2)
    low_dimension_data_matrix = pca.fit_transform(explicit_data_matrix)

    plt.figure(figsize=(size, size))
    embed_dat_matrix_two_dimensions(low_dimension_data_matrix, labels=labels, density_colormap=colormap)
    plt.show()

    def Xtest3(self):
        """
        Test from Kate Marvel
        As the following code snippet demonstrates, regridding a
        cdms2.tvariable.TransientVariable instance using regridTool='regrid2' 
        results in a new array that is masked everywhere.  regridTool='esmf' 
        and regridTool='libcf' both work as expected.

        This is similar to the original test but we construct our own 
        uniform grid. This should passes.
        """
        import cdms2 as cdms
        import numpy as np

        filename = cdat_info.get_sampledata_path() + '/clt.nc'
        a=cdms.open(filename)
        data=a('clt')[0,...]

        print data.mask #verify this data is not masked

        GRID = cdms.grid.createUniformGrid(-90.0, 23, 8.0, -180.0, 36, 10.0, order="yx", mask=None)

        test_data=data.regrid(GRID,regridTool='regrid2')

        # check that the mask does not extend everywhere...
        self.assertNotEqual(test_data.mask.sum(), test_data.size)
        
        if PLOT:
            pylab.subplot(2, 1, 1)
            pylab.pcolor(data[...])
            pylab.title('data')
            pylab.subplot(2, 1, 2)
            pylab.pcolor(test_data[...])
            pylab.title('test_data (interpolated data)')
            pylab.show()

def plotear(xi,yi,zi):
    # mask inner circle
    interior1 = sqrt(((xi+1.5)**2) + (yi**2)) < 1.0 
    interior2 = sqrt(((xi-1.5)**2) + (yi**2)) < 1.0
    zi[interior1] = ma.masked
    zi[interior2] = ma.masked
    p.figure(figsize=(16,10))
    pyplot.jet()
    max=2.8
    min=0.4
    steps = 50
    levels=list()
    labels=list()
    for i in range(0,steps):
	levels.append(int((max-min)/steps*100*i)*0.01+min)
    for i in range(0,steps/2):
	labels.append(levels[2*i])
    CSF = p.contourf(xi,yi,zi,levels,norm=colors.LogNorm())
    CS = p.contour(xi,yi,zi,levels, format='%.3f', labelsize='18')
    p.clabel(CS,labels,inline=1,fontsize=9)
    p.title('electrostatic potential of two spherical colloids, R=lambda/3',fontsize=24)
    p.xlabel('z-coordinate (3*lambda)',fontsize=18)
    p.ylabel('radial coordinate r (3*lambda)',fontsize=18)
    # add a vertical bar with the color values
    cbar = p.colorbar(CSF,ticks=labels,format='%.3f')
    cbar.ax.set_ylabel('potential (reduced units)',fontsize=18)
    cbar.add_lines(CS)
    p.show()

def main():
    pylab.ion();
    ind = [0,];
    ldft = [0,];
    lfft = [0,];
    lpfft = [0,]

    # plot a graph Dft vs Fft, lists just support size until 2**9
    for i in range(1, 9, 1):
        t_before = time.clock();
        dsprocessing.dspDft(rand(2**i).tolist());
        dt = time.clock() - t_before;
        ldft.append(dt);
        print ("dft ", 2**i, dt);
        #pylab.plot([2**i,], [time.clock()-t_before,]);
        t_before = time.clock();
        dsprocessing.dspFft(rand(2**i).tolist());
        dt = time.clock() - t_before;
        print ("fft ", 2**i, dt);
        lfft.append(dt);
        #pylab.plot([2**i,], [time.clock()-t_before,]);
        ind.append(2**i);
        # python fft just to compare
        t_before = time.clock();
        pylab.fft(rand(2**i).tolist());
        dt = time.clock() - t_before;
        lpfft.append(dt);

    pylab.plot(ind, ldft);
    pylab.plot(ind, lfft);
    pylab.plot(ind, lpfft);
    pylab.show();
    return [ind, ldft, lfft, lpfft];

def param_set_averages_plot(results):
    averages_ocr = [
        a[1] for a in sorted(
            param_set_averages(results, metric='ocr').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_q = [
        a[1] for a in sorted(
            param_set_averages(results, metric='q').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_mse = [
        a[1] for a in sorted(
            param_set_averages(results, metric='mse').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    fig = plt.figure(figsize=(6, 4))
    # plt.tight_layout()
    plt.plot(averages_ocr, label='OCR', linewidth=2.0)
    plt.plot(averages_q, label='Q', linewidth=2.0)
    plt.plot(averages_mse, label='MSE', linewidth=2.0)
    plt.ylim([0, 1])
    plt.xlabel(u'Paslėptų neuronų skaičius')
    plt.ylabel(u'Vidurinė Q įverčio pokyčio reikšmė')
    plt.grid(True)
    plt.tight_layout()
    plt.legend(loc='lower right')
    plt.show()

def plotHistogram(data, preTime):
    pylab.figure(1)
    pylab.hist(data, bins=10)
    pylab.xlabel("Virus Population At End of Simulation")
    pylab.ylabel("Number of Trials")
    pylab.title("{0} Time Steps Before Treatment Simulation".format(preTime))
    pylab.show()

    def test_mask_LUT(self):
        """
        The masked image has a masked ring around 1.5deg with value -10
        without mask the pixels should be at -10 ; with mask they are at 0
        """
        x1 = self.ai.xrpd_LUT(self.data, 1000)
#        print self.ai._lut_integrator.lut_checksum
        x2 = self.ai.xrpd_LUT(self.data, 1000, mask=self.mask)
#        print self.ai._lut_integrator.lut_checksum
        x3 = self.ai.xrpd_LUT(self.data, 1000, mask=numpy.zeros(shape=self.mask.shape, dtype="uint8"), dummy= -20.0, delta_dummy=19.5)
#        print self.ai._lut_integrator.lut_checksum
        res1 = numpy.interp(1.5, *x1)
        res2 = numpy.interp(1.5, *x2)
        res3 = numpy.interp(1.5, *x3)
        if logger.getEffectiveLevel() == logging.DEBUG:
            pylab.plot(*x1, label="nomask")
            pylab.plot(*x2, label="mask")
            pylab.plot(*x3, label="dummy")
            pylab.legend()
            pylab.show()
            raw_input()

        self.assertAlmostEqual(res1, -10., 1, msg="Without mask the bad pixels are around -10 (got %.4f)" % res1)
        self.assertAlmostEqual(res2, 0., 4, msg="With mask the bad pixels are actually at 0 (got %.4f)" % res2)
        self.assertAlmostEqual(res3, -20., 4, msg="Without mask but dummy=-20 the dummy pixels are actually at -20 (got % .4f)" % res3)

def cmap_plot(cmdLine):

    pylab.figure(figsize=[5,10])
    a=outer(ones(10),arange(0,1,0.01))
    subplots_adjust(top=0.99,bottom=0.00,left=0.01,right=0.8)
    maps=[m for m in cm.datad if not m.endswith("_r")]
    maps.sort()
    l=len(maps)+1
    for i, m in enumerate(maps):
        print m
        subplot(l,1,i+1)
        pylab.setp(pylab.gca(),xticklabels=[],xticks=[],yticklabels=[],yticks=[])
        imshow(a,aspect='auto',cmap=get_cmap(m),origin="lower")
        pylab.text(100.85,0.5,m,fontsize=10)

# render plot

    if cmdLine: 
        pylab.show(block=True)
    else: 
        pylab.ion()
        pylab.plot([])
        pylab.ioff()
	
    status = 1
    return status

def plotAllWarmJumps():
    jumpAddrs = np.array(getAllWarmJumpsAddr()).reshape((8, 18))
    figure()
    pcolor(jumpAddrs)
    for (x, y), v in np.ndenumerate(jumpAddrs):
        text(y + 0.125, x + 0.5, "0x%03x" % v)
    show()

def plot_heatingrate(data_dict, filename, do_show=True):
    pl.figure(201)
    color_list = ['b','r','g','k','y','r','g','b','k','y','r',]
    fmtlist = ['s','d','o','s','d','o','s','d','o','s','d','o']
    result_dict = {}
    for key in data_dict.keys():
        x = data_dict[key][0]
        y = data_dict[key][1][:,0]
        y_err = data_dict[key][1][:,1]

        p0 = np.polyfit(x,y,1)
        fit = LinFit(np.array([x,y,y_err]).transpose(), show_graph=False)
        p1 = [0,0]
        p1[0] = fit.param_dict[0]['Slope'][0]
        p1[1] = fit.param_dict[0]['Offset'][0]
        print fit
        x0 = np.linspace(0,max(x))
        cstr = color_list.pop(0)
        fstr = fmtlist.pop(0)
        lstr = key + " heating: {0:.2f} ph/ms".format((p1[0]*1e3)) 
        pl.errorbar(x/1e3,y,y_err,fmt=fstr + cstr,label=lstr)
        pl.plot(x0/1e3,np.polyval(p0,x0),cstr)
        pl.plot(x0/1e3,np.polyval(p1,x0),cstr)
        result_dict[key] = 1e3*np.array(fit.param_dict[0]['Slope'])
    pl.xlabel('Heating time (ms)')
    pl.ylabel('nbar')
    if do_show:
        pl.legend()
        pl.show()
    if filename != None:
        pl.savefig(filename)
    return result_dict

def evaluate_result(data, target, result):
	assert(data.shape[0] == target.shape[0])
	assert(target.shape[0] == result.shape[0])
	
	correct = np.where( result == target )
	miss 	= np.where( result != target )
	
	class_rate = float(correct[0].shape[0]) / target.shape[0]

	print "Correct classification rate:", class_rate 
	#get the 3s
	mask 			= np.where(target == wanted[0])
	data_3_correct 	= data[np.intersect1d(mask[0],correct[0])]
	data_3_miss	 	= data[np.intersect1d(mask[0],miss[0])]
	#get the 8s
	mask = np.where(target == wanted[1])
	data_8_correct 	= data[np.intersect1d(mask[0],correct[0])]
	data_8_miss	 	= data[np.intersect1d(mask[0],miss[0])]
	#plot
	plot.title("Scatter")
	plot.xlabel("x_0")
	plot.ylabel("x_1")
	size = 20
	plot.scatter(data_3_correct[:,0], data_3_correct[:,1], marker = "x", c = "r", s = size )
	plot.scatter(   data_3_miss[:,0],    data_3_miss[:,1], marker = "x", c = "b", s = size )
	plot.scatter(data_8_correct[:,0], data_8_correct[:,1], marker = "o", c = "r", s = size )
	plot.scatter(   data_8_miss[:,0],    data_8_miss[:,1], marker = "o", c = "b", s = size )
	plot.show()

def _test_graph():
    i = 10000
    x = np.linspace(0,3.7*pi,i)
    y = (0.3*np.sin(x) + np.sin(1.3 * x) + 0.9 * np.sin(4.2 * x) + 0.06 *
    np.random.randn(i))
    y *= -1
    x = range(i)

    _max, _min = peakdetect(y,x,750, 0.30)
    xm = [p[0] for p in _max]
    ym = [p[1] for p in _max]
    xn = [p[0] for p in _min]
    yn = [p[1] for p in _min]

    plot = pylab.plot(x,y)
    pylab.hold(True)
    pylab.plot(xm, ym, 'r+')
    pylab.plot(xn, yn, 'g+')

    _max, _min = peak_det_bad.peakdetect(y, 0.7, x)
    xm = [p[0] for p in _max]
    ym = [p[1] for p in _max]
    xn = [p[0] for p in _min]
    yn = [p[1] for p in _min]
    pylab.plot(xm, ym, 'y*')
    pylab.plot(xn, yn, 'k*')
    pylab.show()

def createPlot(dataY, dataX, ticksX, annotations, axisY, axisX, dostep, doannotate):
    if not ticksX:
        ticksX = dataX
    
    if dostep:
        py.step(dataX, dataY, where='post', linestyle='-', label=axisY) # where=post steps after point
    else:
        py.plot(dataX, dataY, marker='o', ms=5.0, linestyle='-', label=axisY)
    
    if annotations and doannotate:
        for note, x, y in zip(annotations, dataX, dataY):
            py.annotate(note, (x, y), xytext=(2,2), xycoords='data', textcoords='offset points')

    py.xticks(np.arange(1, len(dataX)+1), ticksX, horizontalalignment='left', rotation=30)
    leg = py.legend()
    leg.draggable()
    py.xlabel(axisX)
    py.ylabel('time (s)')

    # Set X axis tick labels as rungs
    #print zip(dataX, dataY)
  
    py.draw()
    py.show()
    
    return

	def plot(self, GRID=10, ax=None, **kwargs):

		NEWAXIS = False
		if ax == None:
			from pylab import figure, show
			fig = figure()
			ax = fig.add_subplot(111)
			NEWAXIS = True

		if "period" in kwargs:
			period = kwargs.pop("period")
		else:
			period = tl.PI2


		phase = tl.PI2*np.arange(GRID)/float(GRID-1)
		phase += phase[1]/2.
		UV = np.asarray([ [self([phase[i], phase[j]]) for i in xrange(GRID)]	# Spalten
								for j in xrange(GRID)])	# Zeilen

		X, Y = np.meshgrid(phase, phase)
		U, V = UV[:, :, 0], UV[:, :, 1]
	
		Q = ax.quiver(period*X/tl.PI2, period*Y/tl.PI2, U, V, units='width')

		for fp in self.fixedPoints:
			fp.plot(axis=ax, period=period)
			
		if NEWAXIS:
			ax.set_xlim(0., period)
			ax.set_ylim(0., period)
			fig.tight_layout()
			show()

		return Q

def generateKineticsModel(reaction, tunneling='', plot=False):
    
    logging.info('Calculating rate coefficient for {0}...'.format(reaction))
    
    if len(reaction.reactants) == 1:
        kunits = 's^-1'
    elif len(reaction.reactants) == 2:
        kunits = 'm^3/(mol*s)'
    elif len(reaction.reactants) == 3:
        kunits = 'm^6/(mol^2*s)'
    else:
        kunits = ''
    
    Tlist = 1000.0/numpy.arange(0.4, 3.35, 0.05)
    klist = reaction.calculateTSTRateCoefficients(Tlist, tunneling)
    arrhenius = Arrhenius().fitToData(Tlist, klist, kunits)
    klist2 = arrhenius.getRateCoefficients(Tlist)
    
    reaction.kinetics = arrhenius
    
    if plot:
        logging.info('Plotting kinetics model for {0}...'.format(reaction))
        import pylab
        pylab.semilogy(1000.0 / Tlist, klist  * reaction.degeneracy, 'ok')
        pylab.semilogy(1000.0 / Tlist, klist2 * reaction.degeneracy, '-k')
        pylab.xlabel('1000 / Temperature (1000/K)')
        pylab.ylabel('Rate coefficient (SI units)')
        pylab.show()

def qinit(x,y):
    """
    Gaussian hump:
    """
    from numpy import where
    from pylab import plot,show
    nxpoints = 2001
    nypoints = 2001
    #import sympy
        #def sech(x):
        #return sympy.cosh(x)**(-1)
    zmin = 1000.0
    zmin1 = 1000
    g = 9.81
    A = 10.0 # wave height
    k = sqrt(3*A/(4*zmin**3))
    z = zeros(shape=shape(y))
    u = zeros(shape=shape(y))
    hu = zeros(shape=shape(y))
    y0 = 180000
    c = sqrt(g*(A+zmin))
    rho = 1.0e3
    for i in range(nxpoints):
        for j in range(nypoints):
            z[i,j] = A*cosh(k*(y[i,j]-y0))**(-2)
            u[i,j] = -sqrt(g/zmin)*z[i,j]
            hu[i,j] =(z[i,j]+zmin)*u[i,j]
    plot(y,u,'-')
    show()
    #ze = -((y+0e0)**2)/10.
    #z = where(ze>-400000., 400.e0*exp(ze), 0.)
    return hu

def blobFinder(img):
    shape = img.shape[0:2]

    scales = numpy.arange(0.5,7.0,0.5)
    nScales = len(scales)

    bank  = numpy.zeros(shape+(nScales,))


    for si,scale in enumerate(scales):
        print scale
        log         = numpy.abs(scale*vigra.filters.laplacianOfGaussian(img,float(scale)))
        log         = numpy.sum(log,axis=2)
        bank[:,:,si]=log

        #f = pylab.figure()
        #for n, iterImg in enumerate([log,img]):
        #    f.add_subplot(1,2, n)  # this line outputs images side-by-side
        #    pylab.imshow(numpy.swapaxes(norm01(iterImg),0,1))
        #pylab.show()


    maxScale = numpy.argmax(bank,axis=2)

    print "mima",maxScale.min(),maxScale.max()
    print maxScale

    f = pylab.figure()
    for n, iterImg in enumerate([maxScale]):
        f.add_subplot(1,1, n)  # this line outputs images side-by-side
        pylab.imshow(numpy.swapaxes(norm01(iterImg),0,1))
    pylab.show()

def plotEventFlop(library, num, eventNames, sizes, times, events, filename = None):
  from pylab import legend, plot, savefig, semilogy, show, title, xlabel, ylabel
  import numpy as np

  arches = sizes.keys()
  bs     = events[arches[0]].keys()[0]
  data   = []
  names  = []
  for event, color in zip(eventNames, ['b', 'g', 'r', 'y']):
    for arch, style in zip(arches, ['-', ':']):
      if event in events[arch][bs]:
        names.append(arch+'-'+str(bs)+' '+event)
        data.append(sizes[arch][bs])
        data.append(1e-3*np.array(events[arch][bs][event])[:,1])
        data.append(color+style)
      else:
        print 'Could not find %s in %s-%d events' % (event, arch, bs)
  semilogy(*data)
  title('Performance on '+library+' Example '+str(num))
  xlabel('Number of Dof')
  ylabel('Computation Rate (GF/s)')
  legend(names, 'upper left', shadow = True)
  if filename is None:
    show()
  else:
    savefig(filename)
  return