http://takatakamanbou.hatenablog.com/entry/2015/07/18/233000 http://takatakamanbou.hatenablog.com/entry/2015/07/19/233000

ex150718-1.py

import numpy as np

import nnet150718 as nnet
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def genDat3D( N, sig ):

    xy = np.random.random_sample( ( N, 2 ) )
    z = 0.5 * xy[:, 0] + 0.5 * xy[:, 1]
    Xorg   = np.hstack( ( xy, z[:, np.newaxis] ) )
    Xnoisy = np.copy( Xorg )
    Xnoisy[:, 2] += sig * np.random.standard_normal( N )

    return Xorg, Xnoisy
    

def linearAE( D, H ):

    L1 = nnet.Layer( D, H, 'linear', withBias = False, Wini = 0.01 )
    L2 = nnet.Layer( H, D, 'linear', withBias = False, Wini = 0.01 )
    mlp = nnet.MLP( [ L1, L2 ], cost = 'squared_error' )

    return mlp


if __name__ == '__main__':

    D, H = 3, 3
    ae = linearAE( D, H )
    N = 1000
    Xorg, Xnoisy = genDat3D( N, 0.1 )
    Xm = np.mean( Xnoisy, axis = 0 )
    Xorg -= Xm
    Xnoisy -= Xm

    #Xin, Zteacher = Xorg, Xorg  # conventional AE
    #Xin, Zteacher = Xnoisy, Xnoisy  # conventional AE
    Xin, Zteacher = Xnoisy, Xorg    # denoising AE

    eta, mu, lam = 0.01, 0.9, 0.0
    nepoch = 5000
    print '# eta = ', eta, ' mu = ', mu, ' lam = ', lam
    print '### training:   N = ', N, ' D = ', D, ' H = ', H

    i = 0
    tmp, Z = ae.output( Xin )
    msqe = np.mean( ae.cost( Z, Zteacher ) )
    print i, msqe
    
    for i in range( nepoch ):
        msqe = np.mean( ae.train( Xin, Zteacher, eta, mu, lam ) )
        print i+1, msqe

    W = ae.Layers[0].getParams()

    nx = ny = nz = 8
    Nd = nx * ny * nz
    Xgrid = np.empty( ( Nd, 3 ) )
    i = 0
    for iz, z in enumerate( np.arange( 0, 1, 1.0/nz ) ):
        for iy, y in enumerate( np.arange( 0, 1, 1.0/ny ) ):
            for ix, x in enumerate( np.arange( 0, 1, 1.0/nx ) ):
                Xgrid[i, :] = [ x, y, z ]
                i += 1
    Xgrid -= Xm


    #X = Xorg
    #X = Xnoisy
    X = Xgrid
    tmp, Z = ae.output( X )

    fig = plt.figure()
    ax = fig.add_subplot( 111, projection = '3d' )
    ax.set_xlim( -1, 1 )
    ax.set_ylim( -1, 1 )
    ax.set_zlim( -1, 1 )
    ax.set_xlabel( 'x' )
    ax.set_ylabel( 'y' )
    ax.set_zlabel( 'z' )
    
    ax.scatter( X[:, 0], X[:, 1], X[:, 2], color = 'red' )
    ax.scatter( Z[:, 0], Z[:, 1], Z[:, 2], color = 'blue' )
    fig.show()

ex150718-2.py

import numpy as np
import cv2

import mnist0118 as mnist
import nnet150718 as nnet


def genDatMNIST_flip( N, p ):

    mn = mnist.MNIST( 'L' )
    Xorg = mn.getImage()[:N] / 256
    D = Xorg.shape[1]
    Xnoisy = np.copy( Xorg )
    idx = np.random.random_sample( N * D ) < p
    Xnoisy2 = Xnoisy.reshape( -1 )
    Xnoisy2[idx] = 1.0 - Xnoisy2[idx] 

    return Xorg, Xnoisy


def linearAE( D, H ):

    L1 = nnet.Layer( D, H, 'linear', withBias = False, Wini = 0.01 )
    L2 = nnet.Layer( H, D, 'linear', withBias = False, Wini = 0.01 )
    mlp = nnet.MLP( [ L1, L2 ], cost = 'squared_error' )

    return mlp


def makeMNISTImage( X, nx, ny, nrow = 28, ncol = 28, gap = 4 ):

    X2 = X[:nx*ny, :].reshape( ( nx*ny, 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] = X2[iy*nx+ix, :, :]

    return img
    

if __name__ == '__main__':

    np.random.seed( 0 )
    
    NL = NT = 2000
    XorgLT, XnoisyLT = genDatMNIST_flip( NL + NT, 0.1 )
    XorgL = XorgLT[:NL, :]
    XnoisyL = XnoisyLT[:NL, :]
    XorgT = XorgLT[NL:, :]
    XnoisyT = XnoisyLT[NL:, :]
    D = XorgL.shape[1]

    
    H = 1000
    ae = linearAE( D, H )
    Xm = np.mean( XorgL, axis = 0 )
    XorgL -= Xm
    XnoisyL -= Xm
    XorgT -= Xm
    XnoisyT -= Xm

    #Xin, Zteacher = XorgL, XorgL  # conventional AE
    #Xin, Zteacher = XnoisyL, XnoisyL  # conventional AE
    Xin, Zteacher = XnoisyL, XorgL    # denoising AE

    eta, mu, lam = 0.01, 0.9, 0.0
    nepoch = 5000
    print '# eta = ', eta, ' mu = ', mu, ' lam = ', lam
    print '### training:   NL = ', NL, ' D = ', D, ' H = ', H

    i = 0
    tmp, Z = ae.output( Xin )
    msqe = np.mean( ae.cost( Z, Zteacher ) )
    print i, msqe
    
    for i in range( nepoch ):
        msqe = np.mean( ae.train( Xin, Zteacher, eta, mu, lam ) )
        if i % 100 == 99:
            print i+1, msqe

    W = ae.Layers[0].getParams()

    img = makeMNISTImage( ( XorgT + Xm ) * 256, 10, 5 )
    cv2.imwrite( 'hogeXorg.png', img )

    img = makeMNISTImage( ( XnoisyT + Xm ) * 256, 10, 5 )
    cv2.imwrite( 'hogeXnoisy.png', img )

    tmp, Z = ae.output( XnoisyT )
    #tmp, Z = ae.output( XorgT )
    img = makeMNISTImage( ( Z + Xm ) * 256, 10, 5 )
    cv2.imwrite( 'hogeZ.png', img )

nnet150718.py

import numpy as np
import theano
import theano.tensor as T


# activation functions
d_afunc = { 'linear': lambda Y: Y,
            'sigmoid': T.nnet.sigmoid,
            'softmax': T.nnet.softmax,
            'ReLu': lambda Y: T.switch( Y > 0, Y, 0 ) }

### uniform random numbers for weight initialization
#
def randomU( shape, a ):

    # [ -a, a )
    return 2 * a * ( np.random.random_sample( shape ) - 0.5 ) 


### Gaussian random numbers for weight initialization
#
def randomN( shape, sig ):

    # N(0,sig)
    return sig * np.random.standard_normal( shape )



########## Layer ##########

class Layer( object ):

    def __init__( self, Din, Nunit, afunc, withBias = True, Wini = 0.01, floatX = theano.config.floatX ):

        self.Din   = Din
        self.Nunit = Nunit
        self.afunc = d_afunc[afunc]
        self.withBias = withBias
        
        # theano shared variables for weights & biases
        self.W = theano.shared( np.array( randomN( ( Nunit, Din ), Wini ), dtype = floatX ) )
        self.dW = theano.shared( np.zeros( ( Nunit, Din ), dtype = floatX ) )
        if withBias:
            self.b = theano.shared( np.zeros( Nunit, dtype = floatX ) )
            #self.b = theano.shared( np.ones( Nunit, dtype = floatX ) )
            self.db = theano.shared( np.zeros( Nunit, dtype = floatX ) )


    def output( self, X ):

        if self.withBias:
            Y = T.dot( X, self.W.T ) + self.b  # Ndat x Nunit
        else:
            Y = T.dot( X, self.W.T )

        Z = self.afunc( Y )

        return Y, Z


    def getParams( self ):

        if self.withBias:
            return [ self.W.get_value(), self.b.get_value() ]
        else:
            return self.W.get_value()
            
        

########## MLP ##########

class MLP( object ):

    def __init__( self, Layers, cost = 'cross_entropy' ):

        # layers - list of Layer instances
        self.Layers = Layers

        # theano functions
        self.output = self._Tfunc_output()
        self.cost   = self._Tfunc_cost( cost )
        self.train  = self._Tfunc_train( cost )


    ### theano function for output computation
    #
    def _Tfunc_output( self ):

        X = T.matrix()  # N x D
        Y, Z = _T_output( self.Layers, X )

        return theano.function( [ X ], [ Y, Z ] )


    ### theano function for cost computation (cross-entropy)
    #
    def _Tfunc_cost( self, costfunc ):

        if costfunc == 'cross_entropy':
            Z = T.matrix()     # N x K
            lab = T.ivector()  # N-dim
            teacher = lab
            cost    = _T_cost_ce( Z, lab )
        else:
            Z  = T.matrix()  # N x Dout
            ZT = T.matrix()  # N x Dout
            teacher = ZT
            cost    = _T_cost_sqe( Z, ZT )

        return theano.function( [ Z, teacher ], cost )



    ### theano function for gradient descent learning
    #
    def _Tfunc_train( self, costfunc ):

        X = T.matrix( 'X' )  # N x D
        if costfunc == 'cross_entropy':
            Y, Z = _T_output( self.Layers, X )
            lab  = T.ivector( 'lab' )  # N-dim
            teacher = lab
            cost    = T.mean( _T_cost_ce( Z, lab ) )
        else:
            Y, Z = _T_output( self.Layers, X )
            ZT   = T.matrix( 'ZT' )  # N x Dout
            teacher = ZT
            cost    = T.mean( _T_cost_sqe( Z, ZT ) )

        eta  = T.scalar( 'eta' )
        mu   = T.scalar( 'mu' )
        lam  = T.scalar( 'lambda' )

        updatesList = []
        for layer in self.Layers:
            gradW = T.grad( cost, layer.W )
            #dWnew = -eta * gradW + mu * layer.dW
            dWnew = -eta * ( gradW + lam * layer.W ) + mu * layer.dW
            Wnew  = layer.W + dWnew
            updatesList.append( ( layer.W, Wnew ) )
            updatesList.append( ( layer.dW, dWnew ) )
            if layer.withBias:
                gradb = T.grad( cost, layer.b )
                # no weight decay for bias
                dbnew = -eta * gradb + mu * layer.db
                bnew  = layer.b + dbnew
                updatesList.append( ( layer.b, bnew ) )
                updatesList.append( ( layer.db, dbnew ) )

        return theano.function( [ X, teacher, eta, mu, lam ], cost, updates = updatesList )


def _T_output( Layers, X ):

    Zprev = X
    for layer in Layers:
        Y, Z = layer.output( Zprev )
        Zprev = Z

    return Y, Z


def _T_cost_ce( Z, lab ):

    return T.nnet.categorical_crossentropy( Z, lab )

def _T_cost_sqe( Z, ZT ):

    return T.sum( T.sqr( Z - ZT ), axis = 1 )