def nms(dets,proba, T):
     
     dets = dets.astype("float")
     if len(dets) == 0:
         return []
     
     x1 = dets[:, 0]
     y1 = dets[:, 1]
     x2 = dets[:, 2]
     y2 = dets[:, 3]
     scores = proba
     
     areas = (x2 - x1 + 1) * (y2 - y1 + 1)
     order = scores.argsort()[::-1]
     
     keep = []
     while order.size > 0:
         i = order[0]
         keep.append(i)
         xx1 = sp.maximum(x1[i], x1[order[1:]])
         yy1 = sp.maximum(y1[i], y1[order[1:]])
         xx2 = sp.minimum(x2[i], x2[order[1:]])
         yy2 = sp.minimum(y2[i], y2[order[1:]])
     
         w = sp.maximum(0.0, xx2 - xx1 + 1)
         h = sp.maximum(0.0, yy2 - yy1 + 1)
         inter = w * h
         ovr = inter / (areas[i] + areas[order[1:]] - inter)
         inds = sp.where(ovr <= T)[0]
         order = order[inds + 1]
     
     return keep

def fitPairwiseModel(Y,XX=None,S_XX=None,U_XX=None,verbose=False):
    N,P = Y.shape
    """ initilizes parameters """
    RV = fitSingleTraitModel(Y,XX=XX,S_XX=S_XX,U_XX=U_XX,verbose=verbose)
    Cg = covariance.freeform(2)
    Cn = covariance.freeform(2)
    gp = gp2kronSum(mean(Y[:,0:2]),Cg,Cn,XX=XX,S_XX=S_XX,U_XX=U_XX)
    conv2 = SP.ones((P,P),dtype=bool)
    rho_g = SP.ones((P,P))
    rho_n = SP.ones((P,P))
    for p1 in range(P):
        for p2 in range(p1):
            if verbose:
                print '.. fitting correlation (%d,%d)'%(p1,p2)
            gp.setY(Y[:,[p1,p2]])
            Cg_params0 = SP.array([SP.sqrt(RV['varST'][p1,0]),1e-6*SP.randn(),SP.sqrt(RV['varST'][p2,0])])
            Cn_params0 = SP.array([SP.sqrt(RV['varST'][p1,1]),1e-6*SP.randn(),SP.sqrt(RV['varST'][p2,1])])
            params0 = {'Cg':Cg_params0,'Cn':Cn_params0}
            conv2[p1,p2],info = OPT.opt_hyper(gp,params0,factr=1e3)
            rho_g[p1,p2] = Cg.K()[0,1]/SP.sqrt(Cg.K().diagonal().prod())
            rho_n[p1,p2] = Cn.K()[0,1]/SP.sqrt(Cn.K().diagonal().prod())
            conv2[p2,p1] = conv2[p1,p2]; rho_g[p2,p1] = rho_g[p1,p2]; rho_n[p2,p1] = rho_n[p1,p2]
    RV['Cg0'] = rho_g*SP.dot(SP.sqrt(RV['varST'][:,0:1]),SP.sqrt(RV['varST'][:,0:1].T))
    RV['Cn0'] = rho_n*SP.dot(SP.sqrt(RV['varST'][:,1:2]),SP.sqrt(RV['varST'][:,1:2].T))
    RV['conv2'] = conv2
    #3. regularizes covariance matrices
    offset_g = abs(SP.minimum(LA.eigh(RV['Cg0'])[0].min(),0))+1e-4
    offset_n = abs(SP.minimum(LA.eigh(RV['Cn0'])[0].min(),0))+1e-4
    RV['Cg0_reg'] = RV['Cg0']+offset_g*SP.eye(P)
    RV['Cn0_reg'] = RV['Cn0']+offset_n*SP.eye(P)
    RV['params0_Cg']=LA.cholesky(RV['Cg0_reg'])[SP.tril_indices(P)]
    RV['params0_Cn']=LA.cholesky(RV['Cn0_reg'])[SP.tril_indices(P)]
    return RV

    def test_Wells_and_Coppersmith_94(self):
        Wells_and_Copper_94 = conversion_functions['Wells_and_Coppersmith_94']
        calc_max_width_in_slab = conversion_functions['calc_max_width_in_slab']
        fault_type = 'reverse'
        max_width = 15
        max_length = 20
        Mw = array([4.1, 5.2, 6.3, 5.3, 9.5, 6.0, 3.3, 7])
        (area, width) = Wells_and_Copper_94(fault_type, Mw, max_width)
        length = minimum((area / width), max_length)
        msg = ('Unexpected width or length values')
        self.failUnlessAlmostEqual(length[2], 16.25548756, 7, msg)
        self.failUnlessAlmostEqual(width[2], 9.39723311, 7, msg)

        # print Mw
        # print length
        # print width
        slab_width = 8
        out_of_dip = array([0.0, 10, 90, 30, 110.0, 0.0, 111.0, 90.0])
        (area, width) = Wells_and_Copper_94(fault_type, Mw, max_width)
        width = minimum(
            calc_max_width_in_slab(out_of_dip, slab_width, max_width),
            width)
        length = minimum((area / width), max_length)
        # print Mw
        # print length
        # print width

        self.failUnlessAlmostEqual(length[2], 19.09457573, 1, msg)
        self.failUnlessAlmostEqual(width[2], 8., 1, msg)

 def nms(boxes, T = 0.5):
     if len(boxes) == 0:
         return []
     boxes = boxes.astype("float")
     pick = []
     x1 = boxes[:,0]
     y1 = boxes[:,1]
     x2 = boxes[:,2]
     y2 = boxes[:,3]    
     area = (x2 - x1 + 1) * (y2 - y1 + 1)
     idxs = sp.argsort(y2)    
     while len(idxs) > 0:
         last = len(idxs) - 1
         i = idxs[last]
         pick.append(i)
         xx1 = sp.maximum(x1[i], x1[idxs[:last]])
         yy1 = sp.maximum(y1[i], y1[idxs[:last]])
         xx2 = sp.minimum(x2[i], x2[idxs[:last]])
         yy2 = sp.minimum(y2[i], y2[idxs[:last]])
         w = sp.maximum(0, xx2 - xx1 + 1)
         h = sp.maximum(0, yy2 - yy1 + 1)
         I = w * h
         #overlap_ratio = I / area[idxs[:last]]
         overlap_ratio = I /(area[i] +  area[idxs[:last]] - I)
         idxs = sp.delete(idxs, sp.concatenate(([last], sp.where(overlap_ratio > T)[0])))
     return boxes[pick].astype("int")

    def _computeBGDiff(self):
        self._flow.update( self._imageBuffer.getLast() )
        
        n = len(self._imageBuffer)        
        prev_im = self._imageBuffer[0]
        forward = None
        for i in range(0,n/2):
            if forward == None:
                forward = self._imageBuffer[i].to_next
            else:
                forward = forward * self._imageBuffer[i].to_next
                
        w,h = size = prev_im.size
        mask = cv.CreateImage(size,cv.IPL_DEPTH_8U,1)
        cv.Set(mask,0)
        interior = cv.GetSubRect(mask, pv.Rect(2,2,w-4,h-4).asOpenCV()) 
        cv.Set(interior,255)
        mask = pv.Image(mask)

        prev_im = forward(prev_im)
        prev_mask = forward(mask)
        

        next_im = self._imageBuffer[n-1]
        back = None
        for i in range(n-1,n/2,-1):
            if back == None:
                back = self._imageBuffer[i].to_prev
            else:
                back = back * self._imageBuffer[i].to_prev
        
        next_im = back(next_im)
        next_mask = back(mask)
        
        curr_im = self._imageBuffer[n/2]

                
        prevImg = prev_im.asMatrix2D()
        curImg  = curr_im.asMatrix2D()
        nextImg = next_im.asMatrix2D()
        prevMask = prev_mask.asMatrix2D()
        nextMask = next_mask.asMatrix2D()

        # Compute transformed images
        delta1 = sp.absolute(curImg - prevImg)   #frame diff 1
        delta2 = sp.absolute(nextImg - curImg)   #frame diff 2
        
        delta1 = sp.minimum(delta1,prevMask)
        delta2 = sp.minimum(delta2,nextMask)
        
        #use element-wise minimum of the two difference images, which is what
        # gets compared to threshold to yield foreground mask
        return sp.minimum(delta1, delta2)

def logloss(Y_true, Y_pred):
    epsilon = 1e-15
    pred = sp.maximum(epsilon, Y_pred)
    pred = sp.minimum(1-epsilon, Y_pred)
    ll = sum(Y_true*sp.log(pred) + sp.subtract(1,Y_true)*sp.log(sp.subtract(1,Y_pred)))
    ll = ll * -1.0/len(Y_true)
    return ll

def benjamini_hochberg_yekutieli(p_values=None,q_value=0.05,sort_idx=None,return_sort_idx=False):
    p_values = p_values.ravel()
    if sort_idx is None:
        sort_idx = sp.argsort(p_values)
        p_values = p_values[sort_idx]
    else:
        sort_idx = sort_idx.ravel()
        p_values = p_values[sort_idx]
    m = p_values.shape[0]
    idx_line = sp.arange(1,m+1)
    cV = (1.0/idx_line).sum()
    thr_line = (idx_line*q_value*cV)/float(m);
    thr_ind = sp.where(p_values<=thr_line)[0]
    if thr_ind.shape[0]==0:
        thr = 0.0;
    else:
        thr = p_values[thr_ind.max()]
    #adjust p_values
    p_values_adjusted = sp.ones(m)
    prev = 1.0
    for i in range(m,0,-1):
        p_values_adjusted[i-1] = sp.minimum(prev,p_values[i-1]*float(m)*cV/float(i))
        if p_values_adjusted[i-1]>1:
            p_values_adjusted[i-1]=1
        prev = p_values_adjusted[i-1]
    #resort pvalues
    p_tmp = p_values_adjusted.copy()
    p_values_adjusted[sort_idx] = p_tmp
    if return_sort_idx==True:
        return [thr,p_values_adjusted,sort_idx]        
    else:
        return [thr,p_values_adjusted]

def watershed(i_d, logger):
    # compute neighbours
    logger.info('Computing differences...')
    nbs = compute_neighbours(i_d)

    # compute min altitude map
    logger.info('Computing minimal altitude map...')
    minaltitude = nbs[0]
    for nb in nbs[1:]:
        minaltitude = scipy.minimum(minaltitude, nb)
    
    # watershed
    logger.info('Watershed \w minaltitude pre-computation...')
    result = scipy.zeros(i_d.shape, dtype=scipy.uint16)
    nb_labs = 0
    for x in range(result.shape[0]):
        for y in range(result.shape[1]):
            for z in range(result.shape[2]):
                if result[x,y,z] != 0: continue # 10%
                L, lab = stream_set(i_d, result, minaltitude, (x, y, z)) # 90%
                if -1 == lab:
                    nb_labs += 1
                    for p in L:
                        result[p] = nb_labs
                else:
                    for p in L: 
                        result[p] = lab
                        
    return result

def logloss(p, y):
    epsilon = 1e-15
    p = sp.maximum(epsilon, p)
    p = sp.minimum(1-epsilon, p)
    ll = sum(y*sp.log(p) + sp.subtract(1,y)*sp.log(sp.subtract(1,p)))
    ll = ll * -1.0/len(y)
    return ll

def find_largest_hole(parameters,ar):
    
    minimal_distances = []
    all_pixels = sp.array(range(len(ar)))
    empty_pixels = all_pixels[(ar[all_pixels] == parameters.badval)]
    if len(empty_pixels) == 0:
        print "no empty pixels"
        return 2*sp.pi/(6*parameters.nside)
        
        
    print "the number of empty pixels is", len(empty_pixels)
    nonempty_pixels = all_pixels[(ar[all_pixels] != parameters.badval)\
                                & (ar[all_pixels] != parameters.unseen)]
        
    for p in empty_pixels:

        minimal_distance = 3.14
        theta,phi = hp.pix2ang(parameters.nside,p)
        
        for p_i in nonempty_pixels:
            
            theta_i, phi_i = hp.pix2ang(parameters.nside,p_i)
            angular_distance = hp.rotator.angdist([theta,phi],[theta_i,phi_i])
            minimal_distance = sp.minimum(minimal_distance,angular_distance)
            
        minimal_distances.append(minimal_distance)
        
         
    radius_of_largest_hole = sp.amax(minimal_distances)

    print "The angular radius of largest hole = ", radius_of_largest_hole, "rad."
    
    return radius_of_largest_hole

def generateThumbnail(inputFile, thumbSize):
    global size
    # logging.debug('Input File: %s\n' % inputFile)
    # logging.debug('Ouput File: %s\n' % outputFile)
    # logging.debug('Thumb Size: %s\n' % thumbSize)

    h5f = tables.openFile(inputFile)

    dataSource = HDFDataSource.DataSource(inputFile, None)

    md = MetaData.genMetaDataFromSourceAndMDH(dataSource, MetaDataHandler.HDFMDHandler(h5f))

    xsize = h5f.root.ImageData.shape[1]
    ysize = h5f.root.ImageData.shape[2]

    if xsize > ysize:
        zoom = float(thumbSize) / xsize
    else:
        zoom = float(thumbSize) / ysize

    size = (int(xsize * zoom), int(ysize * zoom))

    im = h5f.root.ImageData[min(md.EstimatedLaserOnFrameNo + 10, (h5f.root.ImageData.shape[0] - 1)), :, :].astype("f")

    im = im.T - min(md.Camera.ADOffset, im.min())

    h5f.close()

    im = maximum(minimum(1 * (255 * im) / im.max(), 255), 0)

    return im.astype("uint8")

def llfun(act, pred):
    epsilon = 1e-15
    pred = sp.maximum(epsilon, pred)
    pred = sp.minimum(1 - epsilon, pred)
    ll = sum(act * sp.log(pred) + sp.subtract(1, act) * sp.log(sp.subtract(1, pred)))
    ll = ll * -1.0 / len(act)
    return ll

def binary_logloss(p, y):
    epsilon = 1e-15
    p = sp.maximum(epsilon, p)
    p = sp.minimum(1-epsilon, p)
    res = sum(y * sp.log(p) + sp.subtract(1, y) * sp.log(sp.subtract(1, p)))
    res *= -1.0/len(y)
    return res

def evaluate_ll(y, yhat):
    epsilon = 1e-15
    yhat = sp.maximum(epsilon, yhat)
    yhat = sp.minimum(1-epsilon, yhat)
    ll = sum(y*sp.log(yhat) + sp.subtract(1,y)*sp.log(sp.subtract(1,yhat)))
    ll = ll * -1.0/len(y)
    return ll

def entropyloss(act, pred):
    epsilon = 1e-15
    pred = sp.maximum(epsilon, pred)
    pred = sp.minimum(1-epsilon, pred)
    el = sum(act*sp.log10(pred) + sp.subtract(1,act)*sp.log10(sp.subtract(1,pred)))
    el = el * -1.0/len(act)
    return el

def logloss(act, pred):
    epsilon = 1e-4
    pred = sp.maximum(epsilon, pred)
    pred = sp.minimum(1-epsilon, pred)
    ll = -1.0/len(act) * sum(act*sp.log(pred) +
            sp.subtract(1,act)*sp.log(sp.subtract(1,pred)))
    return ll

def plotgpsonly(TEClist,gpslist,plotdir,m,ax,fig,latlim,lonlim):
    """ Makes a set of plots when only gps data is avalible."""
    maxplot = len(gpslist)
    strlen = int(sp.ceil(sp.log10(maxplot))+1)
    fmstr = '{0:0>'+str(strlen)+'}_'
    plotnum=0
    for gps_cur in gpslist:
        gpshands = []
        gpsmin = sp.inf
        gpsmax = -sp.inf
        for igpsn, (igps,igpslist) in enumerate(zip(TEClist,gps_cur)):
            print('Plotting GPS data from rec {0} of {1}'.format(igpsn,len(gps_cur)))
            # check if there's anything to plot
            if len(igpslist)==0:
                continue

            (sctter,scatercb) = scatterGD(igps,'alt',3.5e5,vbounds=[0,20],time = igpslist,gkey = 'vTEC',cmap='plasma',fig=fig,
                  ax=ax,title='',cbar=True,err=.1,m=m)
            gpsmin = sp.minimum(igps.times[igpslist,0].min(),gpsmin)
            gpsmax = sp.maximum(igps.times[igpslist,0].max(),gpsmax)

            gpshands.append(sctter)
        scatercb.set_label('vTEC in TECu')
        #change he z order

        print('Ploting {0} of {1} plots'.format(plotnum,maxplot))
        plt.savefig(os.path.join(plotdir,fmstr.format(plotnum)+'GPSonly.png'))
        plotnum+=1
        for i in reversed(gpshands):
            i.set_zorder(i.get_zorder()+1)

def histogram_intersection(h1, h2):  # 6 us @array, 30 us @list \w 100 bins
    """
    Calculate the common part of two histograms.
    The histogram intersection between two histograms \f$H\f$ and \f$H'\f$ of size \f$m\f$ is
    defined as
    \f[
        d_{\cap}(H, H') = \sum_{m=1}^M\min(H_m, H'_m)  
    \f]
    
    Attributes:
    - a real metric
    
    Attributes for normalized histograms:
    - \f$d(H, H')\in[0, 1]\f$
    - \f$d(H, H) = 1\f$
    - \f$d(H, H') = d(H', H)\f$
    
    Attributes for not-normalized histograms:
    - not applicable
    
    Attributes for not-equal histograms:
    - not applicable
    
    @param h1 the first histogram, normalized,
    @type h1 array-like sequence
    @param h2 the second histogram, normalized, same bins as h1
    @type h2 array-like sequence
    
    @return histogram intersection
    @rtype float
    """
    h1, h2 = __prepare_histogram(h1, h2)
    return scipy.sum(scipy.minimum(h1, h2))

def _set_reach_dist(setofobjects, point_index, epsilon):

    # Assumes that the query returns ordered (smallest distance first)
    # entries. This is the case for the balltree query...

    dists, indices = setofobjects.query(setofobjects.data[point_index],
                                        setofobjects._nneighbors[point_index])

    # Checks to see if there more than one member in the neighborhood ##
    if sp.iterable(dists):

        # Masking processed values ##
        # n_pr is 'not processed'
        n_pr = indices[(setofobjects._processed[indices] < 1)[0].T]
        rdists = sp.maximum(dists[(setofobjects._processed[indices] < 1)[0].T],
                            setofobjects.core_dists_[point_index])

        new_reach = sp.minimum(setofobjects.reachability_[n_pr], rdists)
        setofobjects.reachability_[n_pr] = new_reach

        # Checks to see if everything is already processed;
        # if so, return control to main loop ##
        if n_pr.size > 0:
            # Define return order based on reachability distance ###
            return n_pr[sp.argmin(setofobjects.reachability_[n_pr])]
        else:
            return point_index

def set_reach_dist(SetOfObjects, point_index, epsilon):

    """
    Sets reachability distance and ordering. This function is the primary workhorse of
    the OPTICS algorithm.
    
    SetofObjects: Instantiated and prepped instance of 'setOfObjects' class
    epsilon: Determines maximum object size that can be extracted. Smaller epsilons
        reduce run time. (float)

    """
    
    row = [SetOfObjects.data[point_index,:]]
    indices = np.argsort(row)
    distances = np.sort(row)

    if scipy.iterable(distances):

        unprocessed = indices[(SetOfObjects._processed[indices] < 1)[0].T]
        rdistances = scipy.maximum(distances[(SetOfObjects._processed[indices] < 1)[0].T],
            SetOfObjects._core_dist[point_index])
        SetOfObjects._reachability[unprocessed] = scipy.minimum(
            SetOfObjects._reachability[unprocessed], rdistances)

        if unprocessed.size > 0:
            return unprocessed[np.argsort(np.array(SetOfObjects._reachability[
                unprocessed]))[0]]
        else:
            return point_index
    else:
        return point_index