Commit: 3d3bd6336937f3c21376f7376bf57e1364e1f943 - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 3d3bd6336937f3c21376f7376bf57e1364e1f943 (tree)
• Time: 2010-07-16 19:51:49
• Author: lorenzo
• Commiter: lorenzo

Log Message

I updated the script; now it is on par with the un-branched version.
Be careful: I replaced a "<" with a "<=" in the role in the definition of time_list
in the main code (loop executed at the end of the code, NOT any function definition).
This is wrong as some slices would be counted twice, but I did that just as a test.

Change Summary

• modified: Python-codes/iterate_tag_statistics_from_new_output_branch_1.py

Incremental Difference

diff -r 3b3765ab3ee7 -r 3d3bd6336937 Python-codes/iterate_tag_statistics_from_new_output_branch_1.py
--- a/Python-codes/iterate_tag_statistics_from_new_output_branch_1.py	Tue Jul 13 16:35:39 2010 +0000
+++ b/Python-codes/iterate_tag_statistics_from_new_output_branch_1.py	Fri Jul 16 10:51:49 2010 +0000
@@ -11,10 +11,18 @@
 #only with contact protocols
 
 def generate_iteration_grid(boot_time_non_unique,\
-                            number_intervals,interval_duration,ini_gap):
+                            number_intervals,interval_duration,t_ini,t_end):
     time_unique=s.unique1d(boot_time_non_unique)
 
-    time_ini=time_unique[0]+ini_gap
+    # time_ini=time_unique[0]+ini_gap
+
+    
+    if (t_ini>0):
+        time_ini=t_ini
+    else:
+        time_ini=time_unique[0]
+
+    
     print "time_ini is, ", time_ini
 
     if (number_intervals<=0 and  interval_duration<=0):
@@ -25,40 +33,22 @@
         print "Error in calling the function"
         #break
     if (number_intervals>0 and interval_duration<=0):
-        time_grid=s.linspace(time_ini,time_unique[-1],number_intervals+1)
+        if (t_end<=0):
+            time_grid=s.linspace(time_ini,time_unique[-1],number_intervals)
+        elif (t_end>0):
+            time_grid=s.linspace(time_ini,t_end,number_intervals)
+            
 
     if (number_intervals<=0 and interval_duration>0):
-        time_grid=s.arange(time_ini,time_unique[-1],interval_duration)
+        if (t_end<=0):
+            time_grid=s.arange(time_ini,time_unique[-1],interval_duration)
+        if (t_end>0):
+            time_grid=s.arange(time_ini,t_end,interval_duration)
+        
 
     return (time_grid.astype("int64"))
 
 
-def decompose_hash_64(my_hashed_number,my_len):
-    #this function takes a 64-bit number and decomposes it into 4 16-bit
-    #integers, namely tag_id A, boot A, tag_id B, boot B.
-
-    ini_arr=s.zeros(4).astype("int64")
-    arr_long=s.zeros(4*my_len).reshape(my_len,4).astype("int64")
-
-    first_part =  my_hashed_number >> 32 #to get back the 1st argument
-    second_part =  my_hashed_number &  0xFFFFFFFF #to get back the second argument 
-
-    tag_A = first_part >> 16
-    boot_A= first_part  &  0xFFFF
-
-    tag_B=second_part >> 16
-    boot_B= second_part &  0xFFFF
-
-    ini_arr[0]=tag_A
-    ini_arr[1]=boot_A
-    ini_arr[2]=tag_B
-    ini_arr[3]=boot_B
-
-    for i in xrange(my_len):
-        arr_long[i, :]=ini_arr
-
-    return (arr_long)
-    
 
 def contact_duration_and_interval_single_couple(time_list, delta_slice):
     #this function is modelled on contact_duration_and_interval_single_tag.
@@ -73,10 +63,6 @@
 
 
     time_list=(time_list-time_list[0])/delta_slice
-
-
-
-    
     gaps=s.diff(time_list) #a bit more efficient than the line above
 
     #print "gaps is, ", gaps
@@ -89,9 +75,7 @@
 
     #p.save("gaps.dat",gaps, fmt='%d')
 
-    # find_gap=s.where(gaps != 1)[0]
-
-    find_gap=s.where(gaps > 1)[0]
+    find_gap=s.where(gaps != 1)[0]
 
     gap_distr=(gaps[find_gap]-1)*delta_slice  #so, this is really the list of the
     #time interval between two contacts for my tag. After the discussion with Ciro,
@@ -139,6 +123,9 @@
 
     res=res*delta_slice
 
+        #Now I add some extra calculations. I also want to save the beginning and
+    #end time of each contact
+
     contact_begin=s.zeros(1)
     contact_end=s.zeros(1)
 
@@ -168,11 +155,23 @@
         contact_begin=s.hstack((tag_times[0],tag_times[find_gap+1]))
         contact_end=s.hstack((tag_times[find_gap],tag_times[-1]))
 
+    # if (0 in contact_begin or 0 in contact_end):
+    #     print "error with a zero!"
+
+    # print "len(find_gap) is, ", len(find_gap)
+    # print "len(res) is, ", len(res)
+
+    # print "res is, ", res
+    # print "contact_begin is, ", contact_begin
+    # print "contact_end is, ", contact_end
+    # print "contact_end-contact_begin is, ",contact_end-contact_begin
+    # print "tag_times is, ", tag_times
+
 
 
     #print "the sum of all the durations is, ", res.sum()
 
-    return [res,gap_distr,contact_begin,contact_end] 
+    return [res,gap_distr, contact_begin, contact_end]
 
 
 
@@ -239,12 +238,11 @@
         print "The chosen tag does not exist hence no analysis can be performed on it"
         return
 
+    tag_times=s.copy(single_tag_no_rep)
 
 
     # delta_slice=int(delta_slice) #I do not need floating point arithmetic
 
-    tag_times=s.copy(single_tag_no_rep)
-
     single_tag_no_rep=(single_tag_no_rep-single_tag_no_rep[0])/delta_slice
     gaps=s.diff(single_tag_no_rep) #a bit more efficient than the line above
 
@@ -266,6 +264,7 @@
     #from all the others when I see an increment larger than  one in the
     #rescaled time.
 
+
     gap_distr=(gaps[find_gap]-1)*delta_slice  #so, this is really the list of the
     #time interval between two contacts for my tag. After the discussion with Ciro,
     #I modified slightly the definition (now there is a -1) in the definition.
@@ -314,7 +313,9 @@
 
     res=res*delta_slice
 
-    #Now I add some extra calculations. I also want to save the beginning and
+
+
+        #Now I add some extra calculations. I also want to save the beginning and
     #end time of each contact
 
     contact_begin=s.zeros(1)
@@ -357,8 +358,12 @@
     # print "contact_end is, ", contact_end
     # print "contact_end-contact_begin is, ",contact_end-contact_begin
     # print "tag_times is, ", tag_times
-    
 
+
+
+    #print "the sum of all the durations is, ", res.sum()
+
+    # return [res,gap_distr]
     return [res,gap_distr, contact_begin, contact_end] 
 
 def contact_duration_and_interval_many_tags(sliced_interactions,\
@@ -375,6 +380,13 @@
     tag_ids= n.unique1d(s.ravel(sliced_interactions[:,1:3])) #to get a list of
     #all tag ID`s, which appear (repeated) on two rows of the matrix output by 
     # time_binned_interaction
+
+    n.savetxt("list_interacting_tag_ids.dat", tag_ids ,  fmt='%d')
+
+
+
+    #n.savetxt("tag_IDs.dat", tag_ids ,  fmt='%d')
+
     
     # tag_ids=tag_ids.astype('int')
 
@@ -395,7 +407,6 @@
     overall_end_contacts=s.zeros(0)
 
 
-
     for i in xrange(len(tag_ids)): 
         track_tag_id=tag_ids[i]  #i.e. iterate on all tags
         
@@ -404,25 +415,19 @@
 
         #print "contact_times is, ", contact_times
 
-        # if (i==0):
-        #     p.save("contact_times_single_tag.dat",contact_times, fmt="%d")
-        #     p.save("contact_times_single_tag_remapped.dat",\
-        #            (contact_times-contact_times[0])/delta_slice, fmt="%d")
-            
-
         results=contact_duration_and_interval_single_tag(contact_times, delta_slice)
 
         tag_duration=results[0]
-        # if (i==0):
-        #     p.save("duration_single_tag.dat",tag_duration, fmt="%d")
 
         
         tag_intervals=results[1] #get
         #an array with the time intervals between two contacts for a given tag
 
+
         tag_contact_begin=results[2]
         tag_contact_end=results[3]
-        
+ 
+
         
         #print "tag_intervals is, ", tag_intervals
         
@@ -434,42 +439,27 @@
         overall_end_contacts= \
                    s.hstack((overall_end_contacts,tag_contact_end))
         
+
+        
+        
         #print "overall_gaps is, ", overall_gaps
 
         overall_duration=s.hstack((overall_duration,tag_duration))
 
-
-    # print "overall_begin_contacts.shape is, ",overall_begin_contacts.shape
-    # print "overall_end_contacts.shape is, ",overall_end_contacts.shape
-    # print "overall_duration.shape is, ",overall_duration.shape
-
     #overall_gaps=overall_gaps[s.where(overall_gaps !=0)]
     #overall_duration=overall_duration[s.where(overall_duration !=0)]
     filename="many_tags_contact_interval_distr2_%01d"%(counter+1)
     filename=filename+"_.dat"
 
-    p.save(filename, overall_gaps ,  fmt='%d')
+    n.savetxt(filename, overall_gaps ,  fmt='%d')
 
     filename="many_tags_contact_duration_distr2_%01d"%(counter+1)
     filename=filename+"_.dat"
 
-    p.save(filename, overall_duration ,  fmt='%d')
-
-    # print "overall_duration,overall_begin_contacts,overall_end_contacts are, ",
-    # overall_duration,overall_begin_contacts, \
-    #          overall_end_contacts
     
-    overall_duration_and_start_end= \
-        s.vstack((overall_duration,overall_begin_contacts, \
-             overall_end_contacts))
+    n.savetxt(filename, overall_duration ,  fmt='%d')
 
-    filename="many_tags_contact_duration_distr_extended_%01d"%(counter+1)
-    filename=filename+"_.dat"
-
-    p.save(filename, overall_duration_and_start_end.transpose() ,  fmt='%d')
-    
-
-    return 0 #I have no need to return anythin else at this point
+    return overall_duration, overall_gaps
 
 # def my_hash(arr):
 #     my_hash=hash((arr[0], arr[1]))
@@ -477,17 +467,20 @@
 #     return my_hash
 
 
+
+
 def my_hash(arr):  #this will operate on a similar function which
     #has already combined tag_id and bootcount (both taking 16 bites each)
 
     # my_hash(arr) >> 32 to get back the 1st argument
     # my_hash(arr) &  0xFFFFFFFF to get back the second argument
 
+    arr=s.sort(arr) #maybe this is not needed at all; the rows of sliced data (which are arr)
+    #are already sorted out elsewhere.
+
     return (arr[0] << 32) | arr[1]
 
 
-
-
 def couple_hash_table(sliced_data,counter):
     hash_list=s.arange(len(sliced_data))
 
@@ -501,13 +494,13 @@
     filename=filename+"_.dat"
 
 
-    p.save(filename, hash_and_time,fmt='%d')
+    n.savetxt(filename, hash_and_time,fmt='%d')
 
     #the aim of this function is to identify each reported binary contact with a single number.
     #then, instead of looking for couples of interacting tags, I will simply look for repeated
     #occurrences of a number
 
-    return hash_and_time
+    return hash_and_time.astype("int64")
 
 def detect_binary_contacts_among_same_two_tags(hash_and_time,\
                                          delta_slice, counter):
@@ -521,28 +514,30 @@
     time_arr=hash_and_time[:,0] #important! Measure everything in units
     #of delta_slice
     
-    couple_interaction_duration=s.zeros(0) #this array will store the times at which the SAME
+    couple_interaction_duration=s.zeros(0).astype("int64") #this array will store the times at which the SAME
     #i-j tags interact.
 
-    couple_interaction_interval=s.zeros(0)
+    couple_interaction_interval=s.zeros(0).astype("int64")
+
+    record_hash=s.zeros(0).astype("int64")
+
+    hash_list_overall=s.zeros(0).astype("int64")
 
 
-    overall_begin_contacts=s.zeros(0)
-
-    overall_end_contacts=s.zeros(0)
+    overall_begin_contacts=s.zeros(0).astype("int64")
 
-    overall_extended_info=-s.ones(4).astype("int64")
-    
+    overall_end_contacts=s.zeros(0).astype("int64")
 
-    
+    # overall_extended_info=-s.ones(4).astype("int64")
+
+    overall_extended_info=-s.ones(2).astype("int64")
+
 
     for i in xrange(len(couple_list_unique)):
         chosen_couple=couple_list_unique[i]
         sel_chosen_couple_times=s.where(couple_list==chosen_couple)
         chosen_couple_times=time_arr[sel_chosen_couple_times]
-
-        
-        # print "chosen_couple is, ", chosen_couple
+        #print "i is, ", i
         #print "chosen_couple_times is, ", chosen_couple_times
         chosen_couple_contact_duration= \
              contact_duration_and_interval_single_couple(chosen_couple_times, delta_slice)[0]
@@ -550,44 +545,72 @@
         chosen_couple_contact_interval= \
              contact_duration_and_interval_single_couple(chosen_couple_times, delta_slice)[1]
 
-        
+        couple_interaction_duration=s.hstack((couple_interaction_duration,\
+                                              chosen_couple_contact_duration))
+
+        couple_interaction_interval=s.hstack((couple_interaction_interval,\
+                                              chosen_couple_contact_interval))
+
+        hash_in_iter=s.ones(len(chosen_couple_contact_duration)).astype("int64")\
+                      *chosen_couple
+
+        record_hash=s.hstack((record_hash,hash_in_iter))
+
+        interacting_couple_arr=s.ones(len(chosen_couple_contact_duration)).astype("int64")*chosen_couple
+
+        hash_list_overall=s.hstack((hash_list_overall,interacting_couple_arr))
+
+
         tag_contact_begin=\
          contact_duration_and_interval_single_couple(chosen_couple_times, delta_slice)[2]
         tag_contact_end=\
               contact_duration_and_interval_single_couple(chosen_couple_times, delta_slice)[3]
 
 
-        couple_interaction_duration=s.hstack((couple_interaction_duration,\
-                                              chosen_couple_contact_duration))
 
-        couple_interaction_interval=s.hstack((couple_interaction_interval,\
-                                              chosen_couple_contact_interval))
 
         overall_begin_contacts= \
                s.hstack((overall_begin_contacts,tag_contact_begin))
         overall_end_contacts= \
                    s.hstack((overall_end_contacts,tag_contact_end))
 
-        extended_info=decompose_hash_64(chosen_couple,len(tag_contact_begin))
+        # extended_info=decompose_hash_64(chosen_couple,len(tag_contact_begin))
+
+        extended_info=decompose_hash_32_vec(chosen_couple,len(tag_contact_begin))
+
+
+        # print "extended_info is, ", extended_info
+
+
+
+
+        
+
 
         overall_extended_info=s.vstack((overall_extended_info,extended_info))
 
 
 
+    filename="hashed_couples_corresponding_to_contact_duration_times_%01d"%(counter+1)
+    filename=filename+"_.dat"
+
+
     filename="couple_interaction_duration_times_%01d"%(counter+1)
     filename=filename+"_.dat"
 
 
-    p.save(filename, couple_interaction_duration,fmt='%d')
+    n.savetxt(filename, couple_interaction_duration,fmt='%d')
 
     filename="couple_interaction_interval_times_%01d"%(counter+1)
     filename=filename+"_.dat"
 
 
-    p.save(filename, couple_interaction_interval,fmt='%d')
+    n.savetxt(filename, couple_interaction_interval,fmt='%d')
 
-    # print "couple_interaction_duration.shape is, ",\
-    #       couple_interaction_duration.shape
+    filename="hashed_couples_in_contact_durations_%01d"%(counter+1)
+    filename=filename+"_.dat"
+
+    n.savetxt(filename, record_hash,fmt='%d')
 
 
 
@@ -614,11 +637,98 @@
         filename="couple_interaction_duration_times_extended_%01d"%(counter+1)
         filename=filename+"_.dat"
 
-        p.save(filename, overall_duration_and_start_end ,  fmt='%d')
+        n.savetxt(filename, overall_duration_and_start_end ,  fmt='%d')
 
 
     
-    return couple_interaction_duration
+    return [couple_interaction_duration,hash_list_overall]
+
+
+
+def decompose_hash_64(my_hashed_number,my_len):
+    #this function takes a 64-bit number and decomposes it into 4 16-bit
+    #integers, namely tag_id A, boot A, tag_id B, boot B.
+
+    ini_arr=s.zeros(4).astype("int64")
+    arr_long=s.zeros(4*my_len).reshape(my_len,4).astype("int64")
+
+    first_part =  my_hashed_number >> 32 #to get back the 1st argument
+    second_part =  my_hashed_number &  0xFFFFFFFF #to get back the second argument 
+
+    tag_A = first_part >> 16
+    boot_A= first_part  &  0xFFFF
+
+    tag_B=second_part >> 16
+    boot_B= second_part &  0xFFFF
+
+    ini_arr[0]=tag_A
+    ini_arr[1]=boot_A
+    ini_arr[2]=tag_B
+    ini_arr[3]=boot_B
+
+    for i in xrange(my_len):
+        arr_long[i, :]=ini_arr
+
+    return (arr_long)
+    
+
+
+
+def decompose_hash_32(my_hashed_number):
+    #this function takes a 64-bit number and decomposes it into 4 16-bit
+    #integers, namely tag_id A, boot A, tag_id B, boot B.
+
+    ini_arr=s.zeros(2).astype("int64")
+
+    tag_A = my_hashed_number >> 32
+    tag_B= my_hashed_number  &  0xFFFFFFFF
+
+    ini_arr[0]=tag_A
+    ini_arr[1]=tag_B
+    
+    return (ini_arr)
+
+
+
+
+def decompose_hash_32_vec(my_hashed_number,my_len):
+    #this function takes a 64-bit number and decomposes it into 4 16-bit
+    #integers, namely tag_id A, boot A, tag_id B, boot B.
+
+    ini_arr=s.zeros(2).astype("int64")
+    arr_long=s.zeros(2*my_len).reshape(my_len,2).astype("int64")
+
+    # first_part =  my_hashed_number >> 32 #to get back the 1st argument
+    # second_part =  my_hashed_number &  0xFFFFFFFF #to get back the second argument 
+
+    # tag_A = first_part >> 16
+    # boot_A= first_part  &  0xFFFF
+
+    # tag_B=second_part >> 16
+    # boot_B= second_part &  0xFFFF
+
+
+    tag_A = my_hashed_number >> 32
+    tag_B= my_hashed_number  &  0xFFFFFFFF
+
+    ini_arr[0]=tag_A
+    ini_arr[1]=tag_B
+
+
+    # ini_arr[0]=tag_A
+    # ini_arr[1]=boot_A
+    # ini_arr[2]=tag_B
+    # ini_arr[3]=boot_B
+
+    for i in xrange(my_len):
+        arr_long[i, :]=ini_arr
+
+    # print "arr_long is, ", arr_long
+
+    return (arr_long)
+    
+
+
 
 
 def remove_self_loops_and_save(slice_lump):
@@ -647,10 +757,13 @@
 #sliced_interactions=p.load("sliced_dynamics.dat")
 #sliced_interactions=p.load("edge_list.dat")
 
-ini_gap=1240876800-1240092599   # 86400 #i.e. one day in seconds
+t_ini=1254384000  # 86400 #i.e. one day in seconds
 
-number_intervals=-12
-interval_duration= 86400 #i.e. one day in seconds
+number_intervals=1
+interval_duration= -86400 #i.e. 1 day #604800 #i.e. one week in seconds
+t_end= 1254391200      #i.e. I choose when to end the time grid
+#whereas with a negative value of time_end everything is as it used to be.
+
 
 f = open(sys.argv[1])
 sliced_interactions = [map(int, string.split(line)) for line in f.readlines()]
@@ -665,9 +778,12 @@
 time_non_unique=sliced_interactions[:,0]
 
 time_grid= generate_iteration_grid(time_non_unique,\
-                            number_intervals,interval_duration,ini_gap)
+                            (number_intervals+1),interval_duration,t_ini,t_end)
 
-p.save("time_grid.dat", time_grid, fmt="%d")
+print "time_grid is, ", time_grid
+
+n.savetxt("time_grid.dat", time_grid, fmt="%d")
+
    
 file_list=-s.ones(len(time_grid))
 
@@ -676,7 +792,7 @@
     print "m+1 is, ", m+1
 
     time_sel=s.where((time_non_unique>= time_grid[m]) & \
-                     (time_non_unique< time_grid[m+1]))[0]
+                     (time_non_unique<= time_grid[m+1]))[0]
     
     sliced_interactions_sel=sliced_interactions[time_sel,:]
 
@@ -684,13 +800,31 @@
         file_list[m]=m+1
 
 
+    if (len(time_grid==2)): #then I just have a single interval I am going to save
+        n.savetxt("raw_interactions_before_processing.dat", sliced_interactions_sel,fmt="%d")
+
+
     contact_duration_and_interval_many_tags(sliced_interactions_sel,\
                                             delta_slice,m)
 
     hash_and_time=couple_hash_table(sliced_interactions_sel,m)
 
 
-    detect_binary_contacts_among_same_two_tags(hash_and_time, delta_slice,m)
+    binary_statistics=detect_binary_contacts_among_same_two_tags(hash_and_time, delta_slice,m)
+
+    tag_couples=s.arange(2*len(binary_statistics[1])).reshape((len(binary_statistics[1]),2)).astype("int64")
+
+    for q in xrange(len(binary_statistics[1])):
+        tag_couples[q,:]=decompose_hash_32(binary_statistics[1][q])
+    
+    filename="decomposed_hashed_couples_corresponding_to_contact_duration_times_%01d"%(m+1)
+    filename=filename+"_.dat"
+
+
+    n.savetxt(filename,tag_couples ,fmt='%d')
+                                                              
+
+                                                              
 
 
 sel=s.where(file_list>0)
@@ -700,7 +834,7 @@
 #print "file_list is, ", file_list
 
 
-p.save("file_list.dat", file_list, fmt='%d')
+n.savetxt("file_list.dat", file_list, fmt='%d')