--- a/python/indicators.py	Fri May 25 18:15:18 2018 -0400
+++ b/python/indicators.py	Sun May 27 23:22:48 2018 -0400
@@ -6,7 +6,7 @@
 from matplotlib.pyplot import plot, ylim
 from matplotlib.pylab import find
 from numpy import array, arange, mean, floor, mean
-
+from scipy import percentile
 
 def multivariateName(indicatorNames):
     return '_'.join(indicatorNames)
@@ -35,7 +35,7 @@
             assert len(values) == timeInterval.length()
             self.timeInterval = timeInterval
             self.values = {}
-            for i in xrange(int(round(self.timeInterval.length()))):
+            for i in range(int(round(self.timeInterval.length()))):
                 self.values[self.timeInterval[i]] = values[i]
         self.maxValue = maxValue
 
@@ -60,7 +60,7 @@
         self.iterInstantNum = 0 # index in the interval or keys of the dict
         return self
 
-    def next(self):
+    def __next__(self):
         if self.iterInstantNum >= len(self.values):#(self.timeInterval and self.iterInstantNum>=self.timeInterval.length())\
            #     or (self.iterInstantNum >= self.values)
             raise StopIteration
@@ -160,14 +160,25 @@
         TemporalIndicator.__init__(self, name, values, timeInterval, maxValue)
         self.mostSevereIsMax = mostSevereIsMax
 
-    def getMostSevereValue(self, minNInstants=1): # TODO use np.percentile
-        values = array(self.values.values())
-        indices = range(len(values))
-        if len(indices) >= minNInstants:
-            values = sorted(values[indices], reverse = self.mostSevereIsMax) # inverted if most severe is max -> take the first values
-            return mean(values[:minNInstants])
+    def getMostSevereValue(self, minNInstants=1, centile=15.):
+        '''if there are more than minNInstants observations, 
+        returns either the average of these maximum values 
+        or if centile is not None the n% centile from the most severe value
+
+        eg for TTC, 15 returns the 15th centile (value such that 15% of observations are lower)'''
+        if self.__len__() < minNInstants:
+            return None
         else:
-            return None
+            values = list(self.values.values())
+            if centile is not None:
+                if self.mostSevereIsMax:
+                    c = 100-centile
+                else:
+                    c = centile
+                return percentile(values, c)
+            else:
+                values = sorted(values, reverse = self.mostSevereIsMax) # inverted if most severe is max -> take the first values
+                return mean(values[:minNInstants])
 
     def getInstantOfMostSevereValue(self):
         '''Returns the instant at which the indicator reaches its most severe value'''
@@ -189,15 +200,15 @@
 
     assert len(indicatorValues) == trajectory.length()
     indicatorMap = {}
-    for k in xrange(trajectory.length()):
+    for k in range(trajectory.length()):
         p = trajectory[k]
         i = floor(p.x/squareSize)
         j = floor(p.y/squareSize)
-        if indicatorMap.has_key((i,j)):
+        if (i,j) in indicatorMap:
             indicatorMap[(i,j)].append(indicatorValues[k])
         else:
             indicatorMap[(i,j)] = [indicatorValues[k]]
-    for k in indicatorMap.keys():
+    for k in indicatorMap:
         indicatorMap[k] = mean(indicatorMap[k])
     return indicatorMap
 
@@ -210,7 +221,7 @@
 #             points.append([x,y])
 #     inside = nx.points_inside_poly(array(points), polygon)
 #     indicatorMap = {}
-#     for i in xrange(len(inside)):
+#     for i in range(len(inside)):
 #         if inside[i]:
 #             indicatorMap[(floor(points[i][0]/squareSize), floor(points[i][1]/squareSize))] = 0
 #     return indicatorMap
@@ -229,12 +240,12 @@
     if more than one maps has a value)'''
     indicatorMap = {}
     for m in maps:
-        for k,v in m.iteritems():
-            if indicatorMap.has_key(k):
+        for k,v in m.items():
+            if k in indicatorMap:
                 indicatorMap[k].append(v)
             else:
                 indicatorMap[k] = [v]
-    for k in indicatorMap.keys():
+    for k in indicatorMap:
         indicatorMap[k] = combinationFunction(indicatorMap[k])
     return indicatorMap