def logpdf(x, mu, sigma):
     if type(x) == np.ndarray:
         if sigma > 0:
             small = np.log(0.5 + 0.5 * special.erf((np.log(1e-6) - mu) /
                            (np.sqrt(2.0) * sigma)))
             I = pp.find(x > 1e-6)
             log_x = np.log(x[I])
             lp = small * np.ones(x.shape)
             lp[I] = -log_x - 0.5 * np.log(2.0 * np.pi) - np.log(sigma) - \
                 0.5 * ((log_x - mu) ** 2) / (sigma ** 2)
         else:
             I = pp.find(x == mu)
             lp = -np.inf * np.ones(x.shape)
             lp[I] = 0
     else:
         if sigma > 0:
             if x > 1e-6:
                 log_x = np.log(x)
                 lp = - log_x - 0.5 * np.log(2.0 * np.pi) - \
                     np.log(sigma) - \
                     0.5 * ((log_x - mu) ** 2) / (sigma ** 2)
             else:
                 lp = np.log(0.5 + 0.5 * special.erf((np.log(1e-6) - mu) /
                             (np.sqrt(2.0) * sigma)))
         else:
             if x == mu:
                 lp = 0
             else:
                 lp = -np.inf
     return lp

def ellipse(width, height):

    horpow = 2
    verpow = horpow
    # to produce a Ellipse with horizontal axis == ceil(2*hor semi axis)
    width = width + 0.5
    #to produce a Ellipse with vertical axis == ceil(2*vert semi axis)
    height = height + 0.5
    [x, y] = np.meshgrid(np.arange(-width, width + 1),
                        np.arange(-height, height + 1))
    #print width, height, x.shape, y.shape
    bll = (abs(x / width) ** horpow + abs(y / height) ** verpow) < 1
    xs = plt.find(bll.sum(0) == 0)
    ys = plt.find(bll.sum(1) == 0)
    bll = bll[ys[0] + 1:ys[1], :]
    bll = bll[:, xs[0] + 1:xs[1]]
    bll = bll.T
    mask_inds = np.where(bll > 0)
    mask_inds = np.array(mask_inds).T
    inv_mask_inds = np.where(bll == 0)
    inv_mask_inds = np.array(inv_mask_inds).T
    box_inds = np.array(np.where(bll != 5555)).T
    box_edges = np.array(
                [[box_inds[:, 0].min(), box_inds[:, 1].min()],
                [box_inds[:, 0].max(), box_inds[:, 1].max()]])
    return mask_inds, inv_mask_inds, box_edges

    def extract_time_series(self, filename, node_yxz):
        dis = self.parent.get_package('DIS')
        if (dis != None):
            yy, xx, zz = dis.get_node_coordinates()

            times, heads, strings = self.read_all(filename)

            if self.silent == 0: print ('Please wait. Extracting time series.')
            ht = np.empty((1, len(node_yxz)))
            rr = 0
            for h in heads:
                if (rr > 0):
                    ht = np.vstack((ht, np.empty(len(node_yxz))))
                cc = 0
                for yxz in node_yxz:
                    y = yxz[0]
                    x = yxz[1]
                    z = yxz[2]
                    r = find(abs(yy - y) < abs(y) / 1e6)[0]
                    c = find(abs(xx - x) < abs(x) / 1e6)[0]
                    l = find(abs(zz[r, c] - z) < abs(z) / 1e6)[0]
                    ht[rr, cc] = h[r, c, l]
                    cc = cc + 1
                rr = rr + 1
            return times, ht

    def energy_spectrum(self, psi, energy_grid):
	"""
	energy_grid, spectrum = energy_spectrum(psi, energy_grid)

	Creates the energy spectrum on the <energy_grid>.

	"""
	all_energies = diag(self.input_function.H_0)
	all_population = abs(psi)**2

	spectrum = zeros(len(energy_grid))
    
	#Density of states.
	for q in range(max(self.input_function.index_array[:,1])+1):
	    q_indices = find(self.input_function.index_array[:,1] == q)
	    energies = all_energies[q_indices]
	    population = all_population[q_indices]
	    dos = 1 / diff(energies)
	    if True:
		#Choose an appropriate output grid.
		start_index = find(energy_grid > energies[0])[0]
		end_index = find(energy_grid < energies[-1])[-1] + 1 

		spectrum[start_index:end_index] += UnivariateSpline(energies[1:], 
		    population[1:] * dos, s=0)(energy_grid[start_index:end_index])	
	    else:
		spectrum += UnivariateSpline(energies[1:], 
		    population[1:] * dos, s=0)(energy_grid)
	    #raise
	return energy_grid, spectrum

    def make_continuum_list(self):
	"""
	Creates a list of the states that are in the continuum, 
	based on the density of states. 
	"""
	energies = diag(self.input_function.H_0)
	self.energies = energies
	continuum = []
	border_value_list = []

	#Looping the q's.
	for i in range(max(self.input_function.index_array[:,1])+1):

	    #Get the indices corresponding to  given q.
	    temp_ind = find(self.input_function.index_array[:,1] == i)
	    
	    #Corresponding energies.
	    temp_E = energies[temp_ind]
	    
	    #The continuum starts when gap between energies increases.
	    border_value = temp_ind[argmin(diff(temp_E))]
	    border_value_list.append(energies[border_value])

	    continuum.append(temp_ind[find((temp_ind > border_value))])
	
	self.continuum = [item for sublist in continuum for item in sublist]
	self.continuum_limit_dos = mean(border_value_list)

def lognormal_cdf( x, mu, sigma ):
  if sigma > 0:
    small = 0.5 + 0.5*special.erf( (np.log(1e-6)-mu)/(np.sqrt(2.0)*sigma))
    if x.__class__ == np.ndarray:
      lp = np.zeros( len(x) )
      I = pp.find( x > 1e-6 )
      J = pp.find( x <= 1e-6)
      lp[I] = 0.5 + 0.5*special.erf( (np.log(x)-mu)/(np.sqrt(2.0)*sigma))
      lp[J] = small
      return lp
    else:
      if x > 1e-6:
        return 0.5 + 0.5*special.erf( (np.log(x)-mu)/(np.sqrt(2.0)*sigma))
      else:
        return small
  else:
    if x.__class__ == np.ndarray:
      logx = np.log(x+1e-6)
      lp = 0.5*np.ones( len(x))
      I1 = pp.find( logx < mu )
      I2 = pp.find( logx > mu )
      
      lp[I1] = 0
      lp[I2] = 1
      return lp
    else:
      if np.log(x) < mu:
        return 0
      elif np.log(x) > mu:
        return 1
      else:
        return 0.5

def plt_annual_WBCt(fign):


	fig = plt.figure(num=fign,figsize=(8,5),facecolor='w')
	ax = plt.gca()


	argo = np.ones(soest_gs.latitude.shape[0])*np.nan
	wind = np.ones(scow_gs.latitude.shape[0])*np.nan


	d = py.find(soest_gs.depth == 1200)[0]
	for i in xrange(argo.shape[0]):
		gs = soestgs_annual_Ipsi[d,i,:]
		good  = py.find(~np.isnan(gs))

		if len(good) > 0:
			argo[i] = gs[good[0]]

		del gs, good


	for i in xrange(wind.shape[0]):
		gs = scowgs_annual_psi[i,:]
		good  = py.find(~np.isnan(gs))

		if len(good) > 0:		
			wind[i] = gs[good[0]]

		del gs, good

	ax.plot(soest_gs.latitude,argo,'b',lw=3.5,label='ARGO')
	ax.plot(scow_gs.latitude,wind,'r',lw=3.5,label='Sverdrup - Ekman')

	ax.legend(numpoints=1,loc=0,prop={'size':20},frameon=False)
	ax.plot([-60,0],[0,0],'k--',lw=3.5)

	ax.set_ylim(-40,40)
	ax.set_xlim(-5,-50)

	ax.set_yticks(range(-40,50,10))
	ax.set_xticks(range(-50,0,5))
	ax.tick_params(labelsize=18)
	rstyle(ax)

	labels = [ur'5$^{\circ}$S',ur'10$^{\circ}$S',ur'15$^{\circ}$S',ur'20$^{\circ}$S',
	ur'25$^{\circ}$S',ur'30$^{\circ}$S',ur'35$^{\circ}$S',ur'40$^{\circ}$S',ur'45$^{\circ}$S',
	ur'50$^{\circ}$S']

	labels.reverse()

	ax.set_xticklabels(labels)

	ax.set_ylabel('Transport at X$_{w}$ (Sv)',fontsize=30,fontweight='demibold')
	ax.set_xlabel(ur'Latitude',fontsize=30,fontweight='demibold')

	plt.draw()
	plt.show(block=False)

	return fig

def get_affine_inliers_RANSAC(num_m, xy1_m, xy2_m,\
                              acd1_m, acd2_m, xy_thresh_sqrd, sigma_thresh_sqrd=None):
    '''Computes initial inliers by iteratively computing affine transformations
    between matched keypoints'''
    aff_inliers = []
    # Enumerate All Hypothesis (Match transformations)
    for mx in xrange(num_m): 
        xy1  = xy1_m[:,mx].reshape(2,1) #  XY Positions
        xy2  = xy2_m[:,mx].reshape(2,1) 
        A1   = matrix(insert(acd1_m[:,mx], [1.], 0.)).reshape(2,2)
        A2   = matrix(insert(acd2_m[:,mx], [1.], 0.)).reshape(2,2)
        # Compute Affine Tranform 
        # from img1 to img2 = (E2\E1) 
        Aff  = linalg.inv(A2).dot(A1)
        #
        # Transform XY-Positions
        xy1_mAt = xy2 + Aff.dot( (xy1_m - xy1) ) 
        xy_err_sqrd = sum( power(xy1_mAt - xy2_m, 2) , 0)
        _inliers = find(xy_err_sqrd < xy_thresh_sqrd)
        #
        # Transform Ellipse Geometry (solved on paper)
        if not sigma_thresh_sqrd is None:
            scale1_mAt = (acd1_m[0]*Aff[0,0]) *\
                         (acd1_m[1]*Aff[1,0]+acd1_m[2]*Aff[1,1])
            scale2_m   = acd2_m[0] * acd2_m[2]
            scale_err  = np.abs(scale1_mAt - scale2_m)
            _inliers_scale = find(scale_err < sigma_thresh_sqrd)
            _inliers = np.bitwise_and(_inliers, _inliers_scale)
        #If this hypothesis transformation is better than the ones we have
        #previously seen then set it as the best
        if len(_inliers) > len(aff_inliers):
            aff_inliers = _inliers
            #bst_xy_err  = xy_err_sqrd 
    return aff_inliers

    def res2_name_consistency(em, res):
        '''Score a single query for name consistency
        Input: 
            res - query result
        Returns: Dict
            error_chip - degree of chip error
            name_error - degree of name error
            gt_pos_name - 
            gt_pos_chip - 
        '''
        # Defaults to -1 if no ground truth is in the top results
        cm, nm = em.hs.get_managers('cm','nm')
        qcx  = res.rr.qcx
        qnid = res.rr.qnid
        qnx   = nm.nid2_nx[qnid]
        ret = {'name_error':-1, 'chip_error':-1,
               'gt_pos_chip':-1, 'gt_pos_name':-1, 
               'chip_precision': -1, 'chip_recall':-1}
        if qnid == nm.UNIDEN_NID: exec('return ret')
        # ----
        # Score Top Chips
        top_cx = res.cx_sort()
        gt_pos_chip_list = (1+pylab.find(qnid == cm.cx2_nid(top_cx)))
        # If a correct chip was in the top results
        # Reward more chips for being in the top X
        if len(gt_pos_chip_list) > 0:
            # Use summation formula sum_i^n i = n(n+1)/2
            ret['gt_pos_chip'] = gt_pos_chip_list.min()
            _N = len(gt_pos_chip_list)
            _SUM_DENOM = float(_N * (_N + 1)) / 2.0
            ret['chip_error'] = float(gt_pos_chip_list.sum())/_SUM_DENOM
        # Calculate Precision / Recall (depends on the # threshold/max_results)
        ground_truth_cxs = np.setdiff1d(np.array(nm.nx2_cx_list[qnx]), np.array([qcx]))
        true_positives  = top_cx[gt_pos_chip_list-1]
        false_positives = np.setdiff1d(top_cx, true_positives)
        false_negatives = np.setdiff1d(ground_truth_cxs, top_cx)

        nTP = float(len(true_positives)) # Correct result
        nFP = float(len(false_positives)) # Unexpected result
        nFN = float(len(false_negatives)) # Missing result
        #nTN = float( # Correct absence of result

        ret['chip_precision'] = nTP / (nTP + nFP)
        ret['chip_recall']    = nTP / (nTP + nFN)
        #ret['true_negative_rate'] = nTN / (nTN + nFP)
        #ret['accuracy'] = (nTP + nFP) / (nTP + nTN + nFP + nFN)
        # ----
        # Score Top Names
        (top_nx, _) = res.nxcx_sort()
        gt_pos_name_list = (1+pylab.find(qnid == nm.nx2_nid[top_nx]))
        # If a correct name was in the top results
        if len(gt_pos_name_list) > 0: 
            ret['gt_pos_name'] = gt_pos_name_list.min() 
            # N should always be 1
            _N = len(gt_pos_name_list)
            _SUM_DENOM = float(_N * (_N + 1)) / 2.0
            ret['name_error'] = float(gt_pos_name_list.sum())/_SUM_DENOM
        # ---- 
        # RETURN RESULTS
        return ret

def detect_signals():
    vector, label = weeklydataset_sg_ndata(
        "/media/4AC0AB31C0AB21E5/Documents and Settings/Claudio/Documenti/Thesis/Workloads/MSClaudio/ews/access_log-20110805.csv",
        [],
    )
    x, target = aggregatebymins_sg_ndata(vector[1])

    starttime = time.time()
    y = array(target)
    t = array(x)
    thr = max(y) * 2 / 3
    print thr
    I = pylab.find(y > thr)
    #    print I
    #    pylab.plot(t,y, 'b',label='signal')
    #    pylab.plot(t[I], y[I],'ro',label='detections')
    #    pylab.plot([0, t[len(t)-1]], [thr,thr], 'g--')

    J = pylab.find(diff(I) > 1)
    argpeak = []
    targetpeak = []
    for K in split(I, J + 1):
        ytag = y[K]
        peak = pylab.find(ytag == max(ytag))
        #        pylab.plot(peak+K[0],ytag[peak],'sg',ms=7)
        argpeak.append(peak + K[0])
        targetpeak.append(ytag[peak])

    eta = time.time() - starttime
    print "time elapsed %f" % eta
    return list(itertools.chain(*argpeak)), list(itertools.chain(*targetpeak))

def detect_saccades(orientations, timestep):
    orientations = orientations[~numpy.isnan(orientations)]

    unwrappedOrientations = numpy.unwrap(orientations, 180)
    SMOOTH_TIME = 1
    smoothWindow = int(round(SMOOTH_TIME / timestep))
    smoothedOrientations = tools.smooth(unwrappedOrientations, smoothWindow)
    angSpds = abs(numpy.diff(smoothedOrientations)) / timestep
    MIN_PEAK_SPD = 30
    sacInds = tools.local_minima(-angSpds) & (angSpds > MIN_PEAK_SPD)
    sacSpds = angSpds[sacInds]

    MIN_SAC_SPD = 5
    inSac = angSpds > MIN_SAC_SPD
    startInds = pylab.find(numpy.diff(inSac.view(dtype=int8)) == 1)
    stopInds = pylab.find(numpy.diff(inSac.view(dtype=int8)) == -1)
    if stopInds[-1] < startInds[-1]:
        l = stopInds.tolist()
        l.append(len(inSac) - 1)
        stopInds = numpy.array(l)
    if startInds[0] > stopInds[0]:
        l = startInds.tolist()
        l.insert(0, 0)
        startInds = numpy.array(l)

    angDisp = numpy.zeros(len(angSpds))
    for i, startInd in enumerate(startInds):
        stopInd = stopInds[i]
        angDisp[startInd:stopInd] = smoothedOrientations[stopInd] - smoothedOrientations[startInd]

    sacAmps = angDisp[sacInds]
    return sacAmps

def vl_test_hikmeans():
	""" VL_TEST_HIKMEANS Test VL_HIKMEANS function
	"""
	K = 3
	nleaves = 100
	data = numpy.array(numpy.random.rand(2, 1000) * 255, 'uint8')
	datat = numpy.array(numpy.random.rand(2, 10000) * 255, 'uint8')
	
	[tree, A] = vlfeat.vl_hikmeans(data, K, nleaves, verb=1)
	AT = vlfeat.vl_hikmeanspush(tree, datat)
		
	pylab.figure()
	plottree(tree)
	pylab.xlim(0, 255)
	pylab.ylim(0, 255)
	print('hikmeans-tree') ;
	
	pylab.figure()
	gen = color_gen()
	for k in range(K*K):
		color = next(gen)
		sel = pylab.find(A[-1, :] == k)
		pylab.plot(data[0, sel], data[1, sel], '.', color=color)
		sel = pylab.find(AT[-1, :] == k)
		pylab.plot(datat[0, sel], datat[1, sel], '+', color=color)
	
	plottree(tree, linewidth=4)
	pylab.xlim(0, 255)
	pylab.ylim(0, 255)
	print('hikmeans-clusters') ;

    def __init__(self, paramfile):
        for row in paramfile:

            if row[0]=='BIAS_COMBINE': self.BIAS_COMBINE = row[1]
            if row[0]=='COMP_COMBINE': self.COMP_COMBINE = row[1]
            if row[0]=='FLAT_COMBINE': self.FLAT_COMBINE = row[1]
            if row[0]=='FIBER_TRACE': self.FIBER_TRACE = row[1]
            if row[0]=='BIAS_SUBTRACT': self.BIAS_SUBTRACT = row[1]
            if row[0]=='LAMBDA_SOLVE': self.LAMBDA_SOLVE = row[1]
            if row[0]=='COSMIC_RAYS': self.COSMIC_RAYS = row[1]
            if row[0]=='WRITE_SPECTRA': self.WRITE_SPECTRA = row[1]
            
            if row[0]=='FIBERS_TOTAL': self.FIBERS_TOTAL = int(row[1])
            if row[0]=='FIBERS_EXCLUDE':
                fibs = pl.array( [ int(i) for i in row[1].split(',') ] )
                fibs.sort()
                fibs = fibs[ pl.find( 0<fibs ) ]
                fibs = fibs[ pl.find( fibs<self.FIBERS_TOTAL+1 ) ]
                self.FIBERS_EXCLUDE = fibs
            if row[0]=='FIBER_WIDTH': self.FIBER_WIDTH = int(row[1])
            if row[0]=='FIBER_SEP': self.FIBER_SEP = int(row[1])
            if row[0]=='LO_BUFFER': self.LO_BUFFER = int(row[1])
            if row[0]=='HI_BUFFER': self.HI_BUFFER = int(row[1])
            if row[0]=='COMP_HEADER': self.COMP_HEADER = row[1]
            if row[0]=='COMP_NAME': self.COMP_NAME = row[1]
            if row[0]=='POLYFIT_ORDER': self.POLYFIT_ORDER = int(row[1])
            if row[0]=='CR_SCAN_DX': self.CR_SCAN_DX = int(row[1])

    def find_continuum(self):
	"""
	Creates a list of el states that are in the continuum.
	Stores the list as an instance variable.
	"""
	f = tables.openFile(self.name_el_couplings)
	
	#Electronic basis state energies.
	r_grid = f.root.R_grid[:]
	index_middle = argmin(abs(r_grid - 4))
	energies = f.root.E[:,index_middle]
	f.close()
	
	continuum = []
	
	#Looping the q's.
	for i in range(max(self.index_array[:,1])+1):

	    #Get the indices corresponding to  given q.
	    temp_ind = find(self.index_array[:,1] == i)
	    
	    #Corresponding energies.
	    temp_E = energies[temp_ind]
	    
	    #Checks if there are positive energies.
	    if temp_E[-1] > 0:
		#The continuum starts when gap between energies increases.
		#OBS: TODO DEBUG
		border_value = temp_ind[argmin(diff(temp_E))] + 0
		continuum.append(temp_ind[find((temp_ind > border_value))])
	
	self.continuum = [item for sublist in continuum for item in sublist]

    def named_state_population(self, psi, name_list):
	"""
	populations = named_state_population(psi, name_list)

	Returns the population in a set of named states.

	Parameters
	----------
	psi : 2D complex array. The wavefunction to be analyzed.
	name_list : string list, e.g. ["2p", "3d"].

	Returns
	-------
	populations : 1D float array, the population of the basis states. 
	"""

	angular_names = {"s":0, "p":1, "d":2, "f":3}

	el_indices = []

	for i in name_list:
	    q_matches = find(self.index_array[:,1] == (angular_names[i[1]]))
	    n_matches = find(self.index_array[:,2] == int(i[0]) - angular_names[i[1]] -1)
	    el_indices.append(intersect1d(q_matches, n_matches))
	
	el_indices = list(ravel(array(el_indices)))
	populations = self.state_population(psi, el_indices)

	return populations

 def _CalculateMaxActivity(self,
                           culture_levels,
                           activities): 
     """Calculates the maximal promoter activity.
     
     Args:
         culture_levels: the culture levels.
         activities: the activities computed.
     
     Returns:
         The maximal (smoothed) reporter level within the
         allowed range of culture levels.
     """
     abovemin_i = pylab.find(culture_levels >= self.min_culture_level)
     abovemax_i = pylab.find(culture_levels >= self.max_culture_level)
     if not abovemin_i.any() or not abovemax_i.any():
         return numpy.NAN
     
     lower_i = numpy.min(abovemin_i)
     upper_i = numpy.min(abovemax_i)
     if lower_i >= upper_i:
         return numpy.NAN
     
     window_size = upper_i - lower_i
     
     activities_window = activities[lower_i:lower_i+window_size]
     diff_activities = self._MovingAverage(activities_window)        
     max_activity = numpy.max(diff_activities)
     if max_activity < 0:
         return numpy.NAN        
     
     return max_activity   

 def ConditionalNRLHist(self, nrls, labels):
     """
      Analyze method
     :param nrls: 
     :type nrls: 
     :param labels:
     :type labels:
     :rtype: Resultlist of figures 
     """
     figures = []
     for m in range(0, self.scen.M):
         q = labels[m, 1, :]
         ind = find(q >= 0.5)
         ind2 = find(q < 0.5)
         r = nrls[ind]
         xmin, xmax = min(nrls), max(nrls)
         lnspc = linspace(xmin, xmax, 100)
         m, s = stats.norm.fit(r)
         pdf_g = stats.norm.pdf(lnspc, m, s)
         r = nrls[ind2]
         # xmin, xmax = min(r), max(r)
         lnspc2 = linspace(xmin, xmax, 100)  # len(r)
         m, s = stats.norm.fit(r)
         pdf_g2 = stats.norm.pdf(lnspc2, m, s)
         fg = figure()
         plot(lnspc, pdf_g / len(pdf_g), label="Norm")
         hold(True)
         plot(lnspc2, pdf_g2 / len(pdf_g2), 'k', label="Norm")
         legend(['Posterior: Activated', 'Posterior: Non Activated'])
         # xmin, xmax = min(xt), max(xt)
         # ind2 = find(q <= 0.5)
         figures.append(fg)
     if self.shower:
         show()
     return figures

def climatology(t, z, w=None, result='year'):
    """Returns monthly climatology of a time-series.

    PARAMETERS
        t (array like) :
            Time in matplotlib time format.
        z (array like) :
            The data to calculate climatology.
        w (array like) :
            Data weight, should have same dimensions as 'z'.
        result (string) :
            If set to 'year', returns only one year of data. If set to
            'full', returns the climatology for every time t.

    RETURN
        z_clim (array like) :
            Climatological averages.

    """
    # Checks for proper dimensions
    shape = z.shape
    if len(shape) == 1:
        z = z[:, None, None]
    c, b, a = z.shape

    if w != None:
        if w.shape != z.shape:
            raise Warning, 'Data and weight arrays are not the same.'

    # Starts converting time to datetime format. Determines the start and end
    # of the relevant dataset to ensure that only whole years are used.
    # Initializes climatology variable and calculates the averages.
    Time = common.num2ymd
    try:
        start = pylab.find(Time[:, 1] == 1)[0]
    except:
        start = None
    try:
        end = pylab.find(Time[:, 1] == 12)[-1]
    except:
        end = None
    time_clim = numpy.arange(12) + 1
    z_clim = numpy.ma.zeros([12, b, a])
    for i in range(12):
        selt = pylab.find(Time[start:end+1, 1] == i + 1) + start
        if w == None:
            z_clim[i, :, :] = z[selt, :, :].mean(axis=0)
        else:
            z_clim[i, :, :] = ((z[selt, :, :] * w[selt, :, :]).sum(axis=0) /
                w[selt, :, :].sum(axis=0))
    #
    z_clim.mask = z_clim.mask | numpy.isnan(z_clim.data)

    if result == 'full':
        z_clim = z_clim[Time[:, 1]-1]

    return z_clim, Time

def GetAverageRange(X, lower_lim, domainheight):
#
# get range of X for averaging such that lower_lim<X<h0
  try: start_val = pylab.find(numpy.array(X)>lower_lim)[0]
  except IndexError: print "not enough values of X, \n lower lim = " + str(lower_lim) + "\n X = " + str(X); sys.exit(1)
  try: end_val = pylab.find(numpy.array(X)>0.4-domainheight)[0]
  except IndexError: end_val = len(X)
#
  return start_val, end_val

def pdf_from_cv(cv, origin, roundoff=1e-9):
    """
    Returns a gaussian pdf from a given covariance matrix
    
    @param cv (n-by-n): covariance matrix
    @param origin (n-by-1): origin
    @param roundoff (float): ratio of largest to smallest eigenvalue 
        s.th. smallest eigenvalue is not considered 0 (degenerate case)
        
    @return Returns a pdf function        
    """
    
    if not np.allclose(cv - cv.T, 0):
        raise ValueError("argument must be a positive definite symmetric "
                "covariance matrix")
        
    dim = cv.shape[0]
    x_ref = np.array(origin)
    if len(x_ref) != dim:
        raise ValueError("x_ref dimension must match covariance dimension")
    
    # this is numerically stable here since cv is symmetric
    # -> eigenvectors are normal
    eigvals, eigvecs = np.linalg.eig(cv) 
    if not all(eigvals > 0):
        raise ValueError("argument must be a positive definite "
                "symmetric covariance matrix")
        
    max_ev = max(eigvals)
    small_ev = find(eigvals < roundoff*max_ev)
    large_ev = find(eigvals >= roundoff*max_ev)
    eigvals_i = eigvals.copy()
    eigvals_i[small_ev] = 0
    eigvals_i[large_ev] = 1./eigvals_i[large_ev]
    
    # compute the pseudo-inverse using eigenvalues
    # (numerically stable here, see above)
    cv_i = np.dot(eigvecs, 
               np.dot(np.diag(eigvals_i), eigvecs.T))
    
    p_det = reduce(lambda x, y: x*y, eigvals[large_ev]) # pseudo-determinant
    
    scale = (1. / np.sqrt((2*np.pi)**dim * p_det))
    projector = np.dot(cv, cv_i)
    def fun(xp, threshold=1e-10, debug=False):
        x = np.dot(projector, xp)
        if np.linalg.norm(x - xp) > threshold:
            if debug:
                print "debug:", svd(projector, 0, 0)
                print "threshold exceeded:", x, xp
                print ("(norm of", x - xp, ":", np.linalg.norm(x - xp), " > ",
                        threshold, ") ")

            return 0
        return scale * np.exp(-0.5 * np.dot(x - x_ref, np.dot(cv_i, x - x_ref)))
    return fun