mdp.py 32.9 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, ndarray, ones, zeros
Steven Cordwell's avatar
Steven Cordwell committed
36
from numpy.random import rand
Steven Cordwell's avatar
Steven Cordwell committed
37
from numpy.random import randint as randi
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
"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 "
    "greater than 1."
}

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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    
    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
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
    """
    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
177
           probability), optional (SxS) (default, random)
Steven Cordwell's avatar
Steven Cordwell committed
178
179
180
181
182
    
    Returns
    ----------
    P : transition probability matrix (SxSxA)
    R : reward matrix (SxSxA)
Steven Cordwell's avatar
Steven Cordwell committed
183
184
185
186
187
188

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

Steven Cordwell's avatar
Steven Cordwell committed
189
190
191
192
193
194
195
196
197
198
199
    """
    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
200
201
202
203
204
        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
205
206
207
    
    if is_sparse:
        # definition of transition matrix : square stochastic matrix
Steven Cordwell's avatar
Steven Cordwell committed
208
        P = zeros((A, ), dtype=object)
Steven Cordwell's avatar
Steven Cordwell committed
209
        # definition of reward matrix (values between -1 and +1)
Steven Cordwell's avatar
Steven Cordwell committed
210
        R = zeros((A, ), dtype=object)
Steven Cordwell's avatar
Steven Cordwell committed
211
        for a in range(A):
Steven Cordwell's avatar
Steven Cordwell committed
212
            PP = mask[a] * rand(S, S)
Steven Cordwell's avatar
Steven Cordwell committed
213
            for s in range(S):
Steven Cordwell's avatar
Steven Cordwell committed
214
215
                if (mask[a, s, :].sum() == 0):
                    PP[s, randint(0, S - 1)] = 1
Steven Cordwell's avatar
Steven Cordwell committed
216
217
218
219
220
221
222
223
224
                PP[s, :] = PP[s, :] / PP[s, :].sum()
            P[a] = sparse(PP)
            R[a] = sparse(mask * (2 * rand(S, S) - ones((S, S))))
    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
225
            P[a, :, :] = mask[a] * rand(S, S)
Steven Cordwell's avatar
Steven Cordwell committed
226
            for s in range(S):
Steven Cordwell's avatar
Steven Cordwell committed
227
228
                if (mask[a, s, :].sum() == 0):
                    P[a, s, randint(0, S - 1)] = 1
Steven Cordwell's avatar
Steven Cordwell committed
229
                P[a, s, :] = P[a, s, :] / P[a, s, :].sum()
Steven Cordwell's avatar
Steven Cordwell committed
230
            R[a, :, :] = mask[a] * (2 * rand(S, S) - ones((S, S), dtype=int))
Steven Cordwell's avatar
Steven Cordwell committed
231
232
233
    
    return (P, R)

Steven Cordwell's avatar
Steven Cordwell committed
234
class MDP(object):
Steven Cordwell's avatar
Steven Cordwell committed
235
    """The Markov Decision Problem Toolbox."""
Steven Cordwell's avatar
Steven Cordwell committed
236
237
238
239
240
    
    def __init__(self):
        """"""
        self.verbose = False
    
Steven Cordwell's avatar
Steven Cordwell committed
241
    def bellmanOperator(self):
Steven Cordwell's avatar
Steven Cordwell committed
242
243
        """
        Applies the Bellman operator on the value function.
Steven Cordwell's avatar
Steven Cordwell committed
244
        
Steven Cordwell's avatar
Steven Cordwell committed
245
246
        Updates the value function and the Vprev-improving policy.
        
247
        Returns
Steven Cordwell's avatar
Steven Cordwell committed
248
        -------
249
        (policy, value) : tuple of new policy and its value
Steven Cordwell's avatar
Steven Cordwell committed
250
251
252
        """
        Q = matrix(zeros((self.S, self.A)))
        for aa in range(self.A):
Steven Cordwell's avatar
Steven Cordwell committed
253
            Q[:, aa] = self.R[:, aa] + (self.discount * self.P[aa] * self.value)
Steven Cordwell's avatar
Steven Cordwell committed
254
        
255
256
257
258
259
        # Which way is better? if choose the first way, then the classes that
        # call this function must be changed
        # 1. Return, (policy, value)
        # return (Q.argmax(axis=1), Q.max(axis=1))
        # 2. change self.policy and self.value directly
Steven Cordwell's avatar
Steven Cordwell committed
260
261
262
263
        self.value = Q.max(axis=1)
        self.policy = Q.argmax(axis=1)
    
    def check(self, P, R):
Steven Cordwell's avatar
Steven Cordwell committed
264
        """Checks if the matrices P and R define a Markov Decision Process.
Steven Cordwell's avatar
Steven Cordwell committed
265
        
Steven Cordwell's avatar
Steven Cordwell committed
266
267
268
269
270
        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.
Steven Cordwell's avatar
Steven Cordwell committed
271
272
273
274
275
276
277
278
279
280
        
        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  
Steven Cordwell's avatar
Steven Cordwell committed
281
282
283
284
285
286
287
288
289
        """
        
        # 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):
Steven Cordwell's avatar
Steven Cordwell committed
290
            raise TypeError(mdperr["P_type"])
Steven Cordwell's avatar
Steven Cordwell committed
291
        
Steven Cordwell's avatar
Steven Cordwell committed
292
        if (not type(R) is ndarray):
Steven Cordwell's avatar
Steven Cordwell committed
293
            raise TypeError(mdperr["R_type"])
Steven Cordwell's avatar
Steven Cordwell committed
294
295
296
297
298
299
            
        if (P.dtype == object):
            if (P.ndim > 1):
                raise ValueError(mdperr["obj_shape"])
            else:
                P_is_object = True
Steven Cordwell's avatar
Steven Cordwell committed
300
        else:
Steven Cordwell's avatar
Steven Cordwell committed
301
302
303
304
            if (P.ndim != 3):
                raise ValueError(mdperr["P_shape"])
            else:
                P_is_object = False
Steven Cordwell's avatar
Steven Cordwell committed
305
            
Steven Cordwell's avatar
Steven Cordwell committed
306
307
308
309
310
        if (R.dtype == object):
            if (R.ndim > 1):
                raise ValueError(mdperr["obj_shape"])
            else:
                R_is_object = True
Steven Cordwell's avatar
Steven Cordwell committed
311
        else:
Steven Cordwell's avatar
Steven Cordwell committed
312
313
314
315
            if (not R.ndim in (2, 3)):
                raise ValueError(mdperr["R_shape"])
            else:
                R_is_object = False
Steven Cordwell's avatar
Steven Cordwell committed
316
317
318
319
320
        
        if P_is_object:
            aP = P.shape[0]
            sP0 = P[0].shape[0]
            sP1 = P[0].shape[1]
Steven Cordwell's avatar
Steven Cordwell committed
321
            # check to see that the other object array elements are the same shape
Steven Cordwell's avatar
Steven Cordwell committed
322
323
324
325
            for aa in range(1, aP):
                sP0aa = P[aa].shape[0]
                sP1aa = P[aa].shape[1]
                if ((sP0aa != sP0) or (sP1aa != sP1)):
Steven Cordwell's avatar
Steven Cordwell committed
326
                    raise ValueError(mdperr["obj_square"])
Steven Cordwell's avatar
Steven Cordwell committed
327
328
329
330
        else:
            aP, sP0, sP1 = P.shape
        
        if ((sP0 < 1) or (aP < 1) or (sP0 != sP1)):
Steven Cordwell's avatar
Steven Cordwell committed
331
332
333
334
335
            raise ValueError(mdperr["P_shape"])
        
        for aa in range(aP):
            if P_is_object:
                self.checkSquareStochastic(P[aa])
Steven Cordwell's avatar
Steven Cordwell committed
336
            else:
Steven Cordwell's avatar
Steven Cordwell committed
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
                self.checkSquareStochastic(P[aa, :, :])
            aa = aa + 1
        
        if R_is_object:
            aR = R.shape[0]
            sR0 = R[0].shape[0]
            sR1 = R[0].shape[1]
            # check to see that the other object array elements are the same shape
            for aa in range(1, aR):
                sR0aa = R[aa].shape[0]
                sR1aa = R[aa].shape[1]
                if ((sR0aa != sR0) or (sR1aa != sR1)):
                    raise ValueError(mdperr["obj_square"])
        elif (R.ndim == 3):
            aR, sR0, sR1 = R.shape
        else:
            sR0, aR = R.shape
            sR1 = sR0
Steven Cordwell's avatar
Steven Cordwell committed
355
        
Steven Cordwell's avatar
Steven Cordwell committed
356
357
358
359
360
        if ((sR0 < 1) or (aR < 1) or (sR0 != sR1)):
            raise ValueError(mdperr["R_shape"])
                
        if (sP0 != sR0) or (aP != aR):
            raise ValueError(mdperr["PR_incompat"])
Steven Cordwell's avatar
Steven Cordwell committed
361
362
363
364
    
    def checkSquareStochastic(self, Z):
        """Check if Z is a square stochastic matrix
        
Steven Cordwell's avatar
Steven Cordwell committed
365
366
        Arguments
        --------------------------------------------------------------
Steven Cordwell's avatar
Steven Cordwell committed
367
            Z = a numpy ndarray SxS, possibly sparse (csr_matrix)
Steven Cordwell's avatar
Steven Cordwell committed
368
369
        Evaluation
        -------------------------------------------------------------
Steven Cordwell's avatar
Steven Cordwell committed
370
371
372
373
            error_msg = error message or None if correct
        """
        s1, s2 = Z.shape
        if (s1 != s2):
Steven Cordwell's avatar
Steven Cordwell committed
374
           raise ValueError(mdperr["mat_square"])
Steven Cordwell's avatar
Steven Cordwell committed
375
        elif (absolute(Z.sum(axis=1) - ones(s2))).max() > 10**(-12):
Steven Cordwell's avatar
Steven Cordwell committed
376
377
378
379
380
            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"]) 
Steven Cordwell's avatar
Steven Cordwell committed
381
382
        else:
            return(None)
Steven Cordwell's avatar
Steven Cordwell committed
383
    
Steven Cordwell's avatar
Steven Cordwell committed
384
385
386
387
    def computePpolicyPRpolicy(self):
        """Computes the transition matrix and the reward matrix for a policy.
        """
        pass    
Steven Cordwell's avatar
Steven Cordwell committed
388
389
390
391
    
    def computePR(self, P, R):
        """Computes the reward for the system in one state chosing an action
        
Steven Cordwell's avatar
Steven Cordwell committed
392
        Arguments
Steven Cordwell's avatar
Steven Cordwell committed
393
        ---------
Steven Cordwell's avatar
Steven Cordwell committed
394
395
396
397
398
399
400
401
        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
402
        Evaluation
Steven Cordwell's avatar
Steven Cordwell committed
403
        ----------
Steven Cordwell's avatar
Steven Cordwell committed
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
            PR(SxA)   = reward matrix
        """
        # make P be an object array with (S, S) shaped array elements
        if (P.dtype == object):
            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)
        if ((R.ndim == 2) and (not R.dtype is object)):
            # R already has shape (S, A)
            self.R = R
        else: 
            # R has shape (A, S, S) or object shaped (A,) with each element
            # shaped (S, S)
            self.R = zeros((self.S, self.A))
            if (R.dtype is object):
                for aa in range(self.A):
                    self.R[:, aa] = sum(P[aa] * R[aa], 2)
            else:
                for aa in range(self.A):
                    self.R[:, aa] = sum(P[aa] * R[aa, :, :], 2)
        
        # 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)
    
Steven Cordwell's avatar
Steven Cordwell committed
440
    def getSpan(self, W):
Steven Cordwell's avatar
Steven Cordwell committed
441
442
443
444
445
        """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
446
    
Steven Cordwell's avatar
Steven Cordwell committed
447
448
449
450
451
452
453
    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
454
455
456
        """Ask for running resolution functions of the MDP Toolbox in verbose
        mode.
        """
Steven Cordwell's avatar
Steven Cordwell committed
457
        self.verbose = True
Steven Cordwell's avatar
Steven Cordwell committed
458
459
460
461
462
463
464
465
466
467
468
469
470

class FiniteHorizon(MDP):
    """Resolution of finite-horizon MDP with backwards induction.
    """
    pass

class LP(MDP):
    """Resolution of discounted MDP with linear programming.
    """
    pass

class PolicyIteration(MDP):
    """Resolution of discounted MDP with policy iteration algorithm.
Steven Cordwell's avatar
Steven Cordwell committed
471
472
473
474
475
476
477
478
    
    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
479
    """
Steven Cordwell's avatar
Steven Cordwell committed
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
    
    def __init__(self, transitions, reward, discount, epsilon=0.01, max_iter=1000, initial_value=0):
        """"""
        MDP.__init__(self)
        
        self.check(transitions, reward)
        
        self.S = transitions.shape[1]
        
        self.A = transitions.shape[0]
        
        self.P = transitions
        
        self.R = reward
        
        #self.computePR(transitions, reward)
        
        if (initial_value == 0):
            self.value = zeros((self.S))
            #self.value = matrix(zeros((self.S, 1)))
        else:
            if (len(initial_value) != self.S):
                raise TypeError("The initial value must be length S")
            
            self.value = matrix(initial_value)
        
        self.policy = randi(0, self.A, self.S)

        self.discount = discount
        if (discount < 1):
            # compute a bound for the number of iterations
            #self.max_iter = self.boundIter(epsilon)
            self.max_iter = 5000
            # computation of threshold of variation for V for an epsilon-optimal policy
            self.thresh = epsilon * (1 - self.discount) / self.discount
        else: # discount == 1
            # bound for the number of iterations
            self.max_iter = max_iter
            # threshold of variation for V for an epsilon-optimal policy
            self.thresh = epsilon 
        
        self.iter = 0
        
        self.time = None
    
    def iterate(self):
        """"""
        done = False
        stop_criterion = 0.01
        
        while not done:
            stop = False
            while not stop:
                change = 0
                for s in range(self.S):
                    v = self.value[s]
                    a = self.policy[s]
                    self.value[s] = (self.P[a, s, :] * (self.R[a, s, :] +
                        (self.discount * self.value))).sum()
                    change = max(change, abs(v - self.value[s]))
                
                if change < stop_criterion:
                    stop = True
            
            policy_stable = True
            for s in range(self.S):
                b = self.policy[s]
                self.policy[s] = (self.P[:, s, :] * (self.R[:, s, :] +
                        (self.discount * self.value))).sum(1).argmax()
                if b !=  self.policy[s]:
                    policy_stable = False
            
            if policy_stable:
                done = True
        
        # store value and policy as tuples
        self.value = tuple(array(self.value).reshape(self.S).tolist())
        self.policy = tuple(array(self.policy).reshape(self.S).tolist())
Steven Cordwell's avatar
Steven Cordwell committed
558
559
560
561
562
563
564

class PolicyIterationModified(MDP):
    """Resolution of discounted MDP with modified policy iteration algorithm.
    """
    pass

class QLearning(MDP):
Steven Cordwell's avatar
Steven Cordwell committed
565
566
567
568
569
570
571
572
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
    """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) 
    
    value : learned value function (S).
    
    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).

    Examples
    ---------
    >>> 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
603
604
           [  0.01062959,   0.79870231],
           [ 10.08191776,   0.35309404]])
Steven Cordwell's avatar
Steven Cordwell committed
605
    >>> ql.value
Steven Cordwell's avatar
Steven Cordwell committed
606
    array([  0.        ,   0.79870231,  10.08191776])
Steven Cordwell's avatar
Steven Cordwell committed
607
608
609
610
611
612
613
614
615
616
    >>> 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
617
618
    array([[ 94.99525115,  99.99999007],
           [ 53.92930199,   5.57331205]])
Steven Cordwell's avatar
Steven Cordwell committed
619
    >>> ql.value
Steven Cordwell's avatar
Steven Cordwell committed
620
    array([ 99.99999007,  53.92930199])
Steven Cordwell's avatar
Steven Cordwell committed
621
    >>> ql.policy
Steven Cordwell's avatar
Steven Cordwell committed
622
623
624
    array([1, 0])
    >>> ql.time
    0.6501460075378418
Steven Cordwell's avatar
Steven Cordwell committed
625
    
Steven Cordwell's avatar
Steven Cordwell committed
626
627
628
629
630
631
    """
    
    def __init__(self, transitions, reward, discount, n_iter=10000):
        """Evaluation of the matrix Q, using the Q learning algorithm
        """
        
Steven Cordwell's avatar
Steven Cordwell committed
632
        MDP.__init__(self)
Steven Cordwell's avatar
Steven Cordwell committed
633
        
Steven Cordwell's avatar
Steven Cordwell committed
634
635
636
637
638
639
        # Check of arguments
        if (discount <= 0) or (discount >= 1):
            raise ValueError("MDP Toolbox Error: Discount rate must be in ]0,1[")
        elif (n_iter < 10000):
            raise ValueError("MDP Toolbox Error: n_iter must be greater than 10000")
        
Steven Cordwell's avatar
Steven Cordwell committed
640
        self.check(transitions, reward)
Steven Cordwell's avatar
Steven Cordwell committed
641
642
643
644
645
646
647
648
649
650
651
        
        self.computePR(transitions, reward)
        
        self.discount = discount
        
        self.n_iter = n_iter
        
        # 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
652
        self.discrepancy = []
Steven Cordwell's avatar
Steven Cordwell committed
653
654
655
656
657
658
659
660
661
        
        self.time = None
        
    def iterate(self):
        """
        """
        self.time = time()
        
        # initial state choice
Steven Cordwell's avatar
Steven Cordwell committed
662
        # s = randint(0, self.S - 1)
Steven Cordwell's avatar
Steven Cordwell committed
663
664
665
666
667
        
        for n in range(self.n_iter):
            
            # Reinitialisation of trajectories every 100 transitions
            if ((n % 100) == 0):
Steven Cordwell's avatar
Steven Cordwell committed
668
                s = randint(0, self.S - 1)
Steven Cordwell's avatar
Steven Cordwell committed
669
670
671
672
673
            
            # 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
674
675
                # optimal_action = self.Q[s, :].max()
                a = self.Q[s, :].argmax()
Steven Cordwell's avatar
Steven Cordwell committed
676
677
            else:
                a = randint(0, self.A - 1)
Steven Cordwell's avatar
Steven Cordwell committed
678
            
Steven Cordwell's avatar
Steven Cordwell committed
679
680
681
            # 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
682
            s_new = -1
Steven Cordwell's avatar
Steven Cordwell committed
683
684
            while ((p < p_s_new) and (s_new < s)):
                s_new = s_new + 1
Steven Cordwell's avatar
Steven Cordwell committed
685
686
687
                p = p + self.P[a][s, s_new]
            
            if (self.R.dtype == object):
Steven Cordwell's avatar
Steven Cordwell committed
688
689
                r = self.R[a][s, s_new]
            elif (self.R.ndim == 3):
Steven Cordwell's avatar
Steven Cordwell committed
690
                r = self.R[a, s, s_new]
Steven Cordwell's avatar
Steven Cordwell committed
691
            else:
Steven Cordwell's avatar
Steven Cordwell committed
692
693
                r = self.R[s, a]
            
Steven Cordwell's avatar
Steven Cordwell committed
694
695
            # Updating the value of Q   
            # Decaying update coefficient (1/sqrt(n+2)) can be changed
Steven Cordwell's avatar
Steven Cordwell committed
696
            delta = r + self.discount * self.Q[s_new, :].max() - self.Q[s, a]
Steven Cordwell's avatar
Steven Cordwell committed
697
698
699
700
701
702
703
            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
704
            self.discrepancy.append(absolute(dQ))
Steven Cordwell's avatar
Steven Cordwell committed
705
706
707
            
            # Computing means all over maximal Q variations values
            if ((n % 100) == 99):
Steven Cordwell's avatar
Steven Cordwell committed
708
709
                self.mean_discrepancy.append(mean(self.discrepancy))
                self.discrepancy = []
Steven Cordwell's avatar
Steven Cordwell committed
710
711
712
713
714
            
            # compute the value function and the policy
            self.value = self.Q.max(axis=1)
            self.policy = self.Q.argmax(axis=1)
            
Steven Cordwell's avatar
Steven Cordwell committed
715
        self.time = time() - self.time
Steven Cordwell's avatar
Steven Cordwell committed
716
717
718

class RelativeValueIteration(MDP):
    """Resolution of MDP with average reward with relative value iteration
Steven Cordwell's avatar
Steven Cordwell committed
719
720
    algorithm.
    """
Steven Cordwell's avatar
Steven Cordwell committed
721
722
723
724
    pass

class ValueIteration(MDP):
    """
Steven Cordwell's avatar
Steven Cordwell committed
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
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
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
    Solves discounted MDP with the value iteration algorithm.
    
    Description
    -----------
    mdp_value_iteration applies the value iteration algorithm to solve
    discounted MDP. The algorithm consists in solving Bellman's equation
    iteratively. 
    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
    ---------------
    value : value function
        A vector which stores the optimal value function. It exists only after
        the iterate() method has been called. Shape is (S, ).
    policy : epsilon-optimal policy
        A vector which stores the optimal policy. It exists only after
        the iterate() method has been called. Shape is (S, ).
    iter : number of done iterations
        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()
    >>> vi.value
    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()
    >>> vi.value
    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()
    >>> vi.value
    array([ 40.04862539,  33.65371176])
    >>> vi.policy
    array([1, 0])
    
Steven Cordwell's avatar
Steven Cordwell committed
837
    """
Steven Cordwell's avatar
Steven Cordwell committed
838
    
Steven Cordwell's avatar
Steven Cordwell committed
839
    def __init__(self, transitions, reward, discount, epsilon=0.01, max_iter=1000, initial_value=0):
Steven Cordwell's avatar
Steven Cordwell committed
840
        """Resolution of discounted MDP with value iteration algorithm."""
Steven Cordwell's avatar
Steven Cordwell committed
841
        
Steven Cordwell's avatar
Steven Cordwell committed
842
        MDP.__init__(self)
Steven Cordwell's avatar
Steven Cordwell committed
843
        
Steven Cordwell's avatar
Steven Cordwell committed
844
        self.check(transitions, reward)
Steven Cordwell's avatar
Steven Cordwell committed
845
        
Steven Cordwell's avatar
Steven Cordwell committed
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
        self.computePR(transitions, reward)
        
        # initialization of optional arguments
        if (initial_value == 0):
            self.value = matrix(zeros((self.S, 1)))
        else:
            if (initial_value.size != self.S):
                raise TypeError("The initial value must be length S")
            
            self.value = matrix(initial_value)
        
        self.discount = discount
        if (discount < 1):
            # compute a bound for the number of iterations
            #self.max_iter = self.boundIter(epsilon)
            self.max_iter = 5000
            # computation of threshold of variation for V for an epsilon-optimal policy
            self.thresh = epsilon * (1 - self.discount) / self.discount
        else: # discount == 1
            # bound for the number of iterations
            self.max_iter = max_iter
            # threshold of variation for V for an epsilon-optimal policy
            self.thresh = epsilon 
        
Steven Cordwell's avatar
Steven Cordwell committed
870
        self.iter = 0
Steven Cordwell's avatar
Steven Cordwell committed
871
872
        
        self.time = None
Steven Cordwell's avatar
Steven Cordwell committed
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
    
    def boundIter(self, epsilon):
        """Computes a bound for the number of iterations for the value iteration
        algorithm to find an epsilon-optimal policy with use of span for the 
        stopping criterion
        
        Arguments --------------------------------------------------------------
        Let S = number of states, A = number of actions
            epsilon   = |V - V*| < epsilon,  upper than 0,
                optional (default : 0.01)
        Evaluation -------------------------------------------------------------
            max_iter  = bound of the number of iterations for the value 
            iteration algorithm to find an epsilon-optimal policy with use of
            span for the stopping criterion
            cpu_time  = used CPU time
        """
        # See Markov Decision Processes, M. L. Puterman, 
        # Wiley-Interscience Publication, 1994 
        # p 202, Theorem 6.6.6
        # k =    max     [1 - S min[ P(j|s,a), p(j|s',a')] ]
        #     s,a,s',a'       j
        k = 0
        h = zeros(self.S)
        
        for ss in range(self.S):
            PP = zeros((self.S, self.A))
            for aa in range(self.A):
                PP[:, aa] = self.P[aa][:, ss]
            # the function "min()" without any arguments finds the
            # minimum of the entire array.
            h[ss] = PP.min()
        
        k = 1 - h.sum()
        V1 = self.bellmanOperator(self.value)
        # p 201, Proposition 6.6.5
Steven Cordwell's avatar
Steven Cordwell committed
908
        max_iter = log( (epsilon * (1 - self.discount) / self.discount) / self.getSpan(V1 - self.value) ) / log(self.discount * k)
Steven Cordwell's avatar
Steven Cordwell committed
909
910
911
912
913
914
915
916
        return ceil(max_iter)
    
    def iterate(self):
        """
        """
        self.time = time()
        done = False
        while not done:
Steven Cordwell's avatar
Steven Cordwell committed
917
            self.iter = self.iter + 1
Steven Cordwell's avatar
Steven Cordwell committed
918
919
920
921
            
            Vprev = self.value
            
            # Bellman Operator: updates "self.value" and "self.policy"
Steven Cordwell's avatar
Steven Cordwell committed
922
            self.bellmanOperator()
Steven Cordwell's avatar
Steven Cordwell committed
923
            
Steven Cordwell's avatar
Steven Cordwell committed
924
925
            # The values, based on Q. For the function "max()": the option
            # "axis" means the axis along which to operate. In this case it
Steven Cordwell's avatar
Steven Cordwell committed
926
            # finds the maximum of the the rows. (Operates along the columns?)
Steven Cordwell's avatar
Steven Cordwell committed
927
            variation = self.getSpan(self.value - Vprev)
Steven Cordwell's avatar
Steven Cordwell committed
928
929
            
            if self.verbose:
Steven Cordwell's avatar
Steven Cordwell committed
930
                print("      %s         %s" % (self.iter, variation))
Steven Cordwell's avatar
Steven Cordwell committed
931
932
933
934
935
            
            if variation < self.thresh:
                done = True
                if self.verbose:
                    print("...iterations stopped, epsilon-optimal policy found")
Steven Cordwell's avatar
Steven Cordwell committed
936
            elif (self.iter == self.max_iter):
Steven Cordwell's avatar
Steven Cordwell committed
937
938
939
                done = True 
                if self.verbose:
                    print("...iterations stopped by maximum number of iteration condition")
Steven Cordwell's avatar
Steven Cordwell committed
940
941
942
943
944
        
        # store value and policy as tuples
        self.value = tuple(array(self.value).reshape(self.S).tolist())
        self.policy = tuple(array(self.policy).reshape(self.S).tolist())
        
Steven Cordwell's avatar
Steven Cordwell committed
945
946
        self.time = time() - self.time

Steven Cordwell's avatar
Steven Cordwell committed
947
948
949
950
class ValueIterationGS(MDP):
    """Resolution of discounted MDP with value iteration Gauss-Seidel algorithm.
    """
    pass