dan mackinlay avatar dan mackinlay committed 17ab71f

now we have MI estimates.

Comments (0)

Files changed (1)

 import numpy as np
 import scipy as sp
 from scipy.spatial import distance
-from copy import copy
+import pyentropy
+from random import sample
 
 # we expect params to contain the following keys:
 # num_agents, dt, noise, radius, steps
         self.params = params.copy()        
     def clone(self):
         """make a new copy of me with arrays pointing to the same memory locations
-        but with params dict cloned."""
+        but with params dict cloned so that branching thing works without trampling state."""
         clone = self.__class__(self.params.copy())
         for name in self.statevars:
             setattr(clone, name, getattr(self, name))
         return clone
 
-def basic_flock_factory(params):
-    """Initialise a sim with basic random state."""
-    sim = SimState(params=params)
-    sim.locs = np.random.uniform(size=(params['num_agents'], 2))
-    sim.vels = random_unit_vectors(params)
-    return sim    
-
 def do_flocking_sim(simstate, steps=100, seed=1):
     """
     Starting from the state in *simstate*, iterate the simulation *steps* times.
         for i in xrange(params['num_branches'])
     ]
 
+def analyse_branched_experiment(branches, n_pairs=10):
+    """Take a list of stacked trials and analyse them."""
+    #I'm punting I can ignore most of the state vector. Anyway, they'll be fiddly.
+    # stacked_vels = np.dstack([b.vels for b in branches])
+    #In fact, do we really need multiple axes of info? Nah.
+    stacked_vels = np.vstack([b.vels[:,0] for b in branches])
+    num_agents = stacked_vels.shape[0]
+    # Now, we estimate distributions of projected velocity in order to, in turn, estimate mean
+    # mutual information.
+    ests = np.zeros(n_pairs)
+    for i in xrange(n_pairs):
+        #this is not sliced efficiently. Should restack. Anyway...
+        pair = stacked_vels[sample(xrange(num_agents), 2), :]
+        ests[i] = ince_mi_dist_cont(pair[0], pair[1])
+    return ests
+
+def basic_flock_factory(params):
+    """Initialise a sim with basic random state."""
+    sim = SimState(params=params)
+    sim.locs = np.random.uniform(size=(params['num_agents'], 2))
+    sim.vels = random_unit_vectors(params)
+    return sim
+
 def normalise_lengths(row_vectors):
     """normalise the lengths of some stacked row vectors.
     """
 def fract(a):
     """shorthand for 'fractional part of'"""
     return np.remainder(a,1.0)
+
+def ince_mi_dist_cont(X, Y, n_bins=None, **kwargs):
+    """wraps pyentropy's grandiose object initialisation and quantisation
+    and just tells us what the bloody MI is, with no intermediate fiddling
+    about, from continuous input.
+    
+    If no n_bins is given, we choose the maximimum statistically valid number
+    by the cochrane criterion, and leave it be. This will give good results in
+    all except particularly pathological cases. (bins will be (marginwise)
+    approximately equal occupancy)
+    """
+    # kwargs.setdefault('method', 'qe')
+    # kwargs.setdefault('sampling', 'kt')
+    kwargs.setdefault('method', 'plugin')
+    kwargs.setdefault('sampling', 'naive')
+    X = np.asanyarray(X).ravel()
+    Y = np.asanyarray(Y).ravel()
+    if n_bins is None:
+        n_bins = choose_n_bins(X.size)
+    X, _ = pyentropy.quantise(X, n_bins, centers=False)
+    Y, _ = pyentropy.quantise(Y, n_bins, centers=False)
+    del(_) #Go away, extraneous data!
+    ds = pyentropy.DiscreteSystem(
+      X, (1, n_bins),
+      Y, (1, n_bins)
+    )
+    ds.calculate_entropies(**kwargs)
+    return ds.I()
+
+def choose_n_bins(n_data_points, test=True):
+    """according to Cellucci et al (2005), the maximal number of bins is given
+    thus"""
+    if n_data_points<20 and test:
+        raise ValueError("%d is too small a number to bin" % n_data_points)
+    return int(float(n_data_points/5.)**(0.5))
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.