def make_grid(self):
     amax = np.amax([np.amax(coor, axis=0) for coor in self.coordinates], axis=0)
     amin = np.amin([np.amin(coor, axis=0) for coor in self.coordinates], axis=0)
     self.x = np.arange(amin[0], amax[0] + self.blocksize, self.blocksize)
     self.y = np.arange(amin[1], amax[1] + self.blocksize, self.blocksize)
     self.z = np.arange(amin[2], amax[2] + self.blocksize, self.blocksize)
     self.values = np.zeros((self.x.shape[0] - 1, self.y.shape[0] - 1, self.z.shape[0] - 1))

def extract_active_site(protein_file, ligand_file, cutoff=4):
  """Extracts a box for the active site."""
  protein_coords = rdkit_util.load_molecule(
      protein_file, add_hydrogens=False)[0]
  ligand_coords = rdkit_util.load_molecule(
      ligand_file, add_hydrogens=True, calc_charges=True)[0]
  num_ligand_atoms = len(ligand_coords)
  num_protein_atoms = len(protein_coords)
  pocket_inds = []
  pocket_atoms = set([])
  for lig_atom_ind in range(num_ligand_atoms):
    lig_atom = ligand_coords[lig_atom_ind]
    for protein_atom_ind in range(num_protein_atoms):
      protein_atom = protein_coords[protein_atom_ind]
      if np.linalg.norm(lig_atom - protein_atom) < cutoff:
        if protein_atom_ind not in pocket_atoms:
          pocket_atoms = pocket_atoms.union(set([protein_atom_ind]))
  # Should be an array of size (n_pocket_atoms, 3)
  pocket_atoms = list(pocket_atoms)
  n_pocket_atoms = len(pocket_atoms)
  pocket_coords = np.zeros((n_pocket_atoms, 3))
  for ind, pocket_ind in enumerate(pocket_atoms):
    pocket_coords[ind] = protein_coords[pocket_ind]

  x_min = int(np.floor(np.amin(pocket_coords[:, 0])))
  x_max = int(np.ceil(np.amax(pocket_coords[:, 0])))
  y_min = int(np.floor(np.amin(pocket_coords[:, 1])))
  y_max = int(np.ceil(np.amax(pocket_coords[:, 1])))
  z_min = int(np.floor(np.amin(pocket_coords[:, 2])))
  z_max = int(np.ceil(np.amax(pocket_coords[:, 2])))
  return (((x_min, x_max), (y_min, y_max), (z_min, z_max)), pocket_atoms,
          pocket_coords)

 def set_data(self, zname, zdata, zcolor):
     if zdata!=None:
         if self.overall_plot_type=="polygon":
            if zname not in self.clts: #plottables['plotted']:#self.pd.list_data():
                clt=PolyCollection(zdata, alpha=0.5, antialiased=True)#, rasterized=False, antialiased=False)
                clt.set_color(colorConverter.to_rgba(zcolor))                
                self.clts[zname]=clt
                self.axe.add_collection(self.clts[zname], autolim=True)
            else:                
                self.clts[zname].set_verts(zdata)
         if self.overall_plot_type=="XY":
             if zname not in self.clts:
                 clt = LineCollection(zdata)#, offsets=offs)
                 clt.set_color(colors)
                 #print dir(clt)
                 self.clts[zname]=clt
                 self.axe.add_collection(self.clts[zname], autolim=True)
                 self.axe.autoscale_view()
             else:
                 self.clts[zname].set_segments(zdata)
         if self.overall_plot_type=="img":
             if zname not in self.clts:
                 axeimg=self.axe.imshow( Magvec, 
                                        vmin=amin(Magvec),
                                        vmax=0.001, #amax(Magvec), 
                                        aspect="auto", origin="lower",
                                 extent=[amin(yoko),amax(yoko), amin(freq),amax(freq)],
                                 #cmap='RdBu'
                                 )
                 self.fig.colorbar(axeimg)

def construct(phi1, phi2, nomod = 0, amp1 =[], amp2=[], eta = 0, ampout= 0): #does ampout need to be there?
	if len(amp1) > 0 or len(amp2) > 0:
		tempshape = phi1.shape
		w = tempshape[1]
		h = tempshape[0]
		if len(amp1) == 0:
			temp1 = np.ones(w)
			temp2 = np.ones(h)
			for r in temp2:
				amp1 += [temp1]
		if len(amp2) == 0:
			temp1 = np.ones(w)
			temp2 = np.ones(h)
			for r in temp2:
				amp2 += [temp1]
		psi1 = amp1 * np.exp(1j*phi1)
		psi2 = amp2 * np.exp(1j*phi2)
		psi = psi1 * psi2
		psi = np.array(psi)
		apsi = abs(psi)
		psi = psi/(np.amax(abs(psi)))
		phi = np.arctan2(sp.real(psi),sp.imag(psi))
		phi -= np.amin(phi)
		phi = phi % (2.*np.pi)
		eta = 2*np.median(abs(psi))
		randarray = np.array([[random.random() for i in range(w)] for j in range(h)])
		shape = (abs(psi) >= (eta*randarray)) 
		index = np.where(shape == False)
		phi[index] = 0
		ampout = abs(psi)
	else:
		phi = phi1 + phi2
		phi = phi - np.amin(phi)
		phi = phi % (2.*np.pi)
	return phi

def fit_min_max(args,p,max_iter,proj_list,proj_axis):
	mins = numpy.array([])
	maxs = numpy.array([])
	
	kind = [item for item in args.kind.split(' ')]
	
	for i in xrange(int(args.fmin)+int(args.step),max_iter+1,int(args.step)):
		args.proj = proj_list[p]
		axis = proj_axis[args.proj]
		dat = load_map(args,p,i)
		unit_l, unit_d, unit_t, unit_m = load_units(i, args)

		if kind[p] == 'dens':
			dat *= unit_d	# in g/cc
		if kind[p] in ['vx','vy','vz']:
			dat *= (unit_l/unit_t)/1e5 # in km/s
		if kind[p] in ['stars','dm']:
			dat += 1e-12
		
		if args.logscale:
			mins = numpy.append(mins,numpy.log10(numpy.amin(dat)))
			maxs = numpy.append(maxs,numpy.log10(numpy.amax(dat)))
		else:
			mins = numpy.append(mins,numpy.amin(dat))
			maxs = numpy.append(maxs,numpy.amax(dat))
		
	ii = range(int(args.fmin)+int(args.step),max_iter+1,int(args.step))
	cmin = polyfit(ii,mins,args.poly)	
	cmax = polyfit(ii,maxs,args.poly)

	return p, cmin, cmax

 def produce_video(self, interval=200, repeat_delay=2000, filename='video_output.gif', override_min_max=None):
     """
     Finalize and save the video of the data.
     
     interval and repeat_delay are the interval between frames and the repeat
         delay before restart, both in milliseconds.
     filename is the name of the file to save in the present working 
         directory. At present, only .gifs will implement reliably without
         tweaking Python's PATHs.
     override_min_max allows the user to set their own maximum and minimum
         for the scale on the plot. Use a len-2 tuple, (min, max).
     """
     #find the limits for the plot:
     if not override_min_max:
         self.min_limit = np.amin(self.data_list[0])
         self.max_limit = np.amax(self.data_list[0])
         assert len(self.data_list) > 1, 'You must include at least two frames to make an animation!'
         for i in self.data_list[1:]: #assumes there is more than one frame in the loop
             self.min_limit = min((self.min_limit, np.amin(i)))
             self.max_limit = max((self.max_limit, np.amax(i)))
     else:
         self.min_limit=override_min_max[0]
         self.max_limit=override_min_max[1]
         
     self.fig.colorbar(self.plotfunc(self.grid, self.data_list[0],limits=(self.min_limit,self.max_limit),allow_colorbar=False, **self.kwds))
     ani = animation.FuncAnimation(self.fig, _make_image, frames=self._yield_image, interval=interval, blit=True, repeat_delay=repeat_delay)
     ani.save(filename, fps=1000./interval)

def _write_data(lock, im, index, outfile, outshape, outtype, rescale_factor, logfilename, cputime, itime):    	      

	lock.acquire()
	try:        
		t0 = time() 			
		f_out = getHDF5(outfile, 'a')					 
		f_out_dset = f_out.require_dataset('exchange/data', outshape, outtype, chunks=tdf.get_dset_chunks(outshape[0])) 
		im = im * rescale_factor
		tdf.write_tomo(f_out_dset,index,im.astype(outtype))
					
		# Set minimum and maximum:
		if (amin(im[:]) < float(f_out_dset.attrs['min'])):
			f_out_dset.attrs['min'] = str(amin(im[:]))
		if (amax(im[:]) > float(f_out_dset.attrs['max'])):
			f_out_dset.attrs['max'] = str(amax(im[:]))		
		f_out.close()			
		t1 = time() 

		# Print out execution time:
		log = open(logfilename,"a")
		log.write(linesep + "\ttomo_%s processed (CPU: %0.3f sec - I/O: %0.3f sec)." % (str(index).zfill(4), cputime, t1 - t0 + itime))
		log.close()	

	finally:
		lock.release()	

def statprint(host_per_pg, pg_per_host):
    val = pg_per_host.values()  # sets val to a list of the values in pg_per_host
    mean = numpy.mean(val)
    maxvalue = numpy.amax(val)
    minvalue = numpy.amin(val)
    std = numpy.std(val)
    median = numpy.median(val)
    variance = numpy.var(val)
    print("for placement groups on hosts: ")
    print( "the mean is: ", mean)
    print( "the max value is: ", maxvalue)
    print( "the min value is: ", minvalue)
    print( "the standard deviation is: ", std)
    print( "the median is: ", median)
    print( "the variance is: ", variance)
    # prints statements for stats
    host_mean = numpy.mean(host_per_pg)
    host_max = numpy.amax(host_per_pg)
    host_min = numpy.amin(host_per_pg)
    host_std = numpy.std(host_per_pg)
    host_median = numpy.median(host_per_pg)
    host_variance = numpy.var(host_per_pg)
    # these are the variables for hosts/pgs
    print("hosts per placement group: ")
    print("the mean is: ", host_mean)
    print("the max value is: ", host_max)
    print("the min value is: ", host_min)
    print("the standard deviation is: ", host_std)
    print("the median is: ", host_median)
    print("the variance is: ", host_variance)

    def get_crange(self, crange=[None,None], 
                         xrange=[None,None], 
                         yrange=[None,None], scale = 'linear'):

        x, y, z = self._eval_xyz()
        if ( x is None or
             y is None) : 
             x = np.arange(z.shape[1])
             y = np.arange(z.shape[0])
        if (self.getvar('offset') is not None and
            (self.getvar('zdir') == 'x' or 
             self.getvar('zdir') == 'y')):
            crange = self._update_range(crange, (np.amin(z), np.amax(z)))
        elif (xrange[0] is not None and
            xrange[1] is not None and
            yrange[0] is not None and
            yrange[1] is not None):
            zt = np.ma.masked_array(z)
            if y.ndim == 1:
               zt[(y < yrange[0]) | (y > yrange[1]),:] = np.ma.masked
            else:
               zt[(y < yrange[0]) | (y > yrange[1])] = np.ma.masked
            if x.ndim == 1:
               zt[:,(x < xrange[0]) | (x > xrange[1])] = np.ma.masked
            else:
               zt[(x < xrange[0]) | (x > xrange[1])] = np.ma.masked
            if scale == 'log': zt[z <= 0] = np.ma.masked
            crange = self._update_range(crange, (np.amin(zt), np.amax(zt)))


        return crange

 def createCube(self, cellSize_xy):
     cellNumber_x = round((self.extent.XMax - self.extent.XMin) / cellSize_xy)
     cellNumber_y = round((self.extent.YMax - self.extent.YMin) / cellSize_xy)
     X = self.ssdo.xyCoords[:,0]
     Y = self.ssdo.xyCoords[:,1]
     time = self.ssdo.fields[self.timeField].data
     time = NUM.array([i for i in time], NUM.datetime64)
     startDateTime = NUM.datetime64('1970-01-01 00:00:00')
     T = time - startDateTime
     self.startTime = NUM.amin(T) + NUM.datetime64('1970-01-01 00:00:00')
     T = NUM.array([i.item().days for i in T], int)
     startT = NUM.amin(T)
     endT = NUM.amax(T)
     cellNumber_t = round((endT - startT) / self.cellSize_t) + 1
     X = (X - self.extent.XMin) / self.cellSize_xy
     Y = (self.extent.YMax - Y) / self.cellSize_xy
     T = (T - startT) / self.cellSize_t
     X = NUM.floor(X)
     Y = NUM.floor(Y)
     T = NUM.floor(T)
     CellIdList = (cellNumber_x * cellNumber_y * T) + (cellNumber_x * Y) + X
     BothEnds = NUM.array([0, (cellNumber_t * cellNumber_x * cellNumber_y -1)])
     CellIdList = NUM.concatenate((CellIdList, BothEnds), axis=0)
     CellIdList = NUM.array(CellIdList, dtype = 'int32')
     counts = NUM.bincount(CellIdList)
     counts[BothEnds[0]] = counts[BothEnds[0]] - 1
     counts[BothEnds[1]] = counts[BothEnds[1]] - 1
     return counts.reshape(cellNumber_t, cellNumber_x, cellNumber_y)

    def __init__(self, template=None, off_pulse_bins=None, 
                       off_pulse_auto_margin=0., 
                       correlation_harmonics=None):

        self.template = template
        self.off_pulse_auto_margin = off_pulse_auto_margin

        if template is not None:
            self.name = "Minimum estimator using cross-correlator"
        elif off_pulse_bins is not None:
            self.name = "Minimum estimator using known minimum"
        else:
            self.name = "Minimum estimator finding minimum"
        self.RMS = False

        if off_pulse_bins is None:
            if template is not None:
                self.off_pulse_bins = list(where(
                    template-amin(template) <= off_pulse_auto_margin*(amax(template)-amin(template))
                    )[0])
            else:
                self.off_pulse_bins = None
        else:
            self.off_pulse_bins = list(off_pulse_bins)
            if not self.off_pulse_bins:
                raise ValueError, "No off-pulse bins specified!"
        self.correlation_harmonics = correlation_harmonics 

def make_plottable(self, method="cubic", nk=50):
   ''' 
   :param method: cubic|nearest|linear
   :param nk: number of k points
   :return: x,y,z, zmin, zmax
   '''
   pl = np.zeros((len(self.mesh), 2))
   ik=0
   for k in self.mesh:
    pl[ik, 0]=k[0]
    pl[ik, 1]=k[1]
    ik+=1
   x = pl[:,0]
   y = pl[:,1]
   xi = np.linspace(min(x), max(x),nk)
   yi = np.linspace(min(y), max(y),nk)
   zmin,zmax=np.zeros((self.data.shape[1], self.data.shape[2]), np.complex64), np.zeros((self.data.shape[1], self.data.shape[2]), np.complex64)
   zi=[]
   for ind_x in range(self.data.shape[1]):
    zi.append([])
    for ind_y in range(self.data.shape[2]):

     z = self.data[:,ind_x,ind_y]
     zmin[ind_x,ind_y]=np.amin(z.real)+np.amin(z.imag)*1j
     zmax[ind_x,ind_y]=np.amax(z.real)+np.amax(z.imag)*1j
     zi[ind_x].append(griddata((x, y), z, (xi[None,:], yi[:,None]), method=method))
   return xi,yi,np.array(zi),zmin,zmax

def normalize_images(images, normalize='max'):
    '''
    Max/min normalizes a set of images

    args:
    	@a images shape = (N,W,H,C) where C is number of channels
    '''
    N, Pw, Ph = images.shape
    images_norm = np.zeros((N,Pw,Ph))

    if normalize=='max':
        maxs = np.amax(images, axis=(1,2))
        mins = np.amin(images,axis=(1,2))

        for i in range(0,N):
            images_norm[i,:] = (images[i]-mins[i])/(maxs[i]-mins[i]+1e-6)
        return images_norm

    if normalize == 'global_max':
        max_ = np.amax(images)
        min_ = np.amin(images)

        images_norm = (images-min_)/(max_-min_)
        return images_norm
    if normalize=='mean':
        pass

def plot_field(X, Y, U, V, filename):
    '''
    Function to plot the potential field.

    Args:
        X (numpy.ndarray): X component of the sample points.
        Y (numpy.ndarray): Y component of the sample points.
        U (numpy.ndarray): X component of field at sample points.
        V (numpy.ndarray): Y component of field at sample points.
    '''
    # Generate plot.
    padding = 0.5
    plt.figure()
    plt.quiver(X, Y, U, V,
               color='#007ce8',
               units='x',
               pivot='tail')
    plt.axis('equal')
    plt.axis([np.amin(X) - padding,
              np.amax(X) + padding,
              np.amin(Y) - padding,
              np.amax(Y) + padding])
    # plt.savefig("potential_field_back1.svg", format='svg')
    plt.savefig(filename + ".svg", format='svg')
    plt.show()

    def get_data_extent(self):
        if (self._data_extent is not None and
            self._data_extent_checked):
            return self._data_extent
        if (self._data_extent is not None and
            self.isempty()):
            return self._data_extent
        x, y = self._eval_xy()
        if self.isempty():
            if x is None:
                self._data_extent=[0, len(y), np.min(y), np.max(y)]
            else:               
                self._data_extent=[np.min(x), np.max(x), np.min(y), np.max(y)]
        else:
            xr = (np.inf, -np.inf)
            yr = (np.inf, -np.inf)
            for a in self._artists:
                xdata = a.get_xdata()
                ydata = a.get_ydata()
                xr = (min((xr[0],np.amin(xdata))), 
                      max((xr[1],np.amax(xdata))))
                yr = (min((yr[0],np.amin(ydata))), 
                      max((yr[1],np.amax(ydata))))

            self._data_extent=[min((min(xr), np.amin(x))),
                               max((max(xr), np.amax(x))),
                               min((min(yr), np.amin(y))),
                               max((max(yr), np.amax(y)))]
            self._data_extent_checked = True
        return self._data_extent

    def _eucl_max(self, nii1, nii2):
        origdata1 = nii1.get_data()
        origdata1 = np.logical_not(
            np.logical_or(origdata1 == 0, np.isnan(origdata1)))
        origdata2 = nii2.get_data()
        origdata2 = np.logical_not(
            np.logical_or(origdata2 == 0, np.isnan(origdata2)))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(
                np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        if origdata1.max() == 0 or origdata2.max() == 0:
            return np.NaN

        border1 = self._find_border(origdata1)
        border2 = self._find_border(origdata2)

        set1_coordinates = self._get_coordinates(border1, nii1.affine)
        set2_coordinates = self._get_coordinates(border2, nii2.affine)
        distances = cdist(set1_coordinates.T, set2_coordinates.T)
        mins = np.concatenate(
            (np.amin(distances, axis=0), np.amin(distances, axis=1)))

        return np.max(mins)

def test():

    filepath = '../pygview_test_data/se_ne5np4_test.cam.h0.0000-01-01-00000.nc'

    dslice = dataslice( filepath, 'T', time_ndx=0, level_ndx=8 )

    print '         filepath : '+filepath
    print 'dslice.data.shape : ',dslice.data.shape
    print 'dslice.structured : ',dslice.structured
    print '     dslice.units : ',dslice.units
    print '  dslice.lon.size : ',dslice.lon.size
    print '  dslice.lat.size : ',dslice.lat.size
    print '  dslice.lon : ',dslice.lon
    print '  dslice.lat : ',dslice.lat

    print '     dslice.lon.min, dslice.lon.max   : ',dslice.lon.min(), dslice.lon.max()

    print '     dslice.lon.min, dslice.lon.max   : ',numpy.amin(dslice.lon), numpy.amax(dslice.lon)
    print '     dslice.lat.min, dslice.lat.max   : ',numpy.amin(dslice.lat), numpy.amax(dslice.lat)
    print '     dslice.data.min,dslice.data.max  : ',numpy.amin(dslice.data), numpy.amax(dslice.data)

    filepath = '../pygview_test_data/fv_10x15_test.cam.h0.0000-01-01-00000.nc'

    dslice = dataslice( filepath, 'V', time_ndx=0, level_ndx=18 )
    
    print '         filepath : '+filepath
    print 'dslice.data.shape : ',dslice.data.shape
    print 'dslice.structured : ',dslice.structured
    print '     dslice.units : ',dslice.units
    print '     dslice.lon   : ',dslice.lon
    print '     dslice.lat   : ',dslice.lat
    print '     dslice.data  : ',dslice.data.min(),dslice.data.max()

    def add_boundary_pores(self, labels, spacing):
        r"""
        Add boundary pores to the specified faces of the network

        Pores are offset from the faces by 1/2 of the given ``spacing``, such
        that they lie directly on the boundaries.

        Parameters
        ----------
        labels : string or list of strings
            The labels indicating the pores defining each face where boundary
            pores are to be added (e.g. 'left' or ['left', 'right'])

        spacing : scalar or array_like
            The spacing of the network (e.g. [1, 1, 1]).  This must be given
            since it can be quite difficult to infer from the network,
            for instance if boundary pores have already added to other faces.

        """
        spacing = np.array(spacing)
        if spacing.size == 1:
            spacing = np.ones(3)*spacing
        for item in labels:
            Ps = self.pores(item)
            coords = np.absolute(self['pore.coords'][Ps])
            axis = np.count_nonzero(np.diff(coords, axis=0), axis=0) == 0
            offset = np.array(axis, dtype=int)/2
            if np.amin(coords) == np.amin(coords[:, np.where(axis)[0]]):
                offset = -1*offset
            topotools.add_boundary_pores(network=self, pores=Ps, offset=offset,
                                         apply_label=item + '_boundary')

def basemap_raster_mercator(lon, lat, grid, cmap = None):
  """
  Render a raster in mercator projection.  Locations with no values are
  rendered transparent.
  """
  # longitude/latitude extent
  lons = (np.amin(lon), np.amax(lon))
  lats = (np.amin(lat), np.amax(lat))

  if cmap is None:
    cmap = mpl.cm.jet
    cmap.set_bad('w', 1.0)

  # construct spherical mercator projection for region of interest
  m = Basemap(projection='merc',llcrnrlat=lats[0], urcrnrlat=lats[1],
              llcrnrlon=lons[0],urcrnrlon=lons[1])

  vmin,vmax = np.nanmin(grid),np.nanmax(grid)
  masked_grid = np.ma.array(grid,mask=np.isnan(grid))
  fig = plt.figure(frameon=False)
  plt.axis('off')
  m.pcolormesh(lon,lat,masked_grid,latlon=True,cmap=cmap,vmin=vmin,vmax=vmax)

  str_io = StringIO.StringIO()
  plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True)
  bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ]

  return str_io.getvalue(), bounds

def _lowest_index(dm):
    """Return the index of the lowest value in the input distance matrix.

    If there are ties for the lowest value, the index of top-left most
    occurrence of that value will be returned.

    This should be ultimately be replaced with a new DistanceMatrix object
    method (#228).

    """
    # get the positions of the lowest value
    results = np.vstack(np.where(dm.data == np.amin(dm.condensed_form()))).T
    # select results in the bottom-left of the array
    results = results[results[:, 0] > results[:, 1]]
    # calculate the distances of the results to [0, 0]
    res_distances = np.sqrt(results[:, 0]**2 + results[:, 1]**2)
    # detect distance ties & return the point which would have
    # been produced by the original function
    if np.count_nonzero(res_distances == np.amin(res_distances)) > 1:
        eqdistres = results[res_distances == np.amin(res_distances)]
        res_coords = eqdistres[np.argmin([r[0] for r in eqdistres])]
    else:
        res_coords = results[np.argmin(res_distances)]

    return tuple([res_coords[0], res_coords[1]])