Adjusted the share folder to run analyse script. Added some analysis of perlin map spatial scale

scripts/spatial_maps/simplex_algorithm_deconstruction.py

@@ -0,0 +1,284 @@
+import numpy as np
+import matplotlib.pyplot as pl
+
+
+# Based on: https://gist.github.com/KdotJPG/b1270127455a94ac5d19
+
+import sys
+from ctypes import c_int64
+from math import floor as _floor
+
+
+if sys.version_info[0] < 3:
+    def floor(num):
+        return int(_floor(num))
+else:
+    floor = _floor
+
+
+STRETCH_CONSTANT_2D = -0.211324865405187    # (1/Math.sqrt(2+1)-1)/2
+SQUISH_CONSTANT_2D = 0.366025403784439      # (Math.sqrt(2+1)-1)/2
+STRETCH_CONSTANT_3D = -1.0 / 6              # (1/Math.sqrt(3+1)-1)/3
+SQUISH_CONSTANT_3D = 1.0 / 3                # (Math.sqrt(3+1)-1)/3
+STRETCH_CONSTANT_4D = -0.138196601125011    # (1/Math.sqrt(4+1)-1)/4
+SQUISH_CONSTANT_4D = 0.309016994374947      # (Math.sqrt(4+1)-1)/4
+
+NORM_CONSTANT_2D = 47
+NORM_CONSTANT_3D = 103
+NORM_CONSTANT_4D = 30
+
+DEFAULT_SEED = 0
+
+
+# Gradients for 2D. They approximate the directions to the
+# vertices of an octagon from the center.
+GRADIENTS_2D = (
+     5,  2,    2,  5,
+    -5,  2,   -2,  5,
+     5, -2,    2, -5,
+    -5, -2,   -2, -5,
+)
+
+# Gradients for 3D. They approximate the directions to the
+# vertices of a rhombicuboctahedron from the center, skewed so
+# that the triangular and square facets can be inscribed inside
+# circles of the same radius.
+GRADIENTS_3D = (
+    -11,  4,  4,     -4,  11,  4,    -4,  4,  11,
+     11,  4,  4,      4,  11,  4,     4,  4,  11,
+    -11, -4,  4,     -4, -11,  4,    -4, -4,  11,
+     11, -4,  4,      4, -11,  4,     4, -4,  11,
+    -11,  4, -4,     -4,  11, -4,    -4,  4, -11,
+     11,  4, -4,      4,  11, -4,     4,  4, -11,
+    -11, -4, -4,     -4, -11, -4,    -4, -4, -11,
+     11, -4, -4,      4, -11, -4,     4, -4, -11,
+)
+
+# Gradients for 4D. They approximate the directions to the
+# vertices of a disprismatotesseractihexadecachoron from the center,
+# skewed so that the tetrahedral and cubic facets can be inscribed inside
+# spheres of the same radius.
+GRADIENTS_4D = (
+     3,  1,  1,  1,      1,  3,  1,  1,      1,  1,  3,  1,      1,  1,  1,  3,
+    -3,  1,  1,  1,     -1,  3,  1,  1,     -1,  1,  3,  1,     -1,  1,  1,  3,
+     3, -1,  1,  1,      1, -3,  1,  1,      1, -1,  3,  1,      1, -1,  1,  3,
+    -3, -1,  1,  1,     -1, -3,  1,  1,     -1, -1,  3,  1,     -1, -1,  1,  3,
+     3,  1, -1,  1,      1,  3, -1,  1,      1,  1, -3,  1,      1,  1, -1,  3,
+    -3,  1, -1,  1,     -1,  3, -1,  1,     -1,  1, -3,  1,     -1,  1, -1,  3,
+     3, -1, -1,  1,      1, -3, -1,  1,      1, -1, -3,  1,      1, -1, -1,  3,
+    -3, -1, -1,  1,     -1, -3, -1,  1,     -1, -1, -3,  1,     -1, -1, -1,  3,
+     3,  1,  1, -1,      1,  3,  1, -1,      1,  1,  3, -1,      1,  1,  1, -3,
+    -3,  1,  1, -1,     -1,  3,  1, -1,     -1,  1,  3, -1,     -1,  1,  1, -3,
+     3, -1,  1, -1,      1, -3,  1, -1,      1, -1,  3, -1,      1, -1,  1, -3,
+    -3, -1,  1, -1,     -1, -3,  1, -1,     -1, -1,  3, -1,     -1, -1,  1, -3,
+     3,  1, -1, -1,      1,  3, -1, -1,      1,  1, -3, -1,      1,  1, -1, -3,
+    -3,  1, -1, -1,     -1,  3, -1, -1,     -1,  1, -3, -1,     -1,  1, -1, -3,
+     3, -1, -1, -1,      1, -3, -1, -1,      1, -1, -3, -1,      1, -1, -1, -3,
+    -3, -1, -1, -1,     -1, -3, -1, -1,     -1, -1, -3, -1,     -1, -1, -1, -3,
+)
+
+
+def overflow(x):
+    # Since normal python ints and longs can be quite humongous we have to use
+    # this hack to make them be able to overflow
+    return c_int64(x).value
+
+
+class OpenSimplex(object):
+    """
+    OpenSimplex n-dimensional gradient noise functions.
+    """
+
+    def __init__(self, seed=DEFAULT_SEED):
+        """
+        Initiate the class using a permutation array generated from a 64-bit seed number.
+        """
+        # Generates a proper permutation (i.e. doesn't merely perform N
+        # successive pair swaps on a base array)
+        perm = self._perm = [0] * 256 # Have to zero fill so we can properly loop over it later
+        perm_grad_index_3D = self._perm_grad_index_3D = [0] * 256
+        source = [i for i in range(0, 256)]
+        seed = overflow(seed * 6364136223846793005 + 1442695040888963407)
+        seed = overflow(seed * 6364136223846793005 + 1442695040888963407)
+        seed = overflow(seed * 6364136223846793005 + 1442695040888963407)
+        for i in range(255, -1, -1):
+            seed = overflow(seed * 6364136223846793005 + 1442695040888963407)
+            r = int((seed + 31) % (i + 1))
+            if r < 0:
+                r += i + 1
+            perm[i] = source[r]
+            perm_grad_index_3D[i] = int((perm[i] % (len(GRADIENTS_3D) / 3)) * 3)
+            source[r] = source[i]
+
+    def _extrapolate2d(self, xsb, ysb, dx, dy):
+        perm = self._perm
+        index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E
+
+        g1, g2 = GRADIENTS_2D[index:index + 2]
+        return g1 * dx + g2 * dy
+
+    def _extrapolate3d(self, xsb, ysb, zsb, dx, dy, dz):
+        perm = self._perm
+        index = self._perm_grad_index_3D[
+            (perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF
+            ]
+
+        g1, g2, g3 = GRADIENTS_3D[index:index + 3]
+        return g1 * dx + g2 * dy + g3 * dz
+
+    def _extrapolate4d(self, xsb, ysb, zsb, wsb, dx, dy, dz, dw):
+        perm = self._perm
+        index = perm[(
+            perm[(
+                perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb
+            ) & 0xFF] + wsb
+        ) & 0xFF] & 0xFC
+
+        g1, g2, g3, g4 = GRADIENTS_4D[index:index + 4]
+        return g1 * dx + g2 * dy + g3 * dz + g4 * dw
+
+    def noise2d(self, x, y):
+        """
+        Generate 2D OpenSimplex noise from X,Y coordinates.
+        """
+        # Place input coordinates onto grid.
+        stretch_offset = (x + y) * STRETCH_CONSTANT_2D
+        xs = x + stretch_offset
+        ys = y + stretch_offset
+
+        # Floor to get grid coordinates of rhombus (stretched square) super-cell origin.
+        xsb = floor(xs)
+        ysb = floor(ys)
+
+        # Skew out to get actual coordinates of rhombus origin. We'll need these later.
+        squish_offset = (xsb + ysb) * SQUISH_CONSTANT_2D
+        xb = xsb + squish_offset
+        yb = ysb + squish_offset
+
+        # Compute grid coordinates relative to rhombus origin.
+        xins = xs - xsb
+        yins = ys - ysb
+
+        # Sum those together to get a value that determines which region we're in.
+        in_sum = xins + yins
+
+        # Positions relative to origin point.
+        dx0 = x - xb
+        dy0 = y - yb
+
+        value = 0
+
+        # Contribution (1,0)
+        dx1 = dx0 - 1 - SQUISH_CONSTANT_2D
+        dy1 = dy0 - 0 - SQUISH_CONSTANT_2D
+        attn1 = 2 - dx1 * dx1 - dy1 * dy1
+        extrapolate = self._extrapolate2d
+        # if attn1 > 0:
+        #     attn1 *= attn1
+        #     value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, dx1, dy1)
+
+        # Contribution (0,1)
+        dx2 = dx0 - 0 - SQUISH_CONSTANT_2D
+        dy2 = dy0 - 1 - SQUISH_CONSTANT_2D
+        attn2 = 2 - dx2 * dx2 - dy2 * dy2
+        # if attn2 > 0:
+        #     attn2 *= attn2
+        #     value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, dx2, dy2)
+
+        if in_sum <= 1: # We're inside the triangle (2-Simplex) at (0,0)
+            zins = 1 - in_sum
+            if zins > xins or zins > yins: # (0,0) is one of the closest two triangular vertices
+                if xins > yins:
+                    xsv_ext = xsb + 1
+                    ysv_ext = ysb - 1
+                    dx_ext = dx0 - 1
+                    dy_ext = dy0 + 1
+                else:
+                    xsv_ext = xsb - 1
+                    ysv_ext = ysb + 1
+                    dx_ext = dx0 + 1
+                    dy_ext = dy0 - 1
+            else: # (1,0) and (0,1) are the closest two vertices.
+                xsv_ext = xsb + 1
+                ysv_ext = ysb + 1
+                dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D
+                dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D
+        else: # We're inside the triangle (2-Simplex) at (1,1)
+            zins = 2 - in_sum
+            if zins < xins or zins < yins: # (0,0) is one of the closest two triangular vertices
+                if xins > yins:
+                    xsv_ext = xsb + 2
+                    ysv_ext = ysb + 0
+                    dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D
+                    dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D
+                else:
+                    xsv_ext = xsb + 0
+                    ysv_ext = ysb + 2
+                    dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D
+                    dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D
+            else: # (1,0) and (0,1) are the closest two vertices.
+                dx_ext = dx0
+                dy_ext = dy0
+                xsv_ext = xsb
+                ysv_ext = ysb
+            xsb += 1
+            ysb += 1
+            dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D
+            dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D
+
+        # Contribution (0,0) or (1,1)
+        attn0 = 2 - dx0 * dx0 - dy0 * dy0
+        if attn0 > 0:
+            attn0 *= attn0
+            value += attn0 * attn0 * extrapolate(xsb, ysb, dx0, dy0)
+
+        # Extra Vertex
+        attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext
+        if attn_ext > 0:
+            attn_ext *= attn_ext
+            value += attn_ext * attn_ext * extrapolate(xsv_ext, ysv_ext, dx_ext, dy_ext)
+
+        return value / NORM_CONSTANT_2D
+
+
+noise = OpenSimplex()
+
+nrow = 60
+size = 900
+scale = 200
+
+x = y = np.linspace(0, size, nrow)
+n = [[noise.noise2d(i/scale, j/scale) for j in y] for i in x]
+m = np.concatenate(n)
+sorted_idx = np.argsort(m)
+max_val = nrow * 2
+idx = len(m) // max_val
+for ii, val in enumerate(range(max_val)):
+    m[sorted_idx[ii * idx:(ii + 1) * idx]] = val
+landscape = (m - nrow) / nrow
+
+fig, ax = pl.subplots(1, 1)
+
+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.arange(xmin, xmax + 2 * dx, dx) - dx / 2., np.arange(ymin, ymax + 2 * dy, dy) - dy / 2.)
+
+c = ax.pcolor(X, Y, landscape.reshape(nrow, -1), vmin=-1, vmax=1, cmap='hsv')
+fig.colorbar(c, ax=ax, label="in/out-degree")
+
+# pl.matshow(landscape.reshape(nrow, -1), vmin=-1, vmax=1)
+# print(x)
+# pl.gca().set_xticklabels(x)
+# pl.gca().set_yticklabels(y)
+# pl.colorbar()
+
+print(1/(1+2*STRETCH_CONSTANT_2D))
+print((1/(1+2*STRETCH_CONSTANT_2D))/np.sqrt(3))
+
+pl.show()

scripts/spatial_maps/simplex_input_tuning_correlation/tuning_difference_per_distance_perlin.py

@@ -0,0 +1,98 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import correlate2d
+
+from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
+    create_interneuron_sheet_entropy_max_orientation, get_correct_position_mesh
+from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap
+
+
+def phi_diff(d_phi):
+    return (d_phi + np.pi/2.) % np.pi - (np.pi/2.)
+
+
+dim = 60
+size = 900
+n_exc = dim * dim
+
+corr_len_list = [100, 200, 500]
+
+fig, axes = plt.subplots(3, 4)
+
+for corr_len, axes_line in zip(corr_len_list, axes):
+
+    perlin_map = UniformPerlinMap(dim, dim, corr_len, size, size, 3)
+
+    ex_positions, ex_tunings = create_grid_of_excitatory_neurons(perlin_map.sheet_x, perlin_map.sheet_y, perlin_map.x_dim,
+                                                                 perlin_map.get_tuning)
+
+    X, Y = get_correct_position_mesh(ex_positions)
+
+    # fig, ax = plt.subplots(1, 1)
+
+    head_dir_preference = np.array(ex_tunings).reshape((dim, dim))
+    # TODO: Why was this transposed for plotting? (now changed)
+    c = axes_line[0].pcolor(X, Y, head_dir_preference, vmin=-np.pi, vmax=np.pi, cmap='hsv')
+
+    axes_line[0].set_title('Corr. len.: {}'.format(corr_len))
+    axes_line[0].set_aspect('equal')
+
+
+    dist_list = []
+    phi_list = []
+    d_phi_list = []
+    for i in range(n_exc):
+        phi_list.append(ex_tunings[i])
+        for j in range(n_exc):
+            xi, yi = ex_positions[i]
+            xj, yj = ex_positions[j]
+            dist = np.sqrt((xi - xj)**2 + (yi - yj)**2)
+            dist_list.append(dist)
+            d_phi = np.abs(phi_diff(ex_tunings[i] - ex_tunings[j]))
+            d_phi_list.append(d_phi)
+
+
+    # bins = np.linspace(0.0, np.sqrt(2*size**2), 61)
+    bins = np.linspace(0.0, 500, 61)
+
+    print(bins)
+    binned_ids = np.digitize(dist_list, bins)
+    binned_d_phi = [[] for _ in range(len(bins))]
+    for id, bin_id in enumerate(binned_ids):
+        current_d_phi = d_phi_list[id]
+        # print(bin_id-1)
+        binned_d_phi[bin_id-1].append(current_d_phi)
+
+    binned_mean = [np.mean(bin) for bin in binned_d_phi]
+
+    axes_line[1].plot(bins, binned_mean)
+    axes_line[1].set_title('Mean tuning difference over distance')
+    axes_line[1].set_ylim(0.0, np.pi / 4. + 0.1)
+
+    phi_bins = np.linspace(-np.pi, np.pi, 31)
+    axes_line[2].hist(phi_list, bins=phi_bins)
+    axes_line[2].set_title('# of pyramidal cells per orientation')
+
+    # print(dist_list)
+    # print(d_phi_list)
+    # print(binned_d_phi)
+    # print(phi_diff(np.linspace(0.0, 3*np.pi, 11)))
+
+    near_far_bins = [0.0, 150.0, 10000.0]
+    near_far_ids = np.digitize(dist_list, near_far_bins)
+    near_far_d_phi = [[] for _ in range(len(near_far_bins)-1)]
+
+    for id, bin_id in enumerate(near_far_ids):
+        near_far_d_phi[bin_id-1].append(d_phi_list[id])
+
+    near_far_mean = [np.mean(bin) for bin in near_far_d_phi]
+
+    # fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5))
+
+    axes_line[3].bar([0., 1.], near_far_mean)
+    axes_line[3].set_xticks([0., 1.])
+    axes_line[3].set_xticklabels(['< 150', '> 150'])
+    axes_line[3].set_title('Mean orientation difference near vs far')
+    axes_line[3].set_ylim(0.0, np.pi/4.+0.1)
+
+plt.show()

scripts/spatial_maps/spreizer_test.py

@@ -1,11 +1,14 @@
 import numpy as np
 import pylab as pl
-import noise
+from opensimplex import OpenSimplex
+
+noise = OpenSimplex()
 
 nrow = 100
-size = 4
+size = 10
+scale = 1
 x = y = np.linspace(0, size, nrow)
-n = [[noise.pnoise2(i, j, repeatx=size, repeaty=size) for j in y] for i in x]
+n = [[noise.noise2d(i/scale, j/scale) for j in y] for i in x]
 m = np.concatenate(n)
 sorted_idx = np.argsort(m)
 max_val = nrow * 2

scripts/spatial_network/perlin/run_entropy_maximisation_perlin_map.py

@@ -31,7 +31,7 @@ def get_local_data_folder():
     return data_folder
 
 
-DATA_FOLDER = get_local_data_folder()
+DATA_FOLDER = "/groups/susanne/moritz/share/"
 LOG_FOLDER = "../../../logs/"
 
 TRAJ_NAME = "full_figure_perlin_map"