# This file was *autogenerated* from the file FHMM1.sage.
from sage.all_cmdline import *   # import sage library
_sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_6 = Integer(6); _sage_const_8 = Integer(8); _sage_const_10 = Integer(10); _sage_const_0p5 = RealNumber('0.5')
action = {}
class Queue: #A class for the queueing system
    def __init__(self,lmbda,mu,R,TimeCost = _sage_const_1 ,K = False,alpha = _sage_const_1 ,epsilon = _sage_const_10 **-_sage_const_2 ):
        self.lmbda = lmbda
        self.mu = mu
        self.R = R
        self.TC = TimeCost
        self.queuelength = _sage_const_0 
        self.Selfish =  int(R*self.mu)
        self.alpha = alpha
        self.epsilon = epsilon
        self.A = [_sage_const_0 ,_sage_const_1 ]

        if K == False:
            self.K = self.Selfish
        else:
            self.K = K
        print self.K
        self.findQ()
        self.findP()
    def findQ(self): #Transition rate matrix for the queue
        Q = MatrixSpace(QQ,self.K,self.K)
        rowdata = []

        for e in range(self.K):

            for i in range(self.K):
                if i == e+_sage_const_1 :
                    rowdata.append(self.lmbda)
                elif i == e-_sage_const_1 :
                        rowdata.append(self.mu)
                else:
                    rowdata.append(_sage_const_0 )

        self.Q = Q(rowdata)
        rowsums = [sum(e) for e in self.Q]
        for i in range(self.K):
            self.Q[i,i] = -rowsums[i]

    def findP(self): #Converts Q to P
        Diag = []

        for e in range(self.K):

            Diag.append(abs(self.Q[e][e]))

        t = _sage_const_1 /max(Diag)
        print t
        self.P = self.Q*t + identity_matrix(self.Q.nrows())
        self.P2 = [self.P[_sage_const_0 ][_sage_const_0 ],self.P[_sage_const_0 ][_sage_const_1 ]]
    def obtainpi(self):
        uniteigenvector = [e for e in self.P.eigenvectors_left() if e[_sage_const_0 ]==_sage_const_1 ][_sage_const_0 ][_sage_const_1 ][_sage_const_0 ]

        self.pi = uniteigenvector / sum(uniteigenvector)



    def Pij(self,i):

        return self.P2[i]


    def Cost(self,state,action): # Returns the cost
        if action == _sage_const_0 :
            return _sage_const_0 
        else:
            return (state+_sage_const_1 )*(self.mu)**-_sage_const_1  - self.R

    def SolveNHor(self,T=_sage_const_0 ,N=_sage_const_0 ): #Recursively calculates cost for current state

        values = {}
        values[N] = []
        if not T==_sage_const_0 :

            for e in self.A:
                values[N].append(((self.Cost(N,e)+self.Pij(_sage_const_0 )*self.SolveNHor(T-_sage_const_1 ,max(N+e-_sage_const_1 ,_sage_const_0 )))+self.Pij(_sage_const_1 )*self.SolveNHor(T-_sage_const_1 ,N+e)))
            action[N] = values[N].index(min(values[N]))
            return min(values[N])
        else:
            for e in self.A:
                values[N].append((self.Cost(N,e)))

            return min(values[N])
Simple = Queue(_sage_const_8 ,_sage_const_10 ,_sage_const_0p5 )
B = Simple.SolveNHor(T=_sage_const_6 ,N=_sage_const_0 )
print B
for e in action:
    print e,action[e]

bd3954d3a02acbafb2b4fb87dc378d8111cb03b7
# This file was *autogenerated* from the file FHMM1.sage.
from sage.all_cmdline import *   # import sage library
_sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_6 = Integer(6); _sage_const_8 = Integer(8); _sage_const_10 = Integer(10); _sage_const_0p5 = RealNumber('0.5')
action = {}
class Queue: #A class for the queueing system
    def __init__(self,lmbda,mu,R,TimeCost = _sage_const_1 ,K = False,alpha = _sage_const_1 ,epsilon = _sage_const_10 **-_sage_const_2 ):
        self.lmbda = lmbda
        self.mu = mu
        self.R = R
        self.TC = TimeCost
        self.queuelength = _sage_const_0 
        self.Selfish =  int(R*self.mu)
        self.alpha = alpha
        self.epsilon = epsilon
        self.A = [_sage_const_0 ,_sage_const_1 ]

        if K == False:
            self.K = self.Selfish
        else:
            self.K = K
        print self.K
        self.findQ()
        self.findP()
    def findQ(self): #Transition rate matrix for the queue
        Q = MatrixSpace(QQ,self.K,self.K)
        rowdata = []

        for e in range(self.K):

            for i in range(self.K):
                if i == e+_sage_const_1 :
                    rowdata.append(self.lmbda)
                elif i == e-_sage_const_1 :
                        rowdata.append(self.mu)
                else:
                    rowdata.append(_sage_const_0 )

        self.Q = Q(rowdata)
        rowsums = [sum(e) for e in self.Q]
        for i in range(self.K):
            self.Q[i,i] = -rowsums[i]

    def findP(self): #Converts Q to P
        Diag = []

        for e in range(self.K):

            Diag.append(abs(self.Q[e][e]))

        t = _sage_const_1 /max(Diag)
        print t
        self.P = self.Q*t + identity_matrix(self.Q.nrows())
        self.P2 = [self.P[_sage_const_0 ][_sage_const_0 ],self.P[_sage_const_0 ][_sage_const_1 ]]
    def obtainpi(self):
        uniteigenvector = [e for e in self.P.eigenvectors_left() if e[_sage_const_0 ]==_sage_const_1 ][_sage_const_0 ][_sage_const_1 ][_sage_const_0 ]

        self.pi = uniteigenvector / sum(uniteigenvector)



    def Pij(self,i):

        return self.P2[i]


    def Cost(self,state,action): # Returns the cost
        if action == _sage_const_0 :
            return _sage_const_0 
        else:
            return (state+_sage_const_1 )*(self.mu)**-_sage_const_1  - self.R

    def SolveNHor(self,T=_sage_const_0 ,N=_sage_const_0 ): #Recursively calculates cost for current state

        values = {}
        values[N] = []
        if not T==_sage_const_0 :

            for e in self.A:
                values[N].append(((self.Cost(N,e)+self.Pij(_sage_const_0 )*self.SolveNHor(T-_sage_const_1 ,max(N+e-_sage_const_1 ,_sage_const_0 )))+self.Pij(_sage_const_1 )*self.SolveNHor(T-_sage_const_1 ,N+e)))
            action[N] = values[N].index(min(values[N]))
            return min(values[N])
        else:
            for e in self.A:
                values[N].append((self.Cost(N,e)))

            return min(values[N])
Simple = Queue(_sage_const_8 ,_sage_const_10 ,_sage_const_0p5 )
B = Simple.SolveNHor(T=_sage_const_6 ,N=_sage_const_0 )
print B
for e in action:
    print e,action[e]
