def nms(dets, thresh):
    x1 = dets[:, 0]
    y1 = dets[:, 1]
    x2 = dets[:, 2]
    y2 = dets[:, 3]
    scores = dets[:, 4]

    areas = (x2 - x1 + 1) * (y2 - y1 + 1)
    order = scores.argsort()[::-1]

    keep = []
    while order.size > 0:
        i = order[0]
        keep.append(i)
        xx1 = np.maximum(x1[i], x1[order[1:]])
        yy1 = np.maximum(y1[i], y1[order[1:]])
        xx2 = np.minimum(x2[i], x2[order[1:]])
        yy2 = np.minimum(y2[i], y2[order[1:]])

        w = np.maximum(0.0, xx2 - xx1 + 1)
        h = np.maximum(0.0, yy2 - yy1 + 1)
        inter = w * h
        ovr = inter / ((areas[i] + areas[order[1:]] - inter)+0.0000001)

        inds = np.where(ovr <= thresh)[0]
        order = order[inds + 1]

    return keep

	def ParseMDH (self, n = 0, dry = False):
		''' Parse one MDH entry '''

		self.bm = struct.unpack_from ('2I', self.buf, n +  20)
		if (self.bm[0] & SYNCDATA):
			return n + MDHSIZE
		
		if (self.bm[0] & NOISEADJSCAN):
			self.lc  = struct.unpack_from ('%dH' % 16, self.buf, n + 28) # loop counters
			self.ch  = struct.unpack_from ('H', self.buf, n + 124)
	
			if (dry):
				if (self.noisedims[0] == 0):
					self.noisedims[0] = self.lc[0]
					self.noisedims[1] = struct.unpack_from ('H', self.buf, n + 30)[0]
					self.noisencolb = self.noisedims[md.COL] * self.ndds
					self.noisedims[2:16] = np.array(self.lc[2:16])
				else:
					self.noisedims[2:16] = np.maximum (self.noisedims[2:16],self.lc[2:16])
			return n + MDHSIZE
				
		self.lc  = struct.unpack_from ('%dH' % 16, self.buf, n + 28) # loop counters
		self.ch  = struct.unpack_from ('H', self.buf, n + 124)

		if (dry):
			if (self.dims[0] == 0):
				self.dims[0] = self.lc[0]
				self.dims[1] = struct.unpack_from ('H', self.buf, n + 30)[0]
				self.ncolb = self.dims[md.COL] * self.ds
				self.dims[2:16] = np.array(self.lc[2:16])
			else:
				self.dims[2:16] = np.maximum (self.dims[2:16],self.lc[2:16])

		return n + MDHSIZE

    def autoRange(self,zmin=None,zmax=None,coordList=None):
        """ determine X,Y,Z limits
        """

        if coordList==None:
            coordList = self.coordList

        # autorange in X,Y
        minX = 1.e50
        maxX = -1.e50
        minY = 1.e50
        maxY = -1.e50
        minZ = 1.e50
        maxZ = -1.e50
        for iCoord in coordList:

            #if self.interpPresent:
            #    X,Y = self.interpGrids[iCoord]
            #    Z = self.interpValues[iCoord]
            #else:
            X,Y,Z = self.getXYZpoints(coordList=coordList)

            minX = numpy.minimum(minX,X.min())
            maxX = numpy.maximum(maxX,X.max())
            minY = numpy.minimum(minY,Y.min())
            maxY = numpy.maximum(maxY,Y.max())
            minZ = numpy.minimum(minZ,Z.min())
            maxZ = numpy.maximum(maxZ,Z.max())

        if zmin!=None:
            minZ = zmin
        if zmax!=None:
            maxZ = zmax

        return minX,maxX,minY,maxY,minZ,maxZ

    def update_swe(self):

        #--------------------------------------------------------
        # Note: The Meteorology component uses air temperature
        # to compute P_rain (precip that falls as liquid) and
        # P_snow (precip that falls as snow or ice) separately.
        # P_snow = (self.P * (self.T_air <= 0)) 
        #----------------------------------------------------------
        # Note: This method must be written to work regardless
        # of whether P_rain and T are scalars or grids. (3/14/07)
        #------------------------------------------------------------
        # If P or T_air is a grid, then h_swe and h_snow are grids.
        # This is set up in initialize_computed_vars().
        #------------------------------------------------------------
      
        #------------------------------------------------
        # Increase snow water equivalent due to snowfall
        #------------------------------------------------
        # Meteorology and Channel components may have
        # different time steps, but then self.P_snow
        # will be a time-interpolated value.
        #------------------------------------------------
        dh1_swe  = (self.P_snow * self.dt)
        self.h_swe  += dh1_swe

        #------------------------------------------------
        # Decrease snow water equivalent due to melting
        # Note that SM depends partly on h_snow.
        #------------------------------------------------
        dh2_swe    = self.SM * self.dt
        self.h_swe -= dh2_swe
        np.maximum(self.h_swe, np.float64(0), self.h_swe)  # (in place)

def _joint_probabilities(distances, desired_perplexity, verbose):
    """Compute joint probabilities p_ij from distances.

    Parameters
    ----------
    distances : array, shape (n_samples * (n_samples-1) / 2,)
        Distances of samples are stored as condensed matrices, i.e.
        we omit the diagonal and duplicate entries and store everything
        in a one-dimensional array.

    desired_perplexity : float
        Desired perplexity of the joint probability distributions.

    verbose : int
        Verbosity level.

    Returns
    -------
    P : array, shape (n_samples * (n_samples-1) / 2,)
        Condensed joint probability matrix.
    """
    # Compute conditional probabilities such that they approximately match
    # the desired perplexity
    distances = distances.astype(np.float32, copy=False)
    conditional_P = _utils._binary_search_perplexity(
        distances, None, desired_perplexity, verbose)
    P = conditional_P + conditional_P.T
    sum_P = np.maximum(np.sum(P), MACHINE_EPSILON)
    P = np.maximum(squareform(P) / sum_P, MACHINE_EPSILON)
    return P

def shifted_corr(reference, image, displacement):
    """Calculate the correlation between the reference and the image shifted
    by the given displacement.

    Parameters
    ----------
    reference : np.ndarray
    image : np.ndarray
    displacement : np.ndarray

    Returns
    -------
    correlation : float

    """

    ref_cuts = np.maximum(0, displacement)
    ref = reference[ref_cuts[0]:, ref_cuts[1]:, ref_cuts[2]:]
    im_cuts = np.maximum(0, -displacement)
    im = image[im_cuts[0]:, im_cuts[1]:, im_cuts[2]:]
    s = np.minimum(im.shape, ref.shape)
    ref = ref[:s[0], :s[1], :s[2]]
    im = im[:s[0], :s[1], :s[2]]
    ref -= nanmean(ref.reshape(-1, ref.shape[-1]), axis=0)
    ref = np.nan_to_num(ref)
    im -= nanmean(im.reshape(-1, im.shape[-1]), axis=0)
    im = np.nan_to_num(im)
    assert np.all(np.isfinite(ref)) and np.all(np.isfinite(im))
    corr = nanmean(
        [old_div(np.sum(i * r), np.sqrt(np.sum(i * i) * np.sum(r * r))) for
         i, r in zip(np.rollaxis(im, -1), np.rollaxis(ref, -1))])
    return corr

    def _init_energy(self, pc):
        if pc is self._pc:
            return
        self.set_transform(self._t, pc)
        self._pc = pc
        self._res[:] = self.data[:, self._t] - self.mu[:]
        self._V = np.maximum(self.offset + np.mean(self._res ** 2), SMALL)
        self._res0[:] = self.data[:, self._t] - self.mu0
        self._V0 = np.maximum(self.offset0 + np.mean(self._res0 ** 2), SMALL)

        if self.use_derivatives:
            # linearize the data wrt the transform parameters
            # use the auxiliary array to save the current resampled data
            self._aux[:] = self.data[:, self._t]
            basis = np.eye(6)
            for j in range(pc.size):
                self.set_transform(self._t, pc + self.stepsize * basis[j])
                self.A[:, j] = (self.data[:, self._t] - self._aux)\
                    / self.stepsize
            self.transforms[self._t].param = pc
            self.data[:, self._t] = self._aux[:]
            # pre-compute gradient and hessian of numerator and
            # denominator
            c = 2 / float(self.data.shape[0])
            self._dV = c * np.dot(self.A.T, self._res)
            self._dV0 = c * np.dot(self.A.T, self._res0)
            self._H = c * np.dot(self.A.T, self.A)

 def _logpmf(self, x, mu, alpha, p):
     mu_p = mu ** (p - 1.)
     a1 = np.maximum(np.nextafter(0, 1), 1 + alpha * mu_p)
     a2 = np.maximum(np.nextafter(0, 1), mu + (a1 - 1.) * x)
     logpmf_ = np.log(mu) + (x - 1.) * np.log(a2)
     logpmf_ -=  x * np.log(a1) + gammaln(x + 1.) + a2 / a1
     return logpmf_

def nms2d(boxes, overlap=0.3):
    """Compute the nms given a set of scored boxes,
    as numpy array with 5 columns <x1> <y1> <x2> <y2> <score>
    return the indices of the tubelets to keep
    """

    if boxes.size == 0:
        return np.array([],dtype=np.int32)

    x1 = boxes[:, 0]
    y1 = boxes[:, 1]
    x2 = boxes[:, 2]
    y2 = boxes[:, 3]

    scores = boxes[:, 4]
    areas = (x2-x1+1) * (y2-y1+1)
    I = np.argsort(scores)
    indices = np.zeros(scores.shape, dtype=np.int32)

    counter = 0
    while I.size > 0:
        i = I[-1]
        indices[counter] = i
        counter += 1

        xx1 = np.maximum(x1[i],x1[I[:-1]])
        yy1 = np.maximum(y1[i],y1[I[:-1]])
        xx2 = np.minimum(x2[i],x2[I[:-1]])
        yy2 = np.minimum(y2[i],y2[I[:-1]])

        inter = np.maximum(0.0, xx2 - xx1 + 1) * np.maximum(0.0, yy2 - yy1 + 1)
        iou = inter / (areas[i] + areas[I[:-1]] - inter)
        I = I[np.where(iou <= overlap)[0]]

    return indices[:counter]

    def pdf(self, t):
        """
        Implementing the code distributed with Logan et al. 2013 using
        the reparametrization given above.

        Also, the constant theta is added as usual.
        """
        t=np.maximum(t-self.theta, 1e-5) # absorbed into pdf
        sqrt_t=np.sqrt(t)

        # reparametrization
        a=self.A/2.0
        k=self.alpha-self.A/2.0
        l=self.gamma

        if self.A<1e-10: # this is the solution without starting-point variability
            r=self.alpha/(np.sqrt(2*np.pi*(t**3)))*np.exp(- ((self.alpha-self.gamma*t)**2)/(2*t))
        elif self.gamma<1e-10:
            r=np.exp( -.5*( np.log(2)+np.log(np.pi)+np.log(t))
                      + np.log( np.exp(-( (k-a)**2/(2*t)))-np.exp(-( (k+a)**2/(2*t) )) )
                      - np.log(2) - np.log(a) )
        else:
            r=np.exp( np.log( (np.exp(- (a-k+t*l)**2/(2*t) )-np.exp(- (a+k-t*l)**2/(2*t) ))/np.sqrt(2*np.pi*t)
                              + np.exp(np.log(.5)+np.log(l))*( 2*pnormP( (-k+a)/sqrt_t + sqrt_t*l)-1
                                                               + 2*pnormP( (k+a)/sqrt_t - sqrt_t*l)-1) )
                      - np.log(2) - np.log(a))

        return np.maximum(0.0, np.where( np.isnan(r), 0, r))

    def cdf(self,t):
        t=np.maximum(t-self.theta, 1e-5) # absorbed into cdf

        sqrt_t=np.sqrt(t)

        # reparametrization
        a=self.A/2.0
        k=self.alpha-self.A/2.0
        l=self.gamma

        if self.A<1e-10: # this is the solution without starting-point variability
            r=pnormP( (self.gamma*t-self.alpha)/sqrt_t)+np.exp(2*self.alpha*self.gamma)*pnormP(-(self.gamma*t+self.alpha)/(sqrt_t))
        elif self.gamma<1e-10:
            r=(( -(k+a)*(2*pnormP( (k+a)/sqrt_t)-1)
                 -(k-a)*(2*pnormP(-(k-a)/sqrt_t)-1))/(2*a)) \
              + (1 + np.exp(-.5*(k-a)**2/t - .5*np.log(2) - .5*np.log(np.pi) + .5*np.log(t) - np.log(a))
                 -   np.exp(-.5*(k+a)**2/t - .5*np.log(2) - .5*np.log(np.pi) + .5*np.log(t) - np.log(a)))
        else:
            t1=np.exp( .5*np.log(t)-.5*np.log(2*np.pi) ) * (  np.exp( -((k-a-t*l)**2/t)/2.)
                                                            - np.exp( -((k+a-t*l)**2/t)/2.) ) # ok
            t2=a+(   np.exp(2*l*(k+a)+np.log(pnormP(-(k+a+t*l)/sqrt_t)))
                   - np.exp(2*l*(k-a)+np.log(pnormP(-(k-a+t*l)/sqrt_t))) )/(2*l) # ok
            t4= (.5*(t*l-a-k+.5/l)) * ( 2*pnormP((k+a)/sqrt_t-sqrt_t*l)-1) \
               + .5*(k-a-t*l-.5/l)*( 2*pnormP((k-a)/sqrt_t-sqrt_t*l)-1)
            r=(t4+t2+t1)/(2*a)


        return np.minimum( np.maximum( 0., np.where( np.isnan(r), 0, r) ), 1.)

def nms(boxes, threshold, method):
    if boxes.size == 0:
        return np.empty((0, 3))
    x1 = boxes[:, 0]
    y1 = boxes[:, 1]
    x2 = boxes[:, 2]
    y2 = boxes[:, 3]
    s = boxes[:, 4]
    area = (x2 - x1 + 1) * (y2 - y1 + 1)
    I = np.argsort(s)
    pick = np.zeros_like(s, dtype=np.int16)
    counter = 0
    while I.size > 0:
        i = I[-1]
        pick[counter] = i
        counter += 1
        idx = I[0:-1]
        xx1 = np.maximum(x1[i], x1[idx])
        yy1 = np.maximum(y1[i], y1[idx])
        xx2 = np.minimum(x2[i], x2[idx])
        yy2 = np.minimum(y2[i], y2[idx])
        w = np.maximum(0.0, xx2 - xx1 + 1)
        h = np.maximum(0.0, yy2 - yy1 + 1)
        inter = w * h
        if method is 'Min':
            o = inter / np.minimum(area[i], area[idx])
        else:
            o = inter / (area[i] + area[idx] - inter)
        I = I[np.where(o <= threshold)]
    pick = pick[0:counter]
    return pick

    def test_elementwise_max_grad(self, n, m, d, gc, dc):
        go = np.random.rand(n, m, d).astype(np.float32)
        X = np.random.rand(n, m, d).astype(np.float32)
        Y = np.random.rand(n, m, d).astype(np.float32)
        Z = np.random.rand(n, m, d).astype(np.float32)
        mx = np.maximum(np.maximum(X, Y), Z)
        inputs = [mx, go, X, Y, Z]

        def max_grad_op(mx, go, X, Y, Z):
            def mx_grad(a):
                return go * (mx == a)

            return [mx_grad(a) for a in [X, Y, Z]]

        op = core.CreateOperator(
            "MaxGradient",
            ["mx", "go", "X", "Y", "Z"],
            ["gX", "gY", "gZ"]
        )

        self.assertReferenceChecks(
            device_option=gc,
            op=op,
            inputs=inputs,
            reference=max_grad_op,
        )
        self.assertDeviceChecks(dc, op, inputs, [0, 1, 2])

def prox_l1(Y, alpha, n_orient):
    """proximity operator for l1 norm with multiple orientation support

    L2 over orientation and L1 over position (space + time)

    Example
    -------
    >>> Y = np.tile(np.array([1, 2, 3, 2, 0], dtype=np.float), (2, 1))
    >>> Y = np.r_[Y, np.zeros_like(Y)]
    >>> print Y
    [[ 1.  2.  3.  2.  0.]
     [ 1.  2.  3.  2.  0.]
     [ 0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.]]
    >>> Yp, active_set = prox_l1(Y, 2, 2)
    >>> print Yp
    [[ 0.          0.58578644  1.58578644  0.58578644  0.        ]
     [ 0.          0.58578644  1.58578644  0.58578644  0.        ]]
    >>> print active_set
    [ True  True False False]
    """
    n_positions = Y.shape[0] // n_orient
    norms = np.sqrt(np.sum((np.abs(Y) ** 2).T.reshape(-1, n_orient), axis=1))
    # Ensure shrink is >= 0 while avoiding any division by zero
    shrink = np.maximum(1.0 - alpha / np.maximum(norms, alpha), 0.0)
    shrink = shrink.reshape(-1, n_positions).T
    active_set = np.any(shrink > 0.0, axis=1)
    shrink = shrink[active_set]
    if n_orient > 1:
        active_set = np.tile(active_set[:, None], [1, n_orient]).ravel()
    Y = Y[active_set]
    if len(Y) > 0:
        for o in range(n_orient):
            Y[o::n_orient] *= shrink
    return Y, active_set

def test_maximum_minimum_scalar():
    data1 = mx.symbol.Variable('data')
    shape = (3, 4)
    data_tmp1 = np.random.rand(3,4)
    data_tmp1[:] = 2
 
    arr_data1 = mx.nd.array(data_tmp1)
    arr_grad1 = mx.nd.empty(shape)

    test =  mx.sym.maximum(data1,3) + mx.sym.maximum(9,data1) + mx.sym.minimum(5,data1) + mx.sym.minimum(data1,4)
    exe_test = test.bind(mx.cpu(), args=[arr_data1], args_grad=[arr_grad1])
    exe_test.forward()
    out = exe_test.outputs[0].asnumpy()
    npout =  np.maximum(data_tmp1,3) + np.maximum(9,data_tmp1) + np.minimum(5,data_tmp1) + np.minimum(data_tmp1,4)
    assert reldiff(out, npout) < 1e-6

    out_grad = mx.nd.empty(shape)
    out_grad[:] = 2
    exe_test.backward(out_grad)
    
    npout_grad = np.ones(shape)
    npout_grad[:] = 2
    mask1 = (data_tmp1 > 3).astype('float')
    mask2 = (9 > data_tmp1).astype('float')
    mask3 = (5 < data_tmp1).astype('float')
    mask4 = (data_tmp1 < 4).astype('float')
    npout_grad1 = npout_grad * mask1 + (npout_grad - npout_grad * mask2) + (npout_grad - npout_grad * mask3) + npout_grad * mask4
    
    assert reldiff(arr_grad1.asnumpy(), npout_grad1) < 1e-6

def rp_gumbel_original(p_zero, loc, scale, flvol, max_return_period=1e9):
    """
    Transforms a unique, or array of flood volumes into the belonging return
    periods, according to gumbel parameters (belonging to non-zero part of the
    distribution) and a zero probability
    Inputs:
        p_zero:        probability that flood volume is zero
        loc:           Gumbel location parameter (of non-zero part of distribution)
        scale:         Gumbel scale parameter (of non-zero part of distribution)
        flvol:         Flood volume that will be transformed to return period
        max_return_period: maximum return period considered. This maximum is needed to prevent that floating point
                        precision becomes a problem (default: 1e9)
    This function is copied from: https://repos.deltares.nl/repos/Hydrology/trunk/GLOFRIS/src/rp_bias_corr.py
    """
    
    np.seterr(divide='ignore')
    np.seterr(invalid='ignore')
    max_p = 1-1./max_return_period
    max_p_residual = np.minimum(np.maximum((max_p-np.float64(p_zero))/(1-np.float64(p_zero)), 0), 1)
    max_reduced_variate = -np.log(-np.log(np.float64(max_p_residual)))
    # compute the gumbel reduced variate belonging to the Gumbel distribution (excluding any zero-values)
    # make sure that the reduced variate does not exceed the one, resembling the 1,000,000 year return period
    reduced_variate = np.minimum((flvol-loc)/scale, max_reduced_variate)
    # reduced_variate = (flvol-loc)/scale
    # transform the reduced variate into a probability (residual after removing the zero volume probability)
    p_residual = np.minimum(np.maximum(np.exp(-np.exp(-np.float64(reduced_variate))), 0), 1)
    # tranform from non-zero only distribution to zero-included distribution
    p = np.minimum(np.maximum(p_residual*(1-p_zero) + p_zero, p_zero), max_p)  # Never larger than max_p
    # transform into a return period    
    return_period = 1./(1-p)
    test_p = p == 1    
    return return_period, test_p

def psfplots():
	tpsf = wise.get_psf_model(1, pixpsf=True)
	
	psfp = tpsf.getPointSourcePatch(0, 0)
	psf = psfp.patch
	
	psf /= psf.sum()
	
	plt.clf()
	plt.imshow(np.log10(np.maximum(1e-5, psf)), interpolation='nearest', origin='lower')
	plt.colorbar()
	ps.savefig()
	
	h,w = psf.shape
	cx,cy = w/2, h/2
	
	X,Y = np.meshgrid(np.arange(w), np.arange(h))
	R = np.sqrt((X - cx)**2 + (Y - cy)**2)
	plt.clf()
	plt.semilogy(R.ravel(), psf.ravel(), 'b.')
	plt.xlabel('Radius (pixels)')
	plt.ylabel('PSF value')
	plt.ylim(1e-8, 1.)
	ps.savefig()
	
	plt.clf()
	plt.loglog(R.ravel(), psf.ravel(), 'b.')
	plt.xlabel('Radius (pixels)')
	plt.ylabel('PSF value')
	plt.ylim(1e-8, 1.)
	ps.savefig()
	
	print('PSF norm:', np.sqrt(np.sum(np.maximum(0, psf)**2)))
	print('PSF max:', psf.max())

def _beta_divergence_dense(X, W, H, beta):
    """Compute the beta-divergence of X and W.H for dense array only.

    Used as a reference for testing nmf._beta_divergence.
    """
    if isinstance(X, numbers.Number):
        W = np.array([[W]])
        H = np.array([[H]])
        X = np.array([[X]])

    WH = np.dot(W, H)

    if beta == 2:
        return squared_norm(X - WH) / 2

    WH_Xnonzero = WH[X != 0]
    X_nonzero = X[X != 0]
    np.maximum(WH_Xnonzero, 1e-9, out=WH_Xnonzero)

    if beta == 1:
        res = np.sum(X_nonzero * np.log(X_nonzero / WH_Xnonzero))
        res += WH.sum() - X.sum()

    elif beta == 0:
        div = X_nonzero / WH_Xnonzero
        res = np.sum(div) - X.size - np.sum(np.log(div))
    else:
        res = (X_nonzero ** beta).sum()
        res += (beta - 1) * (WH ** beta).sum()
        res -= beta * (X_nonzero * (WH_Xnonzero ** (beta - 1))).sum()
        res /= beta * (beta - 1)

    return res

def fuzzy_c_means(points, num_centers, m=2., tol=1e-4, max_iter=100,
                  verbose=False):
  '''Uses Fuzzy C-Means to downsample `points`.
  m : aggregation parameter >1, larger implies smoother clusters
  Returns indices of downsampled points.
  '''
  num_points = points.shape[0]
  if num_centers >= num_points:
    return np.arange(num_points)
  # randomly initialize cluster assignments matrix
  assn = np.random.random((points.shape[0], num_centers))
  # iterate assignments until they converge
  for i in range(max_iter):
    # compute centers
    w = assn ** m
    w /= w.sum(axis=0)
    centers = w.T.dot(points)
    # calculate new assignments
    d = pairwise_distances(points, centers)
    d **= 2. / (m - 1)
    np.maximum(d, 1e-10, out=d)
    new_assn = 1. / np.einsum('ik,ij->ik', d, 1./d)
    # check for convergence
    change = np.linalg.norm(new_assn - assn)
    if verbose:
      print('At iteration %d: change = %g' % (i+1, change))
    if change < tol:
      break
    assn = new_assn
  else:
    warnings.warn("fuzzy_c_means didn't converge in %d iterations" % max_iter)
  # find points closest to the selected cluster centers
  return d.argmin(axis=0)

def clip_lower(arr,lower_bound):
    """
    In-place, one-sided version of numpy.clip().

    i.e. numpy.clip(arr,a_min=lower_bound,out=arr) if it existed.
    """
    maximum(arr,lower_bound,arr)