123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207 |
- import random as rng
- import numpy as np
- from matplotlib.patches import Ellipse
- class Interneuron:
- def __init__(self, x, y, a, b, phi):
- self.x = x
- self.y = y
- self.a = a # Semimajor axis
- self.b = b # Semiminor axis
- self.phi = phi
- self.c = 2.0 * a # Sum distance
- self.foc_len = np.sqrt(a ** 2 - b ** 2)
- self.area = np.pi * a * b
- def get_p1(self):
- x1 = self.x - self.foc_len * np.cos(self.phi)
- y1 = self.y - self.foc_len * np.sin(self.phi)
- return x1, y1
- def get_p2(self):
- x2 = self.x + self.foc_len * np.cos(self.phi)
- y2 = self.y + self.foc_len * np.sin(self.phi)
- return x2, y2
- def get_ellipse(self):
- ell = Ellipse((self.x, self.y), 2.0 * self.a, 2.0 * self.b,
- angle=self.phi * 180. / np.pi, linewidth=1, fill=False, zorder=2)
- return ell
- def get_area(self):
- return np.pi * self.a * self.b
- # TODO: move to plotting utils
- def get_position_mesh(ex_positions):
- n_ex = int(np.sqrt(len(ex_positions)))
- x_wrong = np.array([x for x, _ in ex_positions]).reshape((n_ex, n_ex))
- y_wrong = np.array([y for _, y in ex_positions]).reshape((n_ex, n_ex))
- x = x_wrong[0, :]
- y = y_wrong[:, 0]
- xmin = np.min(x)
- xmax = np.max(x)
- dx = (xmax - xmin) / (x.shape[0] - 1)
- ymin = np.min(y)
- ymax = np.max(y)
- dy = (ymax - ymin) / (y.shape[0] - 1)
- X, Y = np.meshgrid(np.linspace(xmin, xmax + dx, n_ex + 1, endpoint=True) - dx / 2, np.linspace(ymin, ymax + dy, n_ex + 1, endpoint=True) - dy / 2)
- return X, Y
- def get_entropy_of_axonal_coverage(ex_neuron_tuning, ie_connections):
- ex_phase_vals_per_interneuron = []
- phase_entropy_per_interneuron = []
- n_hist_bins = 10
- for connected_idxs in ie_connections:
- ex_phase_vals_of_p = [ex_neuron_tuning[idx] for idx in connected_idxs]
- phase_entropy = 0.0
- ex_phase_vals_per_interneuron.append(ex_phase_vals_of_p)
- ex_phase_hist, _ = np.histogram(ex_phase_vals_of_p, n_hist_bins)
- n_ex_of_p = float(len(ex_phase_vals_of_p))
- ex_phase_hist = 1. / n_ex_of_p * ex_phase_hist.astype(float)
- for n_phi in ex_phase_hist:
- if n_phi != 0:
- phase_entropy -= n_phi * np.log(n_phi)
- phase_entropy_per_interneuron.append(phase_entropy)
- return phase_entropy_per_interneuron
- def get_excitatory_neurons_in_inhibitory_axonal_clouds(ex_neuron_positions, inhibitory_axonal_clouds):
- ie_connections = []
- for p in inhibitory_axonal_clouds:
- ie_conn_p = []
- x_p_1, y_p_1 = p.get_p1()
- x_p_2, y_p_2 = p.get_p2()
- for ex_idx, ex_pos in enumerate(ex_neuron_positions):
- x_ex, y_ex = ex_pos
- d_1 = np.sqrt((x_ex - x_p_1) ** 2 + (y_ex - y_p_1) ** 2)
- d_2 = np.sqrt((x_ex - x_p_2) ** 2 + (y_ex - y_p_2) ** 2)
- if d_1 + d_2 < p.c:
- ie_conn_p.append(ex_idx)
- ie_connections.append(ie_conn_p)
- return ie_connections
- def get_excitatory_phases_in_inhibitory_axon(ex_neuron_positions, ex_neuron_tunings, inhibitory_neuron):
- ex_phase_vals = []
- x_p_1, y_p_1 = inhibitory_neuron.get_p1()
- x_p_2, y_p_2 = inhibitory_neuron.get_p2()
- for ex_idx, ex_pos in enumerate(ex_neuron_positions):
- x_ex, y_ex = ex_pos
- d_1 = np.sqrt((x_ex - x_p_1) ** 2 + (y_ex - y_p_1) ** 2)
- d_2 = np.sqrt((x_ex - x_p_2) ** 2 + (y_ex - y_p_2) ** 2)
- if d_1 + d_2 < inhibitory_neuron.c:
- ex_phase_vals.append(ex_neuron_tunings[ex_idx])
- return ex_phase_vals
- # TODO: There must be some methods twice
- def get_interneuron_entropy(ex_phase_vals, entropy_bins):
- phase_entropy = 0.
- ex_phase_hist, _ = np.histogram(ex_phase_vals, entropy_bins)
- ex_phase_prob = ex_phase_hist / len(ex_phase_vals)
- for prob in ex_phase_prob:
- if prob > 0.0:
- phase_entropy += -prob * np.log(prob)
- return phase_entropy
- def create_interneuron_sheet_entropy_max_orientation(ex_positions, ex_tunings, number_of_interneurons, long_axis,
- short_axis, max_x, max_y, trial_orientations=30,
- number_of_bins=30, placement_jitter=0.0, seed=0):
- inhibitory_axonal_clouds = initialize_axonal_clouds_on_grid_with_random_orientation(number_of_interneurons,
- long_axis, short_axis, max_x,
- max_y, placement_jitter, seed)
- interneuron_positions = [(inh.x, inh.y) for inh in inhibitory_axonal_clouds]
- '''
- Determine initial tunings and entropy values
- '''
- entropies_of_interneurons = [0.]
- if trial_orientations > 0:
- entropies_of_interneurons, optimal_orientations = get_optimal_orientations_for_maximum_entropy(ex_positions,
- ex_tunings,
- interneuron_positions,
- long_axis,
- short_axis,
- number_of_bins,
- trial_orientations)
- for opt_orientation, axon in zip(optimal_orientations, inhibitory_axonal_clouds):
- axon.phi = opt_orientation
- return inhibitory_axonal_clouds, np.mean(entropies_of_interneurons)
- def get_optimal_orientations_for_maximum_entropy(ex_positions, ex_tunings, interneuron_positions, long_axis, short_axis,
- number_of_bins, trial_orientations):
- entropy_bins = np.linspace(-np.pi, np.pi, number_of_bins)
- entropies_of_interneurons = []
- optimal_orientations = []
- # for i_neuron in tqdm(inhibitory_axonal_clouds, desc="Maximizing entropy by rotation"):
- for x, y in interneuron_positions:
- orientations = np.linspace(-np.pi, np.pi, trial_orientations)
- entropy_for_orientation = []
- for orientation in orientations:
- perturbed_interneuron = Interneuron(x, y, long_axis, short_axis, orientation)
- ex_phase_vals = get_excitatory_phases_in_inhibitory_axon(ex_positions, ex_tunings, perturbed_interneuron)
- entropy_for_orientation.append(get_interneuron_entropy(ex_phase_vals, entropy_bins))
- optimal_orientations.append(orientations[np.argmax(entropy_for_orientation)])
- entropies_of_interneurons.append(np.max(entropy_for_orientation))
- return entropies_of_interneurons, optimal_orientations
- def initialize_axonal_clouds_on_grid_with_random_orientation(number_of_interneurons, long_axis, short_axis, sheet_x,
- sheet_y, placement_jitter=0.0, seed=0):
- '''
- Start with uniform distribution
- '''
- rng.seed(seed)
- interneuron_a = long_axis
- interneuron_b = short_axis
- interneuron_n = number_of_interneurons
- x_n = int(np.sqrt(interneuron_n))
- y_n = int(np.sqrt(interneuron_n))
- x_step = int(sheet_x / (x_n + 1))
- y_step = int(sheet_y / (y_n + 1))
- inhibitory_axonal_clouds = []
- pj = placement_jitter
- for n in range(interneuron_n):
- theta = rng.random() * 2. * np.pi
- pj_x = pj * np.cos(theta)
- pj_y = pj * np.sin(theta)
- inhibitory_axonal_clouds.append(Interneuron(((n % x_n) + 1) * x_step + pj_x,
- (int(n / y_n) + 1) * y_step + pj_y,
- interneuron_a, interneuron_b, theta))
- return inhibitory_axonal_clouds
- def create_grid_of_excitatory_neurons(sheet_size, n_e_per_row, tuning_map):
- x = np.linspace(0, sheet_size, n_e_per_row)
- y = np.linspace(0, sheet_size, n_e_per_row)
- ex_neuron_positions = []
- ex_neuron_tuning = []
- for y_idx, yy in enumerate(y):
- for x_idx, xx in enumerate(x):
- noise_value = tuning_map(xx, yy)
- ex_neuron_positions.append([xx, yy])
- ex_neuron_tuning.append(noise_value)
- return ex_neuron_positions, ex_neuron_tuning
|