In [1]:
from __future__ import (division, absolute_import, print_function, unicode_literals)
import sys, os, bz2, re, itertools, pickle
from copy import copy
import numpy as np; import scipy.linalg as linalg
import tom

In [2]:
# Setup the comment %%# magic
from IPython.core.magic import register_line_cell_magic
ip = get_ipython()
def comment_cell_magic(line, cell):
    "Place a comment and execute the cell"
    get_ipython().run_cell(cell)
ip.register_magic_function(comment_cell_magic, 'cell', '#')
del comment_cell_magic, ip

In [3]:
# Setup inline plotting for PhD plots (use `plt.savefig(filename, dpi=600)`)
%matplotlib inline
%config InlineBackend.print_figure_kwargs = {'bbox_inches':None}
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
fig_width = 412.56496/72.27; fig_height = fig_width * (np.sqrt(5)-1.0)/2.0; fig_size =  [fig_width,fig_height]
fig_params = {'savefig.dpi': 144, 'savefig.pad_inches': 0.0, 'savefig.format': 'pdf',
              'figure.figsize': fig_size,
              'text.usetex': True,
              'font.family': 'serif', 'font.serif': 'Computer Modern Roman', 'mathtext.fontset': 'cm',
              'lines.scale_dashes': False,
              'font.size': 9, 'axes.titlesize': 9,'axes.labelsize': 9, 'legend.fontsize': 9, 'xtick.labelsize': 8, 'ytick.labelsize': 8,
              'xtick.top': True, 'xtick.direction': 'in', 'ytick.right': True, 'ytick.direction': 'in',
              'axes.linewidth': 0.5, 'xtick.major.width': 0.4, 'xtick.major.size': 2, 'ytick.major.size': 2, 'ytick.major.width': 0.4, 'xtick.minor.size': 0, 'ytick.minor.size': 0,
              'xtick.major.pad': 3, 'ytick.major.pad': 2, 'axes.labelpad': 4,
              'legend.borderpad': 0.3, 'legend.labelspacing': 0.3,
              'figure.subplot.bottom': 0.1
             }
mpl.rcParams.update(fig_params)
import matplotlib.pyplot as plt
import brewer2mpl
bmap = brewer2mpl.get_map('Paired', 'qualitative', 10)
bmap1 = brewer2mpl.get_map('set1', 'qualitative', 9)
colors = {'red': bmap.mpl_colors[5], 'lightred': bmap.mpl_colors[4],
          'blue': bmap.mpl_colors[1], 'lightblue': bmap.mpl_colors[0],
          'green': bmap.mpl_colors[3], 'lightgreen': bmap.mpl_colors[2],
          'purple': bmap.mpl_colors[9], 'lightpurple': bmap.mpl_colors[8],
          'orange': bmap.mpl_colors[7], 'lightorange': bmap.mpl_colors[6],
          'yellow': bmap1.mpl_colors[5], 'brown': bmap1.mpl_colors[6],
          'pink': bmap1.mpl_colors[7], 'grey': bmap1.mpl_colors[8], 'black': 'black'}
colorNos = ['blue', 'red', 'green', 'purple', 'orange', 'pink', 'brown', 'yellow', 'grey', 'lightblue', 'lightred', 'lightgreen', 'lightpurple']

In [ ]:
# Custom line styles
def myStyle(lw = 1, color='black', ls='-'):
    lw = 0.8 * lw
    s = 0.5
    if color in colors:
        color = colors[color]
    elif type(color) is int:
        color = colors[colorNos[color]]
    dash_capstyle = 'butt'
    dashes = (None, None)
    if ls == ':' or ls == 1:
        dash_capstyle='round'
        dashes = ((1+lw) * 0.1, (2+lw)*s)
    elif ls == '--' or ls == 2:
        dashes = (4*s, 2*s)
    elif ls == ';' or ls == 3:
        dashes = (1*s, 2*s, 4*s, 2*s)
    elif ls == '---' or ls == 4:
        dashes = (8*s, 1.5*s)
    return {'lw': lw, 'dashes': dashes, 'dash_capstyle': dash_capstyle, 'color': color}
# Hack for special legend symbol
from matplotlib.legend_handler import HandlerTuple
class HandlerTupleVertical(HandlerTuple):
    def __init__(self, **kwargs):
        HandlerTuple.__init__(self, **kwargs)
    def create_artists(self, legend, orig_handle,
                       xdescent, ydescent, width, height, fontsize, trans):
        # How many lines are there.
        numlines = len(orig_handle)
        handler_map = legend.get_legend_handler_map()
        # divide the vertical space where the lines will go
        # into equal parts based on the number of lines
        height_y = (height / numlines)
        leglines = []
        for i, handle in enumerate(orig_handle):
            handler = legend.get_legend_handler(handler_map, handle)
            legline = handler.create_artists(legend, handle, xdescent, (2*i + 1)*height_y, width, 2*height,
                                            fontsize, trans)
            leglines.extend(legline)
        return leglines
# Legend style
def makeLegend(lines=None, labels=None, title=None, row=4):
    ax = plt.subplot(row,3,3*row)
    ax.tick_params(left='off', right='off', top='off', bottom='off')
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    for direction, lw in zip(['top', 'left', 'bottom', 'right'], BenchData.plot_info['LEGEND']['lw']):
        ax.spines[direction].set_linewidth(lw)
    if row == 1: ax.spines['top'].set_linewidth(0)
    drawLegend(lines, labels, title)
def drawLegend(lines=None, labels=None, title=None, title_size=10, loc=10):
    legend = plt.legend(lines, labels, loc=loc, fancybox=True, title=title,
                        handler_map = {tuple: HandlerTupleVertical(), list: HandlerTupleVertical()})
    legend.get_title().set_fontsize(str(title_size)) 
    frame = legend.get_frame()
    frame.set_facecolor('0.9')
    frame.set_edgecolor('none')
    return legend
def finishPlot():
    plt.tight_layout(pad=0.1)
    plt.subplots_adjust(hspace=0, wspace=0)
def savePlot(filename):
    os.makedirs(research_directory + '/results', exist_ok=True)
    plt.savefig(research_directory + '/results/' + filename + '.pdf', dpi=600, bbox_inches='tight', pad_inches=0)

In [ ]:
def make_missingness_hmm(SPs):
    N = len(SPs) + 1
    hmm = tom.hmm.Hmm(N, 2)
    T = np.zeros((N,N))
    for i, sp in enumerate(SPs):
        T[i, 0] = 1 - sp
        T[i, i+1] = sp
    T[N-1, 0] = 1
    hmm.T(T)
    hmm.E(0, np.array([[0], *((N-1) * [[1]])], dtype=float))    
    hmm.E(1, np.array([[1], *((N-1) * [[0]])], dtype=float))
    hmm.pi(np.array([[1], *((N-1) * [[0]])], dtype=float))
    hmm.init()
    return hmm
def make_misingness_oom(p_missing):
    if p_missing <= 0:
        oom = tom.Oom(1, 2)
        oom.tau(0, np.array([[1]], dtype=float))
        oom.tau(1, np.array([[0]], dtype=float))
        oom.initialize()
        return oom
    elif p_missing >= 1:
        oom = tom.Oom(1, 2)
        oom.tau(0, np.array([[0]], dtype=float))
        oom.tau(1, np.array([[1]], dtype=float))
        oom.initialize()
        return oom
    N = int(np.ceil(1/p_missing))
    SPs = (N-1) * [1]
    SPs[-1] = 1/p_missing - N + 1
    hmm = make_missingness_hmm(SPs)
    oom = tom.Oom(hmm)
    oom.w0(oom.stationaryState())
    return oom

In [4]:
o = 1
class BenchData:
    info = {
        'ECOLI':       {'nO':  4, 'nU':  0, 'dim':  0, 'maxTrainLen': 4588690, 'nTrain': 1, 'nTest': 5, 'ylim': (1.93, 2)},
        'BIBLE27':     {'nO': 27, 'nU':  0, 'dim':  0, 'maxTrainLen': 3846968, 'nTrain': 1, 'nTest': 5, 'ylim': (1.9, 3.5)},
        'RANDOM4_7':   {'nO':  4, 'nU':  0, 'dim':  7, 'maxTrainLen': 10**7, 'lenX': 5, 'lenY': 5,
                        'entropies': [1.54532, 1.43766, 1.31455, 1.16933, 0.95743], 'split': int(10**3), 'ylim': (0.95743, 1.54532)},
        'RANDOM4_32':  {'nO':  4, 'nU':  0, 'dim': 32, 'maxTrainLen': 10**7, 'lenX': 5, 'lenY': 5,
                        'entropies': [1.92694, 1.88487, 1.71605, 1.51314, 1.08093], 'split': int(10**3.5), 'ylim': (1.08093, 1.92694)},
        'RANDOM27_7':  {'nO': 27, 'nU':  0, 'dim':  7, 'maxTrainLen': 10**7, 'lenX': 2, 'lenY': 2,
                        'entropies': [4.23171, 3.77634, 3.47086, 3.32650, 3.21578], 'split': int(10**4.5), 'ylim': (3.21578, 4.23171)},
        'RANDOM27_32': {'nO': 27, 'nU':  0, 'dim': 32, 'maxTrainLen': 10**8, 'lenX': 2, 'lenY': 2, 
                        'entropies': [4.66523, 4.53882, 4.32753, 4.11195, 3.71746], 'split': 10**6, 'ylim': (3.71746, 4.66523)},
        'TIGER':       {'nO':  2, 'nU':  3, 'dim':  2, 'maxTrainLen': 10**7, 'lenX': 3, 'lenY': 3},
        'PAINT':       {'nO':  2, 'nU':  4, 'dim':  2, 'maxTrainLen': 10**7, 'lenX': 3, 'lenY': 3},
        'BRIDGE':      {'nO':  5, 'nU': 12, 'dim':  5, 'maxTrainLen': 10**7, 'lenX': 1, 'lenY': 1},
        'NETWORK':     {'nO':  2, 'nU':  4, 'dim':  7, 'maxTrainLen': 10**7, 'lenX': 3, 'lenY': 3},
        'SHUTTLE':     {'nO':  5, 'nU':  3, 'dim':  7, 'maxTrainLen': 10**7, 'lenX': 2, 'lenY': 2},
        'MAZE4X3':     {'nO':  6, 'nU':  4, 'dim': 10, 'maxTrainLen': 10**7, 'lenX': 2, 'lenY': 2},
        'CHEESE':      {'nO':  7, 'nU':  4, 'dim': 11, 'maxTrainLen': 10**7, 'lenX': 2, 'lenY': 2}}
#    plot_info = {
#        'RANDOM4_7':   {'pos':1, 'lw': [1, 1, 0.5, 0.5]},
#        'RANDOM4_32':  {'pos':2, 'lw': [1, 0.5, 0.5, 1]},
#        'RANDOM27_7':  {'pos':4, 'lw': [0.5, 1, 1, 0.5]},
#        'RANDOM27_32': {'pos':5, 'lw': [0.5, 0.5, 1,1]},
#        'TIGER':       {'pos':3, 'lw': [1,1,0.5,1]},
#        'PAINT':       {'pos':6, 'lw': [0.5,1,0.5,1]},
#        'BRIDGE':      {'pos':7, 'lw': [1,1,0.5,0.5]},
#        'NETWORK':     {'pos':8, 'lw': [1,0.5,0.5,0.5]},
#        'SHUTTLE':     {'pos':9, 'lw': [0.5,0.5,0.5,1]},
#        'MAZE4X3':     {'pos':10, 'lw': [0.5,1,1,0.5]},
#        'CHEESE':      {'pos':11, 'lw': [0.5,0.5,1,0.5]},
#        'LEGEND':      {'pos':12, 'lw': [1,1,0,0]}
#    }
    plot_info = {
        'RANDOM4_7':   {'pos':1, 'lw': [o, o, 0.5, 0.5]},
        'RANDOM4_32':  {'pos':2, 'lw': [o, 0.5, 0.5, 1]},
        'RANDOM27_7':  {'pos':4, 'lw': [0.5, o, 1, 0.5]},
        'RANDOM27_32': {'pos':5, 'lw': [0.5, 0.5, 1,1]},
        'TIGER':       {'pos':3, 'lw': [o,1,0.5,o]},
        'PAINT':       {'pos':6, 'lw': [0.5,1,0.5,o]},
        'BRIDGE':      {'pos':7, 'lw': [1,o,0.5,0.5]},
        'NETWORK':     {'pos':8, 'lw': [1,0.5,0.5,0.5]},
        'SHUTTLE':     {'pos':9, 'lw': [0.5,0.5,0.5,o]},
        'MAZE4X3':     {'pos':10, 'lw': [0.5,o,o,0.5]},
        'CHEESE':      {'pos':11, 'lw': [0.5,0.5,o,0.5]},
        'LEGEND':      {'pos':12, 'lw': [0.5,0.5,0,0]},
        'ECOLI':       {'pos':1,  'lw': [o,o,o,0.5]},
        'BIBLE27':     {'pos':2,  'lw': [o,0.5,o,0.5]}       
    }    
    realWorldData = ['ECOLI', 'BIBLE27']
    OOMs = ['RANDOM4_7', 'RANDOM4_32', 'RANDOM27_7', 'RANDOM27_32']
    IOOOMs = ['TIGER', 'PAINT', 'BRIDGE', 'NETWORK', 'SHUTTLE', 'MAZE4X3', 'CHEESE']
    colors = {'SPEC':'grey', 'Standard':'grey', 'RCW':'green', 'RowColWeighted':'green','WLS':'blue', 'GLS':'purple', 'ES': 'red', 'EM': 'yellow'}
    def initStatic(masterSeed = 123456789, nTrain = 3, nTest = 5):
        np.random.seed(masterSeed)
        BenchData.seeds = {}
        for oomName in BenchData.OOMs + BenchData.IOOOMs:
            BenchData.info[oomName]['nTrain'] = nTrain
            BenchData.info[oomName]['nTest'] = nTest
            BenchData.seeds[oomName] = {'train': [np.random.randint(0, 2**31-1) for t in range(nTrain)],
                                        'test':  [np.random.randint(0, 2**31-1) for t in range(nTest)]}
            BenchData.plot_info[oomName]['xlim'] = (10**3, 10**7)
            BenchData.plot_info[oomName]['xscale'] = 'log'
        for oomName in BenchData.OOMs:
            BenchData.plot_info[oomName]['yscale'] = 'linear'
            BenchData.plot_info[oomName]['ylim'] = (BenchData.info[oomName]['entropies'][4], BenchData.info[oomName]['entropies'][0])
        BenchData.plot_info['RANDOM4_7']['xlim'] = (10**2, 10**6)
        BenchData.plot_info['RANDOM27_32']['xlim'] = (10**4, 10**8)
        for oomName in BenchData.realWorldData:
            BenchData.plot_info[oomName]['xlim'] = (10**3, BenchData.info[oomName]['maxTrainLen'])
            BenchData.plot_info[oomName]['xscale'] = 'log'
            BenchData.plot_info[oomName]['yscale'] = 'linear'
            BenchData.plot_info[oomName]['ylim'] = BenchData.info[oomName]['ylim']
        for oomName in BenchData.IOOOMs:
            for polType in ['real', 'rand']:
                BenchData.seeds[oomName + polType] = {'train': [np.random.randint(0, 2**31-1) for t in range(nTrain)],
                                                      'test':  [np.random.randint(0, 2**31-1) for t in range(nTest)]}
            BenchData.plot_info[oomName]['yscale'] = 'log'
            BenchData.plot_info[oomName]['ylim'] = (1e-6, 1e-1)
        for oomName in BenchData.OOMs + BenchData.realWorldData:
            BenchData.seeds['m' + oomName] = [np.random.randint(0, 2**31-1) for t in range(nTrain)]
        np.random.seed()
        
    def __init__(self, oomName, policy = ''):
        self.oomName = oomName
        self.oomNameType = self.oomName + ('-real' if policy == 'real' else '')
        self.info = BenchData.info[oomName]
        self.plot_info = BenchData.plot_info[oomName]
        self._oom = None
        self._policy = ('real' if policy == 'real' else '')
        self.tests_ = self.nTest() * [None]
        
    def initPlot(self, rows = 4, allOOMs = False):
        pos = self.plot_info['pos']
        if self.isRealWorld():
            if not allOOMs: rows = 1
            else: pos *= 3
        ax = plt.subplot(rows, 3, pos)
        plotName = 'RANDOM%d\_%d' % (self.nO(), self.dim()) if self.oomName in BenchData.OOMs else self.oomName
        ax.text(0.978, 0.978, plotName, fontsize=8, weight="bold", horizontalalignment='right', verticalalignment='top',
                transform=ax.transAxes, bbox={'facecolor':(0.9,0.9,0.9), 'edgecolor':'none', 'boxstyle':'round'})
        ax.set_xscale(self.plot_info['xscale'])
        ax.set_yscale(self.plot_info['yscale'])
        ax.set_xlim(self.plot_info['xlim'])
        if self.oomName in BenchData.OOMs + BenchData.realWorldData: ax.set_yticks(np.arange(0, 5, 0.1))
        ax.set_ylim(self.plot_info['ylim'])
        ax.set_yticklabels([])
        ax.set_xticklabels([])
        ax.grid(**myStyle(lw=0.2,color='brown', ls=1))
        if self.isRealWorld() and allOOMs:
            if self.oomName == 'ECOLI': plot_info_lw = [o,1,0.5,o]
            if self.oomName == 'BIBLE27': plot_info_lw = [0.5,1,o,o]
        else:
            plot_info_lw = self.plot_info['lw']
        if allOOMs: plot_info_lw[0] = 0.5
        for i, direction in enumerate(['top', 'left', 'bottom', 'right']):
            ax.spines[direction].set_linewidth(plot_info_lw[i])
        return ax

    def nO(self):
        return BenchData.info[self.oomName]['nO']
    def nU(self):
        return BenchData.info[self.oomName]['nU']
    def dim(self):
        return BenchData.info[self.oomName]['dim']
    def o_min(self, trainLength):
        Σₒ, Σᵢ = self.nO(), max(1, self.nU())
        Σ = Σₒ * Σᵢ
        l = np.log(trainLength) / (np.log(Σ) + np.log(Σᵢ))
        return (l+1) * Σᵢ**l
    def maxTrainLen(self):
        return BenchData.info[self.oomName]['maxTrainLen']
    def trainLengths(self, minExp2 = None, maxExp2 = None):
        if maxExp2 == None:
            if self.oomName == 'RANDOM27_32': maxExp2 = 16
            else: maxExp2 = 14
        if minExp2 == None:
            if self.oomName == 'RANDOM4_7': minExp2 = 4
            else: minExp2 = 6            
        tl = [int(10**(i/2)) for i in range(minExp2, maxExp2+1) if int(10**(i/2)) <= self.maxTrainLen()]
        if tl[-1] < self.maxTrainLen() and self.maxTrainLen() < int(10**(maxExp2/2)):
            if self.maxTrainLen() / tl[-1] < 1.5:  tl.pop()
            tl.append(self.maxTrainLen())
        return tl
    def nTrain(self):
        return BenchData.info[self.oomName]['nTrain']
    def nTest(self):
        return BenchData.info[self.oomName]['nTest']
    def tests(self):
        for i in range(self.nTest()):
            if self.tests_[i] is None:
                self.tests_[i] = self.getSequence(i, use='test')
        return self.tests_
    def isRealWorld(self):
        return self.oomName in BenchData.realWorldData
    @property
    def oom(self):
        if self._oom is None:
            if self.oomName in BenchData.IOOOMs:
                self._oom = tom.load(research_directory + '/benchmarks/' + self.oomName.lower() + '.pomdp.json')
                self._oom = tom.Oom(self._oom)
            elif self.oomName in BenchData.OOMs:
                self._oom = tom.load(research_directory + '/benchmarks/' + self.oomName + '.oom')
        return self._oom
    @property
    def policy(self):
        if self._policy == 'real':
            self._policy = tom.load(research_directory + '/benchmarks/' + self.oomName + '.policy')
            self._policy.nU_ = self.nU(); self._policy.exploration_ = 1/3
        elif self._policy == '':
            self._policy = tom.hmm.Policy()
        return self._policy
    def getSequence(self, seedId = 0, use = 'train', length = None, missingProb = None, regular = False, nOextra = 0, nU = 1):
        if self.oom is None:
            seq = tom.load(research_directory + '/benchmarks/' + self.oomNameType + '_' + use + str(seedId) + '.seq.bz2')
            if length is None: length = seq.length()
        else:
            if length is None: length = self.maxTrainLen() if use == 'train' else 10**5
            self.oom.reset()
            seed = self.seeds[self.oomName][use][seedId]
            seq = self.oom.sample(length, tom.Random(seed), self.policy)
        if missingProb is None:
            return seq
        if self.nU() != 0: raise RuntimeError('Missing values require output-only sequences')
        mRand = tom.Random(self.seeds['m' + self.oomName][seedId])
        mSeq = tom.Sequence(seq.length(), seq.nOutputSymbols() + nOextra, nU)
        if not regular:
            for t in range(seq.length()):
                if mRand.random() < missingProb:
                    mSeq.u(t, 1); mSeq.o(t, seq.nOutputSymbols())
                else:
                    mSeq.u(t, 0); mSeq.o(t, seq.o(t))
        else:
            missing = make_misingness_oom(missingProb).sample(length, mRand)
            for t in range(seq.length()):
                if missing[t] == 1:
                    mSeq.u(t, 1); mSeq.o(t, seq.nOutputSymbols())
                else:
                    mSeq.u(t, 0); mSeq.o(t, seq.o(t))
        return mSeq
    def evaluate(self, learntOom, method = None, stabilization = None, tl = None):
        if type(learntOom) is list:
            return [self.evaluate(loom, method, stabilization, tl) for loom in learntOom]
        if method is None: method = 'aospe' if learntOom.isIO() else 'l2l'
        if stabilization is None:
            #learntOom.stabilization(0.0001, 0.1, 0, 1e-8)
            if method == 'l2l':
                learntOom.stabilization(0.0002, 0.03, 5, 1e-8)
            else:
                learntOom.stabilization(0, 1, 0, 1e-8)
        else:
            learntOom.stabilization(*stabilization)
        if method == 'l2l':
            return np.average([learntOom.l2l(test) for test in self.tests()])
        else:
            return np.average([learntOom.averageOneStepPredictionError(test, self.oom) for test in self.tests()])

BenchData.initStatic()

In [ ]:
# Manipulating results
def combineResults(res):
    """If the results have the structure [{key1: ... {keyn: [results for trainLengths]}}], then this function returns
    {key1: ... {keyn: np.array([average results for trainLengths])}}."""
    if type(res[0]) is dict:
        avres = {}
        for key in res[0].keys():
            if key not in  ['oom', 'loom']:
                avres[key] = combineResults([r[key] for r in res])
        return avres
    else:
        if len(res[0]) == 0: return []
        return np.average(res, axis=0) #np.exp(np.average(np.log(res), axis=0))
    
def addResults(res, res_more):
    """If the results have the structure [{key1: ... {keyn: [results for trainLengths]}}], then this function adds the
    entries from res_more to res."""
    if type(res[0]) is dict:
        for k in res_more[0].keys():
            if k not in res[0].keys():
                for i in range(3): res[i][k] = res_more[i][k]
            else:
                addResults([r[k] for r in res], [rm[k] for rm in res_more])
    else:
        print("Warning: res_more contains results for the same settings as res.")

In [2]:
# simple execution timer functionality
import time
try:
    clock = time.process_time
except:
    clock = time.time
def tictoc(function):
    """Measure the execution time of the given function.
    
    Returns
    -------
    (t, ret) :
    t : double
        the execution time in seconds
    ret : 
        the return value of the `function`
    
    Notes
    -----
    Typically, one will use this function in the following way:
        (t, ret) = tictoc( lambda: my_function(arg1, arg2) )
    The returned time `t` will be the CPU time for python >= 3.3
    and the wall-clock time otherwise.
    """
    start = clock()
    ret = function()
    elapsed = clock() - start
    return (elapsed, ret)

class Timer:
    def __init__(self, use_wall_clock = False):
        self.clock = time.time if use_wall_clock else clock
        self.tic_ = self.clock()
    def tic(self):
        self.tic_ = self.clock()
    def printtoc(self):
        toc = self.clock() - self.tic_
        print('%2dh:%2dm:%2ds' % (int(toc / 3600), int((toc % 3600) / 60), int((toc % 60))))
    def toc(self):
        return self.clock() - self.tic_

In [1]:
def display_width(width_percent = 100):
    from IPython.core.display import display, HTML
    display(HTML("<style>.container { width:%d%% !important; }</style>" % width_percent))