mdp.py 58.6 KB
Newer Older
Steven Cordwell's avatar
Steven Cordwell committed
1
2
# -*- coding: utf-8 -*-
"""
3
Copyright (c) 2011, 2012, 2013 Steven Cordwell
Steven Cordwell's avatar
Steven Cordwell committed
4
5
6
7
Copyright (c) 2009, Iadine Chadès
Copyright (c) 2009, Marie-Josée Cros
Copyright (c) 2009, Frédérick Garcia
Copyright (c) 2009, Régis Sabbadin
Steven Cordwell's avatar
Steven Cordwell committed
8
9
10

All rights reserved.

Steven Cordwell's avatar
Steven Cordwell committed
11
12
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
Steven Cordwell's avatar
Steven Cordwell committed
13
14
15
16
17
18

  * Redistributions of source code must retain the above copyright notice, this
    list of conditions and the following disclaimer.
  * Redistributions in binary form must reproduce the above copyright notice,
    this list of conditions and the following disclaimer in the documentation
    and/or other materials provided with the distribution.
Steven Cordwell's avatar
Steven Cordwell committed
19
20
21
  * Neither the name of the <ORGANIZATION> nor the names of its contributors
    may be used to endorse or promote products derived from this software
    without specific prior written permission.
Steven Cordwell's avatar
Steven Cordwell committed
22
23
24
25

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
Steven Cordwell's avatar
Steven Cordwell committed
26
27
28
29
30
31
32
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Steven Cordwell's avatar
Steven Cordwell committed
33
34
"""

Steven Cordwell's avatar
Steven Cordwell committed
35
from numpy import absolute, array, diag, matrix, mean, mod, multiply, ndarray
36
from numpy import ones, zeros
Steven Cordwell's avatar
Steven Cordwell committed
37
from numpy.random import rand
Steven Cordwell's avatar
Steven Cordwell committed
38
from math import ceil, log, sqrt
Steven Cordwell's avatar
Steven Cordwell committed
39
40
from random import randint, random
from scipy.sparse import csr_matrix as sparse
Steven Cordwell's avatar
Steven Cordwell committed
41
42
from time import time

Steven Cordwell's avatar
Steven Cordwell committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
mdperr = {
"mat_nonneg" :
    "PyMDPtoolbox: Probabilities must be non-negative.",
"mat_square" :
    "PyMDPtoolbox: The matrix must be square.",
"mat_stoch" :
    "PyMDPtoolbox: Rows of the matrix must sum to one (1).",
"mask_numpy" :
    "PyMDPtoolbox: mask must be a numpy array or matrix; i.e. type(mask) is "
    "ndarray or type(mask) is matrix.", 
"mask_SbyS" : 
    "PyMDPtoolbox: The mask must have shape SxS; i.e. mask.shape = (S, S).",
"obj_shape" :
    "PyMDPtoolbox: Object arrays for transition probabilities and rewards "
    "must have only 1 dimension: the number of actions A. Each element of "
    "the object array contains an SxS ndarray or matrix.",
"obj_square" :
    "PyMDPtoolbox: Each element of an object array for transition "
    "probabilities and rewards must contain an SxS ndarray or matrix; i.e. "
    "P[a].shape = (S, S) or R[a].shape = (S, S).",
Steven Cordwell's avatar
Steven Cordwell committed
63
"P_type" :
Steven Cordwell's avatar
Steven Cordwell committed
64
65
66
67
68
69
70
71
72
73
    "PyMDPtoolbox: The transition probabilities must be in a numpy array; "
    "i.e. type(P) is ndarray.",
"P_shape" :
    "PyMDPtoolbox: The transition probability array must have the shape "
    "(A, S, S)  with S : number of states greater than 0 and A : number of "
    "actions greater than 0. i.e. R.shape = (A, S, S)",
"PR_incompat" :
    "PyMDPtoolbox: Incompatibility between P and R dimensions.",
"prob_in01" :
    "PyMDPtoolbox: Probability p must be in [0; 1].",
Steven Cordwell's avatar
Steven Cordwell committed
74
"R_type" :
Steven Cordwell's avatar
Steven Cordwell committed
75
    "PyMDPtoolbox: The rewards must be in a numpy array; i.e. type(R) is "
Steven Cordwell's avatar
Steven Cordwell committed
76
    "ndarray, or numpy matrix; i.e. type(R) is matrix.",
Steven Cordwell's avatar
Steven Cordwell committed
77
78
79
80
81
82
83
84
85
86
"R_shape" :
    "PyMDPtoolbox: The reward matrix R must be an array of shape (A, S, S) or "
    "(S, A) with S : number of states greater than 0 and A : number of actions "
    "greater than 0. i.e. R.shape = (S, A) or (A, S, S).",
"R_gt_0" :
    "PyMDPtoolbox: The rewards must be greater than 0.",
"S_gt_1" :
    "PyMDPtoolbox: Number of states S must be greater than 1.",
"SA_gt_1" : 
    "PyMDPtoolbox: The number of states S and the number of actions A must be "
87
88
89
90
91
    "greater than 1.",
"discount_rng" : 
    "PyMDPtoolbox: Discount rate must be in ]0; 1]",
"maxi_min" :
    "PyMDPtoolbox: The maximum number of iterations must be greater than 0"
Steven Cordwell's avatar
Steven Cordwell committed
92
93
}

Steven Cordwell's avatar
Steven Cordwell committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
def check(P, R):
    """Checks if the matrices P and R define a Markov Decision Process.
    
    Let S = number of states, A = number of actions.
    The transition matrix P must be on the shape (A, S, S) and P[a,:,:]
    must be stochastic.
    The reward matrix R must be on the shape (A, S, S) or (S, A).
    Raises an error if P and R do not define a MDP.
    
    Parameters
    ---------
    P : transition matrix (A, S, S)
        P could be an array with 3 dimensions or a object array (A, ),
        each cell containing a matrix (S, S) possibly sparse
    R : reward matrix (A, S, S) or (S, A)
        R could be an array with 3 dimensions (SxSxA) or a object array
        (A, ), each cell containing a sparse matrix (S, S) or a 2D
        array(S, A) possibly sparse  
    """
    
    # Check of P
    # tranitions must be a numpy array either an AxSxS ndarray (with any 
    # dtype other than "object"); or, a 1xA ndarray with a "object" dtype, 
    # and each element containing an SxS array. An AxSxS array will be
    # be converted to an object array. A numpy object array is similar to a
    # MATLAB cell array.
    if (not type(P) is ndarray):
        raise TypeError(mdperr["P_type"])
    
    if (not type(R) is ndarray):
        raise TypeError(mdperr["R_type"])
    
    # NumPy has an array type of 'object', which is roughly equivalent to
    # the MATLAB cell array. These are most useful for storing sparse
    # matrices as these can only have two dimensions whereas we want to be
    # able to store a transition matrix for each action. If the dytpe of
    # the transition probability array is object then we store this as
    # P_is_object = True.
    # If it is an object array, then it should only have one dimension
    # otherwise fail with a message expalining why.
    # If it is a normal array then the number of dimensions must be exactly
    # three, otherwise fail with a message explaining why.
    if (P.dtype == object):
        if (P.ndim > 1):
            raise ValueError(mdperr["obj_shape"])
        else:
            P_is_object = True
    else:
        if (P.ndim != 3):
            raise ValueError(mdperr["P_shape"])
        else:
            P_is_object = False
    
    # As above but for the reward array. A difference is that the reward
    # array can have either two or 3 dimensions.
    if (R.dtype == object):
        if (R.ndim > 1):
            raise ValueError(mdperr["obj_shape"])
        else:
            R_is_object = True
    else:
        if (not R.ndim in (2, 3)):
            raise ValueError(mdperr["R_shape"])
        else:
            R_is_object = False
    
    # We want to make sure that the transition probability array and the 
    # reward array are in agreement. This means that both should show that
    # there are the same number of actions and the same number of states.
    # Furthermore the probability of transition matrices must be SxS in
    # shape, so we check for that also.
    if P_is_object:
        # If the user has put their transition matrices into a numpy array
        # with dtype of 'object', then it is possible that they have made a
        # mistake and not all of the matrices are of the same shape. So,
        # here we record the number of actions and states that the first
        # matrix in element zero of the object array says it has. After
        # that we check that every other matrix also reports the same
        # number of actions and states, otherwise fail with an error.
        # aP: the number of actions in the transition array. This
        # corresponds to the number of elements in the object array.
        aP = P.shape[0]
        # sP0: the number of states as reported by the number of rows of
        # the transition matrix
        # sP1: the number of states as reported by the number of columns of
        # the transition matrix
        sP0, sP1 = P[0].shape
        # Now we check to see that every element of the object array holds
        # a matrix of the same shape, otherwise fail.
        for aa in range(1, aP):
            # sp0aa and sp1aa represents the number of states in each
            # subsequent element of the object array. If it doesn't match
            # what was found in the first element, then we need to fail
            # telling the user what needs to be fixed.
            sP0aa, sP1aa = P[aa].shape
            if ((sP0aa != sP0) or (sP1aa != sP1)):
                raise ValueError(mdperr["obj_square"])
    else:
        # if we are using a normal array for this, then the first
        # dimension should be the number of actions, and the second and 
        # third should be the number of states
        aP, sP0, sP1 = P.shape
    
    # the first dimension of the transition matrix must report the same
    # number of states as the second dimension. If not then we are not
    # dealing with a square matrix and it is not a valid transition
    # probability. Also, if the number of actions is less than one, or the
    # number of states is less than one, then it also is not a valid
    # transition probability.
    if ((sP0 < 1) or (aP < 1) or (sP0 != sP1)):
        raise ValueError(mdperr["P_shape"])
    
    # now we check that each transition matrix is square-stochastic. For
    # object arrays this is the matrix held in each element, but for
    # normal arrays this is a matrix formed by taking a slice of the array
    for aa in range(aP):
        if P_is_object:
            checkSquareStochastic(P[aa])
        else:
            checkSquareStochastic(P[aa, :, :])
        # aa = aa + 1 # why was this here?
    
    if R_is_object:
        # if the rewarad array has an object dtype, then we check that
        # each element contains a matrix of the same shape as we did 
        # above with the transition array.
        aR = R.shape[0]
        sR0, sR1 = R[0].shape
        for aa in range(1, aR):
            sR0aa, sR1aa = R[aa].shape
            if ((sR0aa != sR0) or (sR1aa != sR1)):
                raise ValueError(mdperr["obj_square"])
    elif (R.ndim == 3):
        # This indicates that the reward matrices are constructed per 
        # transition, so that the first dimension is the actions and
        # the second two dimensions are the states.
        aR, sR0, sR1 = R.shape
    else:
        # then the reward matrix is per state, so the first dimension is 
        # the states and the second dimension is the actions.
        sR0, aR = R.shape
        # this is added just so that the next check doesn't error out
        # saying that sR1 doesn't exist
        sR1 = sR0
    
    # the number of actions must be more than zero, the number of states
    # must also be more than 0, and the states must agree
    if ((sR0 < 1) or (aR < 1) or (sR0 != sR1)):
        raise ValueError(mdperr["R_shape"])
    
    # now we check to see that what the transition array is reporting and
    # what the reward arrar is reporting agree as to the number of actions
    # and states. If not then fail explaining the situation
    if (sP0 != sR0) or (aP != aR):
        raise ValueError(mdperr["PR_incompat"])
        
    # We are at the end of the checks, so if no exceptions have been raised
    # then that means there are (hopefullly) no errors and we return None
    return None

def checkSquareStochastic(Z):
    """Check if Z is a square stochastic matrix
    
    Parameters
    ----------
    Z : a SxS matrix. It could be a numpy ndarray SxS, or a scipy.sparse 
        csr_matrix
    
    Evaluation
    ----------
    Returns None if no error has been detected
    """
    s1, s2 = Z.shape
    if (s1 != s2):
        raise ValueError(mdperr["mat_square"])
    elif (absolute(Z.sum(axis=1) - ones(s2))).max() > 10**(-12):
        raise ValueError(mdperr["mat_stoch"])
    elif ((type(Z) is ndarray) or (type(Z) is matrix)) and (Z < 0).any():
        raise ValueError(mdperr["mat_nonneg"])
    elif (type(Z) is sparse) and (Z.data < 0).any():
        raise ValueError(mdperr["mat_nonneg"]) 
    else:
        return(None)

Steven Cordwell's avatar
Steven Cordwell committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
def exampleForest(S=3, r1=4, r2=2, p=0.1):
    """
    Generates a Markov Decision Process example based on a simple forest
    management.
    
    See the related documentation for more detail.
    
    Parameters
    ---------
    S : number of states (> 0), optional (default 3)
    r1 : reward when forest is in the oldest state and action Wait is performed,
        optional (default 4)
    r2 : reward when forest is in the oldest state and action Cut is performed, 
        optional (default 2)
    p : probability of wild fire occurence, in ]0, 1[, optional (default 0.1)
    
    Evaluation
    ----------
    P : transition probability matrix (A, S, S)
    R : reward matrix (S, A)
Steven Cordwell's avatar
Steven Cordwell committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
    
    Examples
    --------
    >>> import mdp
    >>> P, R = mdp.exampleForest()
    >>> P
    array([[[ 0.1,  0.9,  0. ],
            [ 0.1,  0. ,  0.9],
            [ 0.1,  0. ,  0.9]],

           [[ 1. ,  0. ,  0. ],
            [ 1. ,  0. ,  0. ],
            [ 1. ,  0. ,  0. ]]])
    >>> R
    array([[ 0.,  0.],
           [ 0.,  1.],
           [ 4.,  2.]])
    
Steven Cordwell's avatar
Steven Cordwell committed
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
    """
    if (S <= 1):
        raise ValueError(mdperr["S_gt_1"])
    if (r1 <= 0) or (r2 <= 0):
        raise ValueError(mdperr["R_gt_0"])
    if (p < 0 or p > 1):
        raise ValueError(mdperr["prob_in01"])
    
    # Definition of Transition matrix P(:,:,1) associated to action Wait (action 1) and
    # P(:,:,2) associated to action Cut (action 2)
    #             | p 1-p 0.......0  |                  | 1 0..........0 |
    #             | .  0 1-p 0....0  |                  | . .          . |
    #  P(:,:,1) = | .  .  0  .       |  and P(:,:,2) =  | . .          . |
    #             | .  .        .    |                  | . .          . |
    #             | .  .         1-p |                  | . .          . |
    #             | p  0  0....0 1-p |                  | 1 0..........0 |
    P = zeros((2, S, S))
    P[0, :, :] = (1 - p) * diag(ones(S - 1), 1)
    P[0, :, 0] = p
    P[0, S - 1, S - 1] = (1 - p)
    P[1, :, :] = zeros((S, S))
    P[1, :, 0] = 1
    
    # Definition of Reward matrix R1 associated to action Wait and 
    # R2 associated to action Cut
    #           | 0  |                   | 0  |
    #           | .  |                   | 1  |
    #  R(:,1) = | .  |  and     R(:,2) = | .  |	
    #           | .  |                   | .  |
    #           | 0  |                   | 1  |                   
    #           | r1 |                   | r2 |
    R = zeros((S, 2))
    R[S - 1, 0] = r1
    R[:, 1] = ones(S)
    R[0, 1] = 0
    R[S - 1, 1] = r2
    
    return (P, R)

def exampleRand(S, A, is_sparse=False, mask=None):
    """Generates a random Markov Decision Process.
    
    Parameters
    ----------
    S : number of states (> 0)
    A : number of actions (> 0)
    is_sparse : false to have matrices in plain format, true to have sparse
        matrices optional (default false).
    mask : matrix with 0 and 1 (0 indicates a place for a zero
Steven Cordwell's avatar
Steven Cordwell committed
365
           probability), optional (SxS) (default, random)
Steven Cordwell's avatar
Steven Cordwell committed
366
367
368
369
370
    
    Returns
    ----------
    P : transition probability matrix (SxSxA)
    R : reward matrix (SxSxA)
Steven Cordwell's avatar
Steven Cordwell committed
371
372
373
374
375
376

    Examples
    --------
    >>> import mdp
    >>> P, R = mdp.exampleRand(5, 3)

Steven Cordwell's avatar
Steven Cordwell committed
377
378
379
380
381
382
383
384
385
386
387
    """
    if (S < 1 or A < 1):
        raise ValueError(mdperr["SA_gt_1"])
    
    try:
        if (mask != None) and ((mask.shape[0] != S) or (mask.shape[1] != S)):
            raise ValueError(mdperr["mask_SbyS"])
    except AttributeError:
        raise TypeError(mdperr["mask_numpy"])
    
    if mask == None:
Steven Cordwell's avatar
Steven Cordwell committed
388
389
390
391
392
        mask = rand(A, S, S)
        for a in range(A):
            r = random()
            mask[a][mask[a] < r] = 0
            mask[a][mask[a] >= r] = 1
Steven Cordwell's avatar
Steven Cordwell committed
393
394
395
    
    if is_sparse:
        # definition of transition matrix : square stochastic matrix
Steven Cordwell's avatar
Steven Cordwell committed
396
        P = zeros((A, ), dtype=object)
Steven Cordwell's avatar
Steven Cordwell committed
397
        # definition of reward matrix (values between -1 and +1)
Steven Cordwell's avatar
Steven Cordwell committed
398
        R = zeros((A, ), dtype=object)
Steven Cordwell's avatar
Steven Cordwell committed
399
        for a in range(A):
Steven Cordwell's avatar
Steven Cordwell committed
400
            PP = mask[a] * rand(S, S)
Steven Cordwell's avatar
Steven Cordwell committed
401
            for s in range(S):
Steven Cordwell's avatar
Steven Cordwell committed
402
403
                if (mask[a, s, :].sum() == 0):
                    PP[s, randint(0, S - 1)] = 1
Steven Cordwell's avatar
Steven Cordwell committed
404
405
                PP[s, :] = PP[s, :] / PP[s, :].sum()
            P[a] = sparse(PP)
Steven Cordwell's avatar
Steven Cordwell committed
406
            R[a] = sparse(mask[a] * (2 * rand(S, S) - ones((S, S))))
Steven Cordwell's avatar
Steven Cordwell committed
407
408
409
410
411
412
    else:
        # definition of transition matrix : square stochastic matrix
        P = zeros((A, S, S))
        # definition of reward matrix (values between -1 and +1)
        R = zeros((A, S, S))
        for a in range(A):
Steven Cordwell's avatar
Steven Cordwell committed
413
            P[a, :, :] = mask[a] * rand(S, S)
Steven Cordwell's avatar
Steven Cordwell committed
414
            for s in range(S):
Steven Cordwell's avatar
Steven Cordwell committed
415
416
                if (mask[a, s, :].sum() == 0):
                    P[a, s, randint(0, S - 1)] = 1
Steven Cordwell's avatar
Steven Cordwell committed
417
                P[a, s, :] = P[a, s, :] / P[a, s, :].sum()
Steven Cordwell's avatar
Steven Cordwell committed
418
            R[a, :, :] = mask[a] * (2 * rand(S, S) - ones((S, S), dtype=int))
Steven Cordwell's avatar
Steven Cordwell committed
419
420
421
    
    return (P, R)

Steven Cordwell's avatar
Steven Cordwell committed
422
def getSpan(W):
423
424
425
426
427
428
    """Returns the span of W
    
    sp(W) = max W(s) - min W(s)
    """
    return (W.max() - W.min())

Steven Cordwell's avatar
Steven Cordwell committed
429
class MDP(object):
Steven Cordwell's avatar
Steven Cordwell committed
430
    """The Markov Decision Problem Toolbox."""
Steven Cordwell's avatar
Steven Cordwell committed
431
    
432
    def __init__(self, transitions, reward, discount, max_iter):
Steven Cordwell's avatar
Steven Cordwell committed
433
        """"""
434
        
Steven Cordwell's avatar
Steven Cordwell committed
435
436
        # if the discount is None then the algorithm is assumed to not use it
        # in its computations
Steven Cordwell's avatar
Steven Cordwell committed
437
438
439
440
        if (type(discount) is int) or (type(discount) is float):
            if (discount <= 0) or (discount > 1):
                raise ValueError(mdperr["discount_rng"])
            else:
Steven Cordwell's avatar
Steven Cordwell committed
441
442
443
                if discount == 1:
                    print("PyMDPtoolbox WARNING: check conditions of convergence."\
                        "With no discount, convergence is not always assumed.")
Steven Cordwell's avatar
Steven Cordwell committed
444
445
                self.discount = discount
        elif not discount is None:
Steven Cordwell's avatar
Steven Cordwell committed
446
447
            raise ValueError("PyMDPtoolbox: the discount must be a positive " \
                "real number less than or equal to one.")
448
        
Steven Cordwell's avatar
Steven Cordwell committed
449
450
        # if the max_iter is None then the algorithm is assumed to not use it
        # in its computations
Steven Cordwell's avatar
Steven Cordwell committed
451
452
453
454
455
456
        if (type(max_iter) is int) or (type(max_iter) is float):
            if (max_iter <= 0):
                raise ValueError(mdperr["maxi_min"])
            else:
                self.max_iter = max_iter
        elif not max_iter is None:
Steven Cordwell's avatar
Steven Cordwell committed
457
458
            raise ValueError("PyMDPtoolbox: max_iter must be a positive real "\
                "number greater than zero.")
459
        
Steven Cordwell's avatar
Steven Cordwell committed
460
461
462
        # we run a check on P and R to make sure they are describing an MDP. If
        # an exception isn't raised then they are assumed to be correct.
        check(transitions, reward)
463
        
Steven Cordwell's avatar
Steven Cordwell committed
464
        # computePR will assign the variables self.S, self.A, self.P and self.R
465
466
        self.computePR(transitions, reward)
        
Steven Cordwell's avatar
Steven Cordwell committed
467
468
469
470
        # the verbosity is by default turned off
        self.verbose = False
        # Initially the time taken to perform the computations is set to None
        self.time = None
471
472
        # set the initial iteration count to zero
        self.iter = 0
Steven Cordwell's avatar
Steven Cordwell committed
473
        
Steven Cordwell's avatar
Steven Cordwell committed
474
        self.V = None
Steven Cordwell's avatar
Steven Cordwell committed
475
        self.policy = None
Steven Cordwell's avatar
Steven Cordwell committed
476
    
Steven Cordwell's avatar
Steven Cordwell committed
477
    def bellmanOperator(self):
Steven Cordwell's avatar
Steven Cordwell committed
478
479
        """
        Applies the Bellman operator on the value function.
Steven Cordwell's avatar
Steven Cordwell committed
480
        
Steven Cordwell's avatar
Steven Cordwell committed
481
482
        Updates the value function and the Vprev-improving policy.
        
483
        Returns
Steven Cordwell's avatar
Steven Cordwell committed
484
        -------
485
        (policy, value) : tuple of new policy and its value
Steven Cordwell's avatar
Steven Cordwell committed
486
487
488
        """
        Q = matrix(zeros((self.S, self.A)))
        for aa in range(self.A):
Steven Cordwell's avatar
Steven Cordwell committed
489
            Q[:, aa] = self.R[:, aa] + (self.discount * self.P[aa] * self.V)
Steven Cordwell's avatar
Steven Cordwell committed
490
        
491
492
493
        # Which way is better? if choose the first way, then the classes that
        # call this function must be changed
        # 1. Return, (policy, value)
494
        return (Q.argmax(axis=1), Q.max(axis=1))
Steven Cordwell's avatar
Steven Cordwell committed
495
496
        # 2. update self.policy and self.V directly
        # self.V = Q.max(axis=1)
497
        # self.policy = Q.argmax(axis=1)
Steven Cordwell's avatar
Steven Cordwell committed
498
499
500
501
    
    def computePR(self, P, R):
        """Computes the reward for the system in one state chosing an action
        
Steven Cordwell's avatar
Steven Cordwell committed
502
        Arguments
Steven Cordwell's avatar
Steven Cordwell committed
503
        ---------
Steven Cordwell's avatar
Steven Cordwell committed
504
505
506
507
508
509
510
511
        Let S = number of states, A = number of actions
            P(SxSxA)  = transition matrix 
                P could be an array with 3 dimensions or  a cell array (1xA), 
                each cell containing a matrix (SxS) possibly sparse
            R(SxSxA) or (SxA) = reward matrix
                R could be an array with 3 dimensions (SxSxA) or  a cell array 
                (1xA), each cell containing a sparse matrix (SxS) or a 2D 
                array(SxA) possibly sparse  
Steven Cordwell's avatar
Steven Cordwell committed
512
        Evaluation
Steven Cordwell's avatar
Steven Cordwell committed
513
        ----------
Steven Cordwell's avatar
Steven Cordwell committed
514
515
            PR(SxA)   = reward matrix
        """
Steven Cordwell's avatar
Steven Cordwell committed
516
517
518
        # we assume that P and R define a MDP i,e. assumption is that
        # check(P, R) has already been run and doesn't fail.
        
Steven Cordwell's avatar
Steven Cordwell committed
519
        # make P be an object array with (S, S) shaped array elements
520
        if (P.dtype is object):
Steven Cordwell's avatar
Steven Cordwell committed
521
522
523
524
525
526
527
528
529
530
531
            self.P = P
            self.A = self.P.shape[0]
            self.S = self.P[0].shape[0]
        else: # convert to an object array
            self.A = P.shape[0]
            self.S = P.shape[1]
            self.P = zeros(self.A, dtype=object)
            for aa in range(self.A):
                self.P[aa] = P[aa, :, :]
        
        # make R have the shape (S, A)
532
533
        if R.dtype is object:
            # R is object shaped (A,) with each element shaped (S, S)
Steven Cordwell's avatar
Steven Cordwell committed
534
            self.R = zeros((self.S, self.A))
535
536
537
538
539
540
            for aa in range(self.A):
                self.R[:, aa] = multiply(P[aa], R[aa]).sum(1)
        else:
            if R.ndim == 2:
                # R already has shape (S, A)
                self.R = R
Steven Cordwell's avatar
Steven Cordwell committed
541
            else:
542
543
                # R has shape (A, S, S)
                self.R = zeros((self.S, self.A))
Steven Cordwell's avatar
Steven Cordwell committed
544
                for aa in range(self.A):
Steven Cordwell's avatar
Steven Cordwell committed
545
                    self.R[:, aa] = multiply(P[aa], R[aa, :, :]).sum(1)
Steven Cordwell's avatar
Steven Cordwell committed
546
547
548
549
550
551
552
        
        # convert the arrays to numpy matrices
        for aa in range(self.A):
            if (type(self.P[aa]) is ndarray):
                self.P[aa] = matrix(self.P[aa])
        if (type(self.R) is ndarray):
            self.R = matrix(self.R)
553
554
555
556
557

    def iterate(self):
        """This is a placeholder method. Child classes should define their own
        iterate() method.
        """
Steven Cordwell's avatar
Steven Cordwell committed
558
        raise NotImplementedError("You should create an iterate() method.")
Steven Cordwell's avatar
Steven Cordwell committed
559
    
Steven Cordwell's avatar
Steven Cordwell committed
560
561
562
563
564
565
566
    def setSilent(self):
        """Ask for running resolution functions of the MDP Toolbox in silent
        mode.
        """
        self.verbose = False
    
    def setVerbose(self):
Steven Cordwell's avatar
Steven Cordwell committed
567
568
569
        """Ask for running resolution functions of the MDP Toolbox in verbose
        mode.
        """
Steven Cordwell's avatar
Steven Cordwell committed
570
        self.verbose = True
Steven Cordwell's avatar
Steven Cordwell committed
571
572

class FiniteHorizon(MDP):
Steven Cordwell's avatar
Steven Cordwell committed
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
    """Reolution of finite-horizon MDP with backwards induction
    
    Arguments
    ---------
    Let S = number of states, A = number of actions
    P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
    R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
    discount = discount factor, in ]0, 1]
    N        = number of periods, upper than 0
    h(S)     = terminal reward, optional (default [0; 0; ... 0] )
    Evaluation
    ----------
    V(S,N+1)     = optimal value function
                 V(:,n) = optimal value function at stage n
                        with stage in 1, ..., N
                        V(:,N+1) = value function for terminal stage 
    policy(S,N)  = optimal policy
                 policy(:,n) = optimal policy at stage n
                        with stage in 1, ...,N
                        policy(:,N) = policy for stage N
    cpu_time = used CPU time
  
    Notes
    -----
    In verbose mode, displays the current stage and policy transpose.

Steven Cordwell's avatar
Steven Cordwell committed
604
    """
Steven Cordwell's avatar
Steven Cordwell committed
605

Steven Cordwell's avatar
Steven Cordwell committed
606
607
    def __init__(self, transitions, reward, discount, N, h=None):
        """"""
Steven Cordwell's avatar
Steven Cordwell committed
608
        if N < 1:
Steven Cordwell's avatar
Steven Cordwell committed
609
            raise ValueError('PyMDPtoolbox: N must be greater than 0')
Steven Cordwell's avatar
Steven Cordwell committed
610
        else:
Steven Cordwell's avatar
Steven Cordwell committed
611
            self.N = N
Steven Cordwell's avatar
Steven Cordwell committed
612
        
Steven Cordwell's avatar
Steven Cordwell committed
613
        MDP.__init__(self, transitions, reward, discount, None)
Steven Cordwell's avatar
Steven Cordwell committed
614
        
Steven Cordwell's avatar
Steven Cordwell committed
615
        self.V = zeros(self.S, N + 1)
Steven Cordwell's avatar
Steven Cordwell committed
616
        
Steven Cordwell's avatar
Steven Cordwell committed
617
        if not h is None:
Steven Cordwell's avatar
Steven Cordwell committed
618
            self.V[:, N + 1] = h
Steven Cordwell's avatar
Steven Cordwell committed
619
        
Steven Cordwell's avatar
Steven Cordwell committed
620
621
    def iterate(self):
        """"""
Steven Cordwell's avatar
Steven Cordwell committed
622
623
        self.time = time()
        
Steven Cordwell's avatar
Steven Cordwell committed
624
        for n in range(self.N - 1):
Steven Cordwell's avatar
Steven Cordwell committed
625
626
            W, X = self.bellmanOperator(self.P, self.R, self.discount, self.V[:, self.N - n + 1])
            self.V[:, self.N - n] = W
Steven Cordwell's avatar
Steven Cordwell committed
627
628
629
            self.policy[:, self.N - n] = X
            if self.verbose:
                print("stage: %s ... policy transpose : %s") % (self.N - n, self.policy[:, self.N - n].T)
Steven Cordwell's avatar
Steven Cordwell committed
630
631
        
        self.time = time() - self.time
Steven Cordwell's avatar
Steven Cordwell committed
632
633

class LP(MDP):
Steven Cordwell's avatar
Steven Cordwell committed
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
    """Resolution of discounted MDP with linear programming

    Arguments
    ---------
    Let S = number of states, A = number of actions
    P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
    R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
    discount = discount rate, in ]0; 1[
    h(S)     = terminal reward, optional (default [0; 0; ... 0] )
    
    Evaluation
    ----------
    V(S)   = optimal values
    policy(S) = optimal policy
    cpu_time = used CPU time
    
    Notes    
    -----
    In verbose mode, displays the current stage and policy transpose.
    
    Examples
    --------
Steven Cordwell's avatar
Steven Cordwell committed
661
    """
Steven Cordwell's avatar
Steven Cordwell committed
662

Steven Cordwell's avatar
Steven Cordwell committed
663
664
    def __init__(self, transitions, reward, discount):
        """"""
Steven Cordwell's avatar
Steven Cordwell committed
665
666
667
        
        try:
            from cvxopt import matrix, solvers
Steven Cordwell's avatar
Steven Cordwell committed
668
            self.linprog = solvers.lp
Steven Cordwell's avatar
Steven Cordwell committed
669
        except ImportError:
Steven Cordwell's avatar
Steven Cordwell committed
670
671
            raise ImportError("The python module cvxopt is required to use " \
                "linear programming functionality.")
Steven Cordwell's avatar
Steven Cordwell committed
672
        
Steven Cordwell's avatar
Steven Cordwell committed
673
        from scipy.sparse import eye as speye
Steven Cordwell's avatar
Steven Cordwell committed
674

Steven Cordwell's avatar
Steven Cordwell committed
675
        MDP.__init__(self, transitions, reward, discount, None)
Steven Cordwell's avatar
Steven Cordwell committed
676
677
678
679
680
681
682
        
        # The objective is to resolve : min V / V >= PR + discount*P*V
        # The function linprog of the optimisation Toolbox of Mathworks resolves :
        # min f'* x / M * x <= b
        # So the objective could be expressed as : min V / (discount*P-I) * V <= - PR
        # To avoid loop on states, the matrix M is structured following actions M(A*S,S)
    
Steven Cordwell's avatar
Steven Cordwell committed
683
        self.f = ones(self.S, 1)
Steven Cordwell's avatar
Steven Cordwell committed
684
    
Steven Cordwell's avatar
Steven Cordwell committed
685
686
687
688
        self.M = zeros((self.A * self.S, self.S))
        for aa in range(self.A):
            pos = (aa + 1) * self.S
            self.M[(pos - self.S):pos, :] = discount * self.P[aa] - speye(self.S, self.S)
Steven Cordwell's avatar
Steven Cordwell committed
689
        
Steven Cordwell's avatar
Steven Cordwell committed
690
691
692
693
        self.M = matrix(self.M)
    
    def iterate(self):
        """"""
Steven Cordwell's avatar
Steven Cordwell committed
694
695
        self.time = time()
        
Steven Cordwell's avatar
Steven Cordwell committed
696
        self.V = self.linprog(self.f, self.M, -self.R)
Steven Cordwell's avatar
Steven Cordwell committed
697
        
Steven Cordwell's avatar
Steven Cordwell committed
698
        self.V, self.policy =  self.bellmanOperator(self.P, self.R, self.discount, self.V)
Steven Cordwell's avatar
Steven Cordwell committed
699
700
        
        self.time = time() - self.time
Steven Cordwell's avatar
Steven Cordwell committed
701
702
703

class PolicyIteration(MDP):
    """Resolution of discounted MDP with policy iteration algorithm.
Steven Cordwell's avatar
Steven Cordwell committed
704
    
Steven Cordwell's avatar
Steven Cordwell committed
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
    Arguments
    ---------
    Let S = number of states, A = number of actions
    P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
    R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
    discount = discount rate, in ]0, 1[
    policy0(S) = starting policy, optional 
    max_iter = maximum number of iteration to be done, upper than 0, 
             optional (default 1000)
    eval_type = type of function used to evaluate policy: 
             0 for mdp_eval_policy_matrix, else mdp_eval_policy_iterative
             optional (default 0)
             
    Evaluation
    ----------
    V(S)   = value function 
    policy(S) = optimal policy
    iter     = number of done iterations
    cpu_time = used CPU time
    
    Notes
    -----
    In verbose mode, at each iteration, displays the number 
    of differents actions between policy n-1 and n
    
Steven Cordwell's avatar
Steven Cordwell committed
735
736
737
738
739
740
741
    Examples
    --------
    >>> import mdp
    >>> P, R = mdp.exampleRand(5, 3)
    >>> pi = mdp.PolicyIteration(P, R, 0.9)
    >>> pi.iterate()
    
Steven Cordwell's avatar
Steven Cordwell committed
742
    """
Steven Cordwell's avatar
Steven Cordwell committed
743
    
Steven Cordwell's avatar
Steven Cordwell committed
744
    def __init__(self, transitions, reward, discount, policy0=None, max_iter=1000, eval_type=0):
Steven Cordwell's avatar
Steven Cordwell committed
745
746
        """"""
        
747
        MDP.__init__(self, transitions, reward, discount, max_iter)
Steven Cordwell's avatar
Steven Cordwell committed
748
        
Steven Cordwell's avatar
Steven Cordwell committed
749
750
751
        if policy0 == None:
            # initialise the policy to the one which maximises the expected
            # immediate reward
Steven Cordwell's avatar
Steven Cordwell committed
752
            self.V = matrix(zeros((self.S, 1)))
753
754
            self.policy, null = self.bellmanOperator()
            del null
Steven Cordwell's avatar
Steven Cordwell committed
755
756
757
758
759
760
761
762
763
        else:
            policy0 = array(policy0)
            
            if not policy0.shape in ((self.S, ), (self.S, 1), (1, self.S)):
                raise ValueError('PyMDPtolbox: policy0 must a vector with length S')
            
            policy0 = matrix(policy0.reshape(self.S, 1))
            
            if mod(policy0, 1).any() or (policy0 < 0).any() or (policy0 >= self.S).any():
Steven Cordwell's avatar
Steven Cordwell committed
764
                raise ValueError('PyMDPtoolbox: policy0 must be a vector of integers between 1 and S')
Steven Cordwell's avatar
Steven Cordwell committed
765
766
767
768
            else:
                self.policy = policy0
        
        # set or reset the initial values to zero
Steven Cordwell's avatar
Steven Cordwell committed
769
        self.V = matrix(zeros((self.S, 1)))
Steven Cordwell's avatar
Steven Cordwell committed
770
        
Steven Cordwell's avatar
Steven Cordwell committed
771
        if eval_type in (0, "matrix"):
772
            from numpy.linalg import solve
773
774
            from scipy.sparse import eye
            self.speye = eye
775
            self.lin_eq = solve
Steven Cordwell's avatar
Steven Cordwell committed
776
777
778
779
780
781
782
783
            self.eval_type = "matrix"
        elif eval_type in (1, "iterative"):
            self.eval_type = "iterative"
        else:
            raise ValueError("PyMDPtoolbox: eval_type should be 0 for matrix "\
                "evaluation or 1 for iterative evaluation. strings 'matrix' " \
                "and 'iterative' can also be used.")
    
Steven Cordwell's avatar
Steven Cordwell committed
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
    def computePpolicyPRpolicy(self):
        """Computes the transition matrix and the reward matrix for a policy
        
        Arguments
        ---------
        Let S = number of states, A = number of actions
        P(SxSxA)  = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
        R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
        policy(S) = a policy
        
        Evaluation
        ----------
        Ppolicy(SxS)  = transition matrix for policy
        PRpolicy(S)   = reward matrix for policy
        """
        Ppolicy = matrix(zeros((self.S, self.S)))
        Rpolicy = matrix(zeros((self.S, 1)))
        for aa in range(self.A): # avoid looping over S
        
            # the rows that use action a. .getA1() is used to make sure that
            # ind is a 1 dimensional vector
810
            ind = (self.policy == aa).nonzero()[0].getA1()
Steven Cordwell's avatar
Steven Cordwell committed
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
            if ind.size > 0: # if no rows use action a, then no point continuing
                Ppolicy[ind, :] = self.P[aa][ind, :]
                
                #PR = self.computePR() # an apparently uneeded line, and
                # perhaps harmful in this implementation c.f.
                # mdp_computePpolicyPRpolicy.m
                Rpolicy[ind] = self.R[ind, aa]
        
        # self.R cannot be sparse with the code in its current condition, but
        # it should be possible in the future. Also, if R is so big that its
        # a good idea to use a sparse matrix for it, then converting PRpolicy
        # from a dense to sparse matrix doesn't seem very memory efficient
        if type(self.R) is sparse:
            Rpolicy = sparse(Rpolicy)
        
        #self.Ppolicy = Ppolicy
        #self.Rpolicy = Rpolicy
        return (Ppolicy, Rpolicy)
    
Steven Cordwell's avatar
Steven Cordwell committed
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
    def evalPolicyIterative(self, V0=0, epsilon=0.0001, max_iter=10000):
        """Policy evaluation using iteration
        
        Arguments
        ---------
        Let S = number of states, A = number of actions
        P(SxSxA)  = transition matrix 
            P could be an array with 3 dimensions or 
            a cell array (1xS), each cell containing a matrix possibly sparse
        R(SxSxA) or (SxA) = reward matrix
            R could be an array with 3 dimensions (SxSxA) or 
            a cell array (1xA), each cell containing a sparse matrix (SxS) or
            a 2D array(SxA) possibly sparse  
        discount  = discount rate in ]0; 1[
        policy(S) = a policy
        V0(S)     = starting value function, optional (default : zeros(S,1))
        epsilon   = epsilon-optimal policy search, upper than 0,
            optional (default : 0.0001)
        max_iter  = maximum number of iteration to be done, upper than 0, 
            optional (default : 10000)
            
        Evaluation
        ----------
        Vpolicy(S) = value function, associated to a specific policy
        
        Notes
        -----
        In verbose mode, at each iteration, displays the condition which stopped iterations:
        epsilon-optimum value function found or maximum number of iterations reached.
        """
        if V0 == 0:
861
            policy_V = zeros((self.S, 1))
Steven Cordwell's avatar
Steven Cordwell committed
862
        else:
863
            raise NotImplementedError("evalPolicyIterative: case V0 != 0 not implemented. Use default (V0=0) instead.")
Steven Cordwell's avatar
Steven Cordwell committed
864
        
865
        policy_P, policy_R = self.computePpolicyPRpolicy()
Steven Cordwell's avatar
Steven Cordwell committed
866
867
868
        
        if self.verbose:
            print('  Iteration    V_variation')
869
        
Steven Cordwell's avatar
Steven Cordwell committed
870
871
872
873
        itr = 0
        done = False
        while not done:
            itr = itr + 1
874
875
876
877
878
            
            Vprev = policy_V
            policy_V = policy_R + self.discount * policy_P * Vprev
            
            variation = absolute(policy_V - Vprev).max()
Steven Cordwell's avatar
Steven Cordwell committed
879
880
            if self.verbose:
                print('      %s         %s') % (itr, variation)
881
                
Steven Cordwell's avatar
Steven Cordwell committed
882
883
884
            if variation < ((1 - self.discount) / self.discount) * epsilon: # to ensure |Vn - Vpolicy| < epsilon
                done = True
                if self.verbose:
885
                    print('PyMDPtoolbox: iterations stopped, epsilon-optimal value function')
Steven Cordwell's avatar
Steven Cordwell committed
886
887
888
            elif itr == max_iter:
                done = True
                if self.verbose:
889
                    print('PyMDPtoolbox: iterations stopped by maximum number of iteration condition')
Steven Cordwell's avatar
Steven Cordwell committed
890
        
Steven Cordwell's avatar
Steven Cordwell committed
891
        self.V = policy_V
892
893
    
    def evalPolicyMatrix(self):
Steven Cordwell's avatar
Steven Cordwell committed
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
        """Evaluation of the value function of a policy
        
        Arguments 
        ---------
        Let S = number of states, A = number of actions
        P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
        R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
        discount = discount rate in ]0; 1[
        policy(S) = a policy
        
        Evaluation
        ----------
        Vpolicy(S) = value function of the policy
        """
        
914
        Ppolicy, Rpolicy = self.computePpolicyPRpolicy()
Steven Cordwell's avatar
Steven Cordwell committed
915
        # V = PR + gPV  => (I-gP)V = PR  => V = inv(I-gP)* PR
Steven Cordwell's avatar
Steven Cordwell committed
916
        self.V = self.lin_eq((self.speye(self.S, self.S) - self.discount * Ppolicy) , Rpolicy)
Steven Cordwell's avatar
Steven Cordwell committed
917
918
    
    def iterate(self):
919
        """Run the policy iteration algorithm."""
Steven Cordwell's avatar
Steven Cordwell committed
920
        
921
922
923
        if self.verbose:
            print('  Iteration  Number_of_different_actions')
        
Steven Cordwell's avatar
Steven Cordwell committed
924
        done = False
925
        self.time = time()
Steven Cordwell's avatar
Steven Cordwell committed
926
927
        
        while not done:
928
929
            self.iter = self.iter + 1
            
Steven Cordwell's avatar
Steven Cordwell committed
930
931
            # these evalPolicy* functions will update the classes value
            # attribute
Steven Cordwell's avatar
Steven Cordwell committed
932
933
934
935
            if self.eval_type == "matrix":
                self.evalPolicyMatrix()
            elif self.eval_type == "iterative":
                self.evalPolicyIterative()
936
            
Steven Cordwell's avatar
Steven Cordwell committed
937
938
            # This should update the classes policy attribute but leave the
            # value alone
939
940
            policy_next, null = self.bellmanOperator()
            del null
Steven Cordwell's avatar
Steven Cordwell committed
941
            
942
            n_different = (policy_next != self.policy).sum()
943
944
945
946
            
            if self.verbose:
                print('       %s                 %s') % (self.iter, n_different)
            
947
            if n_different == 0:
Steven Cordwell's avatar
Steven Cordwell committed
948
                done = True
949
950
951
952
953
954
955
956
                if self.verbose:
                    print("...iterations stopped, unchanging policy found")
            elif (self.iter == self.max_iter):
                done = True 
                if self.verbose:
                    print("...iterations stopped by maximum number of iteration condition")
            else:
                self.policy = policy_next
Steven Cordwell's avatar
Steven Cordwell committed
957
        
958
959
        self.time = time() - self.time
        
Steven Cordwell's avatar
Steven Cordwell committed
960
        # store value and policy as tuples
Steven Cordwell's avatar
Steven Cordwell committed
961
        self.V = tuple(array(self.V).reshape(self.S).tolist())
Steven Cordwell's avatar
Steven Cordwell committed
962
        self.policy = tuple(array(self.policy).reshape(self.S).tolist())
Steven Cordwell's avatar
Steven Cordwell committed
963
964

class PolicyIterationModified(MDP):
Steven Cordwell's avatar
Steven Cordwell committed
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
    """Resolution of discounted MDP  with policy iteration algorithm
    
    Arguments
    ---------
    Let S = number of states, A = number of actions
    P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
    R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
    discount = discount rate, in ]0, 1[
    policy0(S) = starting policy, optional 
    max_iter = maximum number of iteration to be done, upper than 0, 
             optional (default 1000)
    eval_type = type of function used to evaluate policy: 
             0 for mdp_eval_policy_matrix, else mdp_eval_policy_iterative
             optional (default 0)
    
    Data Attributes
    ---------------
    V(S)   = value function 
    policy(S) = optimal policy
    iter     = number of done iterations
    cpu_time = used CPU time
    
    Notes
    -----
    In verbose mode, at each iteration, displays the number 
    of differents actions between policy n-1 and n
    
    Examples
    --------
    >>> import mdp
Steven Cordwell's avatar
Steven Cordwell committed
1000
    """
1001
1002
    
    def __init__(self, transitions, reward, discount, epsilon=0.01, max_iter=10):
Steven Cordwell's avatar
Steven Cordwell committed
1003
1004
        """"""
        
Steven Cordwell's avatar
Steven Cordwell committed
1005
        MDP.__init__(self, transitions, reward, discount, max_iter)
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
        
        if epsilon <= 0:
            raise ValueError("epsilon must be greater than 0")
        
        # computation of threshold of variation for V for an epsilon-optimal policy
        if self.discount != 1:
            self.thresh = epsilon * (1 - self.discount) / self.discount
        else:
            self.thresh = epsilon
        
        if discount == 1:
Steven Cordwell's avatar
Steven Cordwell committed
1017
            self.V = matrix(zeros((self.S, 1)))
1018
1019
        else:
            # min(min()) is not right
Steven Cordwell's avatar
Steven Cordwell committed
1020
            self.V = 1 / (1 - discount) * min(min(self.R)) * ones((self.S, 1))
Steven Cordwell's avatar
Steven Cordwell committed
1021
    
1022
1023
1024
1025
1026
1027
    def iterate(self):
        """"""
        
        if self.verbose:
            print('  Iteration  V_variation')
        
Steven Cordwell's avatar
Steven Cordwell committed
1028
        self.time = time()
Steven Cordwell's avatar
Steven Cordwell committed
1029
        
1030
1031
1032
        done = False
        while not done:
            self.iter = self.iter + 1
1033
            
1034
1035
            Vnext, policy = self.bellmanOperator(self.P, self.PR, self.discount, self.V)
            #[Ppolicy, PRpolicy] = mdp_computePpolicyPRpolicy(P, PR, policy);
1036
            
Steven Cordwell's avatar
Steven Cordwell committed
1037
            variation = getSpan(Vnext - self.V);
1038
1039
1040
            if self.verbose:
                print("      %s         %s" % (self.iter, variation))
            
Steven Cordwell's avatar
Steven Cordwell committed
1041
            self.V = Vnext
Steven Cordwell's avatar
Steven Cordwell committed
1042
            if variation < self.thresh:
1043
1044
1045
1046
                done = True
            else:
                is_verbose = False
                if self.verbose:
Steven Cordwell's avatar
Steven Cordwell committed
1047
                    self.setSilent
1048
1049
                    is_verbose = True
                
Steven Cordwell's avatar
Steven Cordwell committed
1050
                self.V = self.evalPolicyIterative()
1051
1052
                
                if is_verbose:
Steven Cordwell's avatar
Steven Cordwell committed
1053
                    self.setVerbose
1054
        
1055
        self.time = time() - self.time
Steven Cordwell's avatar
Steven Cordwell committed
1056
1057

class QLearning(MDP):
Steven Cordwell's avatar
Steven Cordwell committed
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
    """Evaluates the matrix Q, using the Q learning algorithm.
    
    Let S = number of states, A = number of actions
    
    Parameters
    ----------
    P : transition matrix (SxSxA)
        P could be an array with 3 dimensions or a cell array (1xA), each
        cell containing a sparse matrix (SxS)
    R : reward matrix(SxSxA) or (SxA)
        R could be an array with 3 dimensions (SxSxA) or a cell array
        (1xA), each cell containing a sparse matrix (SxS) or a 2D
        array(SxA) possibly sparse
    discount : discount rate
        in ]0; 1[    
    n_iter : number of iterations to execute (optional).
        Default value = 10000; it is an integer greater than the default value.
    
    Results
    -------
    Q : learned Q matrix (SxA) 
    
Steven Cordwell's avatar
Steven Cordwell committed
1080
    V : learned value function (S).
Steven Cordwell's avatar
Steven Cordwell committed
1081
1082
1083
1084
1085
1086
1087
    
    policy : learned optimal policy (S).
    
    mean_discrepancy : vector of V discrepancy mean over 100 iterations
        Then the length of this vector for the default value of N is 100 
        (N/100).

1088
    ExamplesPP[:, aa] = self.P[aa][:, ss]
Steven Cordwell's avatar
Steven Cordwell committed
1089
1090
1091
1092
1093
1094
1095
    ---------
    >>> import mdp
    >>> P, R = mdp.exampleForest()
    >>> ql = mdp.QLearning(P, R, 0.96)
    >>> ql.iterate()
    >>> ql.Q
    array([[  0.        ,   0.        ],
Steven Cordwell's avatar
Steven Cordwell committed
1096
1097
           [  0.01062959,   0.79870231],
           [ 10.08191776,   0.35309404]])
Steven Cordwell's avatar
Steven Cordwell committed
1098
    >>> ql.V
Steven Cordwell's avatar
Steven Cordwell committed
1099
    array([  0.        ,   0.79870231,  10.08191776])
Steven Cordwell's avatar
Steven Cordwell committed
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
    >>> ql.policy
    array([0, 1, 0])
    
    >>> import mdp
    >>> import numpy as np
    >>> P = np.array([[[0.5, 0.5],[0.8, 0.2]],[[0, 1],[0.1, 0.9]]])
    >>> R = np.array([[5, 10], [-1, 2]])
    >>> ql = mdp.QLearning(P, R, 0.9)
    >>> ql.iterate()
    >>> ql.Q
Steven Cordwell's avatar
Steven Cordwell committed
1110
1111
    array([[ 94.99525115,  99.99999007],
           [ 53.92930199,   5.57331205]])
Steven Cordwell's avatar
Steven Cordwell committed
1112
    >>> ql.V
Steven Cordwell's avatar
Steven Cordwell committed
1113
    array([ 99.99999007,  53.92930199])
Steven Cordwell's avatar
Steven Cordwell committed
1114
    >>> ql.policy
Steven Cordwell's avatar
Steven Cordwell committed
1115
1116
1117
    array([1, 0])
    >>> ql.time
    0.6501460075378418
Steven Cordwell's avatar
Steven Cordwell committed
1118
    
Steven Cordwell's avatar
Steven Cordwell committed
1119
1120
1121
1122
1123
1124
    """
    
    def __init__(self, transitions, reward, discount, n_iter=10000):
        """Evaluation of the matrix Q, using the Q learning algorithm
        """
        
1125
1126
1127
1128
        # The following check won't be done in MDP()'s initialisation, so let's
        # do it here
        if (n_iter < 10000):
            raise ValueError("PyMDPtoolbox: n_iter should be greater than 10000")
Steven Cordwell's avatar
Steven Cordwell committed
1129
        
1130
1131
        # after this n_iter will be known as self.max_iter
        MDP.__init__(self, transitions, reward, discount, n_iter)
Steven Cordwell's avatar
Steven Cordwell committed
1132
1133
1134
1135
1136
        
        # Initialisations
        self.Q = zeros((self.S, self.A))
        #self.dQ = zeros(self.S, self.A)
        self.mean_discrepancy = []
Steven Cordwell's avatar
Steven Cordwell committed
1137
        self.discrepancy = []
Steven Cordwell's avatar
Steven Cordwell committed
1138
1139
1140
1141
1142
1143
1144
        
    def iterate(self):
        """
        """
        self.time = time()
        
        # initial state choice
Steven Cordwell's avatar
Steven Cordwell committed
1145
        # s = randint(0, self.S - 1)
Steven Cordwell's avatar
Steven Cordwell committed
1146
        
1147
        for n in range(self.max_iter):
Steven Cordwell's avatar
Steven Cordwell committed
1148
1149
1150
            
            # Reinitialisation of trajectories every 100 transitions
            if ((n % 100) == 0):
Steven Cordwell's avatar
Steven Cordwell committed
1151
                s = randint(0, self.S - 1)
Steven Cordwell's avatar
Steven Cordwell committed
1152
1153
1154
1155
1156
            
            # Action choice : greedy with increasing probability
            # probability 1-(1/log(n+2)) can be changed
            pn = random()
            if (pn < (1 - (1 / log(n + 2)))):
Steven Cordwell's avatar
Steven Cordwell committed
1157
1158
                # optimal_action = self.Q[s, :].max()
                a = self.Q[s, :].argmax()
Steven Cordwell's avatar
Steven Cordwell committed
1159
1160
            else:
                a = randint(0, self.A - 1)
Steven Cordwell's avatar
Steven Cordwell committed
1161
            
Steven Cordwell's avatar
Steven Cordwell committed
1162
1163
1164
            # Simulating next state s_new and reward associated to <s,s_new,a>
            p_s_new = random()
            p = 0
Steven Cordwell's avatar
Steven Cordwell committed
1165
            s_new = -1
Steven Cordwell's avatar
Steven Cordwell committed
1166
1167
            while ((p < p_s_new) and (s_new < s)):
                s_new = s_new + 1
Steven Cordwell's avatar
Steven Cordwell committed
1168
1169
1170
                p = p + self.P[a][s, s_new]
            
            if (self.R.dtype == object):
Steven Cordwell's avatar
Steven Cordwell committed
1171
1172
                r = self.R[a][s, s_new]
            elif (self.R.ndim == 3):
Steven Cordwell's avatar
Steven Cordwell committed
1173
                r = self.R[a, s, s_new]
Steven Cordwell's avatar
Steven Cordwell committed
1174
            else:
Steven Cordwell's avatar
Steven Cordwell committed
1175
1176
                r = self.R[s, a]
            
Steven Cordwell's avatar
Steven Cordwell committed
1177
1178
            # Updating the value of Q   
            # Decaying update coefficient (1/sqrt(n+2)) can be changed
Steven Cordwell's avatar
Steven Cordwell committed
1179
            delta = r + self.discount * self.Q[s_new, :].max() - self.Q[s, a]
Steven Cordwell's avatar
Steven Cordwell committed
1180
1181
1182
1183
1184
1185
1186
            dQ = (1 / sqrt(n + 2)) * delta
            self.Q[s, a] = self.Q[s, a] + dQ
            
            # current state is updated
            s = s_new
            
            # Computing and saving maximal values of the Q variation
Steven Cordwell's avatar
Steven Cordwell committed
1187
            self.discrepancy.append(absolute(dQ))
Steven Cordwell's avatar
Steven Cordwell committed
1188
1189
1190
            
            # Computing means all over maximal Q variations values
            if ((n % 100) == 99):
Steven Cordwell's avatar
Steven Cordwell committed
1191
1192
                self.mean_discrepancy.append(mean(self.discrepancy))
                self.discrepancy = []
Steven Cordwell's avatar
Steven Cordwell committed
1193
1194
            
            # compute the value function and the policy
Steven Cordwell's avatar
Steven Cordwell committed
1195
            self.V = self.Q.max(axis=1)
Steven Cordwell's avatar
Steven Cordwell committed
1196
1197
            self.policy = self.Q.argmax(axis=1)
            
Steven Cordwell's avatar
Steven Cordwell committed
1198
        self.time = time() - self.time
1199
1200
1201
1202
        
        # rather than report that we have not done any iterations, assign the
        # value of n_iter to self.iter
        self.iter = self.max_iter
Steven Cordwell's avatar
Steven Cordwell committed
1203
1204

class RelativeValueIteration(MDP):
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
    """Resolution of MDP with average reward with relative value iteration 
    algorithm 
    
    Arguments
    ---------
    Let S = number of states, A = number of actions
    P(SxSxA) = transition matrix 
             P could be an array with 3 dimensions or 
             a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
    R(SxSxA) or (SxA) = reward matrix
             R could be an array with 3 dimensions (SxSxA) or 
             a cell array (1xA), each cell containing a sparse matrix (SxS) or
             a 2D array(SxA) possibly sparse  
    epsilon  = epsilon-optimal policy search, upper than 0, 
             optional (default: 0.01)
    max_iter = maximum number of iteration to be done, upper than 0,
             optional (default 1000)
    
    Evaluation
    ----------
    policy(S)       = epsilon-optimal policy
    average_reward  = average reward of the optimal policy
    cpu_time = used CPU time
    
    Notes
    -----
    In verbose mode, at each iteration, displays the span of U variation
    and the condition which stopped iterations : epsilon-optimum policy found
    or maximum number of iterations reached.
    
    Examples
    --------
Steven Cordwell's avatar
Steven Cordwell committed
1237
    """
1238
1239
1240
1241
1242
1243
1244
1245
    
    def __init__(self, transitions, reward, epsilon=0.01, max_iter=1000):
        
        MDP.__init__(self,  transitions, reward, None, max_iter)
        
        if epsilon <= 0:
            print('MDP Toolbox ERROR: epsilon must be upper than 0')
    
Steven Cordwell's avatar
Steven Cordwell committed
1246
1247
        self.U = zeros(self.S, 1)
        self.gain = self.U[self.S]
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
    
    def iterate(self):
        """"""
        
        done = False
        if self.verbose:
            print('  Iteration  U_variation')
        
        self.time = time()
        
        while not done:
            
            self.iter = self.iter + 1;
            
Steven Cordwell's avatar
Steven Cordwell committed
1262
            Unext, policy = self.bellmanOperator(self.P, self.R, 1, self.U)
1263
1264
            Unext = Unext - self.gain
            
1265
            variation = getSpan(Unext - self.U)
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
            
            if self.verbose:
                print("      %s         %s" % (self.iter, variation))

            if variation < self.epsilon:
                 done = True
                 average_reward = self.gain + min(Unext - self.U)
                 if self.verbose:
                     print('MDP Toolbox : iterations stopped, epsilon-optimal policy found')
            elif self.iter == self.max_iter:
                 done = True 
                 average_reward = self.gain + min(Unext - self.U);
                 if self.verbose:
                     print('MDP Toolbox : iterations stopped by maximum number of iteration condition')
            
            self.U = Unext
            self.gain = self.U(self.S)
        
        self.time = time() - self.time
Steven Cordwell's avatar
Steven Cordwell committed
1285
1286
1287

class ValueIteration(MDP):
    """
Steven Cordwell's avatar
Steven Cordwell committed
1288
1289
1290
1291
    Solves discounted MDP with the value iteration algorithm.
    
    Description
    -----------
Steven Cordwell's avatar
Steven Cordwell committed
1292
    mdp.ValueIteration applies the value iteration algorithm to solve
Steven Cordwell's avatar
Steven Cordwell committed
1293
    discounted MDP. The algorithm consists in solving Bellman's equation
Steven Cordwell's avatar
Steven Cordwell committed
1294
    iteratively.
Steven Cordwell's avatar
Steven Cordwell committed
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
    Iterating is stopped when an epsilon-optimal policy is found or after a
    specified number (max_iter) of iterations. 
    This function uses verbose and silent modes. In verbose mode, the function
    displays the variation of V (value function) for each iteration and the
    condition which stopped iterations: epsilon-policy found or maximum number
    of iterations reached.
    
    Let S = number of states, A = number of actions.
    
    Parameters
    ----------
    P : transition matrix 
        P could be a numpy ndarray with 3 dimensions (AxSxS) or a 
        numpy ndarray of dytpe=object with 1 dimenion (1xA), each 
        element containing a numpy ndarray (SxS) or scipy sparse matrix. 
    R : reward matrix
        R could be a numpy ndarray with 3 dimensions (AxSxS) or numpy
        ndarray of dtype=object with 1 dimension (1xA), each element
        containing a sparse matrix (SxS). R also could be a numpy 
        ndarray with 2 dimensions (SxA) possibly sparse.
    discount : discount rate
        Greater than 0, less than or equal to 1. Beware to check conditions of
        convergence for discount = 1.
    epsilon : epsilon-optimal policy search
        Greater than 0, optional (default: 0.01).
    max_iter : maximum number of iterations to be done
        Greater than 0, optional (default: computed)
    initial_value : starting value function
        optional (default: zeros(S,1)).
    
    Data Attributes
    ---------------
Steven Cordwell's avatar
Steven Cordwell committed
1327
    V : value function
Steven Cordwell's avatar
Steven Cordwell committed
1328
1329
        A vector which stores the optimal value function. Prior to calling the
        iterate() method it has a value of None. Shape is (S, ).
Steven Cordwell's avatar
Steven Cordwell committed
1330
    policy : epsilon-optimal policy
Steven Cordwell's avatar
Steven Cordwell committed
1331
1332
1333
        A vector which stores the optimal policy. Prior to calling the
        iterate() method it has a value of None. Shape is (S, ).
    iter : number of iterations taken to complete the computation
Steven Cordwell's avatar
Steven Cordwell committed
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
        An integer
    time : used CPU time
        A float
    
    Methods
    -------
    iterate()
        Starts the loop for the algorithm to be completed.
    setSilent()
        Sets the instance to silent mode.
    setVerbose()
        Sets the instance to verbose mode.
    
    Notes
    -----
    In verbose mode, at each iteration, displays the variation of V
    and the condition which stopped iterations: epsilon-optimum policy found
    or maximum number of iterations reached.
    
    Examples
    --------
    >>> import mdp
    >>> P, R = mdp.exampleForest()
    >>> vi = mdp.ValueIteration(P, R, 0.96)
    >>> vi.verbose
    False
    >>> vi.iterate()
Steven Cordwell's avatar
Steven Cordwell committed
1361
    >>> vi.V
Steven Cordwell's avatar
Steven Cordwell committed
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
    array([  5.93215488,   9.38815488,  13.38815488])
    >>> vi.policy
    array([0, 0, 0])
    >>> vi.iter
    4
    >>> vi.time
    0.002871990203857422
    
    >>> import mdp
    >>> import numpy as np
    >>> P = np.array([[[0.5, 0.5],[0.8, 0.2]],[[0, 1],[0.1, 0.9]]])
    >>> R = np.array([[5, 10], [-1, 2]])
    >>> vi = mdp.ValueIteration(P, R, 0.9)
    >>> vi.iterate()
Steven Cordwell's avatar
Steven Cordwell committed
1376
    >>> vi.V
Steven Cordwell's avatar
Steven Cordwell committed
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
    array([ 40.04862539,  33.65371176])
    >>> vi.policy
    array([1, 0])
    >>> vi.iter
    26
    >>> vi.time
    0.010202884674072266
    
    >>> import mdp
    >>> import numpy as np
    >>> from scipy.sparse import csr_matrix as sparse
    >>> P = np.zeros((2, ), dtype=object)
    >>> P[0] = sparse([[0.5, 0.5],[0.8, 0.2]])
    >>> P[1] = sparse([[0, 1],[0.1, 0.9]])
    >>> R = np.array([[5, 10], [-1, 2]])
    >>> vi = mdp.ValueIteration(P, R, 0.9)
    >>> vi.iterate()
Steven Cordwell's avatar
Steven Cordwell committed
1394
    >>> vi.V
Steven Cordwell's avatar
Steven Cordwell committed
1395
1396
1397
1398
    array([ 40.04862539,  33.65371176])
    >>> vi.policy
    array([1, 0])
    
Steven Cordwell's avatar
Steven Cordwell committed
1399
    """
Steven Cordwell's avatar
Steven Cordwell committed
1400
    
Steven Cordwell's avatar
Steven Cordwell committed
1401
    def __init__(self, transitions, reward, discount, epsilon=0.01, max_iter=1000, initial_value=0):
Steven Cordwell's avatar
Steven Cordwell committed
1402
        """Resolution of discounted MDP with value iteration algorithm."""
Steven Cordwell's avatar
Steven Cordwell committed
1403
        
1404
        MDP.__init__(self, transitions, reward, discount, max_iter)
Steven Cordwell's avatar
Steven Cordwell committed
1405
1406
1407
        
        # initialization of optional arguments
        if (initial_value == 0):