http://takatakamanbou.hatenablog.com/entry/2015/08/08/233000

ex150808.py

import numpy as np
import cv2

import mnist_sub150808 as mnist_sub
import nnet150718 as nnet


def addNoise( X ):

    #Xnoisy = mnist_sub.addSaltPepperNoise( X, 0.1 )

    #Xnoisy = mnist_sub.shift( X, 2 )

    X2     = mnist_sub.shift( X, 2 )
    Xnoisy = mnist_sub.addSaltPepperNoise( X2, 0.1 )

    return Xnoisy


def setData( X, Xm, aetype ):

    if aetype == 'Org-Org':  # conventional AE
        Xin      = X - Xm
        Zteacher = Xin
    elif aetype == 'Noisy-Org': # denoising AE
        Xin      = addNoise( X ) - Xm
        Zteacher = X - Xm
    else:       # Noisy-Noisy
        Xin      = addNoise( X ) - Xm
        Zteacher = Xin

    return Xin, Zteacher



if __name__ == '__main__':

    dirMNIST = '.'
    np.random.seed( 0 )

    ##### setting the training data & the validation data
    #
    XL, labL, XV, labV = mnist_sub.getDataLV( dirMNIST = dirMNIST )
    Xm = np.mean( XL, axis = 0 )
    NL, D = XL.shape
    NV, D = XV.shape

    ##### mini batch indicies for stochastic gradient ascent
    #
    batchsize = 100    
    biL = mnist_sub.makeBatchIndex( NL, batchsize )
    nbatchL = biL.shape[0]
    biV = mnist_sub.makeBatchIndex( NV, batchsize )
    nbatchV = biV.shape[0]


    ##### training
    #
    H = 1000
    #L1 = nnet.Layer( D, H, 'linear', withBias = False, Wini = 0.01 )
    L1 = nnet.Layer( D, H, 'ReLu', withBias = True, Wini = 0.01 )
    L2 = nnet.Layer( H, D, 'linear', withBias = False, Wini = 0.01 )
    ae = nnet.MLP( [ L1, L2 ], cost = 'squared_error' )

    eta, mu, lam = 0.01, 0.9, 0.0
    nepoch = 50
    print '# eta = ', eta, ' mu = ', mu, ' lam = ', lam

    print '### training:   NL = ', NL, 'NV = ', NV, ' batchsize = ', batchsize
    print '#               eta = ', eta, 'mu = ', mu, 'lam = ', lam
    print '#               D = ', D, ' H = ', H


    for i in range( nepoch ):

        sqeL = 0.0
        for ib in np.random.permutation( nbatchL ):
            ii = biL[ib, :]
            #aetype = 'Org-Org'
            #aetype = 'Noisy-Noisy'
            aetype = 'Noisy-Org'
            Xin, Zteacher = setData( XL[ii, :], Xm, aetype = aetype )
            sqeL += ae.train( Xin, Zteacher, eta, mu, lam ) * Xin.shape[0]

        if i < 10 or i % 10 == 0:
            # validation error for Noisy-Org
            sqeV = 0.0
            for ib in range( nbatchV ):
                ii = biV[ib, :]
                Xin, Zteacher = setData( XV[ii, :], Xm, aetype = 'Noisy-Org' )
                tmp, Z = ae.output( Xin )
                sqeV += np.sum( ae.cost( Z, Zteacher ) )
            print i, sqeL / NL, sqeV / NV


    ##### setting the test data
    #
    XT, labT = mnist_sub.getDataT( dirMNIST = dirMNIST )
    NT, D = XT.shape
    biT = mnist_sub.makeBatchIndex( NT, batchsize )
    nbatchT = biT.shape[0]

    # test error for Noisy-Org
    print '### test:   NT = ', NT
    sqeT = 0.0
    for ib in range( nbatchT ):
        ii = biT[ib, :]
        Xin, Zteacher = setData( XT[ii, :], Xm, aetype = 'Noisy-Org' )
        tmp, Z = ae.output( Xin )
        sqeT += np.sum( ae.cost( Z, Zteacher ) )
    print i, sqeL / NL, sqeV / NV, sqeT / NT

    # visualizing the reconstruction for noisy test data
    Xin, Zteacher = setData( XT[:50, :], Xm, aetype = 'Noisy-Org' )
    tmp, Z = ae.output( Xin )
    img = mnist_sub.visualize( XT[:50, :] * 255, 10, 5 )
    cv2.imwrite( 'hogeXorg.png', img )
    img = mnist_sub.visualize( ( Xin + Xm ) * 255, 10, 5 )
    cv2.imwrite( 'hogeXnoisy.png', img )
    img = mnist_sub.visualize( ( Z + Xm ) * 255, 10, 5 )
    cv2.imwrite( 'hogeZ.png', img )

    # visualizing the network weights
    if ae.Layers[0].withBias:
        W0 = ae.Layers[0].getParams()[0]
    else:
        W0 = ae.Layers[0].getParams()
    absmax = np.max( np.abs( W0 ), axis = 1 )
    W0 = W0 / absmax[:, np.newaxis] * 127 + 128
    img = mnist_sub.visualize( W0[:100], 10, 10 )
    cv2.imwrite( 'hogeW0.png', img )
    if ae.Layers[1].withBias:
        W1 = ae.Layers[1].getParams()[0].T
    else:
        W1 = ae.Layers[1].getParams().T
    absmax = np.max( np.abs( W1 ), axis = 1 )
    W1 = W1 / absmax[:, np.newaxis] * 127 + 128
    img = mnist_sub.visualize( W1[:100], 10, 10 )
    cv2.imwrite( 'hogeW1.png', img )

mnist150808.py

import struct
import os
import numpy as np


class MNIST:

    def __init__( self, LT, dirMNIST = '.'  ):

        self.nclass = 10
        self.LT = LT
        
        if self.LT == 'L':
            self.fnLabel = os.path.join( dirMNIST, 'train-labels-idx1-ubyte' )
            self.fnImage = os.path.join( dirMNIST, 'train-images-idx3-ubyte' )
        else:
            self.fnLabel = os.path.join( dirMNIST, 't10k-labels-idx1-ubyte' )
            self.fnImage = os.path.join( dirMNIST, 't10k-images-idx3-ubyte' )

    def getLabel( self ):
        
        return _readLabel( self.fnLabel )

    def getImage( self ):
        
        return _readImage( self.fnImage )


##### reading the label file
#
def _readLabel( fnLabel ):

    f = open( fnLabel, 'r' )

    ### header (two 4B integers, magic number(2049) & number of items)
    #
    header = f.read( 8 )
    mn, num = struct.unpack( '>2i', header )  # MSB first (bigendian)
    assert( mn == 2049 )
    #print mn, num

    ### labels (unsigned byte)
    #
    label = np.array( struct.unpack( '>%dB' % num, f.read() ), dtype = int )

    f.close()

    return label


##### reading the image file
#
def _readImage( fnImage ):

    f = open( fnImage, 'r' )

    ### header (four 4B integers, magic number(2051), #images, #rows, and #cols
    #
    header = f.read( 16 )
    mn, num, nrow, ncol = struct.unpack( '>4i', header ) # MSB first (bigendian)
    assert( mn == 2051 )

    ### pixels (unsigned byte)
    #
    npixel = ncol * nrow
    pixel = np.empty( ( num, npixel ) )
    for i in range( num ):
        buf = struct.unpack( '>%dB' % npixel, f.read( npixel ) )
        pixel[i, :] = np.asarray( buf )

    f.close()

    return pixel



if __name__ == '__main__':

    dirMNIST = '.'

    print '# MNIST training data'
    mnist = MNIST( 'L', dirMNIST = dirMNIST )
    lab = mnist.getLabel()
    dat = mnist.getImage()
    print lab.shape, dat.shape

    print '# MNIST test data'
    mnist = MNIST( 'T', dirMNIST = dirMNIST )
    lab = mnist.getLabel()
    dat = mnist.getImage()
    print lab.shape, dat.shape
    

mnist_sub150808.py

import numpy as np
import mnist150808 as mnist


def getDataLV( dirMNIST = '.' ):

    mn = mnist.MNIST( 'L', dirMNIST = dirMNIST )
    X = np.asarray( mn.getImage() / 255, dtype = np.float32 ) # => in [0,1]
    lab = np.asarray( mn.getLabel(), dtype = np.int32 )
    XL, labL = X[:50000], lab[:50000]
    XV, labV = X[50000:], lab[50000:]

    return XL, labL, XV, labV


def getDataT( dirMNIST = '.' ):

    mn = mnist.MNIST( 'T', dirMNIST = dirMNIST )
    XT = np.asarray( mn.getImage() / 255, dtype = np.float32 ) # => in [0,1]
    labT = np.asarray( mn.getLabel(), dtype = np.int32 )

    return XT, labT


def makeBatchIndex( N, batchsize ):
    
    idx = np.random.permutation( N )
    nbatch = int( np.ceil( float( N ) / batchsize ) )
    idxB = np.zeros( ( nbatch, N ), dtype = bool )
    for ib in range( nbatch - 1 ):
        idxB[ib, idx[ib*batchsize:(ib+1)*batchsize]] = True
    ib = nbatch - 1
    idxB[ib, idx[ib*batchsize:]] = True

    return idxB


def addSaltPepperNoise( Xbatch, p ):

    Xnoisy = np.copy( Xbatch )
    pvec = np.random.random_sample( Xbatch.shape[0] * Xbatch.shape[1] )
    idxB = pvec < p/2
    idxW = ( p/2 <= pvec ) * ( pvec < p )
    Xnoisy2 = Xnoisy.reshape( -1 )
    Xnoisy2[idxB] = 0.0
    Xnoisy2[idxW] = 1.0

    return Xnoisy


def shift( Xbatch, dmax ):

    nrow = ncol = 28
    N = Xbatch.shape[0]
    delta = np.random.randint( -dmax, dmax + 1, ( N, 2 ) )
    Xnoisy = np.zeros_like( Xbatch )
    for i in range( N ):
        src = Xbatch[i, :].reshape( ( nrow, ncol ) )
        dst = np.zeros( ( nrow + 2*dmax, ncol + 2*dmax ) )
        dx, dy = delta[i, 0], delta[i, 1]
        dst[dmax+dy:dmax+dy+nrow, dmax+dx:dmax+dx+ncol] = src
        Xnoisy[i, :] = dst[dmax:dmax+nrow, dmax:dmax+ncol].reshape( -1 )

    return Xnoisy
        


def visualize( src, nx, ny, nrow = 28, ncol = 28, gap = 4 ):

    src2 = src.reshape( ( ny, nx, nrow, ncol ) )
    
    w = nx * ( ncol + gap ) + gap
    h = ny * ( nrow + gap ) + gap
    img = np.zeros( ( h, w ), dtype = int ) + 128

    for iy in range( ny ):
        lty = iy * ( nrow + gap ) + gap
        for ix in range( nx ):
            ltx = ix * ( ncol + gap ) + gap
            img[lty:lty+nrow, ltx:ltx+ncol] = src2[iy, ix, :, :]

    return img