HMM_Spice_SubmissionBaseline-Normalized-V2.1.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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",
    "    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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def createRandomMatrixA(numOfStates):\n",
    "    matrixA=np.zeros(shape=(numOfStates,numOfStates),dtype=float)\n",
    "    prob=1.0/(numOfStates*numOfStates)\n",
    "    for i in np.arange(numOfStates):\n",
    "        for j in np.arange(numOfStates):\n",
    "            matrixA[i][j]=prob\n",
    "    return matrixA\n",
    "def createRandomMatrixB(numOfStates,distinctObservations):\n",
    "    matrixB=np.zeros(shape=(numOfStates,distinctObservations),dtype=float)\n",
    "    prob=1.0/(numOfStates*distinctObservations)\n",
    "    for i in np.arange(numOfStates):\n",
    "        for j in np.arange(distinctObservations):\n",
    "            matrixB[i][j]=prob\n",
    "    return matrixB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeAlpha(observations,a,b,pi,alphaDP):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    alphaDpScaleTime0=0    \n",
    "    for i in np.arange(statesC):\n",
    "        alphaDP[0][i]=pi[i]*b[i][observations[0]]\n",
    "        alphaDpScaleTime0+=alphaDP[0][i]\n",
    "    for i in np.arange(statesC):\n",
    "        alphaDP[0][i]/=alphaDpScaleTime0\n",
    "    #print(\"Initial Alpha \",alphaDP[0])\n",
    "    for t in np.arange(1,timePts):\n",
    "        alphaDpScaleTimeT=0\n",
    "        for i in np.arange(statesC):\n",
    "            for j in np.arange(statesC):\n",
    "                alphaDP[t][i]+=alphaDP[t-1][j]*a[j][i]\n",
    "            alphaDP[t][i]*=b[i][observations[t]]\n",
    "            alphaDpScaleTimeT+=alphaDP[t][i]\n",
    "        #print(\"Scaling \",alphaDpScaleTimeT)\n",
    "        for i in np.arange(statesC):\n",
    "            alphaDP[t][i]/=alphaDpScaleTimeT\n",
    "    #print(\"Next AlphaDP \",alphaDP[1])\n",
    "    #print(\"alphaDP \",alphaDP)\n",
    "def observationsLikelihood(alphaDP):\n",
    "    timePts=alphaDP.shape[0]\n",
    "    stateC=alphaDP.shape[1]\n",
    "    ans=0.0\n",
    "    for i in np.arange(stateC):\n",
    "        ans+=alphaDP[timePts-1][i]\n",
    "        #print(\"Alpha@Timept \",i ,\" : \", alphaDP[i])\n",
    "    #print(\"Observation Likelihood \",ans)\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"β[t][i]: Probability of seeing the observations from t + 1 to end(T), given that we are in state i at time t\n",
    "beta[t][i]= a[i][0]*b[0][o[t+1]]*beta[t+1][0] + a[i][1]*b[1][o[t+1]]*beta[t+1][1] + ......... +\n",
    "                a[i][N-2]*b[N-2][o[t+1]]*beta[t+1][N-2] + a[i][N-1]*b[N-1][o[t+1]]*beta[t+1][N-1]\"\"\"\n",
    "def computeBeta(observations,a,b,pi,betaDP):\n",
    "    statesC=a.shape[0]\n",
    "    timePts=observations.shape[0]\n",
    "    if timePts<1:\n",
    "        return\n",
    "    for state in np.arange(statesC):\n",
    "            betaDP[timePts-1][state]=1\n",
    "    for t in np.arange(timePts-2,-1,-1):\n",
    "        betaDpScaleTimeT=0\n",
    "        for i in np.arange(statesC):\n",
    "            for j in np.arange(statesC):\n",
    "                betaDP[t][i]+=a[i][j]*b[j][observations[t+1]]*betaDP[t+1][j]\n",
    "            betaDpScaleTimeT+=betaDP[t][i]\n",
    "        #if t==timePts-2:\n",
    "            #print(\"Non Scaled T-2 \",betaDP[t])\n",
    "        #print(\"Scale At \",t,\" \",betaDpScaleTimeT)\n",
    "        for i in np.arange(statesC):\n",
    "            betaDP[t][i]/=betaDpScaleTimeT\n",
    "    #print(\"BetaDP t-2 \",betaDP[timePts-2])\n",
    "    #print(\"betaDP \",betaDP)\n",
    "    return betaDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeDiGammaNum(t,i,j,alphaDP,betaDP,a,b,observations):\n",
    "    return alphaDP[t][i]*a[i][j]*b[j][observations[t+1]]*betaDP[t+1][j]\n",
    "def computeDiGammaDP(alphaDP,betaDP,a,b,observations):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    diGammaDP=np.zeros(shape=(statesC,statesC),dtype=float)\n",
    "    diGammaDenom=observationsLikelihood(alphaDP)\n",
    "    #print(\"Observation Likelihood \",diGammaDenom ,\" \",observations)\n",
    "    #print(\"Observation \",observations)\n",
    "    #print(\"DiGammaDenom \",diGammaDenom ,\"Observation : \",observations,\"AlphaDP \",alphaDP)\n",
    "    for i in np.arange(statesC):\n",
    "        for j in np.arange(statesC):\n",
    "                for t in np.arange(observationsC-1):\n",
    "                    #print(\"DiGama \",i,\" \",j,\" \",t,\" \",computeDiGammaNum(t,i,j,alphaDP,betaDP,a,b,observations))\n",
    "                    diGammaDP[i][j]+=computeDiGammaNum(t,i,j,alphaDP,betaDP,a,b,observations)/diGammaDenom\n",
    "    #print(\"diGammaDP \",diGammaDP)\n",
    "    return diGammaDP\n",
    "def computeTransitionProbabilityA(alphaDP,betaDP,a,b,observations):\n",
    "    #print(\"Got A \",a)\n",
    "    statesC=alphaDP.shape[1]\n",
    "    newlyComputedTransitionProbA=np.zeros(shape=(statesC,statesC),dtype=float)\n",
    "    diGammaDP=computeDiGammaDP(alphaDP,betaDP,a,b,observations)\n",
    "    #print(\"DiGamma \",diGammaDP)\n",
    "    diGammaDPSumGrpByJ=np.zeros(shape=(statesC),dtype=float)\n",
    "    #print(\"DiGammaDP \",diGammaDP)\n",
    "    for i in np.arange(statesC):\n",
    "        sumAcrossJ=0.0\n",
    "        for j in np.arange(statesC):\n",
    "            sumAcrossJ+=diGammaDP[i][j]\n",
    "        diGammaDPSumGrpByJ[i]=sumAcrossJ\n",
    "    #print(\"DiGramGrpByJ \",diGammaDPSumGrpByJ)\n",
    "    for i in np.arange(statesC):    \n",
    "        for j in np.arange(statesC):\n",
    "            if (diGammaDPSumGrpByJ[i]==0):\n",
    "                #print(\"Underflow \",i,\" \",diGammaDPSumGrpByJ[i])\n",
    "                newlyComputedTransitionProbA[i][j]=0.0\n",
    "            else:\n",
    "                newlyComputedTransitionProbA[i][j]=diGammaDP[i][j]/diGammaDPSumGrpByJ[i]\n",
    "    #print(\"Updated A \",newlyComputedTransitionProbA)\n",
    "    return newlyComputedTransitionProbA   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeGammaNum(t,j,alphaDP,betaDP):\n",
    "    return alphaDP[t][j]*betaDP[t][j]\n",
    "def computeGammaDP(alphaDP,betaDP):\n",
    "    observationsC=alphaDP.shape[0]\n",
    "    statesC=alphaDP.shape[1]\n",
    "    gammaDP=np.zeros(shape=(statesC,observationsC),dtype=float)\n",
    "    gammaDenom=observationsLikelihood(alphaDP)\n",
    "    for i in np.arange(statesC):\n",
    "        for t in np.arange(observationsC):\n",
    "            gammaDP[i][t]=computeGammaNum(t,i,alphaDP,betaDP)/gammaDenom  \n",
    "    return gammaDP\n",
    "def computeObsrProbNum(gammaDP,i,vk,observations):\n",
    "    observationC=len(observations)\n",
    "    obsrProbNum=0.0\n",
    "    for t in np.arange(observationC):\n",
    "        if observations[t]==vk:\n",
    "            obsrProbNum+=gammaDP[i][t]\n",
    "    return obsrProbNum\n",
    "def computeTransitionProbabilityB(alphaDP,betaDP,a,b,observations,observationDict):\n",
    "    statesC=a.shape[0]\n",
    "    observationsC=b.shape[1]\n",
    "    newlyComputedObsrProbB=np.zeros(shape=(statesC,observationsC),dtype=float)\n",
    "    gammaDP=computeGammaDP(alphaDP,betaDP) \n",
    "    #print(\"gammaDP \",gammaDP)\n",
    "    #print(\"ALPHADP : \",alphaDP)\n",
    "    #print(\"BETADP : \",betaDP)\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": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Change Convergence Criteria to be more reasonable/Useful\n",
    "def isConverged(count):\n",
    "    if count>=7:\n",
    "        return True\n",
    "    return False\n",
    "def Forward_Backward_EM_Algo(observations,A,B,pi,observationDict):\n",
    "    count=0\n",
    "    updatedA=A\n",
    "    updatedB=B\n",
    "    while isConverged(count)==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",
    "        computeAlpha(observations,updatedA,updatedB,pi,alphaDP)\n",
    "        #print(\"AlphaDP Computed ....\")\n",
    "        #print(\"ALPHA-DP \",alphaDP)\n",
    "        computeBeta(observations,updatedA,updatedB,pi,betaDP)\n",
    "        #Maximization(M)-Step\n",
    "        newA=computeTransitionProbabilityA(alphaDP,betaDP,updatedA,updatedB,observations)\n",
    "        newB=computeTransitionProbabilityB(alphaDP,betaDP,updatedA,updatedB,observations,observationDict)\n",
    "        #print(\"New A =================================>\")\n",
    "        #print(newA)\n",
    "        #print(\"New B =================================>\")\n",
    "        #print(newB)\n",
    "        updatedA=newA\n",
    "        updatedB=newB\n",
    "        count=count+1\n",
    "    return (updatedA,updatedB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trainHMM(trainDataFile,A,B,pi,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",
    "    #for n in range(numSequences):\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",
    "        l = line.split(\" \")\n",
    "        #print(\"For Sequence \",n,\" =====================================>\")\n",
    "        if(int(l[0])<=1):\n",
    "            #print(\"Skipping \",l)\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,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",
    "        #print(\"State Transition Matrix (A) ===>\")\n",
    "        #print(learnedParams[0])\n",
    "        #print(\"Observation Probability Matrix (B) ===>\")\n",
    "        #print(learnedParams[1])\n",
    "        #print(\"Aggregated State Transition Matrix (A) ===>\")\n",
    "        #print(updatedA)\n",
    "        #print(\"Aggregated Observation Probability Matrix (B) ===>\")\n",
    "        #print(updatedB)\n",
    "    updatedA=updatedA/actuallyUsedSeqs\n",
    "    updatedB=updatedB/actuallyUsedSeqs\n",
    "    return (updatedA,updatedB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def trainModel(fileLoc,maxNoOfStates,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",
    "    #print(initialProbs)\n",
    "    A=createRandomMatrixA(numOfStates)\n",
    "    B=createRandomMatrixB(numOfStates,distinctObservations)\n",
    "    #print(A)\n",
    "    #print(B)\n",
    "    trainedParams=trainHMM(fileLoc,A,B,pi,maxSequences)\n",
    "    end=time.time()\n",
    "    #print(\"For \",maxSequences,\" Sequences : Total Training Time \",end-start)\n",
    "    return trainedParams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Initial Prob. in  0.16648101806640625 seconds\n"
     ]
    }
   ],
   "source": [
    "(A,B)=trainModel('Data/1.spice.train.txt',20,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Initial Prob. in  0.24859118461608887 seconds\n",
      "Computed Initial Prob. in  0.26802635192871094 seconds\n",
      "Computed Initial Prob. in  0.36754536628723145 seconds\n",
      "Computed Initial Prob. in  0.25893568992614746 seconds\n",
      "1 loops, best of 3: 45.8 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit trainModel('Data/1.spice.train.txt',20,15)\n",
    "#(A,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Initial Prob. in  0.21947813034057617 seconds\n",
      "(20, 20, [0.036700000000000003, 0.0091999999999999998, 0.020549999999999999, 0.081900000000000001, 0.042299999999999997, 0.04845, 0.11260000000000001, 0.0385, 0.018249999999999999, 0.045999999999999999, 0.016750000000000001, 0.042049999999999997, 0.049799999999999997, 0.065449999999999994, 0.06855, 0.1487, 0.028850000000000001, 0.021999999999999999, 0.0591, 0.044299999999999999])\n",
      "For  150  Sequences : Total Training Time  513.4503910541534\n"
     ]
    }
   ],
   "source": [
    "(A,B)=trainModel('Data/1.spice.train.txt',20,150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Initial Prob. in  0.17683792114257812 seconds\n",
      "(20, 20, [0.036700000000000003, 0.0091999999999999998, 0.020549999999999999, 0.081900000000000001, 0.042299999999999997, 0.04845, 0.11260000000000001, 0.0385, 0.018249999999999999, 0.045999999999999999, 0.016750000000000001, 0.042049999999999997, 0.049799999999999997, 0.065449999999999994, 0.06855, 0.1487, 0.028850000000000001, 0.021999999999999999, 0.0591, 0.044299999999999999])\n",
      "For  1500  Sequences : Total Training Time  6131.2775712013245\n"
     ]
    }
   ],
   "source": [
    "(A,B)=trainModel('Data/1.spice.train.txt',20,1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#(A,B)=trainModel('Data/1.spice.train.txt',20,1600)\n",
    "#(A,B)"
   ]
  },
  {
   "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
}