def plotterIterations(X,nX,nIt,y,ny) :
    #On se ramène à des journées
    EltIteration = X[:,nIt]
    Iterations = np.unique(EltIteration)
    print(Iterations)

    #Gérer plusieurs couleurs
    colors = cm.rainbow(np.linspace(0, 1, len(Iterations)))
    
    i = 0
    for k,c in zip(Iterations,colors)  :
        #On prend que le jour
        condition = EltIteration == k #entre 60 et 153
        xPlot = np.extract(condition, X[:,nX])
        yPlot = np.extract(condition, y[:,ny])
        plt.scatter(xPlot,yPlot, color=c,s=2)
        if i==10 :
            break
        i+=1
        #sort
        #print(np.concatenate(xPlot,yPlot))
        #z = np.sort(np.concatenate(xPlot,yPlot),0)
        #plt.plot(z[:,0], z[:,1], color=c,linewidth=2)
        #print(k)
        #print(xPlot)
        #print(yPlot)
        #print(z[:,0])
        #print(z[:,1])
        #break

    #plt.plot(X_test, regr.predict(X_test), color='blue',linewidth=3)
    plt.xticks(())
    plt.yticks(())

    plt.title('')

def reduce_grid_points(mesh_divisors,
                       grid_address,
                       dense_grid_points,
                       dense_grid_weights=None,
                       coarse_mesh_shifts=None):
    divisors = np.array(mesh_divisors, dtype='intc')
    if (divisors == 1).all():
        coarse_grid_points = np.array(dense_grid_points, dtype='intc')
        if dense_grid_weights is not None:
            coarse_grid_weights = np.array(dense_grid_weights, dtype='intc')
    else:
        grid_weights = []
        if coarse_mesh_shifts is None:
            shift = [0, 0, 0]
        else:
            shift = np.where(coarse_mesh_shifts, divisors // 2, [0, 0, 0])
        modulo = grid_address[dense_grid_points] % divisors
        condition = (modulo == shift).all(axis=1)
        coarse_grid_points = np.extract(condition, dense_grid_points)
        if dense_grid_weights is not None:
            coarse_grid_weights = np.extract(condition, dense_grid_weights)

    if dense_grid_weights is None:
        return coarse_grid_points
    else:
        return coarse_grid_points, coarse_grid_weights

def parzen(M, sym=True):
    """Return a Parzen window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    sym : bool, optional
        When True, generates a symmetric window, for use in filter design.
        When False, generates a periodic window, for use in spectral analysis.

    Returns
    -------
    w : ndarray
        The window, with the maximum value normalized to 1 (though the value 1
        does not appear if the number of samples is even and sym is True).

    Examples
    --------
    Plot the window and its frequency response:

    >>> from scipy import signal
    >>> from scipy.fftpack import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> window = signal.parzen(51)
    >>> plt.plot(window)
    >>> plt.title("Parzen window")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 2048) / (len(window)/2.0)
    >>> freq = np.linspace(-0.5, 0.5, len(A))
    >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    >>> plt.plot(freq, response)
    >>> plt.axis([-0.5, 0.5, -120, 0])
    >>> plt.title("Frequency response of the Parzen window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = np.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0)
    na = np.extract(n < -(M - 1) / 4.0, n)
    nb = np.extract(abs(n) <= (M - 1) / 4.0, n)
    wa = 2 * (1 - np.abs(na) / (M / 2.0)) ** 3.0
    wb = (1 - 6 * (np.abs(nb) / (M / 2.0)) ** 2.0 +
          6 * (np.abs(nb) / (M / 2.0)) ** 3.0)
    w = np.r_[wa, wb, wa[::-1]]
    if not sym and not odd:
        w = w[:-1]
    return w

def factorial2(n,exact=0):
    """n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi)  n odd
           = 2**(n) * n!                                 n even

    If exact==0, then floating point precision is used, otherwise
    exact long integer is computed.

    Notes:
      - Array argument accepted only for exact=0 case.
      - If n<0, the return value is 0.
    """
    if exact:
        if n < -1:
            return 0L
        if n <= 0:
            return 1L
        val = 1L
        for k in xrange(n,0,-2):
            val *= k
        return val
    else:
        from scipy import special
        n = asarray(n)
        vals = zeros(n.shape,'d')
        cond1 = (n % 2) & (n >= -1)
        cond2 = (1-(n % 2)) & (n >= -1)
        oddn = extract(cond1,n)
        evenn = extract(cond2,n)
        nd2o = oddn / 2.0
        nd2e = evenn / 2.0
        place(vals,cond1,special.gamma(nd2o+1)/sqrt(pi)*pow(2.0,nd2o+0.5))
        place(vals,cond2,special.gamma(nd2e+1) * pow(2.0,nd2e))
        return vals

def get_linear_idxs_of_diff_bits(SortedDeltaIdxs, dists_vec, cfg):

    # Identifies spots in matrix where the indices of diff bits will be.
    IsPossDistLteActual = calc_IsPossDistLteActual(dists_vec, cfg)

    # Columnar linear indices for the whole mtx.
    UniqSortedDeltaIdxs = uniq_SortedDeltaIdxs(SortedDeltaIdxs, cfg)
    
    # >> Which bits were actually different!  (Pts were on opposite sides of HP.) <<
    #
    # Hack:
    # We +1 to the idxs, so that none has the value 0.
    # Then we replace all the idxs we don't want with 0.
    # For those left, we shift back down one.
    LinearIdxsOfDiffBits = (UniqSortedDeltaIdxs + 1) * IsPossDistLteActual
    my_assert (lambda : LinearIdxsOfDiffBits.shape == cfg.shape)
    # print 'LinearIdxsOfDiffBits (+1)';  print LinearIdxsOfDiffBits

    diff_bits = np.extract(LinearIdxsOfDiffBits > 0, LinearIdxsOfDiffBits)
    diff_bits -= 1
    # --------------------------------
    # extract(): returns non-zero vals by going along ROWs,
    # not down COLs (as in Matlab).  So the ORDER of these may be wacky.
    # Sorting gives us the same order that Matlab would.
    # --------------------------------
    diff_bits.sort()

    my_assert (lambda :
               len(diff_bits) ==
               len(np.extract(IsPossDistLteActual == True, IsPossDistLteActual)))
    
    print 'diff_bits:', diff_bits
    return diff_bits

def getClusterXYZ(RecHits, clusterID):
    '''
    Computes the log-energy-weighted xyzt coordinates of a given cluster. The energy is initially
    reduced by a threshold amount to discard noisy hits, and then log'd to account for energy
    collection fluctuations in the material.
    '''
    extractMask = np.logical_and.reduce((RecHits['t'] > 0, 
                                         RecHits['clusterID'] == clusterID,
                                         RecHits['isIn3x3']))                                       #|Extraction mask for cluster t and xyz processing; decides initially which hits to be counted
    #minLayer = np.min(newRecHits['layerID'])
    #newRecHits = np.compress(newRecHits['layerID'] == minLayer, newRecHits)
    x = np.extract(extractMask, RecHits['x'])
    y = np.extract(extractMask, RecHits['y'])
    z = np.extract(extractMask, RecHits['z'])
    E = np.extract(extractMask, RecHits['en'])

    return (np.mean(x), np.mean(y), np.mean(z))                                                     #|Flat average

    # Energy weight
    w0 = 7                                                                                          #|Arbitrary weighting term that works well, as found by Geoffrey Monet
    Eweight = np.maximum(np.log(E/np.sum(E)) + w0, 0)


    xw, yw, zw = (np.dot(x,Eweight)/np.sum(Eweight), 
                  np.dot(y,Eweight)/np.sum(Eweight), 
                  np.dot(z,Eweight)/np.sum(Eweight))                                                #|Weighted average

    return xw, yw, zw

def find_next_set_images(location_x,location_y,heading,file_database_sorted,picture_name_list):
    
    image_found=0
    
    heading,direction_vector=phase_wrap_heading(heading)
    # Convert heading
    phase_wrap=np.array([3, 0, 1, 2, 3, 0],dtype='u1')
    heading_array=np.array([phase_wrap[heading], phase_wrap[heading+1],phase_wrap[heading+2]])
    # find x values
    matched_x_loc=np.extract(file_database_sorted['x_loc']==location_x,file_database_sorted)
    # Check values found!!!!!
    if matched_x_loc.size<4:
        print "Not enough images at this x location!!, number img=\t", matched_x_loc.size 
        return (0,0,heading,direction_vector,0)

    # find y values
    matched_y_loc=np.extract(matched_x_loc['y_loc']==location_y,matched_x_loc)
    # Check values found!!!!!
    if matched_y_loc.size<4:
        print "Not enough images at this y location!!, number img=\t", matched_y_loc.size 
        return (0,0,heading,direction_vector,0)
 
    images_2_display=matched_y_loc['file_id'][heading_array]
    combined_img = np.concatenate(img_file[images_2_display] , axis=1) #file_database_north_sortx[0,:]
    resized_img = cv2.resize(combined_img, (image_display_width, image_display_height)) 
    image_found=1
    picture_name=picture_name_list[images_2_display[1]]
    return (resized_img,image_found,heading,direction_vector,picture_name)

def boxoverlap(regions_a, region_b, thre):
    # (x1,y1) top-left coord, (x2,y2) bottom-right coord, (w,h) size
    TP=NP=0;
    TP_all=NP_all=0
    N=len(region_b);
    
    for (xb,yb,wb,hb) in region_b:
        x1=np.maximum(regions_a[:,0],xb);
        y1=np.maximum(regions_a[:,1],yb);
        x2=np.minimum((regions_a[:,2]+regions_a[:,0]),(xb+wb));
        y2=np.minimum((regions_a[:,3]+regions_a[:,1]),(yb+hb));
        print x1,y1,x2,y2
        w=x2-x1+1;
        h=y2-y1+1;
        inter=w*h;
        aarea=(regions_a[:,2]+1)*(regions_a[:,3]+1);
        barea=(wb+1)*(hb+1);

        #intersection over union overlap
        o=inter/(aarea+float(barea)-inter);
        
        #set invalid entries to 0 overlap
        o[w<=0]=0
        o[h<=0]=0
        TP=len(np.extract(o>=thre, o))
        NP=len(np.extract(o<thre, o))
        TP_all=TP_all+TP
        
    NP_all=NP-TP_all
    if NP_all<0:
        NP_all=0
        
    return TP_all, NP_all, N; 

def auroc(y_prob, y_true):
    #threshold = [0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 0.99, 0.9999, 0.99999, 0.999999, 1.01]
    threshold = np.sort(np.unique(y_prob))
    threshold[0] = 0
    threshold[threshold.shape[0]-1] = 1
    tpr = np.empty_like(threshold)
    fpr = np.empty_like(threshold)
    for i in range(0, len(threshold)):
        #print threshold[i]
        predicted = np.empty_like(y_prob)
        for j in range(0,y_prob.shape[0]):
            if threshold[len(threshold)-1-i] != 1:
                predicted[j] = int(y_prob[j] >= threshold[len(threshold)-1-i])
            else:
                predicted[j] = 0
        a = np.extract(y_true == 1, predicted-y_true)
        tpr[i] = float(a.shape[0] - np.count_nonzero(a))/np.count_nonzero(y_true)
        b = np.extract(y_true == 0, predicted-y_true)
        fpr[i] = np.count_nonzero(b)/float(y_true.shape[0]-np.count_nonzero(y_true))
    #roc = interp1d(fpr, tpr, kind='linear')
    roc_auc = trapz(tpr, fpr)
    #print 'tpr', tpr
    #print 'fpr', fpr
    #print 'roc_auc', roc_auc
    return tpr, fpr, roc_auc

	def computeMcNemarSignificance(self, truth, predictions1, predictions2):
		condition = (truth == 1)
		truth = numpy.extract(condition, truth)
		predictions1 = numpy.extract(condition, predictions1)
		predictions2 = numpy.extract(condition, predictions2)
	
		evals1 = (predictions1 == truth)
		evals2 = (predictions2 == truth)
		
		# Misclassified by the first model only: c01.
		# Misclassified by the second model only: c10.
		c01, c10 = 0, 0
	
		for i, eval1 in enumerate(evals1):
			eval2 = evals2[i]
			if eval1 == 0 and eval2 == 1:
				c01 += 1
			if eval1 == 1 and eval2 == 0:
				c10 += 1
		
		if c01 + c10 < 20:
			print "Unreliable conclusion:", c01, c10
			return 0.0
		else:
			return math.pow(abs(c01 - c10) - 1, 2) / (c01 + c10)

    def running_ave(self, rad_pix, good_image, edges):
        rad_out = []
        prof_out = []
        proferr_out = []
        aperflux = []
        included_pix = []

        for curr_edge in np.arange(len(edges)-1):
            #print 'edges ', edges[curr_edge],edges[curr_edge+1]
            tmp_im = np.extract(rad_pix < edges[curr_edge+1], good_image)
            tmp_rad = np.extract(rad_pix < edges[curr_edge+1], rad_pix)

            #print 'rads1'
            #print tmp_rad[:10]
            #raw_input()

            if len(tmp_rad)>0:
                aper_tmp = np.sum(tmp_im)
                inc_tmp = float(tmp_im.size)
            
                tmp_im = np.extract(tmp_rad >= edges[curr_edge], tmp_im)
                tmp_rad = np.extract(tmp_rad >= edges[curr_edge], tmp_rad)
                
                #print 'rads2'
                #print tmp_rad[:10]
                #raw_input()

                if len(tmp_rad) > 0:
                    rad_out.append(np.mean(tmp_rad))
                    prof_out.append(np.mean(tmp_im))
                    proferr_out.append(np.std(tmp_im))#/(tmp_im.size-1))
                    aperflux.append(aper_tmp)
                    included_pix.append(inc_tmp)

        return rad_out, prof_out,proferr_out, aperflux, included_pix

def gap_rain(qlwell, channel_num=0, threshold=None, pct_boundary=0.3, gap_size=10000):
    """
    Return the rain in the gaps between non-rain droplets.
    """
    rain, nonrain = rain_split(qlwell,
                               channel_num=channel_num,
                               threshold=threshold,
                               pct_boundary=pct_boundary)
    
    # ok, now identify the gaps in the gates.
    times = peak_times(nonrain)
    if nonrain is None or len(nonrain) < 2:
        return np.ndarray([0],dtype=peak_dtype(2))
    
    intervals = np.ediff1d(times, to_begin=0, to_end=0)
    big_intervals = intervals > gap_size

    # find beginning of gaps with extract
    beginnings = np.extract(big_intervals[1:], times)
    ends = np.extract(big_intervals[:-1], times)

    gap_intervals = zip(beginnings, ends)
    gap_intervals.insert(0, (0, times[0]))
    gap_intervals.append((times[-1], times[-1]*100))
    
    # count the rain in the intervals
    gap_drops = np.extract(reduce(np.logical_or, [np.logical_and(peak_times(rain) > b,
                                                                 peak_times(rain) < e) for b, e in gap_intervals]),
                           rain)
    return gap_drops

def rain_split(qlwell, channel_num=0, threshold=None, pct_boundary=0.3, split_all_peaks=False):
    """
    Splits between rain and non-rain.  If you want the well's auto threshold to be used,
    use None as a threshold parameter (the default).
    If you do not want a threshold to be calculated, use '0'. (little unclear from the spec)

    Returns tuple (rain, non-rain)
    """
    if threshold is None:
        threshold = qlwell.channels[channel_num].statistics.threshold
    
    ok_peaks = accepted_peaks(qlwell)
    prain, rain, nrain, p_thresh, mh_thresh, ml_thresh, l_thresh = \
        rain_pvalues_thresholds(ok_peaks, channel_num=channel_num, threshold=threshold, pct_boundary=pct_boundary)

    if split_all_peaks:
        peaks = qlwell.peaks
    else:
        peaks = ok_peaks
    # this would be useful as a standalone, but for efficiency's sake will cut out for now        
    rain_condition_arr = [channel_amplitudes(peaks, channel_num) > p_thresh]
    if mh_thresh and ml_thresh:
        rain_condition_arr.append(np.logical_and(channel_amplitudes(peaks, channel_num) > ml_thresh,
                                              channel_amplitudes(peaks, channel_num) < mh_thresh))
    rain_condition_arr.append(channel_amplitudes(peaks, channel_num) < l_thresh)
    rain_condition = reduce(np.logical_or, rain_condition_arr)
    nonrain_condition = np.logical_not(rain_condition)

    rain = np.extract(rain_condition, peaks)
    nonrain = np.extract(nonrain_condition, peaks)
    return rain, nonrain

def explore_transformation(dirty_dir, clean_dir, num_stds=0):
	images = load_images(dirty_dir)
	cleaned_images = load_images(clean_dir)
	for key in images.keys():
		image = images[key]
		image_c = third_pass_filter(image)
		cv2.imshow('original', image)
		cv2.moveWindow('original', 0, 0)
		cv2.imshow('cleaned by me', image_c)
		cv2.moveWindow('cleaned by me', 500, 0)
		them = cleaned_images[key]
		cv2.imshow('clean', them)
		cv2.moveWindow('clean', 500, 300)

		consider_white = 200
		important_me = np.extract(image_c.flatten() < consider_white, image_c.flatten())
		important_them = np.extract(them.flatten() < consider_white, them.flatten())

		bins = 20
		plt.hist(important_them, bins, label='them')
		plt.hist(important_me, bins, label='me')
		plt.legend(loc='upper right')
		#plt.show()

		cv2.waitKey(0)
		plt.close()
		cv2.destroyAllWindows()

def sub_arr(array, lim, con_array = None, min=None, max=None, boundaries=True):
    """Purpose: Extract sub array of values between min and max limits
    arguements:
     array          var     array to take subset of
     lim            var     array containing [min, max]
     boundaries     bool    include boundaries ie <= and >=
    keywords:
     con_array      var     condition array to apply min/max check on
    Outputs:
     array of values in array with indices satisfying min < con_array < max
    Call example: 
     function()

    TODO: implement only using max or min
    """
    array = check_array(array) # check array is a numpy array

    assert lim[1] >= lim[0], 'min > max'

    if con_array == None: # If no separate array supplied use same array for min/max
        con_array = array
    else: 
        assert np.size(con_array) != np.size(array), 'WARNING: size(con_array) != size(array)'

    if boundaries == True:
        sub = np.extract( (con_array>=lim[0]) * (con_array<=lim[1]), array)
    else:
        sub = np.extract( (con_array>lim[0]) * (con_array<lim[1]), array)
    return sub

def normalize_player_name(img_name, debug=False):
    img_name_w = matcher.MM_WHITE(sat=(0, 96), visibility=(48, 255))(img_name)

    img_name_x_hist = np.extract(np.sum(img_name_w, axis=0) > 128, np.arange(img_name_w.shape[1]))

    img_name_y_hist = np.extract(np.sum(img_name_w, axis=1) > 128, np.arange(img_name_w.shape[0]))

    if (len(img_name_x_hist) == 0) or (len(img_name_y_hist) == 0):
        # In some cases, we can't find any pixels.
        return None

    img_name_left = np.min(img_name_x_hist)
    img_name_right = np.max(img_name_x_hist)

    img_name_top = np.min(img_name_y_hist)
    img_name_bottom = np.max(img_name_y_hist)

    # Cropping error? should be handled gracefully.
    if not (img_name_left < img_name_right):
        return None

    img_name_w = img_name_w[img_name_top:img_name_bottom, img_name_left:img_name_right]

    img_name_w_norm = np.zeros((15, 250), dtype=np.uint8)
    img_name_w_norm[:, 0 : img_name_w.shape[1]] = cv2.resize(img_name_w, (img_name_w.shape[1], 15))

    if debug:
        print(img_name_w_norm.shape)
        cv2.imshow("name", img_name_w_norm)
        cv2.waitKey(1)

    return img_name_w_norm

def get_Ic(qNow, sFinal, qMax, x, interp, extrap, weighted_transition=True):
    '''return the corrected intensity based on circular symmetry'''
    u = numpy.sqrt(qNow*qNow + x*x) # circular-symmetric

    # divide integrand into different regions
    # interpolate from existing data
    u_in = numpy.extract(u <= sFinal, u)
    Ic_in = numpy.exp(interp(u_in))
    
    condition = numpy.multiply(sFinal < u, u <= qMax)
    u_mid = numpy.extract(condition, u)
    if u_mid.size < 2 or not weighted_transition:
        Ic_mid = numpy.exp(interp(u_mid))
    else:
        # make smooth transition between sFinal < q < qMax
        #weight = (u_mid - u_mid.min()) / (u_mid.max() - u_mid.min())
        weight = numpy.linspace(0, 1.0, u_mid.size)
        Ic_mid_in = numpy.exp(interp(u_mid))
        Ic_mid_ex = extrap.calc(u_mid)
        Ic_mid = (1-weight) * Ic_mid_in + weight * Ic_mid_ex

    # extrapolate from model beyond range of available data
    u_ex = numpy.extract(qMax < u, u)
    Ic_ex = extrap.calc(u_ex)
    
    # join the parts of the integrand
    return numpy.concatenate((Ic_in, Ic_mid, Ic_ex))

def plot_dec_boundary(theta, x, y, mapFeat=False):
    theta = np.matrix(theta).T
    if mapFeat is True:
        u = np.linspace(-1, 1.5, 50)
        v = np.linspace(-1, 1.5, 50)
    else:
        u = np.linspace(np.min(x[:,1]), np.max(x[:,1]), 50)
        v = np.linspace(np.min(x[:,2]), np.max(x[:,2]), 50)

    z = np.zeros((u.size, v.size))

    #% Evaluate z = theta*x over the grid
    for i in range(0,u.size):
        for j in range(0,v.size):
            if mapFeat is True:
                x_m = mapfeat(u[i], v[j]) 
            else:
                #x_m = np.matrix([1,u[i],v[j],u[i]**2,u[i]**4])
                x_m = np.matrix([1,u[i],v[j],u[i]**2])
            z[i,j] = x_m * theta
    z = z.T # important to transpose z before calling contour
    print z, z.shape

    # Plot z = 0
    # Notice you need to specify the range [0, 0]
    plt.figure()
    #plt.scatter(x[:, 1], x[:, 2], c=y, cmap=plt.cm.Paired)
    plt.scatter(np.extract(y==1, x[:, 1]), np.extract(y==1, x[:, 2]), c='b', marker='o', label='admitted')
    plt.scatter(np.extract(y==0, x[:, 1]), np.extract(y==0, x[:, 2]), c='r', marker='o', label='declined')
    #ax = fig.add_subplot(111)
    plt.xlabel('Test 1 scores')
    plt.ylabel('Test 2 scores')
    plt.legend()
    plt.contour(u, v, z, [0,0], linewidth = 2, cmap=plt.cm.Paired)
    plt.show()

def count_Ave_CorrCoef(NeuCountArray):
    """
    FUNCTION DESCRIPTION
        This function takes as input a N-D numpy.array of spike times, and outputs a spike count vector; the spike
        counts are averaged over a user defined interval

    :param NeuCountArray: N-D numpy.array, units are seconds, neuron spike times for each neuron stored in an numpy.array
    """


    #Spike times array turned into a numpy array
    NeuCountArray = np.array(NeuCountArray)
    print(NeuCountArray)

    #Generate array of normalized correlation coefficients (makes a symmetrical matrix)
    CountCorrCoef = np.corrcoef(NeuCountArray, rowvar = True)
    print(CountCorrCoef)

    #Keep the upper triangle of symmetrical corr. coef. matrix and the diagnol of ones (the redunate information is converted to zeros)
    UpTriCorrCoef = np.triu(CountCorrCoef)
    print (UpTriCorrCoef)

    #Eliminate the zeros (redunant info from symmetrical matrix) from the symmetrical corr. coef matrix
    NoZero = np.extract(abs(UpTriCorrCoef) > 0, UpTriCorrCoef)
    print(NoZero)

    #Eliminate the ones from the symmetrical corr. coef matrix (eliminate the variances and keep covariances)
    CovariancesOnly = np.extract(abs(NoZero) < 1, NoZero)
    print(CovariancesOnly)

    #Average the covariances
    AveCov = np.mean(abs(CovariancesOnly))

    return (AveCov)

def dataExtraction(data = 'train', class1 = 1, class0 = 0):
    import pickle, gzip
    # Load the dataset
    f = gzip.open('mnist.pkl.gz', 'rb')
    train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
    f.close()
    
    if data is 'test':
        [data, labels] = test_set
    else:
        data   = np.concatenate((train_set[0], valid_set[0]), axis = 0)
        labels = np.concatenate((train_set[1], valid_set[1]), axis = 0)
        
    y1 = np.extract(labels == class1, labels)
    X1 = data[labels == class1, :]
    
    y0 = np.extract(labels == class0, labels)
    X0 = data[labels == class0, :]

    y = np.concatenate((y1, y0), axis = 0)
    X = np.concatenate((X1, X0), axis = 0)
    
    #X = (X - np.mean(X, axis = 0)) / (1 + np.std(X, axis = 0)) # Data Normalization
    y[y == class1] = 1
    y[y == class0] = 0
    y = np.reshape(y, (np.shape(X)[0], 1))
    return y, X