def grad_EVzxVzxT_by_hyper_exact(self, EVzxVzxT_list_this, Z, A, B, hyperno):

        P = Z.shape[0]
        R = Z.shape[1]
        N = A.shape[0]

        if hyperno != 0:
            return EVzxVzxT_list_this * 0

        alpha = self.length_scale * self.length_scale

        I = np.identity(R)
        S = np.diag(B[0, :] * B[0, :])
        Sinv = np.diag(1 / B[0, :] * B[0, :])
        C = I * alpha
        Cinv = I * (1 / alpha)
        CinvSinv = 2 * Cinv + Sinv
        CinvSinv_inv = np.diag(1 / CinvSinv.diagonal())

        dC = self.length_scale * I
        dCinv = -Cinv.dot(dC).dot(Cinv)
        dCinvSinv = 2 * dCinv
        dCinvSinv_inv = -CinvSinv_inv.dot(dCinvSinv).dot(CinvSinv_inv)

        S1 = (
            dCinv
            - dCinv.dot(CinvSinv_inv).dot(Cinv)
            - Cinv.dot(dCinvSinv_inv).dot(Cinv)
            - Cinv.dot(CinvSinv_inv).dot(dCinv)
        )
        S2 = -Sinv.dot(dCinvSinv_inv).dot(Sinv)
        S3 = Sinv.dot(dCinvSinv_inv).dot(Cinv) + Sinv.dot(CinvSinv_inv).dot(dCinv)
        S4 = dCinv.dot(CinvSinv_inv).dot(Cinv) + Cinv.dot(dCinvSinv_inv).dot(Cinv) + Cinv.dot(CinvSinv_inv).dot(dCinv)

        T1s = np.tile(Z.dot(S1).dot(Z.T).diagonal(), [P, 1])
        T1 = np.tile(T1s, [N, 1, 1])
        T2s = T1s.T
        T2 = np.tile(T2s, [N, 1, 1])
        T3 = np.tile(Z.dot(S4).dot(Z.T), [N, 1, 1])
        T4 = np.tile(A.dot(S2).dot(A.T).diagonal(), [P, 1]).T
        T4 = np.expand_dims(T4, axis=2)
        T4 = np.repeat(T4, P, axis=2)
        T5 = A.dot(S3).dot(Z.T)
        T5 = np.expand_dims(T5, axis=2)
        T5 = np.repeat(T5, P, axis=2)
        T6 = np.swapaxes(T5, 1, 2)

        SCinvI = 2 * Cinv.dot(S) + I
        SCinvI_inv = np.diag(1 / SCinvI.diagonal())
        (temp, logDetSCinvI) = np.linalg.slogdet(SCinvI)
        detSCinvI = np.exp(logDetSCinvI)
        dDetSCinvI = -0.5 * np.power(detSCinvI, -0.5) * SCinvI_inv.dot(2 * dCinv).dot(S).trace()

        expTerm = EVzxVzxT_list_this / np.power(detSCinvI, -0.5)

        res = EVzxVzxT_list_this * (-0.5 * T1 - 0.5 * T2 + T3 - 0.5 * T4 + T5 + T6) + dDetSCinvI * expTerm

        res = np.sum(res, axis=0)

        return res

    def _h_arrows(self, length):
        """ length is in arrow width units """
        # It might be possible to streamline the code
        # and speed it up a bit by using complex (x,y)
        # instead of separate arrays; but any gain would be slight.
        minsh = self.minshaft * self.headlength
        N = len(length)
        length = length.reshape(N, 1)
        # This number is chosen based on when pixel values overflow in Agg
        # causing rendering errors
        # length = np.minimum(length, 2 ** 16)
        np.clip(length, 0, 2 ** 16, out=length)
        # x, y: normal horizontal arrow
        x = np.array([0, -self.headaxislength,
                      -self.headlength, 0],
                     np.float64)
        x = x + np.array([0, 1, 1, 1]) * length
        y = 0.5 * np.array([1, 1, self.headwidth, 0], np.float64)
        y = np.repeat(y[np.newaxis, :], N, axis=0)
        # x0, y0: arrow without shaft, for short vectors
        x0 = np.array([0, minsh - self.headaxislength,
                       minsh - self.headlength, minsh], np.float64)
        y0 = 0.5 * np.array([1, 1, self.headwidth, 0], np.float64)
        ii = [0, 1, 2, 3, 2, 1, 0, 0]
        X = x.take(ii, 1)
        Y = y.take(ii, 1)
        Y[:, 3:-1] *= -1
        X0 = x0.take(ii)
        Y0 = y0.take(ii)
        Y0[3:-1] *= -1
        shrink = length / minsh if minsh != 0. else 0.
        X0 = shrink * X0[np.newaxis, :]
        Y0 = shrink * Y0[np.newaxis, :]
        short = np.repeat(length < minsh, 8, axis=1)
        # Now select X0, Y0 if short, otherwise X, Y
        np.copyto(X, X0, where=short)
        np.copyto(Y, Y0, where=short)
        if self.pivot == 'middle':
            X -= 0.5 * X[:, 3, np.newaxis]
        elif self.pivot == 'tip':
            X = X - X[:, 3, np.newaxis]   # numpy bug? using -= does not
                                          # work here unless we multiply
                                          # by a float first, as with 'mid'.
        elif self.pivot != 'tail':
            raise ValueError(("Quiver.pivot must have value in {{'middle', "
                              "'tip', 'tail'}} not {0}").format(self.pivot))

        tooshort = length < self.minlength
        if tooshort.any():
            # Use a heptagonal dot:
            th = np.arange(0, 8, 1, np.float64) * (np.pi / 3.0)
            x1 = np.cos(th) * self.minlength * 0.5
            y1 = np.sin(th) * self.minlength * 0.5
            X1 = np.repeat(x1[np.newaxis, :], N, axis=0)
            Y1 = np.repeat(y1[np.newaxis, :], N, axis=0)
            tooshort = np.repeat(tooshort, 8, 1)
            np.copyto(X, X1, where=tooshort)
            np.copyto(Y, Y1, where=tooshort)
        # Mask handling is deferred to the caller, _make_verts.
        return X, Y

    def grad_EVzxVzxT_by_c(self, EVzxVzxT_list_this, Z, A, B, C, Kpred, p, r):

        P = Z.shape[0]
        R = Z.shape[1]
        N = A.shape[0]

        ainv = 1 / (self.length_scale * self.length_scale)
        siginv = 1 / (B[0, 0] * B[0, 0])

        dA = np.zeros([N, R])
        dA[:, r] = Kpred[r][:, p]

        AAt = 2 * A[:, r] * dA[:, r]

        res1 = -0.5 * np.tile(AAt, [P, 1]).T * (siginv - siginv * (1 / (siginv + 2 * ainv)) * siginv)
        res1 = np.expand_dims(res1, axis=2)
        res1 = np.repeat(res1, P, axis=2)

        res2 = dA.dot(Z.T) * (ainv * (1 / (siginv + 2 * ainv)) * siginv)
        res2 = np.expand_dims(res2, axis=2)
        res2 = np.repeat(res2, P, axis=2)

        res3 = np.swapaxes(res2, 1, 2)

        res = EVzxVzxT_list_this * (res1 + res2 + res3)

        res = np.sum(res, axis=0)

        return res

def generate_anchors(base_size=16, ratios=None, scales=None):
    """
    Generate anchor (reference) windows by enumerating aspect ratios X
    scales w.r.t. a reference window.
    """

    if ratios is None:
        ratios = np.array([0.5, 1, 2])

    if scales is None:
        scales = np.array([2 ** 0, 2 ** (1.0 / 3.0), 2 ** (2.0 / 3.0)])

    num_anchors = len(ratios) * len(scales)

    # initialize output anchors
    anchors = np.zeros((num_anchors, 4))

    # scale base_size
    anchors[:, 2:] = base_size * np.tile(scales, (2, len(ratios))).T

    # compute areas of anchors
    areas = anchors[:, 2] * anchors[:, 3]

    # correct for ratios
    anchors[:, 2] = np.sqrt(areas / np.repeat(ratios, len(scales)))
    anchors[:, 3] = anchors[:, 2] * np.repeat(ratios, len(scales))

    # transform from (x_ctr, y_ctr, w, h) -> (x1, y1, x2, y2)
    anchors[:, 0::2] -= np.tile(anchors[:, 2] * 0.5, (2, 1)).T
    anchors[:, 1::2] -= np.tile(anchors[:, 3] * 0.5, (2, 1)).T

    return anchors

def fill_between_steps(x, y1, y2=0, h_align='mid'):
    ''' Fills a hole in matplotlib: fill_between for step plots.
    Parameters :
    ------------
    x : array-like
        Array/vector of index values. These are assumed to be equally-spaced.
        If not, the result will probably look weird...
    y1 : array-like
        Array/vector of values to be filled under.
    y2 : array-Like
        Array/vector or bottom values for filled area. Default is 0.
    '''
    # First, duplicate the x values
    xx = np.repeat(x,2)
    # Now: the average x binwidth
    xstep = np.repeat((x[1:] - x[:-1]), 2)
    xstep = np.concatenate(([xstep[0]], xstep, [xstep[-1]]))
    # Now: add one step at end of row.
    #~ xx = np.append(xx, xx.max() + xstep[-1])

    # Make it possible to change step alignment.
    if h_align == 'mid':
        xx -= xstep / 2.
    elif h_align == 'right':
        xx -= xstep

    # Also, duplicate each y coordinate in both arrays
    y1 = np.repeat(y1,2)#[:-1]
    if type(y2) == np.ndarray:
        y2 = np.repeat(y2,2)#[:-1]

    return xx, y1, y2

    def updateImage(self, contingencies, rect, sup_valmax):
        """
        Makes an image of size rect from contingencies. The image is used to update a rect inside the heatmap.
        """
        interval_width = int(rect.width() / contingencies.shape[2])
        interval_height = int(rect.height() / contingencies.shape[1])

        contingencies -= np.min(contingencies)
        contingencies /= np.max(contingencies)
        contingencies = np.nan_to_num(contingencies)
        contingencies_argmax = contingencies.argmax(axis=0)
        rows, cols = np.indices(contingencies_argmax.shape)
        contingencies_valmax = contingencies[contingencies_argmax, rows, cols]

        colors_argmax = np.repeat(np.repeat(contingencies_argmax, interval_width, axis=0),
                                  interval_height, axis=1)
        colors_valmax = np.repeat(np.repeat(contingencies_valmax, interval_width, axis=0),
                                  interval_height, axis=1)

        colors = self.color_array[colors_argmax] + ((255-self.color_array[colors_argmax]) * (1-colors_valmax[:, :, None]))
        if sup_valmax:
            colors += ((255-colors) * (1-sup_valmax))

        if rect.width() == self.image_width and rect.height() == self.image_height:
            self.hmi = Heatmap(colors)
            self.plot.addItem(self.hmi)
            self.hmi.setRect(QtCore.QRectF(self.X_min, self.Y_min, self.X_max-self.X_min, self.Y_max-self.Y_min))
        else:
            self.hmi.updateImage_(colors, rect)

        return contingencies_valmax

def metric_vs_num_baggers(classifier, attack, percent_poisoning,
                         no_attack_base_error, 
                         no_attack_bag_errors, 
                         attack_base_error,
                         attack_bag_errors,
                         N,
                         metric,
                         ):    
    no_attack_base_errors = np.repeat(no_attack_base_error, N)
    attack_base_errors = np.repeat(attack_base_error, N)
    
    X = np.linspace(1, N, num=N, endpoint=True)
    
    title = get_attack_name(attack, percent_poisoning)
    
    plt.title(title, fontsize=18)
    
    plt.xlabel('Number of Baggers')
    plt.ylabel(metric)
    
    no_attack_base = plt.plot(X, no_attack_base_errors, 'b--', 
                              label=get_classifier_name(classifier))
    no_attack_bag = plt.plot(X, no_attack_bag_errors, 'b',
                             label='Bagged')
    attack_base = plt.plot(X, attack_base_errors, 'r--',
                           label=get_classifier_name(classifier, percent_poisoning))
    attack_bag = plt.plot(X, attack_bag_errors, 'r',
                          label='Bagged (poisoned)')
    
    #legend = plt.legend(loc='upper right', shadow=True, prop={'size':12})
    
    return plt

def spe_sen(target, actual):
    """Compute the (specificity,sensitivity) couple and the Matthews correlation
    coefficient for a desired Boolean function called ftar for a neuron
    implementing the Boolean function f.

    Parameters
    ----------
    target : array Bool
        actions taken
    actual : array Bool
        actions expected

    Returns
    -------
    spe : float between 0 and 1
        specificity of the response
    sen : float between 0 and 1
        sensitivity of the response
    """
    # Use the binary of the vector to see the difference between actual and
    # target
    tp = np.array(target)*2 - actual
    TN = len(np.repeat(tp, tp == 0))
    FN = len(np.repeat(tp, tp == 2))
    TP = len(np.repeat(tp, tp == 1))
    FP = len(np.repeat(tp, tp == -1))

    spe = float(TN)/(TN+FP)
    sen = float(TP)/(TP+FN)

    return spe, sen

 def _get_sorted_theta(self):
     '''sorts the integral points by bond in descending order'''
     depsf_arr = np.array([])
     V_f_arr = np.array([])
     E_f_arr = np.array([])
     xi_arr = np.array([])
     stat_weights_arr = np.array([])
     nu_r_arr = np.array([])
     r_arr = np.array([])
     for reinf in self.cont_reinf_lst:
         n_int = len(np.hstack((np.array([]), reinf.depsf_arr)))
         depsf_arr = np.hstack((depsf_arr, reinf.depsf_arr))
         V_f_arr = np.hstack((V_f_arr, np.repeat(reinf.V_f, n_int)))
         E_f_arr = np.hstack((E_f_arr, np.repeat(reinf.E_f, n_int)))
         xi_arr = np.hstack((xi_arr, np.repeat(reinf.xi, n_int)))
         stat_weights_arr = np.hstack((stat_weights_arr, reinf.stat_weights))
         nu_r_arr = np.hstack((nu_r_arr, reinf.nu_r))
         r_arr = np.hstack((r_arr, reinf.r_arr))
     argsort = np.argsort(depsf_arr)[::-1]
     # sorting the masks for the evaluation of F
     idxs = np.array([])
     for i, reinf in enumerate(self.cont_reinf_lst):
         idxs = np.hstack((idxs, i * np.ones_like(reinf.depsf_arr)))
     masks = []
     for i, reinf in enumerate(self.cont_reinf_lst):
         masks.append((idxs == i)[argsort])
     max_depsf = [np.max(reinf.depsf_arr) for reinf in self.cont_reinf_lst]
     masks = [masks[i] for i in np.argsort(max_depsf)[::-1]]
     return depsf_arr[argsort], V_f_arr[argsort], E_f_arr[argsort], \
             xi_arr[argsort], stat_weights_arr[argsort], \
             nu_r_arr[argsort], masks, r_arr[argsort]

def make_dataset1():
    '''Make a dataset of single samples with labels from which distribution they come from'''
    # now lets make some samples 
    lns      = min_max_scale(lognormal(size=bsize)) #log normal
    powers   = min_max_scale(power(0.1,size=bsize)) #power law
    norms    = min_max_scale(normal(size=bsize))    #normal
    uniforms = min_max_scale(uniform(size=bsize))    #uniform
    # add our data together
    data = np.concatenate((lns,powers,norms,uniforms))
    
    # concatenate our labels
    labels = np.concatenate((
        (np.repeat(LOGNORMAL,bsize)),
        (np.repeat(POWER,bsize)),
        (np.repeat(NORM,bsize)),
        (np.repeat(UNIFORM,bsize))))
    tsize = len(labels)
    
    # make sure dimensionality and types are right
    data = data.reshape((len(data),1))
    data = data.astype(np.float32)
    labels = labels.astype(np.int32)
    labels = labels.reshape((len(data),))
    
    return data, labels, tsize

def test_linearsvc_fit_sampleweight():
    # check correct result when sample_weight is 1
    n_samples = len(X)
    unit_weight = np.ones(n_samples)
    clf = svm.LinearSVC(random_state=0).fit(X, Y)
    clf_unitweight = svm.LinearSVC(random_state=0).\
        fit(X, Y, sample_weight=unit_weight)

    # check if same as sample_weight=None
    assert_array_equal(clf_unitweight.predict(T), clf.predict(T))
    assert_allclose(clf.coef_, clf_unitweight.coef_, 1, 0.0001)

    # check that fit(X)  = fit([X1, X2, X3],sample_weight = [n1, n2, n3]) where
    # X = X1 repeated n1 times, X2 repeated n2 times and so forth

    random_state = check_random_state(0)
    random_weight = random_state.randint(0, 10, n_samples)
    lsvc_unflat = svm.LinearSVC(random_state=0).\
        fit(X, Y, sample_weight=random_weight)
    pred1 = lsvc_unflat.predict(T)

    X_flat = np.repeat(X, random_weight, axis=0)
    y_flat = np.repeat(Y, random_weight, axis=0)
    lsvc_flat = svm.LinearSVC(random_state=0).fit(X_flat, y_flat)
    pred2 = lsvc_flat.predict(T)

    assert_array_equal(pred1, pred2)
    assert_allclose(lsvc_unflat.coef_, lsvc_flat.coef_, 1, 0.0001)

def test_linearsvr_fit_sampleweight():
    # check correct result when sample_weight is 1
    # check that SVR(kernel='linear') and LinearSVC() give
    # comparable results
    diabetes = datasets.load_diabetes()
    n_samples = len(diabetes.target)
    unit_weight = np.ones(n_samples)
    lsvr = svm.LinearSVR(C=1e3).fit(diabetes.data, diabetes.target,
                                    sample_weight=unit_weight)
    score1 = lsvr.score(diabetes.data, diabetes.target)

    lsvr_no_weight = svm.LinearSVR(C=1e3).fit(diabetes.data, diabetes.target)
    score2 = lsvr_no_weight.score(diabetes.data, diabetes.target)

    assert_allclose(np.linalg.norm(lsvr.coef_),
                    np.linalg.norm(lsvr_no_weight.coef_), 1, 0.0001)
    assert_almost_equal(score1, score2, 2)

    # check that fit(X)  = fit([X1, X2, X3],sample_weight = [n1, n2, n3]) where
    # X = X1 repeated n1 times, X2 repeated n2 times and so forth
    random_state = check_random_state(0)
    random_weight = random_state.randint(0, 10, n_samples)
    lsvr_unflat = svm.LinearSVR(C=1e3).fit(diabetes.data, diabetes.target,
                                           sample_weight=random_weight)
    score3 = lsvr_unflat.score(diabetes.data, diabetes.target,
                               sample_weight=random_weight)

    X_flat = np.repeat(diabetes.data, random_weight, axis=0)
    y_flat = np.repeat(diabetes.target, random_weight, axis=0)
    lsvr_flat = svm.LinearSVR(C=1e3).fit(X_flat, y_flat)
    score4 = lsvr_flat.score(X_flat, y_flat)

    assert_almost_equal(score3, score4, 2)

	def generate_Smolyak_points(self):
		# Merge all different approximation points
		all_points = np.array(0)
		all_complex = np.array(0)
		for i in xrange(1, self.mu + 1):
			all_points = np.append(all_points, self.extreme[i, :self.len[i]])
			all_complex = np.append(all_complex, np.repeat(i, self.len[i]))
		one_dim_len = all_complex.size
		# print all_points, all_complex
		res = np.array([all_points])
		cur_len = one_dim_len
		sum_weight = all_complex
		for i in xrange(1, self.d):
			res = np.repeat(res, one_dim_len, axis = 1)
			to_add = np.repeat(all_points[np.newaxis, :], cur_len, 0).reshape(-1)
			# print res.shape, to_add.shape
			res = np.vstack((res, to_add))
			sum_weight = np.repeat(sum_weight, one_dim_len)
			# print cur_len, all_complex
			sum_weight = sum_weight + np.repeat(all_complex[np.newaxis, :], cur_len, 0).reshape(-1)
			tokeep = (sum_weight <= self.mu)
			idx = np.arange(cur_len * one_dim_len)[tokeep]
			res = res[:, sum_weight <= self.mu]
			sum_weight = sum_weight[tokeep]
			cur_len = sum_weight.size
			# print (cur_len == res.shape[1])
			# print res.T, sum_weight
		self.grid = res.T

def exp_dead_new(file_num, name_file, imsz, wcs, flat_list, foc_list, asp_solution, dead, cut, flat_idx, step, out_path, return_dict):
    print imsz
    count = np.zeros(imsz)

    x_lim = imsz[0]
    y_lim = imsz[1]

    length = flat_list[0].shape[0]
    half_len = length/2.
    print half_len
    l = imsz[0]/10
    start = foc_list[0,1]-half_len
    print foc_list.shape
    print start.shape

    ox = np.repeat(np.arange(l)+start,length+1000)
    oy = np.tile(np.arange(length+1000)+foc_list[0,0]-half_len-500,l)
    omask = (ox>=0) & (ox<imsz[0]) & (oy>=0) & (oy<imsz[1])
    ox = ox[omask]
    oy = oy[omask]
    gl,gb = wcs.all_pix2world(oy,ox,0)
    c = SkyCoord(gl*u.degree, gb*u.degree, frame='galactic')
    rd = c.transform_to(FK5)
    for i in range(asp_solution.shape[0]):
        hrflat = flat_list[flat_idx[i]]
        foc = foc_list[i,:]#wcs.sip_pix2foc(wcs.wcs_world2pix(coo,1),1)
        if (foc[1]+half_len)>=(start+l):
            print 'update'
            start = foc[1]-half_len
            ox = np.repeat(np.arange(l)+start,length+1000)
            oy = np.tile(np.arange(length+1000)+foc[0]-half_len-500,l)
            omask = (ox>=0) & (ox<imsz[0]) & (oy>=0) & (oy<imsz[1])
            if np.sum(omask)==0:
                break
            ox = ox[omask]
            oy = oy[omask]
            gl,gb = wcs.all_pix2world(oy,ox,0)
            c = SkyCoord(gl*u.degree, gb*u.degree, frame='galactic')
            rd = c.transform_to(FK5)
        fmask = (ox>=(foc[1]-length/2)) & (ox<(foc[1]+length/2)) & (oy>=(foc[0]-length/2)) & (oy<(foc[0]+length/2))
        if np.sum(fmask)==0:
            continue
        x = ox[fmask]
        y = oy[fmask]
        xi, eta = gn.gnomfwd_simple(rd.ra.deg[fmask], rd.dec.deg[fmask], 
                                        asp_solution[i,1], asp_solution[i,2], -asp_solution[i,3],1/36000.,0.)
        px = ((xi/36000.)/(1.25/2.)*(1.25/(800* 0.001666))+1.)/2.*length
        py = ((eta/36000.)/(1.25/2.)*(1.25/(800* 0.001666))+1.)/2.*length
        pmask = (px>=0) & (px<length) & (py>=0) & (py<length)
        if np.sum(pmask)==0:
            continue
        count[x[pmask].astype(int),y[pmask].astype(int)] += \
            hrflat[px[pmask].astype(int),py[pmask].astype(int)]*step*(1-dead[i])*cut[i]
        if i%100==0:
            with open('/scratch/dw1519/galex/fits/scan_map/%s_gal_sec_exp_tmp%d.dat'%(name_file, file_num),'w') as f:
                f.write('%d'%i)
            print i
    print '%d done'%file_num
    #return_dict[file_num] = count
    np.save('%s/%s_gal_sec_exp_tmp%d.npy'%(out_path, name_file, file_num), count)

def setup_xz(nx=128,nz=128,edge=None):
    if edge==None:
        dx = (2*np.pi)/(nx)
        x = np.arange(-1*np.pi,np.pi,dx)
        #print x
        x = np.repeat(x,nz)
        x = x.reshape(nx,nz)
    else:
        #x_b = edge #in rads, nz long
        x =[]   
        #edge  = -edge
        #edge[1:10] = np.pi/2
        for i,xmax in enumerate(edge):
            #xmax = -2.5
            #xmax = np.min(edge)
            x.append(np.linspace(xmax-np.pi/10.0,xmax+np.pi/20.0,nx))
            #x.append(np.linspace(-np.pi,xmax-.4,nx))
        x = np.array(x)
        print x.shape
        x = np.transpose(x)
        
    dz = (2*np.pi)/(nz)
    z = np.arange(0,2*np.pi,dz)
    z = np.repeat(z,nx)
    z = np.transpose(z.reshape(nz,nx))


    return x,z

        def _mean4(data, offset=(0, 0), block_id=None):
            rows, cols = data.shape
            rows2, cols2 = data.shape
            pad = []
            # we assume that the chunks except the first ones are aligned
            if block_id[0] == 0:
                row_offset = offset[0] % 2
            else:
                row_offset = 0
            if block_id[1] == 0:
                col_offset = offset[1] % 2
            else:
                col_offset = 0
            row_after = (row_offset + rows) % 2
            col_after = (col_offset + cols) % 2
            pad = ((row_offset, row_after), (col_offset, col_after))

            rows2 = rows + row_offset + row_after
            cols2 = cols + col_offset + col_after

            av_data = np.pad(data, pad, 'edge')
            new_shape = (int(rows2 / 2.), 2, int(cols2 / 2.), 2)
            data_mean = np.ma.mean(av_data.reshape(new_shape), axis=(1, 3))
            data_mean = np.repeat(np.repeat(data_mean, 2, axis=0), 2, axis=1)
            data_mean = data_mean[row_offset:row_offset + rows,
                                  col_offset:col_offset + cols]
            return data_mean

 def test_arrays(self):
     s = readsav(path.join(DATA_PATH, 'struct_pointer_arrays.sav'), verbose=False)
     assert_array_identical(s.arrays.g[0], np.repeat(np.float32(4.), 2).astype(np.object_))
     assert_array_identical(s.arrays.h[0], np.repeat(np.float32(4.), 3).astype(np.object_))
     assert_(np.all(vect_id(s.arrays.g[0]) == id(s.arrays.g[0][0])))
     assert_(np.all(vect_id(s.arrays.h[0]) == id(s.arrays.h[0][0])))
     assert_(id(s.arrays.g[0][0]) == id(s.arrays.h[0][0]))

def trim_image(imData, bgValue=255):

    numRows, numCols = imData.shape

    tempCols = [np.all(imData[:, col] == np.repeat(bgValue, numRows)) for col in range(numCols)]
    tempRows = [np.all(imData[row, :] == np.repeat(bgValue, numCols)) for row in range(numRows)]

    if False not in tempRows or False not in tempCols:
        print 'The entire image is blank with background %i. Not trimming...' % bgValue
        return imData

    firstCol = tempCols.index(False)
    firstRow = tempRows.index(False)
    lastCol = -(1 + tempCols[::-1].index(False))
    lastRow = -(1 + tempRows[::-1].index(False))

    if lastRow == -1:
        if lastCol == -1:
            return imData[firstRow:, firstCol:]
        else:
            return imData[firstRow:, firstCol:(lastCol+1)]
    else:
        if lastCol == -1:
            return imData[firstRow:(lastRow+1), firstCol:]
        else:
            return imData[firstRow:(lastRow+1), firstCol:(lastCol+1)]

    def test_arrays_replicated_3d(self):
        with warnings.catch_warnings():
            warnings.filterwarnings('ignore',
                    message="warning: multi-dimensional structures")
            s = readsav(path.join(DATA_PATH,
                                  'struct_pointer_arrays_replicated_3d.sav'),
                        verbose=False)

        # Check column types
        assert_(s.arrays_rep.g.dtype.type is np.object_)
        assert_(s.arrays_rep.h.dtype.type is np.object_)

        # Check column shapes
        assert_equal(s.arrays_rep.g.shape, (4, 3, 2))
        assert_equal(s.arrays_rep.h.shape, (4, 3, 2))

        # Check values
        for i in range(4):
            for j in range(3):
                for k in range(2):
                    assert_array_identical(s.arrays_rep.g[i, j, k],
                            np.repeat(np.float32(4.), 2).astype(np.object_))
                    assert_array_identical(s.arrays_rep.h[i, j, k],
                            np.repeat(np.float32(4.), 3).astype(np.object_))
                    assert_(np.all(vect_id(s.arrays_rep.g[i, j, k]) == id(s.arrays_rep.g[0, 0, 0][0])))
                    assert_(np.all(vect_id(s.arrays_rep.h[i, j, k]) == id(s.arrays_rep.h[0, 0, 0][0])))

def test_rstrip_inplace():

    # Incorrect type
    s = np.array([1, 2, 3])
    with pytest.raises(TypeError) as exc:
        _rstrip_inplace(s)
    assert exc.value.args[0] == 'This function can only be used on string arrays'

    # Bytes array
    s = np.array(['a ', ' b', ' c c   '], dtype='S6')
    _rstrip_inplace(s)
    assert_equal(s, np.array(['a', ' b', ' c c'], dtype='S6'))

    # Unicode array
    s = np.array(['a ', ' b', ' c c   '], dtype='U6')
    _rstrip_inplace(s)
    assert_equal(s, np.array(['a', ' b', ' c c'], dtype='U6'))

    # 2-dimensional array
    s = np.array([['a ', ' b'], [' c c   ', ' a ']], dtype='S6')
    _rstrip_inplace(s)
    assert_equal(s, np.array([['a', ' b'], [' c c', ' a']], dtype='S6'))

    # 3-dimensional array
    s = np.repeat(' a a ', 24).reshape((2, 3, 4))
    _rstrip_inplace(s)
    assert_equal(s, ' a a')

    # 3-dimensional non-contiguous array
    s = np.repeat(' a a ', 1000).reshape((10, 10, 10))[:2, :3, :4]
    _rstrip_inplace(s)
    assert_equal(s, ' a a')