--- a/pylearn/algorithms/mcRBM.py	Tue Aug 24 16:00:08 2010 -0400
+++ b/pylearn/algorithms/mcRBM.py	Tue Aug 24 16:51:53 2010 -0400
@@ -206,6 +206,7 @@
 
 #TODO: This should be in the datasets folder
 import pylearn.datasets.config
+import pylearn.dataset_ops.image_patches
 from pylearn.dataset_ops.protocol import TensorFnDataset
 from pylearn.dataset_ops.memo import memo
 import pylearn
@@ -223,15 +224,6 @@
     updates = [(p, p - plr * gp) for (plr, p, gp) in zip(lr, params, grads)]
     return updates
 
-@memo
-def load_mcRBM_demo_patches():
-    d = scipy.io.loadmat(os.path.join(pylearn.datasets.config.data_root(),'image_patches', 'mcRBM', 'training_colorpatches_16x16_demo.mat'))
-    totnumcases = d["whitendata"].shape[0]
-    #d = d["whitendata"][0:np.floor(totnumcases/batch_size)*batch_size,:].copy() 
-    d = d["whitendata"].copy()
-    return d
-
-
 def hidden_cov_units_preactivation_given_v(rbm, v, small=0.5):
     """Return argument to the sigmoid that would give mean of covariance hid units
 
@@ -398,15 +390,13 @@
 
 if __name__ == '__main__':
 
-    print >> sys.stderr, "TODO: use P matrix (aka FH matrix)"
-
     dataset='MAR'
     if dataset == 'MAR':
-        R,C= 21,5
+        n_vis=105
         n_patches=10240
-        demodata = scipy.io.loadmat(os.path.join(pylearn.datasets.config.data_root(),'image_patches', 'mcRBM', 'training_colorpatches_16x16_demo.mat'))
     else:
         R,C= 16,16 # the size of image patches
+        n_vis=R*C
         n_patches=100000
 
     n_train_iters=5000
@@ -423,53 +413,29 @@
     s_lr = TT.scalar()
     s_l1_penalty=TT.scalar()
     n_K=256
-    n_F=256
     n_J=100
 
-    rbm = MeanCovRBM.new_from_dims(n_I=R*C, n_K=n_K, n_J=n_J) 
+    rbm = MeanCovRBM.new_from_dims(n_I=n_vis, n_K=n_K, n_J=n_J) 
 
-    sampler = sampler(rbm, n_particles=batchsize)
+    smplr = sampler(rbm, n_particles=batchsize)
 
     def l2(X):
         return numpy.sqrt((X**2).sum())
-    def tile(X, fname):
-        if dataset == 'MAR':
-            X = np.dot(X, demodata['invpcatransf'].T)
-            R=16
-            C=16
-            #X = X.reshape((X.shape[0], 3, 16, 16)).transpose([0,2,3,1]).copy()
-            X = (X[:,:256], X[:,256:512], X[:,512:], None)
-        _img = image_tiling.tile_raster_images(X,
-                img_shape=(R,C),
-                min_dynamic_range=1e-2)
-        image_tiling.save_tiled_raster_images(_img, fname)
-    #print "Burning in..."
-    #for burnin in xrange(n_burnin_steps):
-        #sampler.simulate()
-
-    if 0:
-        print "Just SAMPLING..."
-        for jj in xrange(n_burnin_steps):
-            if 0 == jj % 100:
-                tile(sampler.positions[0].value, "sampler_%06i.png"%jj)
-                tile(numpy.random.randn(100, 105), "random_%06i.png"%jj)
-                print "burning in... ", jj
-                sys.stdout.flush()
-            sampler.simulate()
-
-        sys.exit()
+    if dataset == 'MAR':
+        tile = pylearn.dataset_ops.image_patches.save_filters_of_ranzato_hinton_2010
+    else:
+        def tile(X, fname):
+            _img = image_tiling.tile_raster_images(X,
+                    img_shape=(R,C),
+                    min_dynamic_range=1e-2)
+            image_tiling.save_tiled_raster_images(_img, fname)
 
     batch_idx = TT.iscalar()
 
     if dataset == 'MAR':
-        op = TensorFnDataset(floatX,
-                bcast=(False,),
-                fn=load_mcRBM_demo_patches,
-                single_shape=(105,))
-        train_batch = op((batch_idx * batchsize + np.arange(batchsize))%n_patches)
+        train_batch = pylearn.dataset_ops.image_patches.ranzato_hinton_2010_op(batch_idx * batchsize + np.arange(batchsize))
     else:
-        from pylearn.dataset_ops import image_patches
-        train_batch = image_patches.image_patches(
+        train_batch = pylearn.dataset_ops.image_patches.image_patches(
                 s_idx = (batch_idx * batchsize + np.arange(batchsize)),
                 dims = (n_patches,R,C),
                 center=True,
@@ -481,14 +447,13 @@
 
     grads = contrastive_gradient(rbm,
             pos_v=train_batch, 
-            neg_v=sampler.positions[0],
+            neg_v=smplr.positions[0],
             U_l1_penalty=s_l1_penalty,
             W_l1_penalty=s_l1_penalty)
     sgd_ups = sgd_updates(
                 rbm.params,
                 grads,
                 lr=[2*s_lr, .2*s_lr, .02*s_lr, .1*s_lr, .02*s_lr ])
-
     learn_fn = function([batch_idx, s_lr, s_l1_penalty], 
             outputs=[ 
                 grads[0].norm(2),
@@ -496,7 +461,6 @@
                 (sgd_ups[1][1] - sgd_ups[1][0]).norm(2),
                 ],
             updates = sgd_ups)
-    theano.printing.pydotprint(function([batch_idx, s_l1_penalty], grads[0]), 'grads0.png')
 
     print "Learning..."
     normVF=1
@@ -512,7 +476,7 @@
 
         if print_jj:
             tile(imgs_fn(jj), "imgs_%06i.png"%jj)
-            tile(sampler.positions[0].value, "sample_%06i.png"%jj)
+            tile(smplr.positions[0].value, "sample_%06i.png"%jj)
             tile(rbm.U.value.T, "U_%06i.png"%jj)
             tile(rbm.W.value.T, "W_%06i.png"%jj)
 
@@ -527,12 +491,12 @@
             print 'b min max', rbm.b.value.min(), rbm.b.value.max(),
             print 'c min max', rbm.c.value.min(), rbm.c.value.max()
 
-            print 'parts min', sampler.positions[0].value.min(), 
-            print 'max',sampler.positions[0].value.max(),
-            print 'HMC step', sampler.stepsize,
-            print 'arate', sampler.avg_acceptance_rate
+            print 'parts min', smplr.positions[0].value.min(), 
+            print 'max',smplr.positions[0].value.max(),
+            print 'HMC step', smplr.stepsize,
+            print 'arate', smplr.avg_acceptance_rate
 
-        sampler.simulate()
+        smplr.simulate()
 
         l2_of_Ugrad = learn_fn(jj, 
                 lr/max(1, jj/(20*epoch_size/batchsize)),
@@ -557,7 +521,7 @@
         # But the matrix itself is re-scaled to have an arbitrary abslute size.
         U = rbm.U.value
         U_norms = np.sqrt((U*U).sum(axis=0))
-        assert len(U_norms) == n_F
+        assert len(U_norms) == n_K
         normVF = .95 * normVF + .05 * np.mean(U_norms)
         rbm.U.value = rbm.U.value * normVF/U_norms
 

--- a/pylearn/dataset_ops/image_patches.py	Tue Aug 24 16:00:08 2010 -0400
+++ b/pylearn/dataset_ops/image_patches.py	Tue Aug 24 16:51:53 2010 -0400
@@ -2,12 +2,16 @@
 import theano
 
 from pylearn.datasets.image_patches import  (
+        data_root,
         olshausen_field_1996_whitened_images,
         extract_random_patches)
 
 from .protocol import TensorFnDataset # protocol.py __init__.py
 from .memo import memo
 
+import scipy.io
+from pylearn.io import image_tiling
+
 @memo
 def get_dataset(N,R,C,dtype,center,unitvar):
     seed=98234
@@ -48,3 +52,60 @@
 
     return x
 
+
+
+@memo
+def ranzato_hinton_2010(path=None):
+    if path is None:
+        path = os.path.join(data_root(), 'image_patches', 'mcRBM',
+                'training_colorpatches_16x16_demo.mat')
+    dct = scipy.io.loadmat(path)
+    return dct
+def ranzato_hinton_2010_whitened_patches(path=None):
+    """Return the pca of the data, which is 10240 x 105
+    """
+    dct = ranzato_hinton_2010(path)
+    return dct['whitendata'].astype('float32')
+
+def undo_pca_filters_of_ranzato_hinton_2010(X, path=None):
+    """Return tuple (R,G,B,None) of matrices for matrix `X` of filters (one per row)
+    
+    Return value can be passed to `image_tiling.tile_raster_images`.
+    """
+    dct = ranzato_hinton_2010(path)
+    X = numpy.dot(X, dct['invpcatransf'].T)
+    return (X[:,:256], X[:,256:512], X[:,512:], None)
+
+def save_filters_of_ranzato_hinton_2010(X, fname, min_dynamic_range=1e-3, data_path=None):
+    _img = image_tiling.tile_raster_images(
+            undo_pca_filters_of_ranzato_hinton_2010(X, path=data_path),
+            img_shape=(16,16),
+            min_dynamic_range=min_dynamic_range)
+    image_tiling.save_tiled_raster_images(_img, fname)
+
+def ranzato_hinton_2010_op(s_idx,
+        split='train', 
+        dtype=theano.config.floatX, rasterized=True,
+        center=True,
+        unitvar=True):
+    N = 10240
+
+    if split != 'train':
+        raise NotImplementedError('train/test/valid splits for randomly sampled image patches?')
+
+    if not rasterized:
+        # the data is provided as PCA-sphered, so rasterizing does not make sense
+        # TODO: add a param to enable/disable 'PCA', and if disabled, then consider
+        # rasterizing or not
+        raise NotImplementedError('only pca data is provided')
+
+    if dtype != 'float32':
+        raise NotImplementedError('dtype not float32')
+
+    op = TensorFnDataset(dtype,
+            bcast=(False,), 
+            fn=ranzato_hinton_2010_whitened_patches,
+            single_shape=(105,))
+    x = op(s_idx%N)
+    return x
+

--- a/pylearn/datasets/image_patches.py	Tue Aug 24 16:00:08 2010 -0400
+++ b/pylearn/datasets/image_patches.py	Tue Aug 24 16:51:53 2010 -0400
@@ -111,3 +111,4 @@
 
     return rval
 
+