first

Author: aleksejs-fomins

.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,19 @@
+Copyright (c) 2018 The Python Packaging Authority
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

.idea/modules.xml

@@ -0,0 +1,5 @@
+# Example Package
+
+This is a simple example package. You can use
+[Github-flavored Markdown](https://guides.github.com/features/mastering-markdown/)
+to write your content.

.idea/workspace.xml

@@ -0,0 +1,20 @@
+# import standard libraries
+import os, sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+# import mesostat
+from mesostat.utils.signals import approx_decay_conv
+
+# Create signal
+DT =  0.001   # s
+TAU = 0.1     # s
+t = np.arange(0, 10, DT)
+y = (np.sin(1 * t)**2 > 0.7).astype(float)
+
+yc = approx_decay_conv(y, TAU, DT)
+
+plt.figure()
+plt.plot(t, y)
+plt.plot(t, yc)
+plt.show()

LICENSE

@@ -0,0 +1,7 @@
+from mesostat.utils.decorators import time_mem_1starg
+
+@time_mem_1starg
+def myfunc(x):
+    return x**2
+
+print(myfunc(10))

README.md

@@ -0,0 +1,75 @@
+'''
+VERY IMPORTANT:
+DO NOT USE SAME PROCEDURE FOR UPSAMPLING AND DOWNSAMPLING
+
+* When we upsample, we want to interpolate data
+* When we downsample, we want to sample the trendline, not individual fluctuations
+
+[] Maybe impl FFT-based resampling, that joins both
+'''
+
+# import standard libraries
+import os, sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+# import special libraries
+from mesostat.utils.signals import resample
+
+##########################
+# Downsampling
+##########################
+
+# Create data
+T  = 10     # s
+DT = 0.001  # s
+t1 = np.arange(0, T, DT)
+y1 = np.random.normal(0, 1, t1.size)
+for i in range(1, t1.size):
+    y1[i] += y1[i-1]
+
+# Resample
+DT2 = 0.1   # s
+t2 = np.arange(0, t1[-1], DT2)
+y2 = resample(t1, y1, t2, {"method" : "smooth", "kind" : "window"})
+y3 = resample(t1, y1, t2, {"method" : "smooth", "kind" : "kernel", "ker_sig2" : DT2**2})
+y4 = resample(t1, y1, t2, {"method" : "smooth", "kind" : "kernel", "ker_sig2" : (DT2/2)**2})
+y5 = resample(t1, y1, t2, {"method" : "smooth", "kind" : "kernel", "ker_sig2" : (DT2/4)**2})
+
+# Plot
+plt.figure()
+plt.title('Downsampling using window and kernel estimators')
+plt.plot(t1, y1, '.-', label='orig')
+plt.plot(t2, y2, '.-', label='window')
+plt.plot(t2, y3, '.-', label="ker, s2=d2^2")
+plt.plot(t2, y4, '.-', label="ker, s2=(d2/2)^2")
+plt.plot(t2, y5, '.-', label="ker, s2=(d2/4)^2")
+plt.legend()
+
+##########################
+# Upsampling
+##########################
+
+# Create Data
+T  = 10     # s
+DT = 0.1    # s
+t1 = np.arange(0, T, DT)
+y1 = np.sin(10*t1) * np.sin(2*t1)
+
+# Resample
+DT2 = 0.01  # s
+t2 = np.arange(0, t1[-1], DT2)
+y2 = resample(t1, y1, t2, {"method" : "interpolative", "kind" : "linear"})
+y3 = resample(t1, y1, t2, {"method" : "interpolative", "kind" : "quadratic"})
+y4 = resample(t1, y1, t2, {"method" : "interpolative", "kind" : "cubic"})
+
+# Plot
+plt.figure()
+plt.title('Upsampling using interpolation')
+plt.plot(t1, y1, 'o', label='orig')
+plt.plot(t2, y2, '.-', label='linear')
+plt.plot(t2, y3, '.-', label="quadratic")
+plt.plot(t2, y4, '.-', label="cubic")
+plt.legend()
+
+plt.show()

examples/example_decay_conv.py

@@ -0,0 +1,12 @@
+import sys
+import numpy as np
+
+from mesostat.utils.system import mem_now_as_str, bytes2str
+
+def heavyFunc(i):
+    x = np.random.normal(0,1, 10**8)
+    print(i, mem_now_as_str(), bytes2str(sys.getsizeof(x)))
+    return i
+
+
+print(list(map(heavyFunc, np.arange(10))))

examples/example_decorator_lib.py

@@ -0,0 +1,67 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mesostat.stat.stat as stat
+
+plt.ion()
+
+############################
+# Random bool permutation
+############################
+N_TEST1 = 1000
+nTot = 100
+freqArr = np.zeros(nTot)
+nTrueArr = np.random.randint(0, 100, N_TEST1)
+
+for nTrue in zip(nTrueArr):
+    perm = stat.rand_bool_perm(nTrue, nTot).astype(int)
+    assert np.sum(perm) == nTrue
+    freqArr += perm
+
+plt.figure()
+plt.title("Testing array permutation for bias")
+plt.xlabel("array element")
+plt.ylabel("number of occurences")
+plt.plot(freqArr, '.')
+plt.show()
+
+############################
+# Discrete PDF and CDF
+############################
+N_VAL_TEST2 = 20
+keys = np.random.randint(0, 1000, N_VAL_TEST2)
+vals = np.random.uniform(0, 1, N_VAL_TEST2)
+vals /= np.sum(vals)
+
+pdf = dict(zip(keys, vals))
+cdf = stat.discrete_distr_to_cdf(pdf)
+
+plt.figure()
+plt.title("Discrete PDF and CDF")
+plt.semilogy(pdf.keys(), pdf.values(), '.')
+plt.semilogy(cdf.keys(), cdf.values())
+plt.show()
+
+############################
+# Sampling distributions
+############################
+N_TEST3 = 1000
+N_VAL_TEST3 = 20
+keys = np.arange(N_VAL_TEST3)
+vals = np.random.uniform(0,1,N_VAL_TEST3)
+vals /= np.sum(vals)
+pdf = dict(zip(keys,vals))
+cdf = stat.discrete_distr_to_cdf(pdf)
+
+resample = stat.discrete_cdf_sample(cdf, N_TEST3)
+emp_pdf = stat.discrete_empirical_pdf_from_sample(resample)
+
+assert emp_pdf.keys() == pdf.keys()
+
+plt.figure()
+plt.title("Resamling random distribution using " + str(N_TEST3) + " trials")
+plt.bar(pdf.keys(), pdf.values(), alpha=0.5, label="true")
+plt.bar(emp_pdf.keys(), emp_pdf.values(), alpha=0.5, label="resampled")
+plt.legend()
+plt.ioff()
+plt.show()

examples/example_downsample.py

@@ -0,0 +1,38 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from mesostat.stat.htests import difference_test
+from mesostat.stat.hteststatistics import cranksum
+
+A = np.random.uniform(0.5,2,100)#np.abs(np.random.normal(0,1,100))
+B = np.random.uniform(0,1,100)#np.abs(np.random.normal(0,1,100))
+
+def CDF(x):
+    return np.sort(x), np.linspace(0,1,len(x))
+
+def permutation_test_distinct(f, x, y, nResample=2000):
+    M, N = len(x), len(y)
+    fTrue = f(x, y)
+    xy = np.hstack([x, y])
+    fTest = []
+    for iTest in range(nResample):
+        xyShuffle = xy[np.random.permutation(M + N)]
+        xTest, yTest = xyShuffle[:M], xyShuffle[M:]
+        fTest += [f(xTest, yTest)]
+
+    plt.figure()
+    plt.hist(fTest, bins='auto')
+    plt.axvline(x=fTrue, color='r')
+
+
+xa, ca = CDF(A)
+xb, cb = CDF(B)
+
+fig, ax = plt.subplots()
+ax.plot(xa, ca)
+ax.plot(xb, cb)
+
+permutation_test_distinct(cranksum, A, B)
+plt.show()
+

examples/example_map_memory_use.py

@@ -0,0 +1,115 @@
+# import standard libraries
+import os, sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+# import mesostat
+from mesostat.metric.corr import cross_corr_3D
+from mesostat.utils.plotting import plot_matrix
+
+'''
+   Test 1:
+     Generate random data, and shift it by fixed steps for each channel
+     Expected outcomes:
+     * If shift <= max_delay, corr ~ 1, delay = shift
+     * If shift > max_delay, corr ~ 0, delay = rand
+     * Delay is the same for all diagonals, because we compare essentially the same data, both cycled by the same amount
+'''
+    
+nNode = 5
+nData = 1000
+settings = {
+    'min_lag_sources' : 1,
+    'max_lag_sources' : 2,
+    'dim_order'       : 'ps'
+}
+
+#data = np.random.uniform(0, 1, nNode*nData).reshape((nNode, nData))
+# Generate progressively more random data
+data = np.zeros((nNode, nData))
+data[0] = np.random.normal(0, 1, nData)
+for i in range(1, nNode):
+    data[i] = np.hstack((data[i-1][1:], data[i-1][0]))
+
+rezCorr = cross_corr_3D(data, settings, est='corr')
+rezSpr = cross_corr_3D(data, settings, est='spr')
+
+compose = lambda lst1, lst2: [a + "_" + b for a in lst1 for b in lst2]
+
+plot_matrix(
+    "Test 1: Channels are shifts of the same data",
+    (2, 3),
+    [*rezCorr, *rezSpr],
+    np.array(compose(["corr", "spr"], ["val", "lag", "p"])),
+    lims = [[-1, 1], [0, settings['max_lag_sources']], [0, 1]]*2,
+    draw = True
+)
+
+
+'''
+   Test 2:
+     Generate random data, all copies of each other, each following one a bit more noisy than prev
+     Expected outcomes:
+     * Correlation decreases with distance between nodes, as they are separated by more noise
+     * Correlation should be approx the same for any two nodes given fixed distance between them
+'''
+
+nNode = 5
+nData = 1000
+settings = {
+    'min_lag_sources' : 0,
+    'max_lag_sources' : 0,
+    'dim_order'       : 'ps'
+}
+alpha = 0.5
+
+data = np.random.normal(0, 1, (nNode, nData))
+for i in range(1, nNode):
+    data[i] = data[i-1] * np.sqrt(1 - alpha) + np.random.normal(0, 1, nData) * np.sqrt(alpha)
+
+rezCorr = cross_corr_3D(data, settings, est='corr')
+rezSpr  = cross_corr_3D(data, settings, est='spr')
+
+plot_matrix(
+    "Test 2: Channels are same, but progressively more noisy",
+    (2, 3),
+    [*rezCorr, *rezSpr],
+    np.array(compose(["corr", "spr"], ["val", "lag", "p"])),
+    lims = [[-1, 1], [0, settings['max_lag_sources']], [0, 1]]*2,
+    draw = True
+)
+
+
+'''
+   Test 3:
+     Random data structured by trials. Two channels (0 -> 3) connected with lag 6, others unrelated
+     Expected outcomes:
+     * No structure, except for (0 -> 3) connection
+'''
+
+nNode = 5
+lagTrue = 6
+settings = {
+    'min_lag_sources' : 1,
+    'max_lag_sources' : 6,
+    'dim_order'       : 'rsp'
+}
+nData = settings['max_lag_sources']+1
+nTrial = 200
+
+data = np.random.normal(0, 1, (nTrial,nData,nNode))
+data[0, lagTrue:, :] = data[3, :-lagTrue, :]
+
+rezCorr = cross_corr_3D(data, settings, est='corr')
+rezSpr  = cross_corr_3D(data, settings, est='spr')
+
+plot_matrix(
+    "Test 3: Random trial-based cross-correlation",
+    (2, 3),
+    [*rezCorr, *rezSpr],
+    np.array(compose(["corr", "spr"], ["val", "lag", "p"])),
+    lims = [[-1, 1], [0, settings['max_lag_sources']], [0, 1]]*2,
+    draw = True
+)
+
+plt.show()