def InvokeMap(coastfile='/media/sda4/map-data/aust-coast-noaa-2000000-1.dat',
		    lllon=80,
		    urlon=166,
		    lllat=-47,
		    urlat=-9,
		    draw_map=True):
    global PYLIB_PATH

    map = Basemap(projection='cyl',
			llcrnrlon=lllon,
			urcrnrlon=urlon,
			llcrnrlat=lllat,
			urcrnrlat=urlat,
			#lat_ts=-35,
			lat_0=-35,
			lon_0=120,
			resolution='l',
			area_thresh=1000.)


    try: 
	coast = p.load(coastfile)
	coast = p.load(coastfile)
	coast_x,coast_y = map(coast[:,0],coast[:,1])
	p.plot(coast_x,coast_y,color='black')    
    except IOError:
	map.drawcoastlines()

    map.drawmapboundary()
    map.drawmeridians(p.arange(0,360,10),labels=[0,0,1,0])
    map.drawparallels(p.arange(-90,0,10),labels=[1,0,0,0])

    return map

def fBM_nd(dims, H, return_mat = False, use_eig_ev = True):
    """
    creates fractional Brownian motion
    parameters: dims is a tuple of the shape of the sample path (nxd); 
                H: Hurst exponent
    this is the slow version of fBM. It might, however, be more precise than
    fBM, however - sometimes, the matrix square root has a problem, which might
    induce inaccuracy    
    use_eig_ev: use eigenvalue decomposition for matrix square root computation
    (faster)
    """
    n = dims[0]
    d = dims[1]
    Gamma = zeros((n,n))
    print ('building ...\n')
    for t in arange(n):
        for s in arange(n):
            Gamma[t,s] = .5*((s+1)**(2.*H) + (t+1)**(2.*H) - abs(t-s)**(2.*H))
    print('rooting ...\n')    
    if use_eig_ev:
        ev,ew = eig(Gamma.real)
        Sigma = dot(ew, dot(diag(sqrt(ev)),ew.T) )
    else:
        Sigma = sqrtm(Gamma)
    if return_mat:
        return Sigma
    v = randn(n,d)
    return dot(Sigma,v)

	def plot(self,title='',include_baseline=False,equal_aspect=True):
		""" Method that generates a plot of the ROC curve 
			Parameters:
				title: Title of the chart
				include_baseline: Add the baseline plot line if it's True
				equal_aspect: Aspects to be equal for all plot
		"""
		
		pylab.clf()
		pylab.plot([x[0] for x in self.derived_points], [y[1] for y in self.derived_points], self.linestyle)
		if include_baseline:
			pylab.plot([0.0,1.0], [0.0,1.0],'k-.')
		pylab.ylim((0,1))
		pylab.xlim((0,1))
		pylab.xticks(pylab.arange(0,1.1,.1))
		pylab.yticks(pylab.arange(0,1.1,.1))
		pylab.grid(True)
		if equal_aspect:
			cax = pylab.gca()
			cax.set_aspect('equal')
		pylab.xlabel('1 - Specificity')
		pylab.ylabel('Sensitivity')
		pylab.title(title)
		
		pylab.show()

def visualize_labeled_z():
	x_batch, label_batch = sample_x_and_label_from_data_distribution(len(dataset), sequential=True)
	z_batch = gen(x_batch, test=True)
	z_batch = z_batch.data
	# if z_batch[0].shape[0] != 2:
	# 	raise Exception("Latent code vector dimension must be 2.")

	fig = pylab.gcf()
	fig.set_size_inches(20.0, 16.0)
	pylab.clf()
	colors = ["#2103c8", "#0e960e", "#e40402","#05aaa8","#ac02ab","#aba808","#151515","#94a169", "#bec9cd", "#6a6551"]
	for n in xrange(z_batch.shape[0]):
		result = pylab.scatter(z_batch[n, 0], z_batch[n, 1], c=colors[label_batch[n]], s=40, marker="o", edgecolors='none')

	classes = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"]
	recs = []
	for i in range(0, len(colors)):
	    recs.append(mpatches.Rectangle((0, 0), 1, 1, fc=colors[i]))

	ax = pylab.subplot(111)
	box = ax.get_position()
	ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
	ax.legend(recs, classes, loc="center left", bbox_to_anchor=(1.1, 0.5))
	pylab.xticks(pylab.arange(-4, 5))
	pylab.yticks(pylab.arange(-4, 5))
	pylab.xlabel("z1")
	pylab.ylabel("z2")
	pylab.savefig("%s/labeled_z.png" % args.visualization_dir)

def makecoords(expr, BL,UR,gridspacing=.1):
    from pylab import arange, array
    X=arange(BL[0],UR[0],gridspacing)
    Y=arange(BL[1],UR[1],gridspacing)
    fn=deriv(expr)
    Z=array([[fn(x=x,y=y) for x in X] for y in Y])
    return X,Y,Z

 def _pvoc2(self, X_hat, Phi_hat=None, R=None):
     """
     ::
       alternate (batch) implementation of phase vocoder - time-stretch
       inputs:
         X_hat - estimate of signal magnitude
         [Phi_hat] - estimate of signal phase
         [R] - resynthesis hop ratio
       output:
         updates self.X_hat with modified complex spectrum
     """
     N, W, H = self.nfft, self.wfft, self.nhop
     R = 1.0 if R is None else R
     dphi = P.atleast_2d((2*P.pi * H * P.arange(N/2+1)) / N).T
     print "Phase Vocoder Resynthesis...", N, W, H, R
     A = P.angle(self.STFT) if Phi_hat is None else Phi_hat
     U = P.diff(A,1) - dphi
     U = U - P.np.round(U/(2*P.pi))*2*P.pi
     t = P.arange(0,n_cols,R)
     tf = t - P.floor(t)
     phs = P.c_[A[:,0], U] 
     phs += U[:,idx[1]] + dphi # Problem, what is idx ?
     Xh = (1-tf)*Xh[:-1] + tf*Xh[1:]
     Xh *= P.exp( 1j * phs)
     self.X_hat = Xh

def redshift():
    """
    Evolution with redshift of matter power spectrum
    """
    zs = M.arange(0.,5.,2.)

    for z in zs:
        print z
        c = pt.Camb(hubble = 70., ombh2 = 0.05*(0.7)**2, omch2 = 0.25*(0.7)**2,transfer_redshift = [z])
        c.run()
        ps = pt.PowerSpectrum(c.cp)
        c.kextend(-10,60) #To ensure accurate sigma(r) -- if it doesn't, a warning will ensue
        pt.normalizePk(c,0.8*ps.d1(z)/ps.d1(0.)) #sigma_8 at redshift z
        
        #Sheth-Tormen
        h = halo.HaloModel(c,st_big_a = 0., st_little_a=0.707, stq = 0.3, k = 10**M.arange(-2,2.01,0.2),massdivsperdex=5)

        h.pmm = halo.getHaloPknl(c,h)
        M.loglog(h.k, h.pmm, label='z='+str(z))
        M.loglog(h.k, h.pk,'k:',label='linear')

        cp_halofit = c.cp
        cp_halofit['do_nonlinear'] = 1 # Halofit (Smith et al) fit
        chalofit = pt.Camb(cambParam=cp_halofit)
        chalofit.run()
        wheretoplot = N.where(chalofit.k > 1e-2)[0]
        M.loglog(chalofit.k[wheretoplot[::10]],chalofit.pk[wheretoplot[::10]],'--',label='halofit')

    M.legend()
    M.show()

def adjuster(file_name, threshold=0.9):
    record = collections.OrderedDict()
    total = 0
    for data in open(file_name).readlines():
        curr = data.strip().split('\t')
        record[curr[0]] = int(curr[1])
        total += int(curr[1])

    curr_count = 0
    curr_list = list()
    ratio_list = list()
    for item in record.items():
        curr_count += item[1]
        curr_list.append(item[0])
        ratio_list.append(curr_count / float(total))
        if curr_count / float(total) >= threshold:
            break

    x = pylab.arange(1, len(ratio_list) + 1, 1)
    plt.plot(x, ratio_list)
    plt.show()

    dx = 1
    dy = diff(ratio_list) / dx
    print len(dy)
    print len(x)
    x = pylab.arange(1, len(ratio_list), 1)
    plt.plot(x, dy)
    plt.show()

    print dy


    return curr_list

def bars(data, figname, fig = None):
    if fig == None:
        P.ioff()
        fig = P.figure()
        ax  = fig.add_subplot(111)
        save = True
    else:
        ax = fig
        fig = ax.get_figure()
        save = False

    ind = P.arange(len(data[0][0]))  # the x locations for the groups
    width = 0.20       # the width of the bars
    #colors = ['r','b','g','y']

    bar_groups = [  ax.bar( ind+width*i, grp[0], width,color=colors[i],\
                            yerr=grp[1], ecolor='k')\
                    for i,grp in enumerate(data)]
    barsLegend = tuple([grp[0] for grp in bar_groups])
    etiquetas = [str(i) for i in ctoa_list]
    P.legend( barsLegend, ctoa_list, shadow=True)
    #ax.set_title('Incremento de RT frente a CTD en distintos CTOAS', font, fontsize=12)
    #ax.set_xlabel('CTD distance (cm)', font)
    #ax.set_ylabel('RT increment (ms)', font)
    ax.set_xticks(ind+width)
    ax.set_xticklabels( ctd_names[:-1] )
    ax.set_yticks(P.arange(-50,50,5))
    ax.xaxis.set_ticks_position("bottom")
    ax.yaxis.set_ticks_position("left")
    ax.set_xlim(-width,len(ind))
    ax.set_ylim(-60,60)
    if save == True:
        fig.savefig(figname+'bar'+graphext,dpi=dpi)
        P.close(fig)

def test_expert_model_level_value():
    d = data.ModelData()
    ages=pl.arange(101)

    # create model with no priors
    vars = {}
    vars.update(age_pattern.age_pattern('test', ages, knots=pl.arange(0,101,5), smoothing=.01))
    vars.update(expert_prior_model.level_constraints('test', {}, vars['mu_age'], ages))

    # fit model
    m = mc.MCMC(vars)
    m.sample(3)


    # create model with expert priors
    parameters = {}
    parameters['level_value'] = dict(value=.1, age_below=15, age_above=95)
    parameters['level_bound'] = dict(upper=.01, lower=.001)
    vars = {}
    vars.update(age_pattern.age_pattern('test', ages, knots=pl.arange(0,101,5), smoothing=.01))
    vars.update(expert_prior_model.level_constraints('test', parameters, vars['mu_age'], ages))

    # fit model
    m = mc.MCMC(vars)
    m.sample(3)

def drawROC(points,zeTitle,zeFilename,visible,show_fig,save_fig=True,
            special_point=None,special_value=None,special_label=None):
    AUC=computeAUC(points)
    import pylab

    pylab.clf()
    pylab.grid(color='#aaaaaa', linestyle='-', linewidth=1,alpha=0.5)

    pylab.plot([x[0] for x in points], [y[1] for y in points], '-', linewidth=3,color="#000088",zorder=3)
    pylab.fill_between([x[0] for x in points], [y[1] for y in points],0,color='0.9')
    pylab.plot([0.0,1.0], [0.0, 1.0], '-',color="#AAAAAA")

    pylab.ylim((-0.01,1.01))
    pylab.xlim((-0.01,1.01))
    pylab.xticks(pylab.arange(0,1.1,.1))
    pylab.yticks(pylab.arange(0,1.1,.1))
    pylab.grid(True)

    ax=pylab.gca()
    r = pylab.Rectangle((0,0), 1, 1, edgecolor='#444444', facecolor='none',zorder=1)
    ax.add_patch(r)
    [spine.set_visible(False) for spine in ax.spines.values()]

    if len(points)<10:
      for i in range(1,len(points)-1):
        pylab.plot(points[i][0],points[i][1],'o',color="#000066",zorder=6)

    pylab.xlabel('False positive rate')
    pylab.ylabel('True positive rate')

    if special_point is not None:
        pylab.plot(special_point[0],special_point[1],'o',color="#DD9999",zorder=6)
        if special_value is not None:
            pylab.text(special_point[0]+0.01,special_point[1]-0.01, special_value,
                       {'color' : '#DD5555', 'fontsize' : 10},
                       horizontalalignment = 'left',
                       verticalalignment = 'top',
                       rotation = 0,
                       clip_on = False)
    if special_label is not None:
        if special_label!="":
            labels=[special_label]
            colors=['#DD9999']
            circles=[pylab.Circle((0, 0), 1, fc=colors[0])]
            legend_location = 'lower right'
            pylab.legend(circles, labels, loc=legend_location)

    pylab.text(0.5, 0.3,'AUC=%f'%AUC,
     horizontalalignment='center',
     verticalalignment='center',
     fontsize=18)

    pylab.title(zeTitle)

    if save_fig:
        pylab.savefig(zeFilename,dpi=300)
        print("\n result in "+zeFilename)

    if show_fig:
        pylab.show()

def plotcp(nelem, chord, thick, u):
    # Geometry
    xnode = geomwing(nelem, chord, thick)
    TEat1 = 1
    dt = 1.0
    # Boundary conditions
    chisrf = bcondvel(xnode, u)
    # Integral equation solution
    B, C = srfmatbc(xnode)
    phisrf = solvephi(B, C, chisrf)
    cpoint = collocation(xnode)
    # Pressure figure
    spl = subplot(111)
    spl.set_aspect("equal", "box")
    plotgeom(xnode)
    # Pressure calculation and plotting
    U = u.reshape((1, 2))
    cp = calccp(xnode, TEat1, dt, U, phisrf, chisrf)
    cl = calchalfcl(xnode, TEat1, dt, U, phisrf, chisrf)
    print cl
    plot(cpoint[: nelem / 2 + 1, 0], cp[: nelem / 2 + 1, 0])
    title(r"Stationary pressure coefficient, $c_L = %9.3f$" % cl[0])
    ylabel(r"$c_p$", size=18)
    xlabel(r"$x/c$", size=18)
    xticks(arange(-0.2, 1.3, 0.2))
    yticks(arange(-0.8, 1.3, 0.2))

def multi_plot_search_on_eval_metrics(roc_search_em, score_thresholds, labels, metrics, line_styles, query_id='Query'):
    import pylab
    
    pylab.clf()
    pylab.xlim((0, 1))
    pylab.xticks(pylab.arange(0,1.1,.1))
    pylab.ylim((-0.5, 1))
    pylab.yticks(pylab.arange(0,1.1,.1))
    pylab.grid(True)
    pylab.xlabel('Score Thresholds')
    
    for iix, score_key in enumerate(metrics): 
        for ix, eval_dict in enumerate(roc_search_em):
            pylab.plot(score_thresholds, eval_dict[score_key], 
                       linewidth=2, label=labels[ix] + '(%s)' % score_key, 
                       color=METRIC_COLORS[ix], linestyle=line_styles[iix])
        
    pylab.ylabel(' '.join(METRICS_DICT[score_key] for score_key in metrics))
    pylab.title(' '.join(METRICS_DICT[score_key] for score_key in metrics))
    
    if labels: pylab.legend(loc='lower left', prop={'size':9})

    eval_file_name = '%s_%s.png' % (query_id, '_'.join(METRICS_DICT[score_key] for score_key in metrics))
    pylab.savefig(eval_file_name, dpi=300, bbox_inches='tight', pad_inches=0.1)

    print 'Saved figure: ', eval_file_name
        
    print '----------------------------------------------------------------------------------'

    def fresnelConvolutionTransform(self,d) :
        # make intensity distribution
        i2 = Intensity2D(self.nx,self.startx,self.endx,
                         self.ny,self.starty,self.endy,
                         self.wl)       

        # FT on inital distribution 
        u1ft = pl.fft2(self.i)

        # 2d convolution kernel
        k = 2*pl.pi/i2.wl
        
        # make spatial frequency matrix
        maxsfx = 2*pl.pi/self.dx
        maxsfy = 2*pl.pi/self.dy
        
        dsfx = 2*maxsfx/(self.nx)
        dsfy = 2*maxsfy/(self.ny)
        
        self.sfx = pl.arange(-maxsfx/2,maxsfx/2+1e-15,dsfx/2)
        self.sfy = pl.arange(-maxsfy/2,maxsfy/2+1e-15,dsfy/2)

        [self.sfxgrid, self.sfygrid] = pl.fftshift(pl.meshgrid(self.sfx,self.sfy))
                
        # make convolution kernel 
        kern = pl.exp(1j*d*(self.sfxgrid**2+self.sfygrid**2)/(2*k))
        
        # apply convolution kernel and invert
        i2.i = pl.ifft2(kern*u1ft) 

        return i2

    def _make_log_freq_map(self):
        """
        ::

            For the given ncoef (bands-per-octave) and nfft, calculate the center frequencies
            and bandwidths of linear and log-scaled frequency axes for a constant-Q transform.
        """
        fp = self.feature_params
        bpo = float(self.nbpo) # Bands per octave
        self._fftN = float(self.nfft)
        hi_edge = float( self.hi )
        lo_edge = float( self.lo )
        f_ratio = 2.0**( 1.0 / bpo ) # Constant-Q bandwidth
        self._cqtN = float( P.floor(P.log(hi_edge/lo_edge)/P.log(f_ratio)) )
        self._dctN = self._cqtN
        self._outN = float(self.nfft/2+1)
        if self._cqtN<1: print "warning: cqtN not positive definite"
        mxnorm = P.empty(self._cqtN) # Normalization coefficients        
        fftfrqs = self._fftfrqs #P.array([i * self.sample_rate / float(self._fftN) for i in P.arange(self._outN)])
        logfrqs=P.array([lo_edge * P.exp(P.log(2.0)*i/bpo) for i in P.arange(self._cqtN)])
        logfbws=P.array([max(logfrqs[i] * (f_ratio - 1.0), self.sample_rate / float(self._fftN)) 
                         for i in P.arange(self._cqtN)])
        #self._fftfrqs = fftfrqs
        self._logfrqs = logfrqs
        self._logfbws = logfbws
        self._make_cqt()

def plot_multiple_roc(rocList,title='',labels=None, include_baseline=False, equal_aspect=True,plot_average =True):
    pylab.figure()
    pylab.clf()
    pylab.ylim((0,1))
    pylab.xlim((0,1))
    pylab.xticks(pylab.arange(0,1.1,.1))
    pylab.yticks(pylab.arange(0,1.1,.1))
    pylab.grid(True)
    if equal_aspect:
        cax = pylab.gca()
        cax.set_aspect('equal')
    pylab.xlabel("1 - Specificity")
    pylab.ylabel("Sensitivity")
    pylab.title(title)
    if not labels:
        labels = [ '' for x in rocList]
    _remove_duplicate_styles(rocList)
    for ix, r in enumerate(rocList):
        pylab.plot([x[0] for x in r.derived_points], [y[1] for y in r.derived_points], r.linestyle, linewidth=1, label=labels[ix])
    if include_baseline:
        pylab.plot([0.0,1.0], [0.0, 1.0], 'k-', label= 'random')
    if labels:
        pylab.legend(loc='lower right')
    if plot_average:
        stepsize=.1
        av = multiple_roc_average(rocList,binstep=stepsize)
        x=np.arange(stepsize/2,1,stepsize)
        pylab.plot(x,[np.average(a) for a in av])

 def _pvoc(self, X_hat, Phi_hat=None, R=None):
     """
     ::
       a phase vocoder - time-stretch
       inputs:
         X_hat - estimate of signal magnitude
         [Phi_hat] - estimate of signal phase
         [R] - resynthesis hop ratio
       output:
         updates self.X_hat with modified complex spectrum
     """
     N = self.nfft
     W = self.wfft
     H = self.nhop
     R = 1.0 if R is None else R
     dphi = (2*P.pi * H * P.arange(N/2+1)) / N
     print "Phase Vocoder Resynthesis...", N, W, H, R
     A = P.angle(self.STFT) if Phi_hat is None else Phi_hat
     phs = A[:,0]
     self.X_hat = []
     n_cols = X_hat.shape[1]
     t = 0
     while P.floor(t) < n_cols:
         tf = t - P.floor(t)            
         idx = P.arange(2)+int(P.floor(t))
         idx[1] = n_cols-1 if t >= n_cols-1 else idx[1]
         Xh = X_hat[:,idx]
         Xh = (1-tf)*Xh[:,0] + tf*Xh[:,1]
         self.X_hat.append(Xh*P.exp( 1j * phs))
         U = A[:,idx[1]] - A[:,idx[0]] - dphi
         U = U - P.np.round(U/(2*P.pi))*2*P.pi
         phs += (U + dphi)
         t += P.randn()*P.sqrt(PVOC_VAR*R) + R # 10% variance
     self.X_hat = P.np.array(self.X_hat).T

def main():
    # Create the grid
    x = arange(-100, 101)
    y = arange(-100, 101)

    # Create the meshgrid
    Y, X = meshgrid(x, y)
    A = 1
    B = 2
    V = 6*pi / 201
    W = 4*pi / 201
    F = A*sin(V*X) + B*cos(W*Y)
    Fx = V*A*cos(V*X)
    Fy = W*B*-sin(W*Y)

    # Show the images
    show_image(F)
    show_image(Fx)
    show_image(Fy)

    # Create the grid for the quivers
    xs = arange(-100, 101, 10)
    ys = arange(-100, 101, 10)

    # Here we determine the direction of the quivers
    Ys, Xs = meshgrid(ys, xs)
    FFx = V*A*cos(V*Xs)
    FFy = W*B*-sin(W*Ys)

    # Draw the quivers and the image
    clf()
    imshow(F, cmap=cm.gray, extent=(-100, 100, -100, 100))
    quiver(ys, xs, -FFy, FFx, color='red')
    show()

def length_stats_chart(path, prefixes, sortby=1):
  stats = []
  for prefix in prefixes:
    med, m,s = length_stats(prefix)
    stats.append((prefix,med,m,s))

  stats.sort(key=operator.itemgetter(sortby))
  prefixes, med_list, mean_list, std_list = zip(*stats)

  blockSize = 8 
  ind = p.arange(0, blockSize*len(prefixes), blockSize) # y location for groups
  height = 3 # bar height 

  p3 = p.barh(ind, std_list, 2   * height, color = 'b', linewidth = 0)
  p2 = p.barh(ind, med_list, height, color = 'g', linewidth = 0)
  p1 = p.barh(ind+height, mean_list, height, color = 'r', linewidth = 0)
  
  p.ylim(-height, len(prefixes) * blockSize)
  yfontprop = FontProperties(size=4)
  xfontprop = FontProperties(size='smaller')
  p.xlabel('Unicode Codepoints')
  p.ylabel('Language Code')
  p.title('Descriptive Statistics for Document Lengths')
  p.gca().yaxis.tick_left()
  p.yticks(ind+height, prefixes, fontproperties = yfontprop)
  xmin, xmax = p.xlim()
  p.xticks( p.arange(xmin,xmax,1000),fontproperties = xfontprop)
  p.gca().xaxis.grid(linestyle = '-', linewidth=0.15)
  p.legend((p1[0], p2[0], p3[0]), ('Mean','Median','Standard Deviation'), prop = xfontprop, loc = 'lower right' )

  p.savefig(path, dpi=300)
  p.close()
  p.clf()

    def plot2d(self,key1='a',key2='b',scale1=0.1,scale2=0.1):
        """
	under development
	plot the least_squared around the minimum
	assuming a quadratic approximation given by the hessian
	"""
	import pylab
	key1,key2=key2,key1
	keys=self.last_variables.keys()
	keys.sort()
        nv=len(self.last_variables)
	for k in range(len(keys)):
	    if keys[k]==key1: i,vi=k, self.last_variables[key1]
	    if keys[k]==key2: j,vj=k, self.last_variables[key2]
	ax=pylab.arange(vi*(1.0-10*scale1),vi*(1.0+12*scale1),vi*scale1)[:21]
	ay=pylab.arange(vj*(1.0-10*scale2),vj*(1.0+12*scale2),vj*scale2)[:21]

        dv=[0.0]*nv
        z=[[0 for ik in ax] for jk in ay]
        i0=0
        for x in ax:
            j0=0
            for y in ay:
               dv[i],dv[j]=x-vi,y-vj
               least_squares=0.0
               for ik in range(nv):
                   for jk in range(nv): 
                       least_squares+=0.5*dv[ik]*dv[jk]*self.last_hessian[ik][jk]
               z[i0][j0]=least_squares-self.last_least_squares
               j0+=1
            i0+=1
        pylab.contour(z,extent=(min(ay),max(ay),min(ax),max(ax)))
        pylab.show()