def sigma_clip(image, sigma_lo=3, sigma_hi=3, max_iter=5, axis=0):
    """Reference implementation in numpy"""
    image = image.copy()
    mask = numpy.logical_not(numpy.isfinite(image))
    dummies = mask.sum()
    image[mask] = numpy.NaN
    mean = numpy.nanmean(image, axis=axis, dtype="float64")
    std = numpy.nanstd(image, axis=axis, dtype="float64")
    for _ in range(max_iter):
        if axis == 0:
            mean2d = as_strided(mean, image.shape, (0, mean.strides[0]))
            std2d = as_strided(std, image.shape, (0, std.strides[0]))
        else:
            mean2d = as_strided(mean, image.shape, (mean.strides[0], 0))
            std2d = as_strided(std, image.shape, (std.strides[0], 0))
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            delta = (image - mean2d) / std2d
            mask = numpy.logical_or(delta > sigma_hi,
                                    delta < -sigma_lo)
        dummies = mask.sum()
        if dummies == 0:
            break
        image[mask] = numpy.NaN
        mean = numpy.nanmean(image, axis=axis, dtype="float64")
        std = numpy.nanstd(image, axis=axis, dtype="float64")
    return mean, std

def checkAbundanceScatterClusters():
    # First read the cluster data
    cldata= read_clusterdata.read_caldata()
    # Read the allStar data to match
    # For each of the calibration open clusters, calculate the offset from the 
    # mean in our FEHTAG and AFETAG
    clusters= ['M71','N2158','N2420','N188','M67','N7789','N6819',
               'N6791']
    fehoffset= []
    afeoffset= []
    for cluster in clusters:
        tdata= cldata[cldata['CLUSTER'] == cluster.upper()]
        tdata= tdata[(tdata['TEFF'] < _TEFFMAX)\
                         *(tdata['TEFF'] > _TEFFMIN)\
                         *(tdata['LOGG'] < 3.5)]
        # Compute the average feh and afe and save the offsets
        medianfeh= numpy.median(tdata['FE_H'])
        medianafe= numpy.median(tdata[define_rcsample._AFETAG])
        fehoffset.extend(tdata['FE_H']-medianfeh)
        afeoffset.extend(tdata[define_rcsample._AFETAG]-medianafe)
        if cluster == 'M67': print medianfeh, medianafe, len(tdata)
    fehoffset= numpy.array(fehoffset)
    afeoffset= numpy.array(afeoffset)
    print 'FE_H scatter %g' % (numpy.nanstd(fehoffset[numpy.fabs(fehoffset) < 0.3]))
    print 'A_FE scatter %g' % (numpy.nanstd(afeoffset[numpy.fabs(afeoffset) < 0.3]))
    gindx= (numpy.fabs(fehoffset) < 0.3)*(numpy.fabs(afeoffset) < 0.3)
    print 'FE_H/A_FE correlation %g' % (numpy.mean(afeoffset[gindx]*fehoffset[gindx])/numpy.nanstd(fehoffset[numpy.fabs(fehoffset) < 0.3])/numpy.nanstd(afeoffset[numpy.fabs(afeoffset) < 0.3]))
    print 'FE_H robust scatter %g' % (1.4826*numpy.median(numpy.fabs(fehoffset)))
    print 'A_FE robust scatter %g' % (1.4826*numpy.median(numpy.fabs(afeoffset)))
    bovy_plot.bovy_print()
    bovy_plot.bovy_hist(fehoffset,range=[-0.3,0.3],bins=31,histtype='step')
    bovy_plot.bovy_hist(afeoffset,range=[-0.3,0.3],bins=31,histtype='step',
                        overplot=True)
    bovy_plot.bovy_end_print('test.png')
    return None

def ddiff_var_wb(lh_ctx_file, rh_ctx_file, sc_file):
    import numpy as np
    from generic_pipelines.moco_eval import wb_to_tss
    tss = wb_to_tss(lh_ctx_file, rh_ctx_file, sc_file)
    out_fname = 'ddiff.npz'
    np.savez(out_fname,diffvar=np.nanstd(tss,0), ddiffvar=np.nanstd(np.diff(tss,1,0),0))
    return out_fname

def make_NRCS_image( nobj, bandname, fn='', dir='.', max=np.nan, min=np.nan,
        **kwargs):
    if not fn:
        if 'reduced' in bandname:
            fn = bandname[:9]+'.png'
        else:
            fn = bandname+'.png'
    resize(nobj)
    try:
        s0 = 10.0*np.log10(nobj[bandname])
    except:
        n_obj.undo()
        raise
    s0[np.where(np.isinf(s0))]=np.nan
    #if nobj.fileName[-2:]=='nc':
    #    s0 = flipdim(nobj,s0)

    caption='dB'
    if np.isnan(min):
        min = np.nanmedian(s0,axis=None)-2.0*np.nanstd(s0,axis=None)
    if np.isnan(max):
        max = np.nanmedian(s0,axis=None)+2.0*np.nanstd(s0,axis=None)
    nansatFigure(s0, min, max, dir, fn)
    nobj.undo()
    return fn

def compare_iter(arr_len, n_iter):
    """
    Use bubble, quick and merge sorts of random arrays of a set length,
    for n_iter times. Then return mean and standard deviations of results
    The arrays are limited to values less than 1000.
    """
    bubble_comps = []
    quick_comps  = []
    merge_comps  = []

    # Perform sorting the required number of times:
    for ind in range(n_iter):
        rand_arr = np.random.randint(1000, size = arr_len)
        bubble_comps.append(bubble_sort(rand_arr, 0))
        quick_comps.append(quick_sort(rand_arr, 0))
        merge_comps.append(merge_sort(rand_arr, 0))

    # Extract the number of comparisons:
    bub_no = np.array([x[0] for x in bubble_comps])
    qck_no = np.array([x[0] for x in quick_comps])
    mrg_no = np.array([x[0] for x in merge_comps])

    # Calculate mean and standard deviations:
    bub_mean   = np.nanmean(bub_no)
    qck_mean   = np.nanmean(qck_no)
    mrg_mean   = np.nanmean(mrg_no)
    bub_stddev = np.nanstd(bub_no)
    qck_stddev = np.nanstd(qck_no)
    mrg_stddev = np.nanstd(mrg_no)

    # Return the means and standard deviations
    return bub_mean, bub_stddev, qck_mean, qck_stddev, mrg_mean, mrg_stddev

def compute_righthind(forecast, hindcast, observation):

    """
    Esta função faz a correção do hincast usando
    a regressão linear (correção de bias). Mesmo
    método utilizado pelo Júnior.
    """

    # Calcula média e desvio padrão para todos os pontos da climatologia do modelo
    clim_mean = np.nanmean(hindcast, axis=0)
    clim_std = np.nanstd(hindcast, axis=0)

    # Calcula média e desvio padrão para todos os pontos da climatologia do
    # dados observado
    obs_mean = np.nanmean(observation, axis=0)
    obs_std = np.nanstd(observation, axis=0)

    # Calcula variável padronizada do modelo e observado
    clim_pad = (hindcast - clim_mean)/clim_std
    # obs_pad = (observation - obs_mean)/obs_std
    fcst_pad = (forecast - clim_mean)/clim_std

    # newhind é o hindcast corrigido
    newhind = clim_pad * obs_std + obs_mean

    newfcst = fcst_pad * obs_std + obs_mean

    return newhind, newfcst

def attenSig( comp , min_var_warn=0 , min_var_fail=0 , sd = False ):
    # comp: list of values to run the test on. Therefore, time interval determined by user
    # min_var_warn: minimum range of data variation to trip a warning
    # min_var_fail: minimum range of data variation to trip a fail
    # sd: if true, comp data stdev will be tested against stdev thresholds set by user
    if min_var_warn==0 and min_var_fail==0:
        raise ValueError ("Please indicate min&max deviance threshold values.")
    n=len(comp)

    if sd==True:
        if np.nanstd(comp) <= min_var_warn:
            flag = 3
        elif np.nanstd(comp) <= min_var_fail:
            flag = 4
        else:
            flag = 1


    if sd==False:
        test = abs(max(comp) - min(comp))
        if test <= min_var_warn:
            flag = 3
        elif test <= min_var_fail:
            flag = 4
        else:
            flag = 1

    flag_list = [flag*(x/x) for x in np.arange(1,(len(comp)+1))]

    return flag_list

 def predict(self, data=None):
     if(self.use_period):
         # decomfreq = freq
         res = sm.tsa.seasonal_decompose(self.data.tolist(), freq=self.freq, model=self.model)
         #     res.plot()
         median_trend = pd.rolling_median(Series(self.data),window=self.freq, center=True, min_periods=1)
         resid = res.observed - res.seasonal - median_trend
     else:
         resid = self.data
     random = Series(resid)
     mean_nan = 0
     std_nan = 0
     # random = res.resid
     if (self.mode == 'average'):
         mean_nan = np.nanmean(random)
         std_nan = np.nanstd(random)
     elif (self.mode == 'median'):
         rolling_median = pd.rolling_median(random,3,center=True, min_periods=1)
         mean_nan = np.nanmean(rolling_median)
         std_nan = np.nanstd(rolling_median)
     min_val = mean_nan - 4 * std_nan
     # max_val = mean(random, na.rm = T) + 4*sd(random, na.rm = T)
     max_val = mean_nan + 4 * std_nan
     position = Series(resid.tolist(), index=np.arange(resid.shape[0]))
     anomaly = position[(position > max_val) | (position < min_val)]
     # anomalyL = position[(position<min_val)]
     # anomaly = anomalyH.append(anomalyL).drop_duplicates()
     point_anomaly_idx = anomaly.index
     self.anomaly_idx = point_anomaly_idx
     points_anomaly = self.data[point_anomaly_idx]
     self.anomalies = points_anomaly
     return points_anomaly

def main():
    os.system('modprobe w1-gpio')
    os.system('modprobe w1-therm')
    print len(sys.argv)
    if len(sys.argv) == 1:
        number_of_meas = 7
    else:
        print sys.argv[1]
        number_of_meas = int(sys.argv[1])
    print "number_of_measurements = " + str(number_of_meas)
    
    print "getting device files and serials..."
    THEDICT = _get_w1_tree_and_serials()
    
    print "reading sensors " + str(number_of_meas) + " times ..."
    for step in range(int(number_of_meas)):
        for sensor_id in THEDICT:
            if sensor_id[0:2] == '28' or sensor_id[0:2] == '10':
                temp = read_sensor_ds18b20(sensor_id,THEDICT[sensor_id]["path"])
                volt = "n.a."
                THEDICT[sensor_id]["temp"].append(temp)
                THEDICT[sensor_id]["volt"].append(0.)
            if sensor_id[0:2] == '26':
                temp,volt = read_sensor_ds2438(sensor_id,THEDICT[sensor_id]["path"])
                THEDICT[sensor_id]["temp"].append(temp)
                THEDICT[sensor_id]["volt"].append(volt)
            print "step " + str(step) + " " + sensor_id + " " + str(temp) + " " + str(volt)
    
    print "calculating individual and total means:"
    MEAN_IND = {}
    for sensor_id in THEDICT:
        MEAN_IND[sensor_id] = [
                                np.nanmean(np.array(THEDICT[sensor_id]["temp"])), 
                                np.nanmean(np.array(THEDICT[sensor_id]["volt"]))
                              ]
    total_temp = []
    total_volt = []
    for sensor_id in MEAN_IND:
        if sensor_id[0:2] == '28' or sensor_id[0:2] == '10':
            total_temp.append(MEAN_IND[sensor_id][0])
        if sensor_id[0:2] == '26':
            total_volt.append(MEAN_IND[sensor_id][1])
    mean_temp = np.nanmean(np.array(total_temp))
    mean_volt = np.nanmean(np.array(total_volt))
    
    print "temp mean: " + str(mean_temp) + " +/- " + str(np.nanstd(np.array(total_temp)))
    print "volt mean: " + str(mean_volt) + " +/- " + str(np.nanstd(np.array(total_temp)))
        
    
    print "calculating offsets..."
    OFFSETS = {}
    for sensor_id in MEAN_IND:
        OFFSETS[sensor_id] = [
                               MEAN_IND[sensor_id][0] - mean_temp, 
                               MEAN_IND[sensor_id][1] - mean_volt
                             ]
    print OFFSETS
            
    print "writing offsets..."
    write_offset(OFFSETS)

    def plotQE(self, save=False):
        '''
        Plot the measured and theoretical QE
        '''
        fig, ax = plt.subplots()

        QE_median = np.array([np.nanmedian(QE.flatten()) for QE in self.QE])
        QE_upper = np.array([QE_median[ind] + np.nanstd(QE.flatten())
                             for ind, QE in enumerate(self.QE)])
        QE_lower = np.array([QE_median[ind] - np.nanstd(QE.flatten())
                             for ind, QE in enumerate(self.QE)])

        ax.plot(self.wavelengths, QE_median * 100, linewidth=3, color='black',
                label=r'Measured')
        ax.fill_between(self.wavelengths, QE_lower * 100, QE_upper * 100,
                        where=QE_upper >= QE_lower, color='green', facecolor='green',
                        interpolate='True', alpha=0.1)
        ax.plot(self.wvl_theory, 100 * self.QE_theory, linestyle='-.', linewidth=2,
                color='black', label=r'Theoretical')

        ax.set_xlabel(r'Wavelength (nm)')
        ax.set_ylabel(r'QE (%)')
        ax.legend()
        ax.set_xlim([min(self.wavelengths), max(self.wavelengths)])
        ax.set_ylim([0, max(QE_upper) * 100 * 1.2])

        if save:
            file_name = os.path.join(self.out_directory, self.log_file + '.pdf')
            plt.savefig(file_name, format='pdf')
            plt.close(fig)
        else:
            plt.show(block=False)

def get_night_shifts(offset):
    """Returns the mean APASS-based shifts for all the offsets in the same night."""
    shifts = {}
    shifts_std = {}
    for band in ['r', 'i']:
        offset_shifts = []
        for sibling in red_groups[offset]:
            try:
                offset_shifts.append(SHIFTS[sibling][band + 'shift'])
            except KeyError:
                pass
        shifts[band] = np.nanmean(offset_shifts)
        shifts_std[band] = np.nanstd(offset_shifts)
    if offset not in blue_groups:
        for band in ['u', 'g', 'r2']:
            shifts[band] = np.nan
    else:
        for band in ['u', 'g', 'r2']:
            offset_shifts = []
            for sibling in blue_groups[offset]:
                try:
                    offset_shifts.append(SHIFTS[sibling][band + 'shift'])
                except KeyError:
                    pass
            shifts[band] = np.nanmean(offset_shifts)
            shifts_std[band] = np.nanstd(offset_shifts)
    return shifts, shifts_std

def get_variance_map2(a_plus_b, a_minus_b, bias_mask2, pix_mask, gain):
    #variance0 = a_minus_b
    #a_minus_b = a-b
    msk = bias_mask2 | pix_mask | ~np.isfinite(a_minus_b)

    from destriper import destriper
    variance0 = destriper.get_destriped(a_minus_b,
                                        msk,
                                        pattern=64,
                                        remove_vertical=False,
                                        hori=False)
    #variance0 = a_minus_b

    # stsci_median cannot be used due to too many array error.
    #ss = stsci_median([m1 for m1 in variance0],)
    dd1 = np.ma.array(variance0, mask=msk)
    ss = np.ma.median(dd1, axis=0)

    variance_ = variance0.copy()
    variance_[msk] = np.nan

    st = np.nanstd(variance_)
    st = np.nanstd(variance_[np.abs(variance_) < 3*st])

    variance_[np.abs(variance_-ss) > 3*st] = np.nan

    import scipy.ndimage as ni
    x_std = ni.median_filter(np.nanstd(variance_, axis=0), 11)

    variance_map0 = np.zeros_like(variance_) + x_std**2

    variance_map = variance_map0 + np.abs(a_plus_b)/gain # add poison noise in ADU
    return variance_map

def bin_fit(x, y, buckets=3):
     
    assert buckets in [3,25]

    xstd=np.nanstd(x)
    
    if buckets==3:
        binlimits=[np.nanmin(x), -xstd/2.0,xstd/2.0 , np.nanmax(x)]
    elif buckets==25:
    
        steps=xstd/4.0
        binlimits=np.arange(-xstd*3.0, xstd*3.0, steps)
    
        binlimits=[np.nanmin(x)]+list(binlimits)+[np.nanmax(x)]
    
    fit_y=[]
    err_y=[]
    x_values_to_plot=[]
    for binidx in range(len(binlimits))[1:]:
        lower_bin_x=binlimits[binidx-1]
        upper_bin_x=binlimits[binidx]

        x_values_to_plot.append(np.mean([lower_bin_x, upper_bin_x]))

        y_in_bin=[y[idx] for idx in range(len(y)) if x[idx]>=lower_bin_x and x[idx]<upper_bin_x]

        fit_y.append(np.nanmedian(y_in_bin))
        err_y.append(np.nanstd(y_in_bin))

    ## no zeros
    

    return (binlimits, x_values_to_plot, fit_y, err_y)

def radial_contrast_flr(image, xc, yc, seps, zw, coron_thrupt, klip_thrupt=None):
    rad_flr_ctc = np.empty((len(seps)))
    assert(len(seps) == len(coron_thrupt))
    if klip_thrupt is not None:
        assert(len(seps) == len(klip_thrupt))
        rad_flr_ctc_ktc = np.empty((len(seps)))
    else:
        rad_flr_ctc_ktc = None

    imh = image.shape[0]
    imw = image.shape[1]

    xs = np.arange(imw) - xc
    ys = np.arange(imh) - yc
    XXs, YYs = np.meshgrid(xs, ys)
    RRs = np.sqrt(XXs**2 + YYs**2)

    for si, sep in enumerate(seps):
        r_in = np.max([seps[0], sep-zw/2.])
        r_out = np.min([seps[-1], sep+zw/2.])
        meas_ann_mask = np.logical_and(np.greater_equal(RRs, r_in),
                                          np.less_equal(RRs, r_out))
        meas_ann_ind = np.nonzero(np.logical_and(np.greater_equal(RRs, r_in).ravel(),
                                                    np.less_equal(RRs, r_out).ravel()))[0]
        meas_ann = np.ravel(image)[meas_ann_ind]
        rad_flr_ctc[si] = np.nanstd(meas_ann)/coron_thrupt[si]
        if rad_flr_ctc_ktc is not None:
            rad_flr_ctc_ktc[si] = np.nanstd(meas_ann)/coron_thrupt[si]/klip_thrupt[si]

    #pdb.set_trace()
    return rad_flr_ctc, rad_flr_ctc_ktc

 def _get_x_0_stats(self):
     x_diff = np.diff(self.x_arr_0, axis=1)
     mu_mm = np.nanmean(x_diff)
     std_mm = np.nanstd(x_diff)
     mu_px_mm = np.nanmean(x_diff / self.aramis_info.n_px_facet_step_x)
     std_px_mm = np.nanstd(x_diff / self.aramis_info.n_px_facet_step_x)
     return mu_mm, std_mm, mu_px_mm, std_px_mm

def main():
	fname = 'cleaned_data/core_male_cleaned.csv'
	# print(filelength(fname))
	# fname = 'cleaned_data/core_patients_cleaned.csv'
	# code1 = URETHRAL_INJURY_CODES[0]
	PENILE_FRACTURE_CODE = '95913'
	if not sys.argv[1]:
		print('Please enter an argument')
		sys.exit(0)
	code1 = str(sys.argv[1])
	code2 = '95913'
	# try:
	# 	data_mat = np.loadtxt('cleaned_data/pxpywt.txt')
	# 	print('Done loading file! Starting analysis.')
	# except:
	print('cannot load data, going to try generating...')
	data_mat = binary_arrays(fname,code1,code2,1)
	true_stat = mi(data_mat)
	print(true_stat)
	# print(leaders(DXLIST))
	print('beginning surrogate_stats')
	surrogate_stats = surrogate_mi(data_mat)
	print('beginning bootstrap_stats')
	bootstrap_stats = bootstrap_mi(data_mat)
	np.savetxt('cleaned_data/{0}surrogate_stats.txt'.format(str(sys.argv[2]),),surrogate_stats,fmt='%f')
	np.savetxt('cleaned_data/{0}bootstrap_stats.txt'.format(str(sys.argv[2]),),bootstrap_stats,fmt='%f')
	# plt.hist(bootstrap_stats,50)
	# plt.show()
	print(wald_test(true_stat, np.nanstd(bootstrap_stats), np.nanmean(surrogate_stats), np.nanstd(surrogate_stats)))

def pickVMIN(dataIn, devs):
    try:
        if np.nanmean(dataIn) - devs * np.nanstd(dataIn) < np.nanmin(dataIn):
            return np.nanmin(dataIn)
        else:
            return np.nanmean(dataIn) - devs * np.nanstd(dataIn)
    except ValueError:
        return 1

 def estimate_f(self, i, samples = 200):
     
     #minibatch = np.random.choice(X,batch_size,replace=False)
     
     LL_X = np.array([self.f['LL'](i) for s in range(samples if samples != None else self.samples)])
     LL_Y = np.array([self.f['LLY'](i) for s in range(samples if samples != None else self.samples)])
     
     return np.nanmean(LL_X, 0)+np.nanmean(LL_Y, 0), np.nanstd(LL_X, 0)+np.nanstd(LL_Y, 0)

def pickVMAX(dataIn, devs):
    try:
        if np.nanmean(dataIn) + devs * np.nanstd(dataIn) > np.nanmax(dataIn):
            return np.nanmax(dataIn)
        else:
            return np.nanmean(dataIn) + devs * np.nanstd(dataIn)
    except ValueError:
        return 1

def ddiff_var_moco(in_file):
    import numpy as np
    import h5py
    ts = h5py.File(in_file,'r')
    tss = np.asarray(ts['FMRI/DATA']).T
    out_fname = 'ddiff.npz'
    np.savez(out_fname,diffvar=np.nanstd(tss,0), ddiffvar=np.nanstd(np.diff(tss,1,0),0))
    return out_fname