--- a/python/moving.py	Mon Jun 04 10:25:30 2018 -0400
+++ b/python/moving.py	Mon Jun 04 11:25:49 2018 -0400
@@ -10,6 +10,7 @@
 from scipy.stats import scoreatpercentile
 from scipy.spatial.distance import cdist
 from scipy.signal import savgol_filter
+import copy
 
 try:
     from shapely.geometry import Polygon, Point as shapelyPoint
@@ -1155,7 +1156,47 @@
         return MovingObject(num = num, timeInterval = timeInterval, positions = positions, velocities = velocities)
 
     @staticmethod
-    def concatenate(obj1, obj2, num = None, computePositions = False):
+    def generateWithFeatures(num, features, userType):
+    	newObjectTMP = MovingObject.concatenateWith(features[0],features[1:],num)
+    	newObject = MovingObject(newObjectTMP.getNum() if num is None else num, newObjectTMP.getTimeInterval(), newObjectTMP.positions, newObjectTMP.velocities, userType = userType)
+    	newObject.features = [copy.deepcopy(f) for f in features]
+    	return newObject
+
+    @staticmethod
+    def concatenateWith(obj1, ObjectList=[], num = None): #Concatenate with a list of MovingObjects
+    	if len(ObjectList)==0:
+    		newObject = MovingObject(obj1.getNum() if num is None else num, obj1.getTimeInterval(), obj1.positions, obj1.velocities, userType = obj1.getUserType())
+    		if obj1.hasFeatures() :
+    			newObject.features = list(obj1.features)
+    	else:
+    		def getKey(item):
+    			return 0 if not obj1.commonTimeInterval(item).empty() else max(obj1.getFirstInstant(),item.getFirstInstant())-min(obj1.getLastInstant(),item.getLastInstant())
+    		Order = min(ObjectList , key=getKey)
+    		newObject = MovingObject.concatenate(obj1, Order, num)
+    		ObjectList.remove(Order)
+    		newObject = MovingObject.concatenateWith(newObject,ObjectList)
+    	return newObject
+
+    def updatePosition(self):
+    	if self.features is not None:
+	    	positions = Trajectory()
+	    	for t in self.getTimeInterval():
+	    		nTotal = 0.
+	    		p = Point(0.,0.)
+	    		for obj in self.features:
+	    			if obj.existsAtInstant(t):
+	    				if obj.hasFeatures():
+	    					n = len([f for f in obj.getFeatures() if f.existsAtInstant(t)])
+	    				else:
+	    					n = 1.
+	    				p += obj.getPositionAtInstant(t).__mul__(n)
+	    				nTotal += n
+	    		assert nTotal>0, 'there should be at least one point for each instant'
+	    		positions.addPosition(p.divide(nTotal))
+	    	self.positions = positions
+
+    @staticmethod
+    def concatenate(obj1, obj2, num = None):
         '''Concatenates two objects supposed to overlap temporally '''
         if num is None:
             newNum = obj1.getNum()
@@ -1181,14 +1222,29 @@
             	velocities.addPositionXY(vitessex,vitessey)
             	px+=vitessex
             	py+=vitessey
-
             newObject = MovingObject(newNum, emptyInterval, positions, velocities, userType = obj1.getUserType())
+            newObject.features = [MovingObject(newNum, emptyInterval, positions, velocities, userType = obj1.getUserType())] #In case there is features to add when we recursively call concatenate 
             return MovingObject.concatenate(MovingObject.concatenate(obj1, newObject),obj2)
         else:
             newTimeInterval = TimeInterval.union(obj1.getTimeInterval(), obj2.getTimeInterval())
             # positions
-            if computePositions: # TODO it would probably be better to recompute from all features, if features are available, in other method
-                positions = Trajectory()
+            positions = Trajectory()
+            for t in newTimeInterval:
+            	nTotal = 0.
+            	p = Point(0.,0.)
+            	for obj in [obj1, obj2]:
+            		if obj.existsAtInstant(t):
+            			if obj.hasFeatures():
+		                	n = len([f for f in obj.getFeatures() if f.existsAtInstant(t)])
+		            	else:
+		                	n = 1.
+		            	p += obj.getPositionAtInstant(t).__mul__(n)
+		            	nTotal += n
+            	assert nTotal>0, 'there should be at least one point for each instant'
+            	positions.addPosition(p.divide(nTotal))
+            # velocities: if any
+            if hasattr(obj1, 'velocities') and hasattr(obj2, 'velocities'):
+                velocities = Trajectory()
                 for t in newTimeInterval:
                     nTotal = 0.
                     p = Point(0.,0.)
@@ -1198,30 +1254,11 @@
                                 n = len([f for f in obj.getFeatures() if f.existsAtInstant(t)])
                             else:
                                 n = 1.
-                            p += obj.getPositionAtInstant(t).__mul__(n)
+                            p += obj.getVelocityAtInstant(t).__mul__(n)
                             nTotal += n
                     assert nTotal>0, 'there should be at least one point for each instant'
-                    positions.addPosition(p.divide(nTotal))
-                # velocities: if any
-                if obj1.getVelocities() is not None and obj2.getVelocities() is not None:
-                    velocities = Trajectory()
-                    for t in newTimeInterval:
-                        nTotal = 0.
-                        p = Point(0.,0.)
-                        for obj in [obj1, obj2]:
-                            if obj.existsAtInstant(t):
-                                if obj.hasFeatures():
-                                    n = len([f for f in obj.getFeatures() if f.existsAtInstant(t)])
-                                else:
-                                    n = 1.
-                                p += obj.getVelocityAtInstant(t).__mul__(n)
-                                nTotal += n
-                        assert n>0, 'there should be at least one point for each instant'
-                        velocities.addPosition(p.divide(nTotal))
-                else:
-                    velocities = None
+                    velocities.addPosition(p.divide(nTotal))
             else:
-                positions = None
                 velocities = None
             # TODO object envelop (polygon)
             # user type
@@ -1229,8 +1266,6 @@
                 print('The two moving objects have different user types: obj1 {} obj2 {}'.format(userTypeNames[obj1.getUserType()], userTypeNames[obj2.getUserType()]))
 
             newObject = MovingObject(newNum, newTimeInterval, positions, velocities, userType = obj1.getUserType())
-            if hasattr(obj1, 'featureNumbers') and hasattr(obj2, 'featureNumbers'):
-                newObject.featureNumbers = obj1.featureNumbers+obj2.featureNumbers
             if obj1.hasFeatures() and obj2.hasFeatures():
                 newObject.features = obj1.getFeatures()+obj2.getFeatures()
             return newObject
@@ -1943,7 +1978,7 @@
     axis('equal')
 
 
-if __name__ == "__main__":
+'''if __name__ == "__main__":
     import doctest
     import unittest
     suite = doctest.DocFileSuite('tests/moving.txt')
@@ -1954,3 +1989,4 @@
     if shapelyAvailable: 
         suite = doctest.DocFileSuite('tests/moving_shapely.txt')
         unittest.TextTestRunner().run(suite)
+            for t in newTimeInterval:'''