def test_autos_l2(corruption=0):
#load data
    dataset='mnist.pkl.gz'
    datasets = load_data(dataset)
    train_set_x, train_set_y = datasets[0]
    valid_set_x,valid_set_y = datasets[1]


#test against validation set
    n_hiddens = [10,25,50,100]
    x = T.matrix('x')
    rng = numpy.random.RandomState(123)
    theano_rng = RandomStreams(rng.randint(2 ** 30))

    for n_hidden in n_hiddens:
    #load model
        da = dA(numpy_rng = rng,
                theano_rng = theano_rng,
                input=x,
                n_visible=28*28,
                n_hidden=n_hidden)
        da.load(open('../data/dA_l2/dA_l2_nhid'+str(n_hidden)+'_corr'+str(corruption)+'.p','r'))
        
        reconstructed = da.get_reconstructed_input(input=valid_set_x)
        image = Image.fromarray(tile_raster_images(X=reconstructed.eval(),img_shape=(28, 28), tile_shape=(10, 10),tile_spacing=(1, 1)))
        image.save('../data/dA_l2/pics/dAs_reconstructed_nhid'+str(da.n_hidden)+'_corr'+str(corruption)+'.png')

    image = Image.fromarray(tile_raster_images(X=valid_set_x.get_value(),img_shape=(28, 28), tile_shape=(10, 10),tile_spacing=(1, 1)))
    image.save('../data/dA_l2/pics/original.png')

    def visualize_hidden(self,threshold,bounds):
        print '\nSaving hidden layer filters...\n'

        #Visualizing 1st hidden layer
        f_name = 'my_filter_layer_0.png'
        im_side = sqrt(self.i_size)
        im_count = int(sqrt(self.h_sizes[0]))
        image = Image.fromarray(tile_raster_images(
        X=self.sa_layers[0].W1.get_value(borrow=True).T,
        img_shape=(im_side, im_side), tile_shape=(im_count, im_count),
        tile_spacing=(1, 1)))
        image.save(f_name)

        index = T.lscalar('index')
        max_inputs =[]
        #Higher level hidden layers
        for i in xrange(1,self.n_layers):
            print "Calculating features for higher layers\n"
            inp = 1e-8 + np.random.random_sample((self.i_size,))*0.05
            inp = np.asarray(inp,dtype=config.floatX)
            input = shared(value=inp, name='input',borrow=True)

            max_ins = self.nlopt_optimization(input,threshold,bounds,i)

            f_name = 'my_filter_layer_'+str(i)+'.png'
            im_side = sqrt(self.i_size)
            im_count = int(sqrt(self.h_sizes[i]))
            image = Image.fromarray(tile_raster_images(
                X=max_ins,
                img_shape=(im_side, im_side), tile_shape=(im_count, im_count),
                tile_spacing=(1, 1)))
            image.save(f_name)

def test(dataset = 'mnist.pkl.gz', output_folder = 'plots'):
    
    datasets = load_data(dataset)
    train_set_x, train_set_y = datasets[0]
    
    if not os.path.isdir(output_folder):
        os.makedirs(output_folder)
    os.chdir(output_folder)
    
    rng = numpy.random.RandomState(123)
    theano_rng = RandomStreams(rng.randint(2 ** 30))
    
    input = T.matrix('input')
    
    output = get_corrupted_input_gaussian(theano_rng = theano_rng, input = input)
    
    corrupt = theano.function([input], output)
    
    
    mnist_noise = corrupt(train_set_x.get_value(borrow = True))
    mnist_noise = theano.shared(value=mnist_noise, name='mnist_noise', borrow = True)
#     print train_set_x.get_value(borrow=True)[0]
#     print mnist_noise.get_value(borrow=True)[0]
    
    image_clean = Image.fromarray(tile_raster_images(X = train_set_x.get_value(borrow = True),
                                               img_shape=(28, 28), tile_shape=(1, 6),
                                               tile_spacing=(1,1)))
    image_clean.save('clean_6.png')
    
    image_noise = Image.fromarray(tile_raster_images(X = mnist_noise.get_value(borrow = True),
                                               img_shape=(28, 28), tile_shape=(1, 6),
                                               tile_spacing=(1,1)))
    image_noise.save('noise_6.png')
    
    print 'Done!'

def test(Weights, counter, ext, channel=1):
	"""this is an utility that takes weights and plot there feature as image"""
	tile_shape = (8, 8)
	image_resize_shape = (10, 10)
	img_shape = (window_size, window_size)
	newimg = None
	if channel == 1:
		img = tile_raster_images(X=Weights.T, img_shape=img_shape, tile_shape=tile_shape, tile_spacing=(1, 1))
		newimg = np.zeros((img.shape[0]*image_resize_shape[0], img.shape[1]*image_resize_shape[1]))
		for i in xrange(img.shape[0]):
			for j in xrange(img.shape[1]):
				newimg[i*image_resize_shape[0]:(i+1)*image_resize_shape[0], j*image_resize_shape[1]:(j+1)*image_resize_shape[1]] = img[i][j] * np.ones(image_resize_shape)
		cv2.imwrite('tmp/'+str(counter)+'_'+ext+'.jpg', newimg)
	elif channel == 3:
		tile = Weights.shape[0] / channel
		i = 0
		temp = (Weights.T[:, tile*i:(i+1)*tile], Weights.T[:, (i+1)*tile:(i+2)*tile], Weights.T[:, (i+2)*tile:tile*(i+3)])
		img = tile_raster_images(X=temp, img_shape=img_shape, tile_shape=tile_shape, tile_spacing=(1, 1))
		newimg = cv2.resize(img, (img.shape[0] * image_resize_shape[0],img.shape[1] * image_resize_shape[1]))
		cv2.imwrite('tmp/'+str(counter)+'_'+ext+'.jpg', newimg)
	else:
		temp = []
		Weights = Weights.reshape((window_size*window_size, 64, 64))
		for k in xrange(Weights.shape[1]):
			img = tile_raster_images(X=Weights[:,k, :].T, img_shape=img_shape, tile_shape=tile_shape, tile_spacing=(1, 1))
			newimg = np.zeros((img.shape[0]*image_resize_shape[0], img.shape[1]*image_resize_shape[1]))
			for i in xrange(img.shape[0]):
				for j in xrange(img.shape[1]):
					newimg[i*image_resize_shape[0]:(i+1)*image_resize_shape[0], j*image_resize_shape[1]:(j+1)*image_resize_shape[1]] = img[i][j] * np.ones(image_resize_shape)
			temp.append(newimg)
		result = np.mean(temp, axis=0)
		cv2.imwrite('tmp/'+str(k)+'_'+str(counter)+'_'+ext+'.jpg', result)

def test(Weights, counter, ext, channel=1):
    """this is an utility that takes weights and plot there feature as image"""
    tile_shape = (10, 10)
    image_resize_shape = (2, 2)
    if channel == 1:
        img = tile_raster_images(
            X=Weights.T, img_shape=(window_size, window_size), tile_shape=tile_shape, tile_spacing=(1, 1)
        )
        newimg = np.zeros((img.shape[0] * image_resize_shape[0], img.shape[1] * image_resize_shape[1]))
        for i in xrange(img.shape[0]):
            for j in xrange(img.shape[1]):
                newimg[
                    i * image_resize_shape[0] : (i + 1) * image_resize_shape[0],
                    j * image_resize_shape[1] : (j + 1) * image_resize_shape[1],
                ] = img[i][j] * np.ones(image_resize_shape)
        cv2.imwrite("tmp/" + str(counter) + "_" + ext + ".jpg", newimg)
    else:
        tile = Weights.shape[0] / channel
        i = 0
        temp = (
            Weights.T[:, tile * i : (i + 1) * tile],
            Weights.T[:, (i + 1) * tile : (i + 2) * tile],
            Weights.T[:, (i + 2) * tile : tile * (i + 3)],
        )
        img = tile_raster_images(
            X=temp, img_shape=(window_size, window_size), tile_shape=tile_shape, tile_spacing=(1, 1)
        )
        newimg = cv2.resize(img, (img.shape[0] * image_resize_shape[0], img.shape[1] * image_resize_shape[1]))
        cv2.imwrite("tmp/" + str(counter) + "_" + ext + ".jpg", newimg)

def Find_cifa_10():
    """
    balabala
    """
    which_layer = 2
    
    ''' -------------- '''
    sam_file = open('SubSet1000.pkl', 'r')
    samples = cPickle.load( sam_file )
    sam_file.close()
    print 'Dimension of Samples', np.shape(samples)
    Net = DeConvNet()
    
    #kernel_list = [ 2,23,60,12,45,9 ]    
    kernel_list = [0,5,10,15]
    
    Heaps = findmaxactivation( Net, samples, 9, kernel_list, which_layer=which_layer)
    bigbigmap = None # what is this?
    for kernel_index in Heaps:        
        print 'dealing with',kernel_index,'th kernel'
        heap = Heaps[kernel_index]
        this_sams = []
        this_Deconv = []
        for pairs in heap:
            this_sam = pairs.sam
            this_sams.append( this_sam.reshape([3,32,32]) )
            this_Deconv.append( Net.DeConv( this_sam, kernel_index, which_layer=which_layer ).reshape([3,32,32]) )
        
        this_sams = np.array( this_sams )
        this_sams = np.transpose( this_sams, [ 1, 0, 2, 3 ])
        this_sams = this_sams.reshape( [ 3, 9, 32*32 ])
        this_sams = tuple( [ this_sams[i] for i in xrange(3)] + [None] )    
        
        this_Deconv = np.array( this_Deconv )
        this_Deconv = np.transpose( this_Deconv, [ 1, 0, 2, 3 ])
        this_Deconv = this_Deconv.reshape( [ 3, 9, 32*32 ])
        this_Deconv = tuple( [ this_Deconv[i] for i in xrange(3)] + [None] )

        this_map = tile_raster_images( this_sams, img_shape = (32,32), tile_shape = (3,3), 
                                   tile_spacing=(1, 1), scale_rows_to_unit_interval=True, 
                                    output_pixel_vals=True)
        this_Deconv = tile_raster_images( this_Deconv, img_shape = (32,32), tile_shape = (3,3), 
                                   tile_spacing=(1, 1), scale_rows_to_unit_interval=True, 
                                    output_pixel_vals=True)
        this_pairmap = np.append( this_map, this_Deconv, axis = 0)

        if bigbigmap == None:
            bigbigmap = this_pairmap
            segment_line = 255*np.ones([bigbigmap.shape[0],1,4],dtype='uint8')
        else:
            bigbigmap = np.append(bigbigmap, segment_line, axis = 1)            
            bigbigmap = np.append(bigbigmap, this_pairmap, axis = 1)
            
            
    plt.imshow(bigbigmap)
    plt.show()

def Find_plankton(model_name="plankton_conv_visualize_model.pkl.params"):
    """
    Find plankton that activates the given layers most
    """
    which_layer = 2
    
    import plankton_vis1
    samples = plankton_vis1.loadSamplePlanktons(numSamples=3000)
    print 'Dimension of Samples', np.shape(samples)
    Net = DeConvNet(model_name)
    
    #kernel_list = [ 2,23,60,12,45,9 ]
    kernel_list = range(48,64)
    
    
    num_of_maximum = 9
    Heaps = findmaxactivation( Net, samples, num_of_maximum, kernel_list, which_layer=which_layer)
    bigbigmap = None # what is this?
    for kernel_index in Heaps:        
        print 'dealing with',kernel_index,'th kernel'
        heap = Heaps[kernel_index]
        this_sams = []
        this_Deconv = []
        for pairs in heap:
            this_sam = pairs.sam
            this_sams.append( this_sam.reshape([NUM_C,MAX_PIXEL,MAX_PIXEL]) )
            this_Deconv.append( Net.DeConv( this_sam, kernel_index, which_layer=which_layer ).reshape([NUM_C,MAX_PIXEL,MAX_PIXEL]) )
        
        this_sams = np.array( this_sams )
        this_sams = np.transpose( this_sams, [ 1, 0, 2, 3 ])
        this_sams = this_sams.reshape( [ NUM_C, 9, MAX_PIXEL*MAX_PIXEL ])
        this_sams = tuple( [ this_sams[0] for i in xrange(3)] + [None] )    
        
        this_Deconv = np.array( this_Deconv )
        this_Deconv = np.transpose( this_Deconv, [ 1, 0, 2, 3 ])
        this_Deconv = this_Deconv.reshape( [ NUM_C, num_of_maximum, MAX_PIXEL*MAX_PIXEL ])
        this_Deconv = tuple( [ this_Deconv[0] for i in xrange(3)] + [None] )

        this_map = tile_raster_images( this_sams, img_shape = (MAX_PIXEL,MAX_PIXEL), tile_shape = (3,3), 
                                   tile_spacing=(1, 1), scale_rows_to_unit_interval=True, 
                                    output_pixel_vals=True)
        this_Deconv = tile_raster_images( this_Deconv, img_shape = (MAX_PIXEL,MAX_PIXEL), tile_shape = (3,3), 
                                   tile_spacing=(1, 1), scale_rows_to_unit_interval=True, 
                                    output_pixel_vals=True)
        this_pairmap = np.append( this_map, this_Deconv, axis = 0)

        if bigbigmap == None:
            bigbigmap = this_pairmap
            segment_line = 255*np.ones([bigbigmap.shape[0],1,4],dtype='uint8')
        else:
            bigbigmap = np.append(bigbigmap, segment_line, axis = 1)            
            bigbigmap = np.append(bigbigmap, this_pairmap, axis = 1)
            
            
    plt.imshow(bigbigmap)
    plt.show()

def main():
    if not os.path.exists('rbm_vis'):
        os.makedirs('rbm_vis')

    batch_size = 20
    train_set, test_set, val_set = load_data()

    test_x, test_y = test_set
    test_idx = np.where(test_y == 2)[0]
    test_sub_set = test_x[test_idx]

    train_x, train_y = train_set
    train_idx = np.where(train_y == 2)[0]
    train_sub_set = train_x[train_idx]
    rbm = RBM(
          n_visible= 28 * 28,
          n_hidden=100,
          input=train_sub_set,
          lr=0.5,
          b_size=batch_size
          )
    epoches = 20
    for i in range(epoches):

        w, b, a = rbm.compute_updates()
        img = Image.fromarray(utils.tile_raster_images(w.T,
                img_shape=(28, 28),
                tile_shape=(10, 10),
                tile_spacing=(1, 1)))
        img.save('rbm_vis/filter_eapoch_{}.png'.format(i))

        # test phase
        random_idx = numpy.random.choice(range(len(test_sub_set)))
        random_smpl = test_sub_set[random_idx]
        corr_random_smpl = copy.deepcopy(random_smpl)
        corr_random_smpl[:392] = 0

        # reshape
        test_itm = test_sub_set[random_idx].reshape((1,-1))
        visible_p = corr_random_smpl.reshape((1,-1))
        v2_sample = rbm.reconstruct_visible([corr_random_smpl])
        before_after = np.concatenate((test_itm, visible_p, v2_sample), axis=0)
        reconstruct = Image.fromarray(utils.tile_raster_images(before_after,
                img_shape=(28, 28),
                tile_shape=(1, 3),
                tile_spacing=(1, 1)))
        reconstruct.save('rbm_vis/reconst_eapoch_{}.png'.format(i))
        print('finished eapoch {}'.format(i))

    def Plot_RF(self, network_Q=None, layer=0, filenum=''):
        if network_Q != None:
            Q = network_Q[layer].get_value()
            filenum = str(filenum)
            function = ''
        else:
            Q = self.network.Q[layer]
            function = self.network.parameters.function
        im_size, num_dict = Q.shape

        side = int(np.round(np.sqrt(im_size)))
        im_rows = int(np.sqrt(num_dict))
        if im_rows**2 < num_dict:
            im_cols = im_rows+1
        else:
            im_cols = im_rows
        OC = num_dict/im_size

        img = tile_raster_images(Q.T, img_shape=(side, side),
                                 tile_shape=(im_rows, im_cols), tile_spacing=(1, 1),
                                 scale_rows_to_unit_interval=True, output_pixel_vals=True)
        fig = plt.figure()
        plt.title('Receptive Fields' + filenum)
        plt.axis('off')
        plt.imsave(self.directory + '/Images/RFs/Receptive_Fields'+function+filenum+'.png', img, cmap=plt.cm.Greys)
        plt.close(fig)

    def visualise_weight(self, rbm, image_name):
            assert rbm.associative
            if rbm.v_n in [784, 625, 2500, 5000]:
                plotting_start = time.clock()  # Measure plotting time

                w = rbm.W.get_value(borrow=True).T
                u = rbm.U.get_value(borrow=True).T

                weight = np.hstack((w, u))

                tile_shape = (rbm.h_n / 10 + 1, 10)

                image = Image.fromarray(
                    utils.tile_raster_images(
                        X=weight,
                        img_shape=(self.img_shape[0] *2, self.img_shape[1]),
                        tile_shape=tile_shape,
                        tile_spacing=(1, 1)
                    )
                )
                image.save(image_name)

                plotting_end = time.clock()
                return plotting_end - plotting_start
            return 0

    def visualize_filters(self):
        W = self.W.eval()
        patchSize = np.sqrt(W.shape[0])

        return tile_raster_images(X=W.T, img_shape=(patchSize, patchSize), tile_shape=(20,20), tile_spacing=(0, 0),
                       scale_rows_to_unit_interval=True,
                       output_pixel_vals=True)

def train_da(da, index, x, train_set_x, fig = 'filters_corruption_30.png', corruption_level = 0.3, learning_rate = 0.1, training_epochs = 15, batch_size=20):

    cost, updates = da.get_cost_updates(corruption_level = corruption_level, learning_rate = learning_rate)
    train_da = theano.function([index], cost, updates = updates,
                givens = {x : train_set_x[index * batch_size : (index + 1) * batch_size]})

    n_train_batches = train_set_x.get_value(borrow=True).shape[0] / batch_size

    start_time = timeit.default_timer()

    for epoch in xrange(training_epochs):
        c = []
        for batch_index in xrange(n_train_batches):
            c.append(train_da(batch_index))
        print 'train epoch %d cost' %epoch, numpy.mean(c)

    end_time = timeit.default_timer()
    print 'train time %.2f m' %((end_time - start_time) / 60)

    img_size = numpy.ceil( numpy.sqrt( da.W.get_value(borrow = True).shape[0]))

    image = Image.fromarray(tile_raster_images(
        X=da.W.get_value(borrow=True).T,
        img_shape=(img_size, img_size), tile_shape=(10, 10),
        tile_spacing=(1, 1)))
    image.save(location + fig)

    return (end_time - start_time)

def test_cA(learning_rate=0.01, training_epochs=20,
            dataset='../datasets/mnist.pkl.gz',
            batch_size=10, output_folder='cA_plots', contraction_level=.1):

    datasets = load_data(dataset)
    
    train_set_x, train_set_y = datasets[0]
    
    n_train_batches = train_set_x.get_value(borrow=True).shape[0] 
    n_train_batches /= batch_size

    index = T.lscalar() 
    x = T.matrix('x')

    if not os.path.isdir(output_folder):
        os.makedirs(output_folder)
    os.chdir(output_folder)

    rng = numpy.random.RandomState(123)

    ca = cA(numpy_rng=rng, input=x,
            n_visible=28 * 28, n_hidden=500, n_batchsize=batch_size)

    cost, updates = ca.get_cost_updates(contraction_level=contraction_level,
                                        learning_rate=learning_rate)

    train_ca = theano.function(
        [index],
        [T.mean(ca.L_rec), ca.L_jacob],
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size]
        }
    )

    start_time = timeit.default_timer()

    for epoch in xrange(training_epochs):
        c = []
        for batch_index in xrange(n_train_batches):
            c.append(train_ca(batch_index))

        c_array = numpy.vstack(c)
        print 'Training epoch %d, reconstruction cost ' % epoch, numpy.mean(
            c_array[0]), ' jacobian norm ', numpy.mean(numpy.sqrt(c_array[1]))

    end_time = timeit.default_timer()

    training_time = (end_time - start_time)

    print >> sys.stderr, ('The code for file ' + os.path.split(__file__)[1] +
                          ' ran for %.2fm' % ((training_time) / 60.))
    image = Image.fromarray(tile_raster_images(
        X=ca.W.get_value(borrow=True).T,
        img_shape=(28, 28), tile_shape=(10, 10),
        tile_spacing=(1, 1)))

    image.save('cae_filters.png')

    os.chdir('../')

def sample_DBN():
    persistent_vis_chain = theano.shared(
        numpy.asarray(test_set_x.get_value(borrow=True)[test_idx : test_idx + n_chains], dtype=theano.config.floatX)
    )

    plot_every = 1000
    ([presig_hids, hid_mfs, hid_samples, presig_vis, vis_mfs, vis_samples], updates) = theano.scan(
        rbm.gibbs_vhv, outputs_info=[None, None, None, None, None, persistent_vis_chain], n_steps=plot_every
    )

    updates.update({persistent_vis_chain: vis_samples[-1]})
    sample_fn = theano.function([], [vis_mfs[-1], vis_samples[-1]], updates=updates, name="sample_fn")

    image_data = numpy.zeros((5 * n_samples + 1, 5 * n_chains - 1), dtype="uint8")
    for idx in xrange(n_samples):

        vis_mf, vis_sample = sample_fn()
        print " ... plotting sample ", idx
        image_data[5 * idx : 5 * idx + 4, :] = tile_raster_images(
            X=vis_mf, img_shape=(4, 4), tile_shape=(1, n_chains), tile_spacing=(1, 1)
        )

    image = Image.fromarray(image_data)
    image.save("samples.png")
    os.chdir("../")

def visualize(best_params, best_samples, alg_name, tile_shape_sample, filter_shape=(10,10), tile_shape_J = (10,10), spacing = (1,1)):
    J = best_params
    J_inv = linalg.pinv(J) 
    n_sample = best_samples.shape[0]
    save_name = '-' + 'nsamples' + str(n_sample) + '-'+ alg_name + '.png'
    receptive_field = utils.tile_raster_images(J.T, filter_shape,tile_shape_J, spacing)
    image_rf = Image.fromarray(receptive_field)
    rf_name = 'J' + save_name
    image_rf.save(rf_name)     
    receptive_field_inv = utils.tile_raster_images(J_inv, filter_shape,tile_shape_J, spacing)
    image1 = Image.fromarray(receptive_field_inv)  
    rf_inv_name = 'J_inv' + save_name
    image1.save(rf_inv_name)     
    samples_vis = utils.tile_raster_images(best_samples, filter_shape,tile_shape_sample, spacing)  
    samples_vis_image = Image.fromarray(samples_vis)
    samples_vis_image.save('representative samples.png')

def train_dA(learning_rate=0.1, training_epochs=250, batch_size=30):
    d = load_data_chunk()
    train_x, train_y = d[0]
    
    n_train_batches = train_x.get_value(borrow=True).shape[0] / batch_size
    index = T.lscalar() # index to a [mini]batch
    x = T.matrix('x')  # the data is presented as rasterized images    

    rng = numpy.random.RandomState(123)
    theano_rng = RandomStreams(rng.randint(2 ** 30))    
    image_size = train_x.get_value(borrow=True).shape[1]
    da = dA(numpy_rng=rng, theano_rng=theano_rng, input=x,
            n_visible=image_size, n_hidden=500)

    cost, updates = da.get_cost_updates(corruption_level=0.3,
                                        learning_rate=learning_rate)

    train_da = theano.function([index], cost, updates=updates,
         givens={x: train_x[index * batch_size:(index + 1) * batch_size]})

    for epoch in xrange(training_epochs):
        c = []
        for batch_index in xrange(n_train_batches):
            c.append(train_da(batch_index))
        print 'Training epoch %d, cost ' % epoch, numpy.mean(c)        
        
    image = PIL.Image.fromarray(tile_raster_images(
        X=da.W.get_value(borrow=True).T,
        img_shape=(IMAGE_WIDTH, IMAGE_HEIGHT), tile_shape=(10, 10),
        tile_spacing=(1, 1)))
    image.save('filters_corruption_30.png')

def test_dA(learning_rate=0.1, training_epochs=15,
			dataset='./data/mnist.pkl.gz',
			batch_size=20, output_folder='dA_plots'):
	
	datasets = load_data(dataset)
	train_set_x, train_set_y = datasets[0]

	n_train_batches = train_set_x.get_value(borrow=True).shape[0] / batch_size

	index = T.lscalar()
	x = T.matrix('x')

	if not os.path.isdir(output_folder):
		os.makedirs(output_folder)
	os.chdir(output_folder)

	rng = numpy.random.RandomState(123)
	theano_rng = RandomStreams(rng.randint(2 ** 30))

	da = dA(numpy_rng=rng, theano_rng=theano_rng, input=x,
			n_visible=28*28, n_hidden=500)

	cost, updates = da.get_cost_updates(corruption_level=0.,
										learning_rate=learning_rate)

	train_da = theano.function([index], cost, updates=updates,
		givens={x: train_set_x[index * batch_size:
							   (index + 1) * batch_size]})

	start_time = time.clock();

	# training
	for epoch in xrange(training_epochs):
		c = []
		for batch_index in xrange(n_train_batches):
			c.append(train_da(batch_index))

		print 'Training epoch %d, cost ' % epoch, numpy.mean(c)

	end_time = time.clock()

	training_time = (end_time - start_time)

	print >> sys.stderr, ('The no corruption code for file ' +
						  os.path.split(__file__)[1] + 
						  ' ran for %.2fm' % ((training_time) / 60.))

	image = PIL.Image.fromarray(
		tile_raster_images(X=da.W.get_value(borrow=True).T,
						   img_shape=(28, 28), tile_shape=(10, 10),
						   tile_spacing=(1, 1)))

	image.save('filters_corruption_0.png')

	# training with corruption_level is 30% ......

	os.chdir('../')

 def saveImage(self, fname):
     z = self.W.T if type(self.W) == np.ndarray else self.W.eval().T
     a = int(np.sqrt(self.inputShape))
     b = int(np.sqrt(self.outputShape))
     image = PIL.Image.fromarray(tile_raster_images(
             X=z,
             img_shape=(a, a), tile_shape=(b, b),
             tile_spacing=(1, 1)))
     image.save(fname)

def plot(data, Urs, Ua, dwhite):
    '''
    plot the pictures results.
    '''
    image = Image.fromarray(
        tile_raster_images(
            X=data[:100],
            img_shape=(28, 28),
            tile_shape=(10, 10),
            tile_spacing=(1, 1)
        )
    )
    image.save('original.png')

    image = Image.fromarray(
        tile_raster_images(
            X=dwhite[:100],
            img_shape=(28, 28),
            tile_shape=(10, 10),
            tile_spacing=(1, 1)
        )
    )
    image.save('whitened.png')

    for i in range(4):
        zimage = Image.fromarray(
            tile_raster_images(
                X=Urs[i][:100],
                img_shape=(28, 28),
                tile_shape=(10, 10),
                tile_spacing=(1, 1)
            )
        )
        zimage.save('reduced_k%i.png' % i)

        uimage = Image.fromarray(
            tile_raster_images(
                X=Ua[i][:100],
                img_shape=(28, 28),
                tile_shape=(10, 10),
                tile_spacing=(1, 1)
            )
        )
        uimage.save('reconstructed_k%i.png' % i)

def test_autoencoder():
    learning_rate = 0.1
    training_epochs = 30
    batch_size = 20

    datasets = load_data('data/mnist.pkl.gz')

    train_set_x = datasets[0][0]

    # ミニバッチの数(教師データをbatch数で割るだけ)
    n_train_batches = train_set_x.get_value(borrow=True).shape[0] / batch_size

    # ミニバッチのindexシンボル
    index = T.lscalar()

    # ミニバッチの学習データシンボル
    x = T.matrix('x')

    rng = np.random.RandomState(123)
    theano_rng = RandomStreams(rng.randint(2 ** 30))

    # autoencoder モデル
    da = dA(numpy_rng=rng, theano_rng=theano_rng, input=x, n_visible=28*28, n_hidden=500)

    # コスト関数と更新式のシンボル
    cost, updates = da.get_cost_updates(corruption_level=0.0, learning_rate=learning_rate)

    # trainingの関数
    train_da = theano.function([index], cost, updates=updates, givens={
            x : train_set_x[index*batch_size : (index+1)*batch_size]
        })

    fp = open("log/ae_cost.txt", "w")

    # training
    start_time = time.clock()
    for epoch in xrange(training_epochs):
        c = []
        for batch_index in xrange(n_train_batches):
            c.append(train_da(batch_index))
        print 'Training epoch %d, cost ' % epoch, np.mean(c)
        fp.write('%d\t%f\n' % (epoch, np.mean(c)))

    end_time = time.clock()

    training_time = (end_time - start_time)

    fp.close()

    print "The no corruption code for file " + os.path.split(__file__)[1] + " ran for %.2fm" % ((training_time / 60.0))
    
    image = Image.fromarray(tile_raster_images(
    X=da.W.get_value(borrow=True).T,
    img_shape=(28, 28), tile_shape=(10, 10),
    tile_spacing=(1, 1)))
    image.save('log/dae_filters_corruption_00.png')