HMM_Spice-V5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import sys\n",
    "from collections import defaultdict\n",
    "from math import expm1\n",
    "import time\n",
    "try:\n",
    "    # Python 2.7\n",
    "    import urllib2 as ur\n",
    "    orl2 = True\n",
    "except:\n",
    "    #Python 3.4\n",
    "    import urllib.request as ur\n",
    "    orl2 = False\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def computeInitialProb1(trainDataFile,numOfStates):\n",
    "    trainFile=open(trainDataFile,\"r\")\n",
    "    metaDataLine = trainFile.readline()\n",
    "    headerLine = metaDataLine.split(\" \")\n",
    "    numSequences = int(headerLine[0])\n",
    "    distinctObservations= int(headerLine[1])#Total Number of Distinct Observations\n",
    "    numOfStates=min(numOfStates,distinctObservations)\n",
    "    empiricalCount=np.zeros(shape=numOfStates)\n",
    "    empiricalFreq=defaultdict(int)\n",
    "    for n in range(numSequences):\n",
    "        line = trainFile.readline()#Reading Sequences 1 by 1\n",
    "        l = line.split(\" \")\n",
    "        startState=int(l[1])\n",
    "        empiricalFreq[startState] = empiricalFreq[startState]+1\n",
    "    totalObservations=0\n",
    "    for i in np.arange(numOfStates):\n",
    "        empiricalCount[i]=empiricalFreq[i]\n",
    "        totalObservations=totalObservations+empiricalCount[i]\n",
    "    initialProb=[count/totalObservations for count in empiricalCount]\n",
    "    return (numOfStates,distinctObservations,initialProb)\n",
    "def computeInitialProb(trainDataFile,numOfStates):\n",
    "    trainFile=open(trainDataFile,\"r\")\n",
    "    metaDataLine = trainFile.readline()\n",
    "    headerLine = metaDataLine.split(\" \")\n",
    "    numSequences = int(headerLine[0])\n",
    "    distinctObservations= int(headerLine[1])#Total Number of Distinct Observations\n",
    "    numOfStates=min(numOfStates,distinctObservations)\n",
    "    empiricalDistr=Counter()\n",
    "    for n in range(numSequences):\n",
    "        line = trainFile.readline()#Reading Sequences 1 by 1\n",
    "        line=line.rstrip(\"\\n\")\n",
    "        l = line.split(\" \")\n",
    "        l=l[1:]\n",
    "        lDistr=Counter(l)\n",
    "        empiricalDistr+=lDistr\n",
    "    totalSymbolsSeen=sum(empiricalDistr.values())\n",
    "    initialProb=[]\n",
    "    for i in np.arange(numOfStates):\n",
    "        initialProb.append((1.0*empiricalDistr[str(i)])/totalSymbolsSeen)\n",
    "    return (numOfStates,distinctObservations,initialProb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def createRandomMatrixA(numOfStates):\n",
    "    matrixA=np.ndarray(shape=(numOfStates,numOfStates),dtype=float)\n",
    "    prob=1.0/(numOfStates*numOfStates)\n",
    "    matrixA.fill(prob)\n",
    "    return matrixA\n",
    "def createRandomMatrixB(numOfStates,distinctObservations):\n",
    "    matrixB=np.ndarray(shape=(numOfStates,distinctObservations),dtype=float)\n",
    "    prob=1.0/(numOfStates*distinctObservations)\n",
    "    matrixB.fill(prob)\n",
    "    return matrixB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeAlpha(observations,a,aTranspose,b,bTranspose,pi,alphaDP,alphaDPScaleT):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    alphaDpScaleTime0=0\n",
    "    #Vec(alphaDP[0])=pi[0]*b[0][ob[0]],pi[1]*b[1][ob[0]],pi[2]*b[2][ob[0]],....,pi[n-1]*b[n-1][ob[0]]\n",
    "    #Vec(alphaDP[0])=pi[0]*bTrans[ob[0]][0],pi[1]*bTrans[ob[0]][1],pi[2]*bTrans[ob[0]][2],....,pi[n-1]*bTrans[ob[0]][n-1]\n",
    "    alphaDP[0]=pi*bTranspose[observations[0]]\n",
    "    alphaDPScaleT[0]=1.0/np.sum(alphaDP[0])\n",
    "    alphaDP[0]*=alphaDPScaleT[0]\n",
    "    for t in np.arange(1,timePts):\n",
    "        for i in np.arange(statesC):\n",
    "            alphaDP[t][i]=b[i][observations[t]]*(np.sum(alphaDP[t-1]*aTranspose[i]))\n",
    "        alphaDPScaleT[t]=1.0/np.sum(alphaDP[t])\n",
    "        alphaDP[t]*=alphaDPScaleT[t]\n",
    "def observationsLikelihood(alphaDPScaleT):\n",
    "    return 1.0/np.prod(alphaDPScaleT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeBeta(observations,a,b,bTranspose,pi,betaDP):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    betaDP[timePts-1].fill(1)\n",
    "    for t in np.arange(timePts-2,-1,-1):\n",
    "        betaDpScaleTimeT=0\n",
    "        for i in np.arange(statesC):\n",
    "            betaDP[t][i]=np.sum(a[i]*bTranspose[observations[t+1]]*betaDP[t+1])\n",
    "        betaDpScaleTimeT=1.0/np.sum(betaDP[t])\n",
    "        betaDP[t]*=betaDpScaleTimeT\n",
    "    return betaDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def validateAlphaDP(alphaDP,betaDP):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    ans=[]\n",
    "    for t in np.arange(observationsC-1):\n",
    "        ans.append(np.sum(alphaDP[t]*betaDP[t]))\n",
    "    print(\"Alpha-BETA Validation \",ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def computeGammaDP(diGammaDP):\n",
    "    return np.apply_along_axis(np.sum,2,diGammaDP)\n",
    "def computeDiGammaDP(alphaDP,alphaDPScaleT,betaDP,a,b,bTranspose,observations):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    diGammaDP=np.zeros(shape=(observationsC,statesC,statesC),dtype=float)\n",
    "    diGammaDenom=observationsLikelihood(alphaDPScaleT)\n",
    "    for i in np.arange(statesC):\n",
    "        for t in np.arange(observationsC-1):\n",
    "            diGammaDP[t][i]=alphaDP[t][i]*a[i]*bTranspose[observations[t+1]]*betaDP[t+1]\n",
    "    diGammaDP/=diGammaDenom\n",
    "    return diGammaDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeTransitionProbabilityA(diGammaDP,gammaDP):\n",
    "    diGammaIJSumMatrix=np.apply_along_axis(np.sum,0,diGammaDP)\n",
    "    gammaDPISumMatrix=np.apply_along_axis(np.sum,0,gammaDP)\n",
    "    timePts=diGammaIJSumMatrix.shape[0]\n",
    "    for i in np.arange(timePts):\n",
    "        diGammaIJSumMatrix[i]/=gammaDPISumMatrix[i]\n",
    "    return diGammaIJSumMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def computeObsrProbNum(gammaDPT,i,vk,observations):\n",
    "    gammaDPi=gammaDPT[i]\n",
    "    return np.sum(gammaDPi[np.where(observations==vk)])\n",
    "def computeTransitionProbabilityB(gammaDP,observations,observationDict):\n",
    "    statesC=gammaDP.shape[1]\n",
    "    observationsC=len(observationDict)\n",
    "    newlyComputedObsrProbB=np.zeros(shape=(statesC,observationsC),dtype=float)#Ideal Shape should be transposed\n",
    "    gammaDPISumMatrix=np.apply_along_axis(np.sum,0,gammaDP)\n",
    "    gammaDPT=gammaDP.transpose()\n",
    "    for i in np.arange(statesC):\n",
    "        for vk in observationDict:\n",
    "            newlyComputedObsrProbB[i][vk]=computeObsrProbNum(gammaDPT,i,vk,observations)/gammaDPISumMatrix[i]\n",
    "    return newlyComputedObsrProbB\n",
    "def computeTransitionProbabilityB1(alphaDP,alphaDPScaleT,betaDP,a,b,observations,observationDict):\n",
    "    statesC=a.shape[0]\n",
    "    observationsC=b.shape[1]\n",
    "    newlyComputedObsrProbB=np.zeros(shape=(observationsC,statesC),dtype=float)#Ideal Shape should be transposed\n",
    "    #gammaDP=computeGammaDP(alphaDP,betaDP,alphaDPScaleT)#[t][state]\n",
    "    gammaDP=gammaDP.transpose()\n",
    "    for i in np.arange(statesC):\n",
    "        obsrProbDenom =np.sum(gammaDP[i])\n",
    "        for vk in observationDict:\n",
    "            newlyComputedObsrProbB[i][vk]=computeObsrProbNum(gammaDP,i,vk,observations)/obsrProbDenom\n",
    "    return newlyComputedObsrProbB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Change Convergence Criteria to be more reasonable/Useful\n",
    "def isConverged(count,convergenceIters):\n",
    "    if count>=convergenceIters:\n",
    "        return True\n",
    "    return False\n",
    "def Forward_Backward_EM_Algo(observations,A,B,pi,convergenceIters,observationDict):\n",
    "    count=0\n",
    "    updatedA=A\n",
    "    updatedB=B\n",
    "    while isConverged(count,convergenceIters)==False:\n",
    "        #Expectation(E)-Step\n",
    "        alphaDP=np.zeros(shape=(observations.shape[0],updatedA.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "        betaDP=np.zeros(shape=(observations.shape[0],updatedA.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "        alphaDPScaleT=np.ones(shape=(observations.shape[0]))\n",
    "        updatedATranspose=updatedA.transpose()\n",
    "        updatedBTranspose=updatedB.transpose()\n",
    "        computeAlpha(observations,updatedA,updatedATranspose,updatedB,updatedBTranspose,pi,alphaDP,alphaDPScaleT)\n",
    "        computeBeta(observations,updatedA,updatedB,updatedBTranspose,pi,betaDP)\n",
    "        #validateAlphaDP(alphaDP,betaDP)\n",
    "        diGammaDP=computeDiGammaDP(alphaDP,alphaDPScaleT,betaDP,updatedA,updatedB,updatedBTranspose,observations)\n",
    "        gammaDP=computeGammaDP(diGammaDP)#[t][state]\n",
    "        #Maximization(M)-Step\n",
    "        newA=computeTransitionProbabilityA(diGammaDP,gammaDP)\n",
    "        #newB=computeTransitionProbabilityB(alphaDP,alphaDPScaleT,betaDP,updatedA,updatedB,observations,observationDict)\n",
    "        newB=computeTransitionProbabilityB(gammaDP,observations,observationDict)\n",
    "        updatedA=newA\n",
    "        updatedB=newB\n",
    "        count=count+1\n",
    "    return (updatedA,updatedB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trainHMM(trainDataFile,A,B,pi,convergenceIters,maxSequences=-1):\n",
    "    trainFile=open(trainDataFile,\"r\")\n",
    "    metaDataLine = trainFile.readline()\n",
    "    headerLine = metaDataLine.split(\" \")\n",
    "    numSequences = int(headerLine[0])\n",
    "    distinctObservations= int(headerLine[1])#Total Number of Distinct Observations\n",
    "    observationDict=np.arange(distinctObservations)\n",
    "    updatedA=np.NaN\n",
    "    updatedB=np.NaN\n",
    "    isAUpdated=False\n",
    "    if(maxSequences==-1):\n",
    "        usedSeqs=numSequences\n",
    "    else:\n",
    "        usedSeqs=min(maxSequences,numSequences)\n",
    "    actuallyUsedSeqs=0\n",
    "    for n in range(usedSeqs):\n",
    "        line = trainFile.readline()#Reading Sequences 1 by 1\n",
    "        line=line.rstrip(\"\\n\")\n",
    "        l = line.split(\" \")\n",
    "        if(int(l[0])<=1):\n",
    "            continue\n",
    "        actuallyUsedSeqs+=1\n",
    "        observations=np.array([int(i) for i in l[1:len(l)]])\n",
    "        learnedParams=Forward_Backward_EM_Algo(observations,A,B,pi,convergenceIters,observationDict)\n",
    "        if isAUpdated==False:\n",
    "            isAUpdated=True\n",
    "            updatedA=learnedParams[0]\n",
    "            updatedB=learnedParams[1]\n",
    "        else:\n",
    "            updatedA+=learnedParams[0]\n",
    "            updatedB+=learnedParams[1]\n",
    "    updatedA=updatedA/actuallyUsedSeqs\n",
    "    updatedB=updatedB/actuallyUsedSeqs\n",
    "    return (updatedA,updatedB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def trainModel(fileLoc,maxNoOfStates,convergenceIters,maxSequences=-1):\n",
    "    start = time.time()\n",
    "    initialProbs=computeInitialProb(fileLoc,maxNoOfStates)\n",
    "    end = time.time()\n",
    "    print(\"Computed Initial Prob. in \", end - start ,\"seconds\")\n",
    "    pi=initialProbs[2]\n",
    "    numOfStates=initialProbs[0]\n",
    "    distinctObservations=initialProbs[1]\n",
    "    A=createRandomMatrixA(numOfStates)\n",
    "    B=createRandomMatrixB(numOfStates,distinctObservations)\n",
    "    trainedParams=trainHMM(fileLoc,A,B,pi,convergenceIters,maxSequences)\n",
    "    trainedParams=trainedParams+(pi,)\n",
    "    end=time.time()\n",
    "    print(\"For \",maxSequences,\" Sequences : Total Training Time \",end-start,\" seconds\")\n",
    "    return trainedParams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Initial Prob. in  0.0057408809661865234 seconds\n"
     ]
    },
    {
     "ename": "FloatingPointError",
     "evalue": "invalid value encountered in true_divide",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFloatingPointError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-60-04ad8f83b84f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mold_settings\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mseterr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mall\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'ignore'\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m#seterr to known value\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mseterr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mall\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'raise'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;31m#{'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore','under': 'ignore'}\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mpi\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtrainModel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Data/0.spice.train1.txt'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mseterr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mold_settings\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# reset to default {'over': 'raise', 'divide': 'ignore', 'invalid': 'ignore', 'under': 'ignore'}\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;31m#(A,B,pi)=trainModel('Data/1.spice.train.txt',20,7,1)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-57-457bc1c7f9df>\u001b[0m in \u001b[0;36mtrainModel\u001b[1;34m(fileLoc, maxNoOfStates, convergenceIters, maxSequences)\u001b[0m\n\u001b[0;32m      9\u001b[0m     \u001b[0mA\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcreateRandomMatrixA\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnumOfStates\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m     \u001b[0mB\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcreateRandomMatrixB\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnumOfStates\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdistinctObservations\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 11\u001b[1;33m     \u001b[0mtrainedParams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtrainHMM\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfileLoc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mpi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mconvergenceIters\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmaxSequences\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     12\u001b[0m     \u001b[0mtrainedParams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtrainedParams\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m     \u001b[0mend\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-56-0c26b6d78bc6>\u001b[0m in \u001b[0;36mtrainHMM\u001b[1;34m(trainDataFile, A, B, pi, convergenceIters, maxSequences)\u001b[0m\n\u001b[0;32m     22\u001b[0m         \u001b[0mactuallyUsedSeqs\u001b[0m\u001b[1;33m+=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     23\u001b[0m         \u001b[0mobservations\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0ml\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ml\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 24\u001b[1;33m         \u001b[0mlearnedParams\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mForward_Backward_EM_Algo\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobservations\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mpi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mconvergenceIters\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mobservationDict\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     25\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0misAUpdated\u001b[0m\u001b[1;33m==\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     26\u001b[0m             \u001b[0misAUpdated\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-44-a32af5615e24>\u001b[0m in \u001b[0;36mForward_Backward_EM_Algo\u001b[1;34m(observations, A, B, pi, convergenceIters, observationDict)\u001b[0m\n\u001b[0;32m     21\u001b[0m         \u001b[0mgammaDP\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcomputeGammaDP\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdiGammaDP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;31m#[t][state]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     22\u001b[0m         \u001b[1;31m#Maximization(M)-Step\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 23\u001b[1;33m         \u001b[0mnewA\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcomputeTransitionProbabilityA\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdiGammaDP\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgammaDP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     24\u001b[0m         \u001b[1;31m#newB=computeTransitionProbabilityB(alphaDP,alphaDPScaleT,betaDP,updatedA,updatedB,observations,observationDict)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     25\u001b[0m         \u001b[0mnewB\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcomputeTransitionProbabilityB\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgammaDP\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mobservations\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mobservationDict\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-42-b04b24358d2d>\u001b[0m in \u001b[0;36mcomputeTransitionProbabilityA\u001b[1;34m(diGammaDP, gammaDP)\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mtimePts\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdiGammaIJSumMatrix\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimePts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m         \u001b[0mdiGammaIJSumMatrix\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m/=\u001b[0m\u001b[0mgammaDPISumMatrix\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mdiGammaIJSumMatrix\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mFloatingPointError\u001b[0m: invalid value encountered in true_divide"
     ]
    }
   ],
   "source": [
    "old_settings = np.seterr(all='ignore')  #seterr to known value\n",
    "np.seterr(all='raise')#{'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore','under': 'ignore'}\n",
    "(A,B,pi)=trainModel('Data/0.spice.train1.txt',4,4,2)\n",
    "np.seterr(**old_settings)  # reset to default {'over': 'raise', 'divide': 'ignore', 'invalid': 'ignore', 'under': 'ignore'}\n",
    "#(A,B,pi)=trainModel('Data/1.spice.train.txt',20,7,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.24748808,  0.24642955,  0.24962722,  0.25645515],\n",
       "       [ 0.24923185,  0.24287502,  0.2405963 ,  0.26729682],\n",
       "       [ 0.25084286,  0.23898292,  0.23128661,  0.27888761],\n",
       "       [ 0.24439687,  0.25169186,  0.26399802,  0.23991324]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10557188,  0.5276764 ,  0.        ,  0.36675172],\n",
       "       [ 0.10957474,  0.54557272,  0.        ,  0.34485254],\n",
       "       [ 0.12045131,  0.593053  ,  0.        ,  0.28649569],\n",
       "       [ 0.08524658,  0.43384218,  0.        ,  0.48091124]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getHmmRank(prefix,A,ATranspose,B,BTranspose,pi,uniqueSymbols):\n",
    "    likelihoods=[]\n",
    "    for i in np.arange(uniqueSymbols):\n",
    "        prefix.append(i)\n",
    "        observations=np.array(prefix)\n",
    "        alphaDP=np.zeros(shape=(observations.shape[0],A.shape[0]))# Count_of_Observations*Count_of_Hidden_States\n",
    "        alphaDPScaleT=np.ones(shape=(observations.shape[0]))\n",
    "        computeAlpha(observations,B,ATranspose,B,BTranspose,pi,alphaDP,alphaDPScaleT)\n",
    "        obsrLikelihood=observationsLikelihood(alphaDPScaleT)\n",
    "        prefix.pop()\n",
    "        likelihoods.append((i,obsrLikelihood))\n",
    "    likelihoods=sorted(likelihoods, key=lambda x: -x[1])\n",
    "    ranks=[i[0] for i in likelihoods]\n",
    "    return ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def list_to_string(l):\n",
    "    s=str(l[0])\n",
    "    for x in l[1:]:\n",
    "        s+= \" \" + str(x)\n",
    "    return(s)\n",
    "def formatString(string_in):\n",
    "    \"\"\" Replace white spaces by %20 \"\"\"\n",
    "    return string_in.strip().replace(\" \", \"%20\")\n",
    "# get the test first prefix: the only element of the test set\n",
    "def get_first_prefix(test_file):\n",
    "    \"\"\" This function is called for the public test file(Which only has 1 line)\n",
    "    \"\"\"\n",
    "    f = open(test_file)\n",
    "    prefix = f.readline()\n",
    "    f.close()\n",
    "    return prefix\n",
    "def predictOnSpicePublicData(problem_number,name):\n",
    "    problem_number = str(problem_number)\n",
    "    user_id = '68'\n",
    "    #name = \"hmm_Baseline\"\n",
    "    #train_file = 'Data/0.spice.train.txt'\n",
    "    prefix_file = 'Data/'+problem_number+'.spice.public.test.txt'\n",
    "    first_prefix = get_first_prefix(prefix_file)\n",
    "    prefix_number=1\n",
    "    # get the next symbol ranking on the first prefix\n",
    "    p=first_prefix.split()\n",
    "    prefix=[int(i) for i in p[1:len(p)]]#prefix holds the sequence of values in the public test file(Note:It has only 1 Seq)\n",
    "    print(\"Prefix \",prefix)\n",
    "    ranking=getHmmRank(prefix,A,A.transpose(),B,B.transpose(),pi,A.shape[0])\n",
    "    print(\"Model Ranking \",ranking)\n",
    "    ranking_string=list_to_string(ranking[:5])\n",
    "    #print(\"Prefix number: \" + str(prefix_number) + \" Ranking: \" + ranking_string + \" Prefix: \" + first_prefix)\n",
    "    first_prefix = formatString(first_prefix)\n",
    "\n",
    "    # transform the ranking to follow submission format\n",
    "    ranking_string=formatString(ranking_string)\n",
    "\n",
    "    # create the url to submit the ranking\n",
    "    #name=name+\"_Ver1.7.2\"\n",
    "    name=name\n",
    "    url_base = 'http://spice.lif.univ-mrs.fr/submit.php?user=' + user_id +\\\n",
    "        '&problem=' + problem_number + '&submission=' + name + '&'\n",
    "    url = url_base + 'prefix=' + first_prefix + '&prefix_number=1' + '&ranking=' +\\\n",
    "        ranking_string\n",
    "    response = ur.urlopen(url)\n",
    "    print(\"URL \",url)\n",
    "    content = response.read()\n",
    "    print(\"Response from SPiCe \",content)#Content is a new Sequence returned from the SPiCe server: We will need to predict for this seq\n",
    "    if not orl2:\n",
    "        # Needed for python 3.4...\n",
    "        content= content.decode('utf-8')\n",
    "    list_element = content.split()\n",
    "    head = str(list_element[0])\n",
    "    return content,url_base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prefix  [3, 0, 3, 0, 1, 3, 3]\n",
      "Model Ranking  [3, 0, 1, 2]\n",
      "URL  http://spice.lif.univ-mrs.fr/submit.php?user=68&problem=0&submission=hmm_baseline_v1&prefix=7%203%200%203%200%201%203%203&prefix_number=1&ranking=3%200%201%202\n",
      "Response from SPiCe  b'2 3 3 \\n'\n"
     ]
    }
   ],
   "source": [
    "spiceContentOnPubFile,url_base=predictOnSpicePublicData(0,\"hmm_baseline_v1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def evaluateOnSpiceTrainDataSet(prevContent,url_base):\n",
    "    prefix_number = 2\n",
    "    head=''\n",
    "    content=prevContent\n",
    "    while(head != '[Error]' and head != '[Success]'):\n",
    "        prefix = content[:-1]#Fetch the Sequence returned from Spice Server and exclude the last '\\n'\n",
    "        # Get the ranking\n",
    "        p=prefix.split()\n",
    "        prefix_list=[int(i) for i in p[1:len(p)]]\n",
    "        ranking = getHmmRank(prefix_list,A,A.transpose(),B,B.transpose(),pi,A.shape[0])\n",
    "        ranking_string=list_to_string(ranking[:5])#Here At least alphabet should be 4: Else may get Runtime error\n",
    "        if prefix_number % 200 == 0:\n",
    "            print(\"Prefix number: \" + str(prefix_number) + \" Ranking: \" + ranking_string + \" Prefix: \" + prefix)\n",
    "        # Format the ranking\n",
    "        ranking_string = formatString(ranking_string)\n",
    "        # create prefix with submission needed format\n",
    "        prefix=formatString(prefix)\n",
    "        # Create the url with your ranking to get the next prefix\n",
    "        url = url_base + 'prefix=' + prefix + '&prefix_number=' +\\\n",
    "            str(prefix_number) + '&ranking=' + ranking_string\n",
    "        # Get the answer of the submission on current prefix\n",
    "        response = ur.urlopen(url)\n",
    "        content = response.read()\n",
    "        if not orl2:\n",
    "            # Needed for Python 3.4...\n",
    "            content= content.decode('utf-8')\n",
    "        list_element = content.split()\n",
    "        # modify head in case it is finished or an erro occured\n",
    "        head = str(list_element[0])\n",
    "        # change prefix number\n",
    "        prefix_number += 1\n",
    "    # Post-treatment\n",
    "    # The score is the last element of content (in case of a public test set)\n",
    "    print(content)\n",
    "    list_element = content.split()\n",
    "    score = (list_element[-1])\n",
    "    print(score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prefix number: 200 Ranking: 3 0 1 2 Prefix: 8 3 0 1 3 3 1 2 1 \n",
      "Prefix number: 400 Ranking: 3 0 1 2 Prefix: 7 3 3 1 0 1 3 3 \n",
      "Prefix number: 600 Ranking: 3 0 1 2 Prefix: 3 3 0 1 \n",
      "Prefix number: 800 Ranking: 3 0 1 2 Prefix: 1 3 \n",
      "Prefix number: 1000 Ranking: 3 0 1 2 Prefix: 2 3 0 \n",
      "[Success] Last prefix of the test set. The score of the submission named hmm_baseline_v1 on problem 0 is 0.85064155226946\n",
      "\n",
      "0.85064155226946\n"
     ]
    }
   ],
   "source": [
    "evaluateOnSpiceTrainDataSet(spiceContentOnPubFile,url_base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prefix number: 200 Ranking: 3 0 1 2 Prefix: 8 3 0 1 3 3 1 2 1 \n",
      "Prefix number: 400 Ranking: 3 0 1 2 Prefix: 7 3 3 1 0 1 3 3 \n",
      "Prefix number: 600 Ranking: 3 0 1 2 Prefix: 3 3 0 1 \n",
      "Prefix number: 800 Ranking: 3 0 1 2 Prefix: 1 3 \n",
      "Prefix number: 1000 Ranking: 3 0 1 2 Prefix: 2 3 0 \n",
      "[Success] Last prefix of the test set. The score of the submission named hmm_baseline_Ver1.5 on problem 0 is 0.85064155226946\n",
      "\n",
      "0.85064155226946\n"
     ]
    }
   ],
   "source": [
    "evaluateOnSpiceTrainDataSet(spiceContentOnPubFile,url_base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}