def biweight(x, cst):
    """
    Computes the biweight average and midvariance for a given 1D array.
    Returns a tuple (biweight mean, biweight variance).

    Parameters
    ----------
    x: {ndarray}
        Input Array
    cst : {float} 
        Parameter controlling how outliers are censored.

    Notes
    -----
    The function is restricted to 1D data only.
    """
    assert (x.ndim == 1, "1D array only !")
    xmed = ma.median(x, 0)
    manom = x - xmed
    mad = ma.median(ma.absolute(manom))
    u_i = (manom/float(cst*mad))
    u_i *= ma.less_equal(ma.absolute(u_i), 1.).astype(float)
    w_i = (1-u_i**2)
    if ma.count(w_i) > 0:
        biw_m = xmed + ma.sum(manom * w_i**2)/ma.sum(w_i**2)
    else:
        biw_m = xmed
    biw_sd = ma.sqrt(ma.count(x)*ma.sum(manom**2 * w_i**4))
    biw_sd *= 1./ma.absolute(ma.sum(w_i * (1-5*u_i**2)))
    return (biw_m, biw_sd.item())

 def test_remove_clip_box(self) : 
     """test that we can remove the clip box once set."""
     self.gca.set_clip_box(-90,-75,-180,-165)
     testvec = self.gca.compress(self.grid)
     self.assertTrue(self.gca.is_masked())
     self.assertEqual(ma.count(testvec), 2)
     
     self.gca.remove_mask()
     self.assertFalse(self.gca.is_masked())
     testvec = self.gca.compress(self.grid)
     self.assertEqual(ma.count(testvec), 4)

def get_dims(squaremask):
    """ return number of unmasked pixels along horizontal and vertical profiles
        through the center of the image (=size of mask)
    """

    dimx,dimy = np.shape(squaremask)
    halfpointval = squaremask[dimx/2,dimy/2]
    
    horvec = squaremask[dimx/2,:]
    vervec = squaremask[:,dimy/2]

    return ma.count(horvec), ma.count(vervec)

 def ReturnIs(marray, xc, yc, rbins, oversamp, total=0):
     """Returns the average quantities at different radius of a
       masked array. total=1 just return a total count within rbins""" 
     SizeY = marray.shape[0] # Y size
     SizeX = marray.shape[1] # X size
     x = np.reshape(np.arange(SizeX * SizeY), (SizeY, SizeX)) % SizeX
     x = x.astype(np.float32)
     x /= oversamp 
     y = np.reshape(np.arange(SizeX * SizeY), (SizeY, SizeX)) / SizeX
     y = y.astype(np.float32)
     y /= oversamp
     rx = (x - xc)* self.co + (y - yc) * self.si
     ry = (xc - x) * self.si + (y - yc) * self.co
     r = np.sqrt(rx**2.0 + ry**2.0 / self.one_minus_eg_sq)
     if total:
         con = (r < rbins)
         TotI = marray[con].sum()
         TotN = ma.count(marray[con]) / (oversamp * oversamp * 1.0)
         return TotI, TotN
     else:
         AvgIAtR = []
         AvgIInR = []
         IInRArr = []
         RArr = []
         NInRArr = []
         letbreak = 0 # this will be used to break the loop if eta is 
                      # less than 0.2 for 20 ri's
         for ri in rbins:
             con = (r > ri - 1/oversamp) & (r < ri + 1/oversamp)
             IAtR = marray[con].sum()
             NAtR = ma.count(marray[con]) * 1.0 
             con = (r < ri)
             IInR = marray[con].sum()
             NInR = ma.count(marray[con]) * 1.0
             if NAtR == 0 or NInR == 0 or ri > 20 and NAtR < 30 or \
                ri > 20 and NInR < 30:
                 pass
             else:
                 AvgIAtR.append(IAtR / NAtR)
                 AvgIInR.append(IInR / NInR)
                 IInRArr.append(IInR)
                 RArr.append(ri)
                 NInRArr.append(NInR)
                 if IAtR * NInR / (NAtR * IInR) < 0.2:
                     letbreak += 1
             if letbreak > 20:
                 break
         AvgIAtR = np.asarray(AvgIAtR)
         AvgIInR = np.asarray(AvgIInR)
         IInRArr = np.asarray(IInRArr)
         RArr = np.asarray(RArr) 
         NInRArr = np.asarray(NInRArr) / (oversamp * oversamp * 1.0)
         return AvgIAtR, AvgIInR, IInRArr, RArr, NInRArr

def calc_bulk_stats(stats_found, num_pts_section):
	
	if (stats_found==1):	
		ice_area = ma.count(elevation2d)*(xy_res**2)
		ridge_area_all = ma.count(elevation2d_ridge_ma)*(xy_res**2)
		mean_ridge_height_all = np.mean(elevation2d_ridge_ma) - level_elev
		mean_ridge_heightL = np.mean(ridge_height_mesh)
		ridge_areaL = ma.count(ridge_height_mesh)*(xy_res**2)
		return [mean_x, mean_y, ice_area, num_ridges, ridge_area_all, ridge_areaL, mean_ridge_height_all, mean_ridge_heightL, mean_alt, mean_pitch, mean_roll, mean_vel, num_pts_section, stats_found]
	elif (stats_found==0):
		#a = ma.masked_all((0))
		#masked_val = mean(a)
		return [mean_x, mean_y, -999, 0,-999, -999, -999, -999, mean_alt, mean_pitch, mean_roll, mean_vel, num_pts_section, stats_found]

    def __init__(self, MetricTable):

        # Create empty ratio table
        nprobs = MetricTable.nprobs
        nsolvs = MetricTable.nsolvs
        self.ratios = ma.masked_array(1.0 * ma.zeros((nprobs+1, nsolvs)))

        # Compute best relative performance ratios across
        # solvers for each problem
        for prob in range(nprobs):
            metrics  = MetricTable.prob_mets(prob)
            best_met = ma.minimum(metrics)
	    if (ma.count(metrics)==nsolvs and
                ma.maximum(metrics)<=opts.minlimit):
                self.ratios[prob+1,:] = 1.0;
	    else:
                self.ratios[prob+1,:] = metrics * (1.0 / best_met)

        # Sort each solvers performance ratios
        for solv in range(nsolvs):
            self.ratios[:,solv] = ma.sort(self.ratios[:,solv])

        # Compute largest ratio and use to replace failures entries
        self.maxrat = ma.maximum(self.ratios)
        self.ratios = ma.filled(self.ratios, 1.01 * self.maxrat)

def _pivot_col(T, tol=1.0E-12, bland=False):
    """
    Given a linear programming simplex tableau, determine the column
    of the variable to enter the basis.

    Parameters
    ----------
    T : 2D ndarray
        The simplex tableau.
    tol : float
        Elements in the objective row larger than -tol will not be considered
        for pivoting.  Nominally this value is zero, but numerical issues
        cause a tolerance about zero to be necessary.
    bland : bool
        If True, use Bland's rule for selection of the column (select the
        first column with a negative coefficient in the objective row, regardless
        of magnitude).

    Returns
    -------
    status: bool
        True if a suitable pivot column was found, otherwise False.  A return
        of False indicates that the linear programming simplex algorithm is complete.
    col: int
        The index of the column of the pivot element.  If status is False, col
        will be returned as nan.
    """
    ma = np.ma.masked_where(T[-1, :-1] >= -tol, T[-1, :-1], copy=False)
    if ma.count() == 0:
        return False, np.nan
    if bland:
        return True, np.where(ma.mask == False)[0][0]
    return True, np.ma.where(ma == ma.min())[0][0]

 def get_array_attributes(self):
     lat_end = 45
     lon_end = 180
     ### 5 day data or daily data
     time_end = 730
     # time_end = 3650
     ### Choose array (decadel mean, annual mean or all data)
     ### All Data
     self.what_data = "AllData"
     self.array = load_cflux_masked.load_file(time_end=time_end, lat_end=lat_end, lon_end=lon_end)
     ### decadel mean
     # ~ self.what_data = 'DecadalMean'
     # ~ self.array = ma.mean(self.array, axis=0)
     # ~ self.array = np.reshape(self.array, (1, lat_end, lon_end))
     ### annual mean
     self.what_data = "AnnualCycle"
     self.array = ma.mean(np.split(self.array, 10, axis=0), axis=0)
     ###
     self.array_shape = np.shape(self.array)
     print self.array_shape
     self.count_non_masked = ma.count(self.array)
     self.time_len = self.array_shape[0]  # need to set interpolated and masked array time_end to be equal NB!!!
     self.lat_len = self.array_shape[1]
     self.lon_len = self.array_shape[2]
     self.string_length = len(bin(self.count_non_masked)[2:])
     for item in itertools.product(range(self.lat_len), range(self.lon_len)):
         self.actual_data_dict[item] = np.std(self.array[:, item[0], item[1]])

 def test_ll_corner(self) : 
     """test that we filter out everything but ll corner"""
     self.gca.set_clip_box(-90,-82,-180,-173)
     testvec = self.gca.compress(self.grid)
     self.assertTrue(self.gca.is_masked())
     testmask = self.gca.get_vec_mask()
     self.assertEqual(np.count_nonzero(testmask), 3)
     self.assertEqual(ma.count(testvec), 1)

 def add_chunk(self, chunk):
     if self.masked:
         ma.sum(chunk, axis=self.axis, out=self.temp)
         self.running_total += self.temp.filled(0)
         self.running_count += ma.count(chunk, axis=self.axis)
     else:
         np.sum(chunk, axis=self.axis, out=self.temp)
         self.running_total += self.temp

 def test_xtestCount(self):
     # Test count
     ott = array([0., 1., 2., 3.], mask=[1, 0, 0, 0])
     assert_(count(ott).dtype.type is np.intp)
     assert_equal(3, count(ott))
     assert_equal(1, count(1))
     assert_(eq(0, array(1, mask=[1])))
     ott = ott.reshape((2, 2))
     assert_(count(ott).dtype.type is np.intp)
     assert_(isinstance(count(ott, 0), np.ndarray))
     assert_(count(ott).dtype.type is np.intp)
     assert_(eq(3, count(ott)))
     assert_(getmask(count(ott, 0)) is nomask)
     assert_(eq([1, 2], count(ott, 0)))

    def test_testAverage2(self):
        # More tests of average.
        w1 = [0, 1, 1, 1, 1, 0]
        w2 = [[0, 1, 1, 1, 1, 0], [1, 0, 0, 0, 0, 1]]
        x = arange(6)
        assert_(allclose(average(x, axis=0), 2.5))
        assert_(allclose(average(x, axis=0, weights=w1), 2.5))
        y = array([arange(6), 2.0 * arange(6)])
        assert_(allclose(average(y, None),
                                 np.add.reduce(np.arange(6)) * 3. / 12.))
        assert_(allclose(average(y, axis=0), np.arange(6) * 3. / 2.))
        assert_(allclose(average(y, axis=1),
                                 [average(x, axis=0), average(x, axis=0)*2.0]))
        assert_(allclose(average(y, None, weights=w2), 20. / 6.))
        assert_(allclose(average(y, axis=0, weights=w2),
                                 [0., 1., 2., 3., 4., 10.]))
        assert_(allclose(average(y, axis=1),
                                 [average(x, axis=0), average(x, axis=0)*2.0]))
        m1 = zeros(6)
        m2 = [0, 0, 1, 1, 0, 0]
        m3 = [[0, 0, 1, 1, 0, 0], [0, 1, 1, 1, 1, 0]]
        m4 = ones(6)
        m5 = [0, 1, 1, 1, 1, 1]
        assert_(allclose(average(masked_array(x, m1), axis=0), 2.5))
        assert_(allclose(average(masked_array(x, m2), axis=0), 2.5))
        assert_(average(masked_array(x, m4), axis=0) is masked)
        assert_equal(average(masked_array(x, m5), axis=0), 0.0)
        assert_equal(count(average(masked_array(x, m4), axis=0)), 0)
        z = masked_array(y, m3)
        assert_(allclose(average(z, None), 20. / 6.))
        assert_(allclose(average(z, axis=0),
                                 [0., 1., 99., 99., 4.0, 7.5]))
        assert_(allclose(average(z, axis=1), [2.5, 5.0]))
        assert_(allclose(average(z, axis=0, weights=w2),
                                 [0., 1., 99., 99., 4.0, 10.0]))

        a = arange(6)
        b = arange(6) * 3
        r1, w1 = average([[a, b], [b, a]], axis=1, returned=1)
        assert_equal(shape(r1), shape(w1))
        assert_equal(r1.shape, w1.shape)
        r2, w2 = average(ones((2, 2, 3)), axis=0, weights=[3, 1], returned=1)
        assert_equal(shape(w2), shape(r2))
        r2, w2 = average(ones((2, 2, 3)), returned=1)
        assert_equal(shape(w2), shape(r2))
        r2, w2 = average(ones((2, 2, 3)), weights=ones((2, 2, 3)), returned=1)
        assert_(shape(w2) == shape(r2))
        a2d = array([[1, 2], [0, 4]], float)
        a2dm = masked_array(a2d, [[0, 0], [1, 0]])
        a2da = average(a2d, axis=0)
        assert_(eq(a2da, [0.5, 3.0]))
        a2dma = average(a2dm, axis=0)
        assert_(eq(a2dma, [1.0, 3.0]))
        a2dma = average(a2dm, axis=None)
        assert_(eq(a2dma, 7. / 3.))
        a2dma = average(a2dm, axis=1)
        assert_(eq(a2dma, [1.5, 4.0]))

 def myfunction(d,mx,mn):
   from numpy.ma import maximum,minimum,absolute,greater,count
   try:
     if count(d)==0 : return mx,mn
     mx=float(maximum(mx,float(maximum(d))))
     mn=float(minimum(mn,float(minimum(d))))
   except:
     for i in d:
       mx,mn=myfunction(i,mx,mn)
   return mx,mn

def getStatVal(imageFile,longitude,latitude,winsize,statistic,site):
    """
    Caculates the statistics on the pixels in the window array
    """

    band1,band2,band3,band4,band5,band6,count = 'None','None','None','None','None','None','None'

    if imageFile != 'None' and imageFile != None:
        imageFile=qvf.changestage(imageFile,'tmp')
        temp = '%s_%s_%spix.tif' % (imageFile.split('.')[0],site.strip(),winsize)

        if not os.path.exists(temp):
            subsetRaster = getWindow(imageFile,longitude,latitude,winsize,site)
        else:
            subsetRaster = temp

        try:
            imgInfo = gdalcommon.info(subsetRaster)
            handle = gdal.Open(subsetRaster)

            for band in [1,2,3,4,5,6]:
                if handle != None:
                    bandHandle = handle.GetRasterBand(band)
                    bandArray = bandHandle.ReadAsArray()

                    maskedBand = ma.masked_values(bandArray, 0)
                    count = ma.count(maskedBand)

                    if statistic == 'mean': statVal = maskedBand.mean()
                    elif statistic == 'std': statVal = maskedBand.std()
                    else:
                        statVal = None

                    if band == 1:
                        band1 = statVal

                    elif band == 2:
                        band2 = statVal

                    elif band == 3:
                        band3 = statVal

                    elif band == 4:
                        band4 = statVal

                    elif band == 5:
                        band5 = statVal

                    elif band == 6:
                        band6 = statVal
        except:
            pass

    return band1, band2, band3, band4, band5, band6, count

 def myfunction(d,mx,mn):
   from numpy.ma import maximum,minimum,masked_where,absolute,greater,count
   try:
     d=masked_where(greater(absolute(d),9.9E19),d)
     if count(d)==0 : return mx,mn
     mx=float(maximum(mx,float(maximum(d))))
     mn=float(minimum(mn,float(minimum(d))))
   except:
     for i in d:
       mx,mn=myfunction(i,mx,mn)
   return mx,mn

    def test_reset_clip_box(self) : 
        """test that we can define a different clip box once set"""
        self.gca.set_clip_box(-90,-82,-180,-173)
        testvec = self.gca.compress(self.grid)
        self.assertTrue(self.gca.is_masked())
        testmask = self.gca.get_vec_mask()
        self.assertEqual(np.count_nonzero(testmask), 3)
        self.assertEqual(ma.count(testvec), 1)

        self.gca.set_clip_box(-90,-75,-180,-165)
        testvec = self.gca.compress(self.grid)
        self.assertTrue(self.gca.is_masked())
        self.assertEqual(ma.count(testvec), 2)
        testmask = self.gca.get_vec_mask()
        self.assertEqual(np.count_nonzero(testmask), 2)
        self.assertEqual(ma.count(testvec), 2)

        
        
        
        

 def test_testBasic1d(self):
     # Test of basic array creation and properties in 1 dimension.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     assert_(not isMaskedArray(x))
     assert_(isMaskedArray(xm))
     assert_equal(shape(xm), s)
     assert_equal(xm.shape, s)
     assert_equal(xm.dtype, x.dtype)
     assert_equal(xm.size, reduce(lambda x, y:x * y, s))
     assert_equal(count(xm), len(m1) - reduce(lambda x, y:x + y, m1))
     assert_(eq(xm, xf))
     assert_(eq(filled(xm, 1.e20), xf))
     assert_(eq(x, xm))

def maskImageStats(mimage):
        n=ma.count(mimage)
        mimagesq=mimage*mimage
        sum1=ma.sum(mimage)
        sum2=ma.sum(sum1)
        sumsq1=ma.sum(mimagesq)
        sumsq2=ma.sum(sumsq1)
        avg=sum2/n
        if (n > 1):
                stdev=math.sqrt((sumsq2-sum2*sum2/n)/(n-1))
        else:
                stdev=2e20
        return n,avg,stdev

 def directionality(self):
     # Create series to store data
     self.binDF['self'] = self.binDF['inter'] = self.binDF['up'] = self.binDF['down'] = self.binDF['log2'] = np.nan
     # Loop through rows of the matrix
     for rowNo, row in enumerate(self.probMatrix):
         # Set none values if bin is entirely masked
         if ma.count(row) == 0:
             continue
         # Else calculate values
         else:
             # Extract self frequency
             selflig = row[rowNo]
             # Extract up frequency
             up = row[:rowNo]
             if ma.count(up) == 0:
                 up = 0.
             else:
                 up = up.sum()
             # Extract down frequency
             down = row[rowNo + 1:]
             if ma.count(down) == 0:
                 down = 0.
             else:
                 down = down.sum()
             # Calculate inter value
             inter = 1 - selflig - up - down
             # Calculate log2 value
             if up == 0:
                 if down == 0:
                     log2 = np.nan
                 else:
                     log2 = -np.inf
             elif down == 0:
                 log2 = np.inf
             else:
                 log2 = np.log2(up/down)
             # Store results
             self.binDF.loc[rowNo,['self','inter','up','down','log2']] = (
                 selflig, inter, up, down, log2)

 def test_set_window(self) : 
     window_data = np.ones( (self.lat_size_win, self.lon_size_win) )
     x = self.w.set_window(window_data)
     
     # check output geometry
     self.assertEqual(x.shape[0], 360)
     self.assertEqual(x.shape[1], 720)
     
     # check output is masked
     self.assertTrue(ma.is_masked(x))
     
     # check that the window is only thing in returned array
     win_masked = ma.count_masked(x)
     win = ma.count(x)
     self.assertEqual(win, window_data.size)
     self.assertEqual(win_masked, x.size - window_data.size)
     self.assertTrue(np.all(x[self.w._window] == window_data))