def loadSamplePlanktons(numSamples=100, rotate=False, dim=28):
    if dim == 28:
        if not rotate:
            from pylearn2_plankton.planktonDataPylearn2 import PlanktonData
            ds = PlanktonData(which_set='train')
            designMatrix = ds.get_data()[0] # index 1 is the label
            print "Shape of Design Matrix", np.shape(designMatrix)
            designMatrix = np.reshape(designMatrix, 
                                      (ds.get_num_examples(), 1, MAX_PIXEL, MAX_PIXEL) )
            if numSamples != 'All':
                return np.array(designMatrix[:numSamples,...], dtype=np.float32)
            else:
                return np.array(designMatrix, dtype=np.float32)
        else:
            print "Loading Rotated Data"
            designMatrix = np.load(open(os.path.join(os.environ['PYLEARN2_DATA_PATH'] ,'planktonTrainRotatedX.p'), 'r'))
            return np.reshape(np.array(designMatrix[:numSamples,...], dtype=np.float32),
                              (numSamples,1,MAX_PIXEL,MAX_PIXEL))
    elif dim == 40:
        from pylearn2_plankton.planktonData40pixels import PlanktonData
        ds = PlanktonData(which_set='train')
        designMatrix = ds.get_data()[0] # index 1 is the label
        print "Shape of Design Matrix", np.shape(designMatrix)
        designMatrix = np.reshape(designMatrix, 
                                  (ds.get_num_examples(), 1, 40, 40) )
        if numSamples != 'All':
            return np.array(designMatrix[:numSamples,...], dtype=np.float32)
        else:
            return np.array(designMatrix, dtype=np.float32)

def discrete(seg, n_classes):
    original_shape = seg.shape
    discrete_seg = seg.argmax(axis=3)
    discrete_seg = np.reshape(discrete_seg, (-1,))
    discrete_seg = np.reshape(one_hot(discrete_seg, n_classes), original_shape)

    return discrete_seg

 def Plot2d(self, fignumStart):        
     # Plot xTrue
     plt.figure(fignumStart)
     plt.imshow(self.Theta, interpolation='none')
     plt.colorbar()    
     plt.title('xTrue')
     
     # Plot the reconstructed result
     plt.figure()
     plt.imshow(np.reshape(self.ThetaEstimated, self.Theta.shape), interpolation='none')
     plt.colorbar()
     plt.title('Reconstructed x')
     
     # Plot yErr and its histogram
     yErr = self.NoisyObs - np.reshape(self._reconstructor.hx, self.NoisyObs.shape)
     plt.figure()
     plt.imshow(yErr, interpolation='none')
     plt.colorbar()
     plt.title('yErr')
             
     plt.figure()
     plt.hist(yErr.flat, 20)
     plt.title('Histogram of yErr')       
 
     plt.show()

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

 def setParams(self, params):
     #Set W1 and W2 using single paramater vector.
     W1_start = 0
     W1_end = self.hiddenLayerSize * self.inputLayerSize
     self.W1 = np.reshape(params[W1_start:W1_end], (self.inputLayerSize , self.hiddenLayerSize))
     W2_end = W1_end + self.hiddenLayerSize*self.outputLayerSize
     self.W2 = np.reshape(params[W1_end:W2_end], (self.hiddenLayerSize, self.outputLayerSize))

    def dot(self, coords_a, coords_b, frac_coords=False):
        """
        Compute the scalar product of vector(s).

        Args:
            coords_a, coords_b: Array-like objects with the coordinates.
            frac_coords (bool): Boolean stating whether the vector
                corresponds to fractional or cartesian coordinates.

        Returns:
            one-dimensional `numpy` array.
        """
        coords_a, coords_b = np.reshape(coords_a, (-1,3)), \
                             np.reshape(coords_b, (-1,3))

        if len(coords_a) != len(coords_b):
            raise ValueError("")

        if np.iscomplexobj(coords_a) or np.iscomplexobj(coords_b):
            raise TypeError("Complex array!")

        if not frac_coords:
            cart_a, cart_b = coords_a, coords_b
        else:
            cart_a = np.reshape([self.get_cartesian_coords(vec)
                                 for vec in coords_a], (-1,3))
            cart_b = np.reshape([self.get_cartesian_coords(vec)
                                 for vec in coords_b], (-1,3))

        return np.array([np.dot(a,b) for a,b in zip(cart_a, cart_b)])

    def update_state(self, time, dtime, temp, dtemp, energy, rho, F0, F,
        stran, d, elec_field, stress, statev, **kwargs):
        """Compute updated stress given strain increment"""
        log = logging.getLogger('matmodlab.mmd.simulator')

        # defaults
        cmname = '{0:8s}'.format('umat')
        dfgrd0 = reshape(F0, (3, 3), order='F')
        dfgrd1 = reshape(F, (3, 3), order='F')
        dstran = d * dtime
        ddsdde = zeros((6, 6), order='F')
        ddsddt = zeros(6, order='F')
        drplde = zeros(6, order='F')
        predef = zeros(1, order='F')
        dpred = zeros(1, order='F')
        coords = zeros(3, order='F')
        drot = eye(3)
        ndi = nshr = 3
        spd = scd = rpl = drpldt = pnewdt = 0.
        noel = npt = layer = kspt = kinc = 1
        sse = mmlabpack.ddot(stress, stran) / rho
        celent = 1.
        kstep = 1
        time = array([time, time])

        self.lib.umat(stress, statev, ddsdde,
            sse, spd, scd, rpl, ddsddt, drplde, drpldt, stran, dstran,
            time, dtime, temp, dtemp, predef, dpred, cmname, ndi, nshr,
            self.num_sdv, self.params, coords, drot, pnewdt, celent, dfgrd0,
            dfgrd1, noel, npt, layer, kspt, kstep, kinc, log.info, log.warn,
            StopFortran)

        return stress, statev, ddsdde

def calc_alm_chisq_fromfits(almfile, clfile):
    alm = pyfits.open(almfile)[0].data
    cls = pyfits.open(clfile)[0].data
    numiter = alm.shape[0]
    numchain = alm.shape[1]
    alm = np.reshape(alm, (numiter*numchain,alm.shape[2], alm.shape[3], alm.shape[4], alm.shape[5]))
    alm = alm[:, :, :, 2:, :]
    cls = cls[1:]
    cls = np.reshape(cls, (numiter * numchain, cls.shape[2], cls.shape[3]))
    if alm.shape[1] == 3:
        cls = np.concatenate((cls[:, 0:1, :], cls[:, 3:4, :], cls[:, 5:6, :]), 1)
    elif alm.shape[1] == 1:
        cls = cls[:, 0:1, :]
    cls = cls[:, :, 2:]
    cls = np.transpose(cls).copy()
    alm = np.transpose(alm).copy()
    chisq = np.zeros(cls.shape)
    for i in range(cls.shape[0]):
        l = i + 2
        for m in range(l):
            if m == 0:
                chisq[i, :, :] += alm[0, i, m, :, :] ** 2
            else:
                chisq[i, :, :] += np.sum(2 * alm[:, i, m, :, :] ** 2, 0)
        chisq[i, :, :] = chisq[i, :, :] / cls[i, :, :] / (2 * l + 1) * (l * (l + 1)) / (2 * np.pi)
    return chisq

def bp(theta, ninput, nhidden, noutput, Lambda, X, y):
    '''反向传播, 求得theta的梯度, 这里有很多计算是和fp重复的, 原因在于迭代函数
    fmin_cg的参数格式要求, 重复的程度很高, 很影响效率
    '''
    theta1 = np.reshape(theta[0:nhidden*(ninput+1)], [nhidden, ninput + 1])
    theta2 = np.reshape(theta[nhidden*(ninput+1):],  [noutput, nhidden + 1])

    m = X.shape[0]

    a1 = np_extend(X, 1)
    z2 = np.dot(a1, theta1.T)
    a2 = np_extend(sigmoid(z2), 1)
    z3 = np.dot(a2, theta2.T)
    a3 = sigmoid(z3)

    yTmp = np.eye(noutput)
    yy = yTmp[y][:]

    delta3 = a3 - yy
    delta2 = np.dot(delta3, theta2[:, 1:]) * a2[:, 1:] * (1-a2[:, 1:])

    theta1_g = np_extend(Lambda / m * theta1[:, 1:])
    theta2_g = np_extend(Lambda / m * theta2[:, 1:])
    theta1_g += 1.0 / m * np.dot(delta2.T, a1)
    theta2_g += 1.0 / m * np.dot(delta3.T, a2)

    grad = np.empty(theta.shape)
    grad[0:nhidden*(ninput+1)] = np.reshape(theta1_g, nhidden * (ninput + 1))
    grad[nhidden*(ninput+1):] = np.reshape(theta2_g, noutput * (nhidden + 1))

    return grad

def kron(a,b):
    """Kronecker product of a and b.

    The result is the block matrix::

        a[0,0]*b    a[0,1]*b  ... a[0,-1]*b
        a[1,0]*b    a[1,1]*b  ... a[1,-1]*b
        ...
        a[-1,0]*b   a[-1,1]*b ... a[-1,-1]*b

    Parameters
    ----------
    a : array, shape (M, N)
    b : array, shape (P, Q)

    Returns
    -------
    A : array, shape (M*P, N*Q)
        Kronecker product of a and b

    Examples
    --------
    >>> from scipy import kron, array
    >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
    array([[1, 1, 1, 2, 2, 2],
           [3, 3, 3, 4, 4, 4]])

    """
    if not a.flags['CONTIGUOUS']:
        a = np.reshape(a, a.shape)
    if not b.flags['CONTIGUOUS']:
        b = np.reshape(b, b.shape)
    o = np.outer(a,b)
    o = o.reshape(a.shape + b.shape)
    return np.concatenate(np.concatenate(o, axis=1), axis=1)

def ReadBPLASMA(file_name,BNORM,Ns):
    #Read the BPLASMA output file from MARS-F
    #Return BM1, BM2, BM3
    BPLASMA = num.loadtxt(open(file_name))
 
    Nm1 = BPLASMA[0,0]
    n = num.round(BPLASMA[0,2])
    Mm = num.round(BPLASMA[1:Nm1+1,0])
    Mm.resize([len(Mm),1])


    BM1 = BPLASMA[Nm1+1:,0] + BPLASMA[Nm1+1:,1]*1j
    BM2 = BPLASMA[Nm1+1:,2] + BPLASMA[Nm1+1:,3]*1j
    BM3 = BPLASMA[Nm1+1:,4] + BPLASMA[Nm1+1:,5]*1j

    BM1 = num.reshape(BM1,[Ns,Nm1],order='F')
    BM2 = num.reshape(BM2,[Ns,Nm1],order='F')
    BM3 = num.reshape(BM3,[Ns,Nm1],order='F')

    BM1 = BM1[0:Ns,:]*BNORM
    BM2 = BM2[0:Ns,:]*BNORM
    BM3 = BM3[0:Ns,:]*BNORM

    #NEED TO KNOW WHY THIS SECTION IS INCLUDED - to do with half grid???!!
    #BM2[1:,:] = BM2[0:-1,:] Needed to comment out to compare with RZPlot3
    #BM3[1:,:] = BM3[0:-1,:]

    return BM1, BM2, BM3,Mm

def load_data(dirname="cifar-10-batches-py", one_hot=False):
    tarpath = maybe_download("cifar-10-python.tar.gz",
                             "http://www.cs.toronto.edu/~kriz/",
                             dirname)
    X_train = []
    Y_train = []

    for i in range(1, 6):
        fpath = os.path.join(dirname, 'data_batch_' + str(i))
        data, labels = load_batch(fpath)
        if i == 1:
            X_train = data
            Y_train = labels
        else:
            X_train = np.concatenate([X_train, data], axis=0)
            Y_train = np.concatenate([Y_train, labels], axis=0)

    fpath = os.path.join(dirname, 'test_batch')
    X_test, Y_test = load_batch(fpath)

    X_train = np.dstack((X_train[:, :1024], X_train[:, 1024:2048],
                         X_train[:, 2048:])) / 255.
    X_train = np.reshape(X_train, [-1, 32, 32, 3])
    X_test = np.dstack((X_test[:, :1024], X_test[:, 1024:2048],
                        X_test[:, 2048:])) / 255.
    X_test = np.reshape(X_test, [-1, 32, 32, 3])

    if one_hot:
        Y_train = to_categorical(Y_train, 10)
        Y_test = to_categorical(Y_test, 10)

    return (X_train, Y_train), (X_test, Y_test)

def randmeanfor(R):

    mean_p = 0
    count_p = 0
    mean_n = 0
    count_n = 0
    shape = np.shape(R)
    R = np.reshape(R, np.size(R))

    for k in np.arange(np.size(R)):
        if R[k] > 0:
            mean_p = mean_p + R[k]
            count_p = count_p + 1
        elif R[k] < 0:
            mean_n = mean_n + R[k]
            count_n = count_n + 1

    mean_p = mean_p / count_p
    mean_n = mean_n / count_n

    for k in np.arange(size(R)):
        if R[k] > 0:
            R[k] = mean_p
        elif R[k] < 0:
            R[k] = mean_n
    R = np.reshape(R, shape)

    return R

def load_data_wrapper():
    """Return a tuple containing ``(training_data, validation_data,
    test_data)``. Based on ``load_data``, but the format is more
    convenient for use in our implementation of neural networks.

    In particular, ``training_data`` is a list containing 50,000
    2-tuples ``(x, y)``.  ``x`` is a 784-dimensional numpy.ndarray
    containing the input image.  ``y`` is a 10-dimensional
    numpy.ndarray representing the unit vector corresponding to the
    correct digit for ``x``.

    ``validation_data`` and ``test_data`` are lists containing 10,000
    2-tuples ``(x, y)``.  In each case, ``x`` is a 784-dimensional
    numpy.ndarry containing the input image, and ``y`` is the
    corresponding classification, i.e., the digit values (integers)
    corresponding to ``x``.

    Obviously, this means we're using slightly different formats for
    the training data and the validation / test data.  These formats
    turn out to be the most convenient for use in our neural network
    code."""
    tr_d, va_d, te_d = load_data()
    training_inputs = [np.reshape(x, (784, 1)) for x in tr_d[0]]
    training_results = [vectorized_result(y) for y in tr_d[1]]
    training_data = zip(training_inputs, training_results)
    validation_inputs = [np.reshape(x, (784, 1)) for x in va_d[0]]
    validation_data = zip(validation_inputs, va_d[1])
    test_inputs = [np.reshape(x, (784, 1)) for x in te_d[0]]
    test_data = zip(test_inputs, te_d[1])
    return (training_data, validation_data, test_data)

    def test_cmac(self):
        input_train = np.reshape(np.linspace(0, 2 * np.pi, 100), (100, 1))
        input_train_before = input_train.copy()
        input_test = np.reshape(np.linspace(np.pi, 2 * np.pi, 50), (50, 1))
        input_test_before = input_test.copy()

        target_train = np.sin(input_train)
        target_train_before = target_train.copy()
        target_test = np.sin(input_test)

        cmac = algorithms.CMAC(
            quantization=100,
            associative_unit_size=32,
            step=0.2,
            verbose=False,
        )
        cmac.train(input_train, target_train, epochs=100)

        predicted_test = cmac.predict(input_test)
        predicted_test = predicted_test.reshape((len(predicted_test), 1))
        error = metrics.mean_absolute_error(target_test, predicted_test)

        self.assertAlmostEqual(error, 0.0024, places=4)

        # Test that algorithm didn't modify data samples
        np.testing.assert_array_equal(input_train, input_train_before)
        np.testing.assert_array_equal(input_train, input_train_before)
        np.testing.assert_array_equal(target_train, target_train_before)

def predict(theta, X, ninput, nhidden, noutput):
    theta1 = np.reshape(theta[0:nhidden*(ninput+1)], [nhidden, ninput + 1])
    theta2 = np.reshape(theta[nhidden*(ninput+1):],  [noutput, nhidden + 1])

    h1 = sigmoid(np.dot(np_extend(X, 1), theta1.T))
    h2 = sigmoid(np.dot(np_extend(h1, 1), theta2.T))    
    return np.argmax(h2, axis = 1)

    def build(self):
        import re

        with open(self.filepath, 'r') as f:
            # Currently these are unused, but they are in the format
            # Could possibly store as metadata?
            # Assume first result for regexes
            re_rows = re.compile(u'([0-9]+) rows')
            n_rows = int(re_rows.findall(f.readline())[0])
            re_cols = re.compile(u'([0-9]+) columns')
            n_cols = int(re_cols.findall(f.readline())[0])

        # This also loads the mask
        #   >>> image_data[:, 0]
        image_data = np.loadtxt(self.filepath, skiprows=3, unpack=True)

        # Replace the lowest value with nan so that we can render properly
        data_view = image_data[:, 1:]
        corrupt_value = np.min(data_view)
        data_view[np.any(np.isclose(data_view, corrupt_value), axis=1)] = np.nan

        return MaskedImage(
            np.rollaxis(np.reshape(data_view, [n_rows, n_cols, 3]), -1),
            np.reshape(image_data[:, 0], [n_rows, n_cols]).astype(np.bool),
            copy=False)

def ClusterObjects(farn, struct_elem):
    magn_img = farn.magnitude_image
    dir_img = farn.direction_image    
    
    bin_img = np.zeros(shape=(magn_img.shape[0], magn_img.shape[1]), dtype=np.uint8)
    bin_img[magn_img < 25] = 0
    bin_img[magn_img >= 25] = 1
    
    bin_img = ndimage.binary_dilation(bin_img, structure=struct_elem, iterations=3).astype(bin_img.dtype)

    labels, nb_labels = Morphology.ConnenctedComponents(bin_img)    
    filt_labels, areas, nb_new_labels = Morphology.FilterArea(bin_img, labels, nb_labels, 480)
    
    
    temp_magn = ndimage.mean(magn_img, filt_labels, range(nb_new_labels + 1))
    temp_dir = ndimage.mean(dir_img, filt_labels, range(nb_new_labels + 1))
    
    data = np.concatenate((np.reshape(temp_magn, (-1,1)), np.reshape(temp_dir, (-1,1))), axis=1)
    
    clusters = -1
    if nb_new_labels >= 1:
        Y = pdist(data, 'euclidean')
        agglo = AgglomerativeClustering.Agglomerative(Y, 50.)
        agglo.AggloClustering(criterion = 'distance', method = 'single', metric = 'euclidean', normalized = False)

        clusters = agglo.clusters
             
    bin_img[filt_labels == 0] = 0
    bin_img[filt_labels >= 1] = 1
    
    
    
    return bin_img, nb_new_labels, temp_magn, temp_dir, data, clusters

def timeit_plot3D(data, xlabel='xlabel', ylabel='ylabel', **kwargs): 
    """3D plot of timeit data, one chart per function. 
    """
    dataT = {}
    figs = []
    series = kwargs.get('series', (0,1))
    cmap = kwargs.get('cmap', cm.coolwarm)
    for k, v in data.items():
        dataT[k] = zip(*v)
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        X, Y, Z = dataT[k][series[0]], dataT[k][series[1]], dataT[k][-1]
        wide, tall = (max(X)-min(X)+1), (max(Y)-min(Y)+1)
        intervalX = max(X) - min(heapq.nlargest(2,set(X)))
        intervalY = max(Y) - min(heapq.nlargest(2,set(Y)))
        wide, tall = 1+wide/intervalX, 1+tall/intervalY
        X = np.reshape(X, [wide, tall])
        Y = np.reshape(Y, [wide, tall])
        # TODO: BUG: fix so that Z transposes with x & y reversed
        Z = np.reshape(Z, [wide, tall])
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cmap, linewidth=0, antialiased=False)
        ax.zaxis.set_major_locator(LinearLocator(10))
        ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(substitute_titles(k,series))
        fig.colorbar(surf, shrink=0.5, aspect=5)
        figs.append(fig)
    return figs

 def test(self):
     # Expected
     in_channels = 3
     in_dim = 11
     out_channels = 5
     out_dim = (in_dim/2 + 1)
     img = np.arange(0,in_dim*in_dim*in_channels*1, dtype=np.float32)
     img = np.reshape(img,[in_dim,in_dim,in_channels,1])
     filter = np.arange(0,3*3*in_channels*out_channels, dtype=np.float32)
     filter = np.reshape(filter,[3,3,in_channels,out_channels])
     bias = np.zeros([5])
     expected = np.zeros([out_dim,out_dim,out_channels])
     for och in range(out_channels):
         tmp = np.zeros([out_dim,out_dim,1])
         for ich in range(in_channels):
             imgslice = np.reshape(img[:,:,ich,0],[in_dim,in_dim])
             filterslice = np.reshape(filter[:,:,ich,och],[3,3])
             tmp += np.reshape(convolve(imgslice,filterslice,mode='constant',cval = 0.0)[::2,::2] , [out_dim, out_dim, 1])
         expected[:,:,och] = np.squeeze(tmp) + bias[och]
         
     # test
     owlimg = owl.from_numpy(np.transpose(img))
     owlfilter = owl.from_numpy(np.transpose(filter))
     owlbias = owl.from_numpy(bias)
     convolver = owl.conv.Convolver(1,1,2,2)   
     test = convolver.ff(owlimg, owlfilter, owlbias)
     
     print 'Expected\n',expected
     print "Actual\n",test.to_numpy()
     self.assertTrue(np.allclose(expected, test))