def downsample(data):
    data["_tr_X"] = np.zeros((len(data["tr_X"]), 14 * 14), dtype="float32")
    data["_va_X"] = np.zeros((len(data["va_X"]), 14 * 14), dtype="float32")
    data["_te_X"] = np.zeros((len(data["te_X"]), 14 * 14), dtype="float32")

    for i in xrange(0, len(data["tr_X"])):
        data["_tr_X"][i] = block_reduce(
            data["tr_X"][i].reshape(data["shape_x"]), block_size=(2, 2), func=np.mean
        ).flatten()

    for i in xrange(0, len(data["va_X"])):
        data["_va_X"][i] = block_reduce(
            data["va_X"][i].reshape(data["shape_x"]), block_size=(2, 2), func=np.mean
        ).flatten()

    for i in xrange(0, len(data["te_X"])):
        data["_te_X"][i] = block_reduce(
            data["te_X"][i].reshape(data["shape_x"]), block_size=(2, 2), func=np.mean
        ).flatten()

    data["tr_X"] = data["_tr_X"]
    data["va_X"] = data["_va_X"]
    data["te_X"] = data["_te_X"]

    data["shape_x"] = (14, 14)
    data["n_x"] = 14 * 14
    return data

def test_invalid_block_size():
    image = np.arange(4 * 6).reshape(4, 6)

    with testing.raises(ValueError):
        block_reduce(image, [1, 2, 3])
    with testing.raises(ValueError):
        block_reduce(image, [1, 0.5])

    def get_recorded_data(self, vec):
        '''Extract recorded voltages and timestamps given the recorded Vector instance.  
        If self.stimulus_sampling_rate is smaller than self.simulation_sampling_rate, 
        resample to self.stimulus_sampling_rate.

        Parameters
        ----------
        vec : neuron.Vector 
           constructed by self.record_values

        Returns
        -------
        dict with two keys: 'v' = numpy.ndarray with voltages, 't' = numpy.ndarray with timestamps

        '''
        junction_potential = self.description.data['fitting'][0]['junction_potential']
        
        v = np.array(vec["v"])
        t = np.array(vec["t"])

        if self.stimulus_sampling_rate < self.simulation_sampling_rate:
            factor = self.simulation_sampling_rate / self.stimulus_sampling_rate
                
            Utils._log.debug("subsampling recorded traces by %dX", factor)
            v = block_reduce(v, (factor,), np.mean)[:len(self.stim_curr)]
            t = block_reduce(t, (factor,), np.min)[:len(self.stim_curr)]

        mV = 1.0e-3
        v = (v - junction_potential) * mV
        
        return { "v": v, "t": t }

def k8_down_hints(x):
    RGB = x[:, :, 0:3].astype(np.float)
    A = x[:, :, 3:4].astype(np.float)
    RGB = RGB * A / 255.0
    RGB = block_reduce(RGB, (8, 8, 1), np.max)
    A = block_reduce(A, (8, 8, 1), np.max)
    y = np.concatenate([RGB, A], axis=2)
    return y

def mycorrelate2d(fixed,moved,skip=1):
    """a 2d correlation function for numpy 2d matrices
    
    arguments
    fixed) is the larger matrix which should stay still 
    moved) is the smaller matrix which should move left/right up/down and sample the correlation
    skip) is the number of positions to skip over when sampling, 
    so if skip =3 it will sample at shift 0,0 skip,0 2*skip,0... skip,0 skip,skip...
    
    returns
    corrmat) the 2d matrix with the corresponding correlation coefficents of the data at that offset
    note the 0,0 entry of corrmat corresponds to moved(0,0) corresponding to fixed(0,0)
    and the 1,1 entry of corrmat corresponds to moved(0,0) corresponding to fixed(skip,skip)
    NOTE) the height of corrmat is given by corrmat.height=ceil((fixed.height-moved.height)/skip)
    and the width in a corresonding manner.
    NOTE)the standard deviation is measured over the entire dataset, so particular c values can be above 1.0
    if the variance in the subsampled region of fixed is lower than the variance of the entire matrix
    
    """
    if skip>1:
        fixed = block_reduce(fixed,block_size = (int(skip),int(skip)),func = np.mean,cval=np.mean(fixed))
        moved = block_reduce(moved,block_size = (int(skip),int(skip)),func = np.mean,cval=np.mean(moved))

    (fh,fw)=fixed.shape
    (mh,mw)=moved.shape
    deltah=(fh-mh)
    deltaw=(fw-mw)
    #if (deltah<1 or deltaw<1):
    #    return
    #fixed=fixed-fixed.mean()
    #fixed=fixed/fixed.std()
    #moved=moved-moved.mean()
    #moved=moved/moved.std()
    # ch=np.ceil(deltah*1.0/skip)
    # cw=np.ceil(deltaw*1.0/skip)
    
    # corrmat=np.zeros((ch,cw))
    
    # #print (fh,fw,mh,mw,ch,cw,skip,deltah,deltaw)
    # for shiftx in range(0,deltaw,skip):
    #     for shifty in range(0,deltah,skip):
    #         fixcut=fixed[shifty:shifty+mh,shiftx:shiftx+mw]
    #         corrmat[shifty/skip,shiftx/skip]=(fixcut*moved).sum()
           
    # corrmat=corrmat/(mh*mw)
    

    corrmatt = norm_xcorr.norm_xcorr(moved,fixed, trim=True, method='fourier')

    print 'corrmatt',corrmatt.shape
    print 'moved',moved.shape
    print 'fixed',fixed.shape

    #image_product = np.fft.fft2(fixed) * np.fft.fft2(moved).conj()
    #corrmat = np.fft.fftshift(np.fft.ifft2(image_product))

    return corrmatt

def align_georasters(raster,alignraster,how=np.mean,cxsize=None,cysize=None):
    '''
    Align two rasters so that data overlaps by geographical location
    Usage: (alignedraster_o, alignedraster_a) = AlignRasters(raster, alignraster, how=np.mean)
    where 
        raster: string with location of raster to be aligned
        alignraster: string with location of raster to which raster will be aligned
        how: function used to aggregate cells (if the rasters have different sizes)
    It is assumed that both rasters have the same size
    '''
    (NDV1, xsize1, ysize1, GeoT1, Projection1, DataType1)=(raster.nodata_value, raster.shape[1], raster.shape[0], raster.geot, raster.projection, raster.datatype)
    (NDV2, xsize2, ysize2, GeoT2, Projection2, DataType2)=(alignraster.nodata_value, alignraster.shape[1], alignraster.shape[0], alignraster.geot, alignraster.projection, alignraster.datatype)
    if Projection1.ExportToMICoordSys()==Projection2.ExportToMICoordSys():
        blocksize=(np.round(max(GeoT2[1]/GeoT1[1],1)),np.round(max(GeoT2[-1]/GeoT1[-1],1)))
        mraster=raster.raster
        mmin=mraster.min()
        if block_reduce!=(1,1):
            mraster=block_reduce(mraster,blocksize,func=how)
        blocksize=(np.round(max(GeoT1[1]/GeoT2[1],1)),np.round(max(GeoT1[-1]/GeoT2[-1],1)))
        araster=alignraster.raster
        amin=araster.min()
        if block_reduce!=(1,1):
            araster=block_reduce(araster,blocksize,func=how)
        if GeoT1[0]<=GeoT2[0]:
            row3,mcol=map_pixel(GeoT2[0], GeoT2[3], GeoT1[1] *blocksize[0],GeoT1[-1]*blocksize[1], GeoT1[0], GeoT1[3])
            acol=0
        else:
            row3,acol=map_pixel(GeoT1[0], GeoT1[3], GeoT2[1],GeoT2[-1], GeoT2[0], GeoT2[3])
            mcol=0
        if GeoT1[3]<=GeoT2[3]:
            arow,col3=map_pixel(GeoT1[0], GeoT1[3], GeoT2[1],GeoT2[-1], GeoT2[0], GeoT2[3])
            mrow=0
        else:
            mrow,col3=map_pixel(GeoT2[0], GeoT2[3], GeoT1[1] *blocksize[0],GeoT1[-1]*blocksize[1], GeoT1[0], GeoT1[3])
            arow=0
        mraster=mraster[mrow:,mcol:]
        araster=araster[arow:,acol:]
        if cxsize and cysize:
            araster=araster[:cysize,:cxsize]
            mraster=mraster[:cysize,:cxsize]
        else:
            rows = min(araster.shape[0],mraster.shape[0])
            cols = min(araster.shape[1],mraster.shape[1])
            araster=araster[:rows,:cols]
            mraster=mraster[:rows,:cols]
        mraster=np.ma.masked_array(mraster,mask=mraster<mmin, fill_value=NDV1)
        araster=np.ma.masked_array(araster,mask=araster<amin, fill_value=NDV2)
        GeoT=(max(GeoT1[0],GeoT2[0]), GeoT1[1]*blocksize[0], GeoT1[2], min(GeoT1[3],GeoT2[3]), GeoT1[4] ,GeoT1[-1]*blocksize[1])
        mraster=GeoRaster(mraster, GeoT, projection=Projection1, nodata_value=NDV1, datatype=DataType1)
        araster=GeoRaster(araster, GeoT, projection=Projection2, nodata_value=NDV2, datatype=DataType2)
        return (mraster,araster)
    else:
        print("Rasters need to be in same projection")
        return (-1,-1)

def test_block_reduce_mean():
    image1 = np.arange(4 * 6).reshape(4, 6)
    out1 = block_reduce(image1, (2, 3), func=np.mean)
    expected1 = np.array([[  4.,   7.],
                          [ 16.,  19.]])
    assert_equal(expected1, out1)

    image2 = np.arange(5 * 8).reshape(5, 8)
    out2 = block_reduce(image2, (4, 5), func=np.mean)
    expected2 = np.array([[14. , 10.8],
                          [ 8.5,  5.7]])
    assert_equal(expected2, out2)

def test_block_reduce_min():
    image1 = np.arange(4 * 6).reshape(4, 6)
    out1 = block_reduce(image1, (2, 3), func=np.min)
    expected1 = np.array([[ 0, 3],
                          [12, 15]])
    assert_equal(expected1, out1)

    image2 = np.arange(5 * 8).reshape(5, 8)
    out2 = block_reduce(image2, (4, 5), func=np.min)
    expected2 = np.array([[0, 0],
                          [0, 0]])
    assert_equal(expected2, out2)

def test_block_reduce_max():
    image1 = np.arange(4 * 6).reshape(4, 6)
    out1 = block_reduce(image1, (2, 3), func=np.max)
    expected1 = np.array([[ 8, 11],
                          [20, 23]])
    assert_equal(expected1, out1)

    image2 = np.arange(5 * 8).reshape(5, 8)
    out2 = block_reduce(image2, (4, 5), func=np.max)
    expected2 = np.array([[28, 31],
                          [36, 39]])
    assert_equal(expected2, out2)

def test_block_reduce_sum():
    image1 = np.arange(4 * 6).reshape(4, 6)
    out1 = block_reduce(image1, (2, 3))
    expected1 = np.array([[ 24,  42],
                          [ 96, 114]])
    assert_equal(expected1, out1)

    image2 = np.arange(5 * 8).reshape(5, 8)
    out2 = block_reduce(image2, (3, 3))
    expected2 = np.array([[ 81, 108,  87],
                          [174, 192, 138]])
    assert_equal(expected2, out2)

def Location_Shape(img,segments,segments_label):
    # 72-D Feature
    row,col = segments.shape
    location_block_row = int(math.ceil(row/6.))
    location_block_col = int(math.ceil(col/6.))

    Location_Shape_Features = []
    for label in range(len(segments_label)):
        # Make mask for each segment
        seg_mask = Segment_Mask(segments, label)

        ### Get Location Features
        # Downsample to 6*6
        try:
            downsample = block_reduce(seg_mask, block_size=(location_block_row, location_block_col), cval = 0, func=np.max)
            # Convert to 36-D Location Features
            Location_Features = downsample.flatten().tolist()
        except:
            Location_Features = [0 for x in range(36)]

        ### Get Shape Features
        # Bounding Box
        left,up,right,down = Image.fromarray(np.uint8(seg_mask)).getbbox()

        # Cropped the mask
        cropped_mask =  seg_mask[up:down,left:right]

        # Downsample to 6*6
        cropped_row,cropped_col = cropped_mask.shape

        ### When the number is too small, there would be a bug
        ### Consider this special situation
        if cropped_row < 26:
            cropped_mask = cropped_mask[:(cropped_row-cropped_row%6),:]
        if cropped_col < 26:
            cropped_mask = cropped_mask[:,:(cropped_col-cropped_col%6)]

        cropped_row,cropped_col = cropped_mask.shape
        cropped_block_row = int(math.ceil(cropped_row/6.))
        cropped_block_col = int(math.ceil(cropped_col/6.))

        try:
            downsample = block_reduce(cropped_mask, block_size=(cropped_block_row, cropped_block_col), cval = 0, func=np.max)
            # Convert to 36-D Shape Features
            Shape_Features = downsample.flatten().tolist()
        except:
            Shape_Features = [0 for x in range(36)]

        Location_Shape_Features.append(Location_Features+Shape_Features)

    return Location_Shape_Features

def test_block_reduce_median():
    image1 = np.arange(4 * 6).reshape(4, 6)
    out1 = block_reduce(image1, (2, 3), func=np.median)
    expected1 = np.array([[  4.,   7.],
                          [ 16.,  19.]])
    assert_equal(expected1, out1)

    image2 = np.arange(5 * 8).reshape(5, 8)
    out2 = block_reduce(image2, (4, 5), func=np.median)
    expected2 = np.array([[ 14.,  6.5],
                          [  0.,  0. ]])
    assert_equal(expected2, out2)

    image3 = np.array([[1, 5, 5, 5], [5, 5, 5, 1000]])
    out3 = block_reduce(image3, (2, 4), func=np.median)
    assert_equal(5, out3)

def block_reduce(data, block_size, func=np.sum):
    
    # Backported from Astropy 1.1 for compatibility 

    from skimage.measure import block_reduce

    data = np.asanyarray(data)

    block_size = np.atleast_1d(block_size)
    if data.ndim > 1 and len(block_size) == 1:
        block_size = np.repeat(block_size, data.ndim)

    if len(block_size) != data.ndim:
        raise ValueError('`block_size` must be a scalar or have the same '
                         'length as `data.shape`')

    block_size = np.array([int(i) for i in block_size])
    size_resampled = np.array(data.shape) // block_size
    size_init = size_resampled * block_size

    # trim data if necessary
    for i in range(data.ndim):
        if data.shape[i] != size_init[i]:
            data = data.swapaxes(0, i)
            data = data[:size_init[i]]
            data = data.swapaxes(0, i)

    return block_reduce(data, tuple(block_size), func=func)

def convolve(img, sigma=4):
    '''
    2D Gaussian convolution
    '''

    if img.sum() == 0:
        return img

    img_pad = np.zeros((3 * img.shape[0], 3 * img.shape[1]))
    img_pad[img.shape[0]:2 * img.shape[0], img.shape[1]:2 * img.shape[1]] = img

    x = np.arange(3 * img.shape[0])
    y = np.arange(3 * img.shape[1])
    g = spinterp.interp2d(y, x, img_pad, kind='linear')

    if img.shape[0] == 16:
        upsample = 4
        offset = -(1 - .625)
    elif img.shape[0] == 8:
        upsample = 8
        offset = -(1 - .5625)
    else:
        raise NotImplementedError
    ZZ_on = g(offset + np.arange(0, img.shape[1] * 3, 1. / upsample), offset + np.arange(0, img.shape[0] * 3, 1. / upsample))
    ZZ_on_f = gaussian_filter(ZZ_on, float(sigma), mode='constant')

    z_on_new = block_reduce(ZZ_on_f, (upsample, upsample))
    z_on_new = z_on_new / z_on_new.sum() * img.sum()
    z_on_new = z_on_new[img.shape[0]:2 * img.shape[0], img.shape[1]:2 * img.shape[1]]

    return z_on_new

    def test_block_reduce_mask_array(self):
        test_array = np.arange(16).reshape((4, 4))
        assert_array_equal(test_array[0], [0, 1, 2, 3], verbose=True)
        assert_array_equal(test_array[1], [4, 5, 6, 7], verbose=True)
        assert_array_equal(test_array[2], [8, 9, 10, 11], verbose=True)
        assert_array_equal(test_array[3], [12, 13, 14, 15], verbose=True)

        mask_array = np.full((4, 4), False, dtype=np.bool)
        mask_array[1:3, 1:3] = True
        assert_array_equal(mask_array[0], [False, False, False, False], verbose=True)
        assert_array_equal(mask_array[1], [False,  True,  True, False], verbose=True)
        assert_array_equal(mask_array[2], [False,  True,  True, False], verbose=True)
        assert_array_equal(mask_array[3], [False, False, False, False], verbose=True)

        masked_array = np.ma.array(test_array, mask=mask_array)
        self.assertTrue(np.ma.is_masked(masked_array))
        assert_array_equal(masked_array[0], [0, 1, 2, 3], verbose=True)
        assert_array_equal(masked_array[1], [4, np.nan, np.nan, 7], verbose=True)
        assert_array_equal(masked_array[2], [8, np.nan, np.nan, 11], verbose=True)
        assert_array_equal(masked_array[3], [12, 13, 14, 15], verbose=True)

        mean_aggregated_array = block_reduce(masked_array, (2, 2), func=np.mean)

        self.assertEqual((2, 2), mean_aggregated_array.shape)
        # The mask is ignored in the block_reduce function
        assert_array_equal(mean_aggregated_array[0], [(0 + 1 + 4 + 5) / 4, (2 + 3 + 6 + 7) / 4], verbose=True)
        assert_array_equal(mean_aggregated_array[1], [(8 + 9 + 12 + 13) / 4, (10 + 11 + 14 + 15) / 4], verbose=True)

def get_patch(image, coords, offset, nodule_list, patch_flag=True):
    xyz = image[int(coords[0] - offset): int(coords[0] + offset), int(coords[1] - offset): int(coords[1] + offset),
          int(coords[2] - offset): int(coords[2] + offset)]

    if patch_flag:
        output = np.expand_dims(xyz, axis=-1)
    else:
        # resize xyz
        """
        xyz = scipy.ndimage.zoom(input=xyz, zoom=1/8, order=1) # nearest
        xyz = np.where(xyz > 0, 1.0, 0.0)
        """
        xyz = block_reduce(xyz, (9, 9, 9), np.max)
        output = np.expand_dims(xyz, axis=-1)

        output = indices_to_one_hot(output.astype(np.int32), 2)
        output = np.reshape(output, (label_size, label_size, label_size, 2))
        output = output.astype(np.float32)

        # print('------------------')
        # print(output)

        # print(output)
        # print(np.shape(output))

    nodule_list.append(output)

    def downsample(self, factor, method=np.nansum):
        """
        Down sample image by a given factor.

        The image is down sampled using `skimage.measure.block_reduce`. If the
        shape of the data is not divisible by the down sampling factor, the image
        must be padded beforehand to the correct shape.

        Parameters
        ----------
        factor : int
            Down sampling factor.
        method : np.ufunc (np.nansum), optional
            Method how to combine the image blocks.

        Returns
        -------
        image : `SkyImage`
            Down sampled image.
        """
        from skimage.measure import block_reduce

        shape = self.data.shape

        if not (np.mod(shape, factor) == 0).all():
            raise ValueError('Data shape {0} is not divisable by {1} in all axes.'
                             'Pad image prior to downsamling to correct'
                             ' shape.'.format(shape, factor))

        data = block_reduce(self.data, (factor, factor), method)

        # Adjust WCS
        wcs = get_resampled_wcs(self.wcs, factor, downsampled=True)
        return SkyImage(data=data, wcs=wcs)

 def downsample(self, factor, preserve_counts=True):
     from skimage.measure import block_reduce
     geom = self.geom.downsample(factor)
     block_size = tuple([factor, factor] + [1] * (self.geom.ndim - 2))
     data = block_reduce(self.data, block_size[::-1], np.nansum)
     if not preserve_counts:
         data /= factor**2
     return self.__class__(geom, data, meta=copy.deepcopy(self.meta))

def downsample_array(a, ds):
    '''
    Downsample ndarray by factors in vector ds (same length as dimension of a).
    '''
    from skimage.measure import block_reduce
    a_downsampled = block_reduce(a, ds)

    return a_downsampled

  def featureExtractor(self):

    screen = np.array(ImageGrab.grab(bbox = (50, 120, 1250, 650)))
    small_screen = block_reduce(screen, (530/224, 1200/224, 1), np.max)
    features = np.array(self.tf.transform(small_screen[:,:,0:3])).flatten()

    feat = dict(zip(range(2000), features))

    return  feat