nmfa.py

# -*- coding: utf-8 -*-
import numpy
import pylab

dstr = """\
----------\
---####---\
--######--\
-##----##-\
-##----##-\
-########-\
-########-\
-##----##-\
-##----##-\
----------\
"""

def nmfa3_iter(s, J, b, h, t):
    dot = numpy.dot
    tanh = numpy.tanh
    for i in xrange(1,s.shape[1]):
        s[:,i] = tanh( 0.5 * b * dot(dot(J, s[:,i-1]), s[:,i-1]) + h * t )

def xi2xxx(xi):
    N = xi.shape[0]
    xxx = numpy.zeros((N,N,N))
    for i in xrange(N):
        for j in xrange(N):
            for k in xrange(N):
                xxx[i,j,k] = xi[i]*xi[j]*xi[k]
    return xxx

def add_noise(x, mean_scale, var):
    noise = numpy.random.randn(*x.shape) * var
    return mean_scale * x + noise

def plot_lcurve(s, xi):
    (N, times) = s.shape
    xit = numpy.repeat([xi], times, axis=0).transpose()
    tm = numpy.arange(times, dtype=int )
    Mo = numpy.mean(xit * s, 0)
    Er = numpy.mean((xit - s)**2, 0)
    #Df = numpy.mean(numpy.abs(s[:,1:])-numpy.abs(s[:,:-1]), 0)
    pylab.clf()
    pylab.plot(tm, Mo, '-o')
    pylab.plot(tm, Er, '-o')
    #pylab.plot(tm[1:], Df)
    #pylab.axis([None,None,-0.1,1.1])
    pylab.legend(('Mo', 'Error'), loc='center right')
    pylab.show()

def plot_code(xi, t, s, i0=None, i1=None, width=10):
    pylab.clf()
    if i1 is None:
        i1 = s.shape[1] - 1
    if i0 is None:
        i0 = i1/2

    pylab.subplot(2,2,1)
    pylab.imshow(xi.reshape(width,width), interpolation='nearest')
    pylab.colorbar()
    pylab.title('(a) Sending message')
    pylab.xticks([])
    pylab.yticks([])
    
    subplot(2,2,2)
    pylab.imshow(t.reshape(width,width), interpolation='nearest')
    pylab.colorbar()
    pylab.title('(b) Message with noise')
    pylab.xticks([])
    pylab.yticks([])
    
    pylab.subplot(2,2,3)
    pylab.imshow(s[:,i0].reshape(width,width), interpolation='nearest')
    pylab.colorbar()
    pylab.title('(c) Restored message (t=%d)' % i0)
    pylab.xticks([])
    pylab.yticks([])
    
    pylab.subplot(2,2,4)
    pylab.imshow(s[:,i1].reshape(width,width), interpolation='nearest')
    pylab.colorbar()
    pylab.title('(d) Restored message (t=%d)' % i1)
    pylab.xticks([])
    pylab.yticks([])

def s2a(s, m='-'):
    a = numpy.ones( len(s), dtype=float )
    for (i,si) in enumerate(s):
        if si == m:
            a[i] = -1
    return a

def a2s(a, p='#', m='-'):
    l = []
    for ai in a:
        if ai > 0:
            l.append(p)
        else:
            l.append(m)
    return "".join(l)

def print_s(s, width=10):
    slen = len(s)
    for i in range(slen/width):
        print s[i*width:(i+1)*width]
    
# print_s(a2s(s2a(dstr)))

def print_shape(**dct):
    for (k,v) in dct.iteritems():
        print "%s: " % k, v.shape

def main():
    width = 10
    N = width**2
    times = 21

    j0 = 1.0
    Jv = 1.0
    t0 = 1.0
    tv = 1.0
    
    xxx_scale = j0*6.0/N**2
    xxx_var = Jv**2*6.0/N**2/2.0
    xi_scale = t0
    xi_var = tv

    b = j0/Jv**2/2.0
    h = t0/tv**2

    xi = s2a(dstr)
    xxx = xi2xxx(xi)
    
    # Communication Channel
    J = add_noise(xxx, xxx_scale, xxx_var)
    t = add_noise(xi, xi_scale, xi_var)

    # NMFA
    s = numpy.zeros((N,times))
    s[:,0] = numpy.random.rand(N)
    #print_shape(xi=xi,xxx=xxx,J=J,t=t,s=s)
    nmfa3_iter(s, J, b, h, t)

    # result
    print_s(a2s(s[:,-1]))
    pylab.figure(1)
    plot_lcurve(s, xi)
    pylab.figure(2)
    plot_code(xi, t, s, i0=2)


if __name__ == '__main__':
    main()