MEUL with a non-parametric homeostasis

In this notebook, we will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ). Compared to other posts, such as this previous post, we improve the code to not depend on any parameter (namely the Cparameter of the rescaling function). For that, we will use a non-parametric approach based on the use of cumulative histograms.

This is joint work with Victor Boutin and Angelo Francisioni. See also the other posts on unsupervised learning.

In [1]:
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
import numpy as np
np.set_printoptions(formatter = dict( float = lambda x: "%.3g" % x ), precision=3, suppress=True, threshold=np.inf)
from shl_scripts.shl_experiments import SHL
import time
%load_ext autoreload
%autoreload 2
In [2]:
DEBUG = True
DEBUG = False
if not DEBUG:
    tag = '2017-12-01_Testing_COMPs'
    DEBUG_DOWNSCALE = 1
else:
    tag = '2017-12-01_Testing_COMPs-DEBUG'
    DEBUG_DOWNSCALE = 10

seed = 42
nb_quant = 512
C = 0.
do_sym = False
verbose = False
i_sample = 42

from shl_scripts.shl_experiments import SHL
matname = tag + '_nohomeo'
shl = SHL(DEBUG_DOWNSCALE=DEBUG_DOWNSCALE,
          datapath='/tmp/database', seed=seed, 
          eta=0.05, verbose=2, record_each=50, n_iter=1000, 
          eta_homeo=0., alpha_homeo=1., 
          do_sym=do_sym, nb_quant=nb_quant)
data = shl.get_data(matname=matname)
loading the data called : data_cache/2017-12-01_Testing_COMPs_nohomeo_data
In [3]:
test_size = shl.batch_size # data.shape[0]//2
data_training = data[:test_size, :]
data_test = data[test_size:,:]   
if DEBUG:
    test_size = data.shape[0]//20
    data_training = data[:(data.shape[0]-test_size),:].copy()
    data_test = data[:test_size, :].copy()
In [4]:
dico_partial_learning = shl.learn_dico(data=data_training, matname=matname)
loading the dico called : data_cache/2017-12-01_Testing_COMPs_nohomeo_dico.pkl

We start off by using a short learning with no homeostasis such that we end up with a unbalanced dictionary:

In [5]:
fig, ax = shl.show_dico(dico_partial_learning, data=data, title=matname)
fig.show()
fig, ax = shl.time_plot(dico_partial_learning, variable='prob_active');
fig.show()
No description has been provided for this image
No description has been provided for this image

MP classique

In [6]:
n_samples, n_pixels = data_test.shape
n_dictionary, n_pixels = dico_partial_learning.dictionary.shape
norm_each_filter = np.sqrt(np.sum(dico_partial_learning.dictionary**2, axis=1))
dico_partial_learning.dictionary /= norm_each_filter[:,np.newaxis]

sparse_code_mp = shl.code(data_test, dico_partial_learning, matname=matname)

from shl_scripts.shl_tools import plot_proba_histogram
fig, ax = plot_proba_histogram(sparse_code_mp)
loading the code called : data_cache/2017-12-01_Testing_COMPs_nohomeo_coding.npy
No description has been provided for this image

COMP : learning rescaling on sparse coefficients - basic principle

In [7]:
n_samples, n_dictionary = sparse_code_mp.shape
N = n_samples * n_dictionary
q = np.zeros_like(sparse_code_mp)
for i in range(i_sample, i_sample+1): #range(n_samples):
    for k in range(n_dictionary):
        if sparse_code_mp[i, k] > 0:
            q[i, k] = np.sum(sparse_code_mp[i, k]>sparse_code_mp.ravel())/N
            if q[i, k]==0:
                print(k, i, sparse_code_mp[i, k])
            if i < 22 and k < 20:
                print(i, k, 'Raw value=', sparse_code_mp[i, k], 
                  ' is transformed into=', q[i, k]) 
In [8]:
def rescaling(code, do_sym=do_sym):
    if do_sym:
        code = np.abs(code)
    else:
        code *= code>0
    n_samples, n_dictionary = code.shape
    N = n_samples * n_dictionary
    q = np.zeros_like(code)
    for i in range(n_samples):
        for k in range(n_dictionary):
            if code[i, k] > 0:
                q[i, k] = np.sum(code[i, k]>code.ravel())/N
    return q

Testing on one sample:

In [9]:
%%time
sparse_code_mp_ = rescaling(sparse_code_mp[i_sample:(i_sample+1), :], do_sym=do_sym)
CPU times: user 275 µs, sys: 10 µs, total: 285 µs
Wall time: 281 µs
In [10]:
def rescaling(code, do_sym=do_sym):
    if do_sym: code = np.abs(code)
    n_samples, n_dictionary = code.shape
    N = n_samples * n_dictionary
    q = np.sum(code[:, :, np.newaxis] > code.ravel()[np.newaxis, np.newaxis, :], axis=-1)
    return q
In [11]:
%%time
sparse_code_mp_ = rescaling(sparse_code_mp[i_sample:(i_sample+1), :], do_sym=do_sym)
CPU times: user 590 µs, sys: 545 µs, total: 1.13 ms
Wall time: 748 µs
In [12]:
def get_P_cum(code, nb_quant=512, do_rescale=True, do_sym=do_sym):
    if do_rescale:
        p_c = rescaling(code, do_sym=do_sym)

    n_samples, nb_filter = code.shape
    code_bins = np.linspace(0., 1, nb_quant, endpoint=True)
    P_cum = np.zeros((nb_filter, nb_quant))
    for i in range(nb_filter):
        p, bins = np.histogram(p_c[:, i], bins=code_bins, density=True)
        p /= p.sum()
        P_cum[i, :] = np.hstack((0, np.cumsum(p)))
    return P_cum
In [13]:
%%time
P_cum = get_P_cum(sparse_code_mp[i_sample:(i_sample+1), :], nb_quant=nb_quant, do_sym=do_sym)
CPU times: user 65.6 ms, sys: 2.27 ms, total: 67.9 ms
Wall time: 68.1 ms
P_cum = get_P_cum(sparse_code_mp[i_sample:(i_sample+1), :], nb_quant=nb_quant, do_sym=do_sym) from shl_scripts.shl_tools import plot_P_cum fig, ax = plot_P_cum(P_cum, verbose=False, alpha=.05) ax.plot([0, 1], [0, 1], 'r--') #ax.set_ylim(0.8, 1.01);ax.set_xlim(0.8, 1.01);

COMP : learning modulations on linear coefficients using ranks

Sparse coefficients are distributed according to:

In [14]:
print('Min=', np.min(sparse_code_mp), 'Max=', np.max(sparse_code_mp), 'Shape=', sparse_code_mp.shape)
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
ax.plot(np.sort(sparse_code_mp.ravel()))
ax.set_xlabel('rank'); ax.set_ylabel('coefficient');
Min= 0.0 Max= 36.576175154 Shape= (81792, 324)
No description has been provided for this image

and linear (rectified) coefficients to:

In [15]:
corr = (data_test @ dico_partial_learning.dictionary.T)
corr *= corr>0

print('Min=', np.min(corr), 'Max=', np.max(corr), 'Shape=', corr.shape)

fig, ax = plt.subplots(1, 1, figsize=(8, 5))
ax.plot(np.sort(corr.ravel()))
ax.set_xlabel('rank'); ax.set_ylabel('coefficient');
Min= -0.0 Max= 36.576175154 Shape= (81792, 324)
No description has been provided for this image

The rescaling functions aims at achieving a better numerical stability, such that the quantification won't affect the learning. But it fundamentally does not change the learning algorithm. The later choice (linear rectified coefficients) seems to be a more robust choice and is similar to the prior choice of an exponential function in a past implementation.

In [16]:
def get_rescaling(code, nb_quant=nb_quant, do_sym=do_sym, verbose=False):
#     if do_sym:
#         code = np.abs(code)
#     else:
#         code *= code>0
    sorted_coeffs = np.sort(code.ravel())
    indices = [int(q*(sorted_coeffs.size-1) ) for q in np.linspace(0, 1, nb_quant, endpoint=True)]
    C = sorted_coeffs[indices]
    if verbose:
        print ('At indices (ranks) ', indices)
        print ('the coefficients are ', C )
    return C

C = get_rescaling(corr, nb_quant=nb_quant, do_sym=do_sym, verbose=True)
At indices (ranks)  [0, 51860, 103720, 155580, 207441, 259301, 311161, 363022, 414882, 466742, 518602, 570463, 622323, 674183, 726044, 777904, 829764, 881624, 933485, 985345, 1037205, 1089066, 1140926, 1192786, 1244646, 1296507, 1348367, 1400227, 1452088, 1503948, 1555808, 1607668, 1659529, 1711389, 1763249, 1815110, 1866970, 1918830, 1970690, 2022551, 2074411, 2126271, 2178132, 2229992, 2281852, 2333712, 2385573, 2437433, 2489293, 2541154, 2593014, 2644874, 2696734, 2748595, 2800455, 2852315, 2904176, 2956036, 3007896, 3059756, 3111617, 3163477, 3215337, 3267198, 3319058, 3370918, 3422778, 3474639, 3526499, 3578359, 3630220, 3682080, 3733940, 3785801, 3837661, 3889521, 3941381, 3993242, 4045102, 4096962, 4148823, 4200683, 4252543, 4304403, 4356264, 4408124, 4459984, 4511845, 4563705, 4615565, 4667425, 4719286, 4771146, 4823006, 4874867, 4926727, 4978587, 5030447, 5082308, 5134168, 5186028, 5237889, 5289749, 5341609, 5393469, 5445330, 5497190, 5549050, 5600911, 5652771, 5704631, 5756491, 5808352, 5860212, 5912072, 5963933, 6015793, 6067653, 6119513, 6171374, 6223234, 6275094, 6326955, 6378815, 6430675, 6482535, 6534396, 6586256, 6638116, 6689977, 6741837, 6793697, 6845557, 6897418, 6949278, 7001138, 7052999, 7104859, 7156719, 7208579, 7260440, 7312300, 7364160, 7416021, 7467881, 7519741, 7571602, 7623462, 7675322, 7727182, 7779043, 7830903, 7882763, 7934624, 7986484, 8038344, 8090204, 8142065, 8193925, 8245785, 8297646, 8349506, 8401366, 8453226, 8505087, 8556947, 8608807, 8660668, 8712528, 8764388, 8816248, 8868109, 8919969, 8971829, 9023690, 9075550, 9127410, 9179270, 9231131, 9282991, 9334851, 9386712, 9438572, 9490432, 9542292, 9594153, 9646013, 9697873, 9749734, 9801594, 9853454, 9905314, 9957175, 10009035, 10060895, 10112756, 10164616, 10216476, 10268336, 10320197, 10372057, 10423917, 10475778, 10527638, 10579498, 10631358, 10683219, 10735079, 10786939, 10838800, 10890660, 10942520, 10994380, 11046241, 11098101, 11149961, 11201822, 11253682, 11305542, 11357403, 11409263, 11461123, 11512983, 11564844, 11616704, 11668564, 11720425, 11772285, 11824145, 11876005, 11927866, 11979726, 12031586, 12083447, 12135307, 12187167, 12239027, 12290888, 12342748, 12394608, 12446469, 12498329, 12550189, 12602049, 12653910, 12705770, 12757630, 12809491, 12861351, 12913211, 12965071, 13016932, 13068792, 13120652, 13172513, 13224373, 13276233, 13328093, 13379954, 13431814, 13483674, 13535535, 13587395, 13639255, 13691115, 13742976, 13794836, 13846696, 13898557, 13950417, 14002277, 14054137, 14105998, 14157858, 14209718, 14261579, 14313439, 14365299, 14417159, 14469020, 14520880, 14572740, 14624601, 14676461, 14728321, 14780181, 14832042, 14883902, 14935762, 14987623, 15039483, 15091343, 15143204, 15195064, 15246924, 15298784, 15350645, 15402505, 15454365, 15506226, 15558086, 15609946, 15661806, 15713667, 15765527, 15817387, 15869248, 15921108, 15972968, 16024828, 16076689, 16128549, 16180409, 16232270, 16284130, 16335990, 16387850, 16439711, 16491571, 16543431, 16595292, 16647152, 16699012, 16750872, 16802733, 16854593, 16906453, 16958314, 17010174, 17062034, 17113894, 17165755, 17217615, 17269475, 17321336, 17373196, 17425056, 17476916, 17528777, 17580637, 17632497, 17684358, 17736218, 17788078, 17839938, 17891799, 17943659, 17995519, 18047380, 18099240, 18151100, 18202960, 18254821, 18306681, 18358541, 18410402, 18462262, 18514122, 18565982, 18617843, 18669703, 18721563, 18773424, 18825284, 18877144, 18929004, 18980865, 19032725, 19084585, 19136446, 19188306, 19240166, 19292027, 19343887, 19395747, 19447607, 19499468, 19551328, 19603188, 19655049, 19706909, 19758769, 19810629, 19862490, 19914350, 19966210, 20018071, 20069931, 20121791, 20173651, 20225512, 20277372, 20329232, 20381093, 20432953, 20484813, 20536673, 20588534, 20640394, 20692254, 20744115, 20795975, 20847835, 20899695, 20951556, 21003416, 21055276, 21107137, 21158997, 21210857, 21262717, 21314578, 21366438, 21418298, 21470159, 21522019, 21573879, 21625739, 21677600, 21729460, 21781320, 21833181, 21885041, 21936901, 21988761, 22040622, 22092482, 22144342, 22196203, 22248063, 22299923, 22351783, 22403644, 22455504, 22507364, 22559225, 22611085, 22662945, 22714806, 22766666, 22818526, 22870386, 22922247, 22974107, 23025967, 23077828, 23129688, 23181548, 23233408, 23285269, 23337129, 23388989, 23440850, 23492710, 23544570, 23596430, 23648291, 23700151, 23752011, 23803872, 23855732, 23907592, 23959452, 24011313, 24063173, 24115033, 24166894, 24218754, 24270614, 24322474, 24374335, 24426195, 24478055, 24529916, 24581776, 24633636, 24685496, 24737357, 24789217, 24841077, 24892938, 24944798, 24996658, 25048518, 25100379, 25152239, 25204099, 25255960, 25307820, 25359680, 25411540, 25463401, 25515261, 25567121, 25618982, 25670842, 25722702, 25774562, 25826423, 25878283, 25930143, 25982004, 26033864, 26085724, 26137584, 26189445, 26241305, 26293165, 26345026, 26396886, 26448746, 26500607]
the coefficients are  [-0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 0.000709 0.00149 0.00227 0.00307 0.00386 0.00461
 0.00537 0.00613 0.00693 0.00772 0.00851 0.00929 0.0101 0.0109 0.0117
 0.0126 0.0134 0.0142 0.0151 0.016 0.0169 0.0178 0.0187 0.0196 0.0206
 0.0215 0.0225 0.0235 0.0245 0.0255 0.0265 0.0275 0.0286 0.0297 0.0308
 0.0321 0.0333 0.0345 0.0358 0.0371 0.0384 0.0398 0.0413 0.0428 0.0443
 0.046 0.0477 0.0494 0.0511 0.0529 0.0547 0.0567 0.0587 0.0607 0.0628 0.065
 0.0671 0.0693 0.0717 0.0741 0.0766 0.0792 0.0819 0.0846 0.0873 0.0901
 0.093 0.0959 0.0989 0.102 0.105 0.108 0.111 0.115 0.118 0.121 0.125 0.128
 0.132 0.135 0.139 0.143 0.146 0.15 0.154 0.158 0.162 0.166 0.17 0.174
 0.178 0.183 0.187 0.191 0.195 0.2 0.204 0.209 0.214 0.218 0.223 0.228
 0.233 0.238 0.243 0.248 0.253 0.258 0.263 0.269 0.274 0.28 0.285 0.291
 0.296 0.302 0.308 0.314 0.32 0.326 0.332 0.338 0.344 0.351 0.357 0.363
 0.37 0.377 0.384 0.39 0.397 0.405 0.412 0.419 0.426 0.434 0.441 0.449
 0.457 0.464 0.472 0.48 0.489 0.497 0.505 0.514 0.522 0.531 0.54 0.549
 0.558 0.567 0.577 0.586 0.596 0.606 0.615 0.626 0.636 0.646 0.657 0.667
 0.678 0.689 0.701 0.712 0.723 0.735 0.747 0.759 0.772 0.784 0.797 0.81
 0.823 0.836 0.85 0.864 0.878 0.893 0.907 0.922 0.937 0.952 0.968 0.984 1
 1.02 1.03 1.05 1.07 1.09 1.11 1.12 1.14 1.16 1.18 1.2 1.22 1.25 1.27 1.29
 1.31 1.34 1.36 1.39 1.41 1.44 1.46 1.49 1.52 1.55 1.58 1.61 1.64 1.67 1.7
 1.74 1.77 1.81 1.85 1.89 1.93 1.97 2.01 2.06 2.1 2.15 2.2 2.26 2.31 2.37
 2.43 2.5 2.57 2.64 2.71 2.8 2.88 2.98 3.08 3.19 3.31 3.44 3.59 3.75 3.94
 4.16 4.42 4.74 5.17 5.8 6.91 36.6]
In [17]:
def rescaling(code, C, do_sym=do_sym):
#     if do_sym:
#         code = np.abs(code)
#     else:
#         code *= code>0
    n_samples, n_dictionary = code.shape
    N = n_samples * n_dictionary
    q = np.zeros_like(code)
    code_bins = np.linspace(0., 1, C.size, endpoint=True)

    for i in range(n_samples):
        for k in range(n_dictionary):
            if code[i, k] > 0:
                q[i, k] = np.interp(code[i, k], C, code_bins)
    return q

rescaled_code_ = rescaling(corr[i_sample:(i_sample+1), :], C=C, do_sym=do_sym)
In [18]:
%%time
rescaled_code_ = rescaling(corr[i_sample:(i_sample+1), :], C=C, do_sym=do_sym)
CPU times: user 1.89 ms, sys: 1.36 ms, total: 3.25 ms
Wall time: 1.98 ms

This can be vectorized:

In [19]:
def rescaling(code, C, do_sym=do_sym):
#     if do_sym:
#         code = np.abs(code)
#     else:
#         code *= code>0
    code_bins = np.linspace(0., 1, C.size, endpoint=True)
    return np.interp(code, C, code_bins) * (code > 0.)

print('on ', corr.shape[0], ' samples')
rescaled_code_ = rescaling(corr, C=C, do_sym=do_sym)
on  81792  samples
In [20]:
%%time
sparse_code_mp_ = rescaling(sparse_code_mp, C=C, do_sym=do_sym)
CPU times: user 444 ms, sys: 117 ms, total: 561 ms
Wall time: 635 ms

We notice also that we use half of the vector for negative coefficients. This could be improved as:

In [21]:
def rescaling(code, C, do_sym=do_sym):
#     if do_sym:
#         code = np.abs(code)
#     else:
#         code *= code>0
    p_c = np.zeros_like(code)
    ind_nz = code>0.
    code_bins = np.linspace(0., 1, C.size, endpoint=True)
    p_c[ind_nz] = np.interp(code[ind_nz], C, code_bins)
    return p_c # * (code > 0.)

rescaled_code__ = rescaling(corr, C=C, do_sym=do_sym)
assert(np.array_equal(rescaled_code_, rescaled_code__))
In [22]:
%%time
sparse_code_mp_ = rescaling(sparse_code_mp, C=C, do_sym=do_sym)
CPU times: user 246 ms, sys: 131 ms, total: 377 ms
Wall time: 492 ms

Finally, one can get the P_cum look-up-tables:

In [23]:
def get_P_cum(code, nb_quant, C=C, do_sym=do_sym):
    p_c = rescaling(code, C=C, do_sym=do_sym)

    n_samples, nb_filter = code.shape
    code_bins = np.linspace(0., 1, nb_quant, endpoint=True)
    P_cum = np.zeros((nb_filter, nb_quant))
    for i in range(nb_filter):
        p, bins = np.histogram(p_c[:, i], bins=code_bins, density=True)
        p /= p.sum()
        P_cum[i, :] = np.hstack((0., np.cumsum(p)))
    return P_cum

computing histograms is rather fast:

In [24]:
%%time
P_cum = get_P_cum(sparse_code_mp, nb_quant=nb_quant, C=C)
CPU times: user 1.23 s, sys: 126 ms, total: 1.36 s
Wall time: 1.44 s
In [25]:
C = get_rescaling(corr, nb_quant=nb_quant, do_sym=do_sym, verbose=True)
C_vec_init = C.copy()
fig=plt.figure(figsize=(13, 8))
ax = plt.subplot(111)
x = np.linspace(0, 30, 100)
ax.plot(x, np.ones_like(x), '--')
ax.plot(x, rescaling(x, C=C, do_sym=do_sym))
ax.set_xlabel('raw coefficient')
ax.set_ylabel('rescaled coefficient')
ax.set_xscale('log');
At indices (ranks)  [0, 51860, 103720, 155580, 207441, 259301, 311161, 363022, 414882, 466742, 518602, 570463, 622323, 674183, 726044, 777904, 829764, 881624, 933485, 985345, 1037205, 1089066, 1140926, 1192786, 1244646, 1296507, 1348367, 1400227, 1452088, 1503948, 1555808, 1607668, 1659529, 1711389, 1763249, 1815110, 1866970, 1918830, 1970690, 2022551, 2074411, 2126271, 2178132, 2229992, 2281852, 2333712, 2385573, 2437433, 2489293, 2541154, 2593014, 2644874, 2696734, 2748595, 2800455, 2852315, 2904176, 2956036, 3007896, 3059756, 3111617, 3163477, 3215337, 3267198, 3319058, 3370918, 3422778, 3474639, 3526499, 3578359, 3630220, 3682080, 3733940, 3785801, 3837661, 3889521, 3941381, 3993242, 4045102, 4096962, 4148823, 4200683, 4252543, 4304403, 4356264, 4408124, 4459984, 4511845, 4563705, 4615565, 4667425, 4719286, 4771146, 4823006, 4874867, 4926727, 4978587, 5030447, 5082308, 5134168, 5186028, 5237889, 5289749, 5341609, 5393469, 5445330, 5497190, 5549050, 5600911, 5652771, 5704631, 5756491, 5808352, 5860212, 5912072, 5963933, 6015793, 6067653, 6119513, 6171374, 6223234, 6275094, 6326955, 6378815, 6430675, 6482535, 6534396, 6586256, 6638116, 6689977, 6741837, 6793697, 6845557, 6897418, 6949278, 7001138, 7052999, 7104859, 7156719, 7208579, 7260440, 7312300, 7364160, 7416021, 7467881, 7519741, 7571602, 7623462, 7675322, 7727182, 7779043, 7830903, 7882763, 7934624, 7986484, 8038344, 8090204, 8142065, 8193925, 8245785, 8297646, 8349506, 8401366, 8453226, 8505087, 8556947, 8608807, 8660668, 8712528, 8764388, 8816248, 8868109, 8919969, 8971829, 9023690, 9075550, 9127410, 9179270, 9231131, 9282991, 9334851, 9386712, 9438572, 9490432, 9542292, 9594153, 9646013, 9697873, 9749734, 9801594, 9853454, 9905314, 9957175, 10009035, 10060895, 10112756, 10164616, 10216476, 10268336, 10320197, 10372057, 10423917, 10475778, 10527638, 10579498, 10631358, 10683219, 10735079, 10786939, 10838800, 10890660, 10942520, 10994380, 11046241, 11098101, 11149961, 11201822, 11253682, 11305542, 11357403, 11409263, 11461123, 11512983, 11564844, 11616704, 11668564, 11720425, 11772285, 11824145, 11876005, 11927866, 11979726, 12031586, 12083447, 12135307, 12187167, 12239027, 12290888, 12342748, 12394608, 12446469, 12498329, 12550189, 12602049, 12653910, 12705770, 12757630, 12809491, 12861351, 12913211, 12965071, 13016932, 13068792, 13120652, 13172513, 13224373, 13276233, 13328093, 13379954, 13431814, 13483674, 13535535, 13587395, 13639255, 13691115, 13742976, 13794836, 13846696, 13898557, 13950417, 14002277, 14054137, 14105998, 14157858, 14209718, 14261579, 14313439, 14365299, 14417159, 14469020, 14520880, 14572740, 14624601, 14676461, 14728321, 14780181, 14832042, 14883902, 14935762, 14987623, 15039483, 15091343, 15143204, 15195064, 15246924, 15298784, 15350645, 15402505, 15454365, 15506226, 15558086, 15609946, 15661806, 15713667, 15765527, 15817387, 15869248, 15921108, 15972968, 16024828, 16076689, 16128549, 16180409, 16232270, 16284130, 16335990, 16387850, 16439711, 16491571, 16543431, 16595292, 16647152, 16699012, 16750872, 16802733, 16854593, 16906453, 16958314, 17010174, 17062034, 17113894, 17165755, 17217615, 17269475, 17321336, 17373196, 17425056, 17476916, 17528777, 17580637, 17632497, 17684358, 17736218, 17788078, 17839938, 17891799, 17943659, 17995519, 18047380, 18099240, 18151100, 18202960, 18254821, 18306681, 18358541, 18410402, 18462262, 18514122, 18565982, 18617843, 18669703, 18721563, 18773424, 18825284, 18877144, 18929004, 18980865, 19032725, 19084585, 19136446, 19188306, 19240166, 19292027, 19343887, 19395747, 19447607, 19499468, 19551328, 19603188, 19655049, 19706909, 19758769, 19810629, 19862490, 19914350, 19966210, 20018071, 20069931, 20121791, 20173651, 20225512, 20277372, 20329232, 20381093, 20432953, 20484813, 20536673, 20588534, 20640394, 20692254, 20744115, 20795975, 20847835, 20899695, 20951556, 21003416, 21055276, 21107137, 21158997, 21210857, 21262717, 21314578, 21366438, 21418298, 21470159, 21522019, 21573879, 21625739, 21677600, 21729460, 21781320, 21833181, 21885041, 21936901, 21988761, 22040622, 22092482, 22144342, 22196203, 22248063, 22299923, 22351783, 22403644, 22455504, 22507364, 22559225, 22611085, 22662945, 22714806, 22766666, 22818526, 22870386, 22922247, 22974107, 23025967, 23077828, 23129688, 23181548, 23233408, 23285269, 23337129, 23388989, 23440850, 23492710, 23544570, 23596430, 23648291, 23700151, 23752011, 23803872, 23855732, 23907592, 23959452, 24011313, 24063173, 24115033, 24166894, 24218754, 24270614, 24322474, 24374335, 24426195, 24478055, 24529916, 24581776, 24633636, 24685496, 24737357, 24789217, 24841077, 24892938, 24944798, 24996658, 25048518, 25100379, 25152239, 25204099, 25255960, 25307820, 25359680, 25411540, 25463401, 25515261, 25567121, 25618982, 25670842, 25722702, 25774562, 25826423, 25878283, 25930143, 25982004, 26033864, 26085724, 26137584, 26189445, 26241305, 26293165, 26345026, 26396886, 26448746, 26500607]
the coefficients are  [-0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
 -0 -0 -0 -0 -0 -0 -0 -0 0.000709 0.00149 0.00227 0.00307 0.00386 0.00461
 0.00537 0.00613 0.00693 0.00772 0.00851 0.00929 0.0101 0.0109 0.0117
 0.0126 0.0134 0.0142 0.0151 0.016 0.0169 0.0178 0.0187 0.0196 0.0206
 0.0215 0.0225 0.0235 0.0245 0.0255 0.0265 0.0275 0.0286 0.0297 0.0308
 0.0321 0.0333 0.0345 0.0358 0.0371 0.0384 0.0398 0.0413 0.0428 0.0443
 0.046 0.0477 0.0494 0.0511 0.0529 0.0547 0.0567 0.0587 0.0607 0.0628 0.065
 0.0671 0.0693 0.0717 0.0741 0.0766 0.0792 0.0819 0.0846 0.0873 0.0901
 0.093 0.0959 0.0989 0.102 0.105 0.108 0.111 0.115 0.118 0.121 0.125 0.128
 0.132 0.135 0.139 0.143 0.146 0.15 0.154 0.158 0.162 0.166 0.17 0.174
 0.178 0.183 0.187 0.191 0.195 0.2 0.204 0.209 0.214 0.218 0.223 0.228
 0.233 0.238 0.243 0.248 0.253 0.258 0.263 0.269 0.274 0.28 0.285 0.291
 0.296 0.302 0.308 0.314 0.32 0.326 0.332 0.338 0.344 0.351 0.357 0.363
 0.37 0.377 0.384 0.39 0.397 0.405 0.412 0.419 0.426 0.434 0.441 0.449
 0.457 0.464 0.472 0.48 0.489 0.497 0.505 0.514 0.522 0.531 0.54 0.549
 0.558 0.567 0.577 0.586 0.596 0.606 0.615 0.626 0.636 0.646 0.657 0.667
 0.678 0.689 0.701 0.712 0.723 0.735 0.747 0.759 0.772 0.784 0.797 0.81
 0.823 0.836 0.85 0.864 0.878 0.893 0.907 0.922 0.937 0.952 0.968 0.984 1
 1.02 1.03 1.05 1.07 1.09 1.11 1.12 1.14 1.16 1.18 1.2 1.22 1.25 1.27 1.29
 1.31 1.34 1.36 1.39 1.41 1.44 1.46 1.49 1.52 1.55 1.58 1.61 1.64 1.67 1.7
 1.74 1.77 1.81 1.85 1.89 1.93 1.97 2.01 2.06 2.1 2.15 2.2 2.26 2.31 2.37
 2.43 2.5 2.57 2.64 2.71 2.8 2.88 2.98 3.08 3.19 3.31 3.44 3.59 3.75 3.94
 4.16 4.42 4.74 5.17 5.8 6.91 36.6]
No description has been provided for this image
In [26]:
P_cum = get_P_cum(rescaling(sparse_code_mp, C=C, do_sym=do_sym), nb_quant=nb_quant, C=C, do_sym=do_sym)
from shl_scripts.shl_tools import plot_P_cum
fig, ax = plot_P_cum(P_cum, verbose=False, alpha=.05)
ax.plot([0, 1], [0, 1], 'r--');
ax.set_ylim(0.85, 1.01);ax.set_xlim(0., 1.);
No description has been provided for this image

Finally, we notice that the C vector is of the same length as each P_cum vector and we will integrate it to this array to ease message passing.

COMP : optimizing the quantile function

These quantile and rescaling functions are implanted in the shl_scripts:

In [27]:
from shl_scripts.shl_encode import quantile, rescaling

corr = (data_test @ dico_partial_learning.dictionary.T)
#corr *= corr>0
print('correlation=', corr[i_sample, :])
print('transformed correlation=', rescaling(corr, C=C, do_sym=do_sym)[i_sample, :])
correlation= [-0.606 0.428 -0.566 -0.457 0.304 0.687 0.55 -1.84 0.43 0.755 -0.286 0.173
 0.365 -0.185 0.615 0.263 0.614 0.819 -1.68 1.78 0.0753 -0.536 0.404 -1.77
 1.14 0.117 0.625 -0.296 -0.0527 0.591 0.175 -0.639 0.812 0.0353 0.825 1.07
 1.97 -0.282 0.187 -0.305 -1.18 0.818 -0.2 1.31 -0.218 -1.13 0.809 0.509
 0.069 -0.673 -0.131 -0.165 -0.0833 -1.9 -0.244 -0.414 -0.0485 0.718 -0.154
 0.374 0.0455 1.17 -0.048 1.44 0.602 0.00801 -1.67 1.16 0.307 -0.152 1.36
 -0.173 0.444 -0.175 0.561 1.35 -0.738 1.36 2.26 0.855 0.666 0.527 -0.456
 1.3 0.556 -0.217 -0.28 1.52 0.744 -0.307 -1.02 0.891 6.32e-05 1.09 -1.94
 0.568 -0.473 1.05 0.245 0.163 1.13 -0.0183 0.686 0.543 0.379 0.6 1.38
 0.456 2.21 -1.09 0.411 -0.872 2.18 -0.301 0.681 1.35 0.568 -0.337 1.66
 -1.09 -0.486 -0.598 -0.241 1.52 -0.126 0.42 -0.0964 -1.19 1.09 0.0617
 0.542 -0.0428 -1.21 0.549 -0.385 0.2 0.0937 -0.74 -0.845 -0.783 0.0349
 0.973 1.13 -0.306 -1.02 0.0624 0.0613 -0.278 0.967 -0.658 0.579 -0.0586
 -0.463 -0.726 0.412 -0.0301 0.929 0.765 -0.3 -1.11 -0.222 -0.213 -0.21
 0.231 -0.907 -0.117 0.0441 0.918 -0.195 0.603 -0.661 0.419 -0.522 -0.182
 1.41 0.414 -0.331 -0.256 0.636 0.407 0.156 0.0501 1.68 0.811 1.4 -0.833
 1.01 1.8 0.263 1.28 -0.559 -0.0339 1.19 0.828 -0.212 2.02 0.708 -0.0837
 -0.791 -1 0.226 -0.353 1.26 1.18 0.64 0.21 0.61 -0.0966 -1.19 0.58 0.0121
 0.283 -1.74 0.911 0.98 -0.381 0.11 0.359 0.502 -0.574 -0.475 0.0498 0.0972
 0.775 -1.33 -0.757 -0.0152 1.11 -0.464 0.0499 -1.77 1.81 1.26 0.227 -0.858
 0.734 -0.228 -0.602 1.15 1.47 1.52 1.4 1.02 -0.039 0.217 0.654 2.4 0.481
 0.161 -0.47 0.626 0.356 0.708 0.68 0.992 -0.364 -0.106 0.948 -0.418 0.712
 0.529 -0.207 1.03 -0.315 -0.741 0.481 -0.457 -1.01 -0.638 0.423 -0.278
 1.09 1.46 1.95 -0.22 0.571 -0.524 -0.645 -1.01 1.1 0.0848 1.58 0.34 0.143
 0.187 0.363 1.24 0.583 -0.577 0.118 -0.611 1.7 -2.21 0.435 0.0115 0.717
 0.984 0.666 -1.57 1.36 -0.272 -0.204 -0.0694 0.885 0.157 0.598 0.687
 -0.432 0.737 -0.863 0.506 -0.223 -1.04 -0.0878 0.246 -0.237 0.823 -0.371
 -0.939 1.31 0.642 -0.593 0.0982 1.84]
transformed correlation= [0 0.768 0 0 0.73 0.823 0.797 0 0.768 0.835 0 0.679 0.75 0 0.81 0.716 0.81
 0.845 0 0.932 0.621 0 0.761 0 0.884 0.649 0.812 0 0 0.805 0.68 0 0.844
 0.578 0.846 0.876 0.941 0 0.685 0 0 0.845 0 0.9 0 0 0.843 0.788 0.616 0 0
 0 0 0 0 0 0 0.829 0 0.753 0.592 0.887 0 0.91 0.808 0.523 0 0.887 0.732 0
 0.904 0 0.772 0 0.799 0.904 0 0.904 0.953 0.85 0.82 0.792 0 0.899 0.798 0
 0 0.916 0.833 0 0 0.855 0.503 0.879 0 0.801 0 0.874 0.709 0.674 0.883 0
 0.823 0.795 0.754 0.807 0.906 0.775 0.951 0 0.763 0 0.95 0 0.822 0.903
 0.801 0 0.925 0 0 0 0 0.916 0 0.765 0 0 0.879 0.61 0.795 0 0 0.796 0 0.691
 0.635 0 0 0 0.578 0.866 0.883 0 0 0.61 0.609 0 0.865 0 0.803 0 0 0 0.763 0
 0.86 0.837 0 0 0 0 0 0.704 0 0 0.591 0.859 0 0.808 0 0.765 0 0 0.908 0.764
 0 0 0.814 0.762 0.67 0.598 0.926 0.844 0.907 0 0.87 0.933 0.716 0.897 0 0
 0.889 0.846 0 0.944 0.827 0 0 0 0.702 0 0.895 0.889 0.815 0.695 0.809 0 0
 0.803 0.533 0.723 0 0.858 0.866 0 0.645 0.748 0.786 0 0 0.597 0.637 0.838
 0 0 0 0.881 0 0.597 0 0.934 0.895 0.702 0 0.832 0 0 0.885 0.912 0.916
 0.908 0.871 0 0.698 0.817 0.958 0.781 0.673 0 0.812 0.747 0.827 0.822
 0.868 0 0 0.862 0 0.828 0.792 0 0.872 0 0 0.781 0 0 0 0.766 0 0.879 0.912
 0.94 0 0.801 0 0 0 0.88 0.628 0.92 0.742 0.664 0.685 0.749 0.894 0.804 0
 0.649 0 0.927 0 0.769 0.532 0.829 0.867 0.82 0 0.904 0 0 0 0.854 0.671
 0.807 0.824 0 0.832 0 0.787 0 0 0 0.71 0 0.845 0 0 0.9 0.815 0 0.638 0.935]
In [28]:
C = get_rescaling(corr, nb_quant=nb_quant, do_sym=do_sym, verbose=True)
fig=plt.figure(figsize=(13, 8))
ax = plt.subplot(111)
x = np.linspace(0, 30, 100)
ax.plot(x, np.ones_like(x), '--')
ax.plot(x, rescaling(x, C=C_vec_init, do_sym=do_sym))
ax.plot(x, rescaling(x, C=.5, do_sym=do_sym))
ax.set_xlabel('raw coefficient')
ax.set_ylabel('rescaled coefficient')
ax.set_xscale('log')
ax.axis('tight');
At indices (ranks)  [0, 51860, 103720, 155580, 207441, 259301, 311161, 363022, 414882, 466742, 518602, 570463, 622323, 674183, 726044, 777904, 829764, 881624, 933485, 985345, 1037205, 1089066, 1140926, 1192786, 1244646, 1296507, 1348367, 1400227, 1452088, 1503948, 1555808, 1607668, 1659529, 1711389, 1763249, 1815110, 1866970, 1918830, 1970690, 2022551, 2074411, 2126271, 2178132, 2229992, 2281852, 2333712, 2385573, 2437433, 2489293, 2541154, 2593014, 2644874, 2696734, 2748595, 2800455, 2852315, 2904176, 2956036, 3007896, 3059756, 3111617, 3163477, 3215337, 3267198, 3319058, 3370918, 3422778, 3474639, 3526499, 3578359, 3630220, 3682080, 3733940, 3785801, 3837661, 3889521, 3941381, 3993242, 4045102, 4096962, 4148823, 4200683, 4252543, 4304403, 4356264, 4408124, 4459984, 4511845, 4563705, 4615565, 4667425, 4719286, 4771146, 4823006, 4874867, 4926727, 4978587, 5030447, 5082308, 5134168, 5186028, 5237889, 5289749, 5341609, 5393469, 5445330, 5497190, 5549050, 5600911, 5652771, 5704631, 5756491, 5808352, 5860212, 5912072, 5963933, 6015793, 6067653, 6119513, 6171374, 6223234, 6275094, 6326955, 6378815, 6430675, 6482535, 6534396, 6586256, 6638116, 6689977, 6741837, 6793697, 6845557, 6897418, 6949278, 7001138, 7052999, 7104859, 7156719, 7208579, 7260440, 7312300, 7364160, 7416021, 7467881, 7519741, 7571602, 7623462, 7675322, 7727182, 7779043, 7830903, 7882763, 7934624, 7986484, 8038344, 8090204, 8142065, 8193925, 8245785, 8297646, 8349506, 8401366, 8453226, 8505087, 8556947, 8608807, 8660668, 8712528, 8764388, 8816248, 8868109, 8919969, 8971829, 9023690, 9075550, 9127410, 9179270, 9231131, 9282991, 9334851, 9386712, 9438572, 9490432, 9542292, 9594153, 9646013, 9697873, 9749734, 9801594, 9853454, 9905314, 9957175, 10009035, 10060895, 10112756, 10164616, 10216476, 10268336, 10320197, 10372057, 10423917, 10475778, 10527638, 10579498, 10631358, 10683219, 10735079, 10786939, 10838800, 10890660, 10942520, 10994380, 11046241, 11098101, 11149961, 11201822, 11253682, 11305542, 11357403, 11409263, 11461123, 11512983, 11564844, 11616704, 11668564, 11720425, 11772285, 11824145, 11876005, 11927866, 11979726, 12031586, 12083447, 12135307, 12187167, 12239027, 12290888, 12342748, 12394608, 12446469, 12498329, 12550189, 12602049, 12653910, 12705770, 12757630, 12809491, 12861351, 12913211, 12965071, 13016932, 13068792, 13120652, 13172513, 13224373, 13276233, 13328093, 13379954, 13431814, 13483674, 13535535, 13587395, 13639255, 13691115, 13742976, 13794836, 13846696, 13898557, 13950417, 14002277, 14054137, 14105998, 14157858, 14209718, 14261579, 14313439, 14365299, 14417159, 14469020, 14520880, 14572740, 14624601, 14676461, 14728321, 14780181, 14832042, 14883902, 14935762, 14987623, 15039483, 15091343, 15143204, 15195064, 15246924, 15298784, 15350645, 15402505, 15454365, 15506226, 15558086, 15609946, 15661806, 15713667, 15765527, 15817387, 15869248, 15921108, 15972968, 16024828, 16076689, 16128549, 16180409, 16232270, 16284130, 16335990, 16387850, 16439711, 16491571, 16543431, 16595292, 16647152, 16699012, 16750872, 16802733, 16854593, 16906453, 16958314, 17010174, 17062034, 17113894, 17165755, 17217615, 17269475, 17321336, 17373196, 17425056, 17476916, 17528777, 17580637, 17632497, 17684358, 17736218, 17788078, 17839938, 17891799, 17943659, 17995519, 18047380, 18099240, 18151100, 18202960, 18254821, 18306681, 18358541, 18410402, 18462262, 18514122, 18565982, 18617843, 18669703, 18721563, 18773424, 18825284, 18877144, 18929004, 18980865, 19032725, 19084585, 19136446, 19188306, 19240166, 19292027, 19343887, 19395747, 19447607, 19499468, 19551328, 19603188, 19655049, 19706909, 19758769, 19810629, 19862490, 19914350, 19966210, 20018071, 20069931, 20121791, 20173651, 20225512, 20277372, 20329232, 20381093, 20432953, 20484813, 20536673, 20588534, 20640394, 20692254, 20744115, 20795975, 20847835, 20899695, 20951556, 21003416, 21055276, 21107137, 21158997, 21210857, 21262717, 21314578, 21366438, 21418298, 21470159, 21522019, 21573879, 21625739, 21677600, 21729460, 21781320, 21833181, 21885041, 21936901, 21988761, 22040622, 22092482, 22144342, 22196203, 22248063, 22299923, 22351783, 22403644, 22455504, 22507364, 22559225, 22611085, 22662945, 22714806, 22766666, 22818526, 22870386, 22922247, 22974107, 23025967, 23077828, 23129688, 23181548, 23233408, 23285269, 23337129, 23388989, 23440850, 23492710, 23544570, 23596430, 23648291, 23700151, 23752011, 23803872, 23855732, 23907592, 23959452, 24011313, 24063173, 24115033, 24166894, 24218754, 24270614, 24322474, 24374335, 24426195, 24478055, 24529916, 24581776, 24633636, 24685496, 24737357, 24789217, 24841077, 24892938, 24944798, 24996658, 25048518, 25100379, 25152239, 25204099, 25255960, 25307820, 25359680, 25411540, 25463401, 25515261, 25567121, 25618982, 25670842, 25722702, 25774562, 25826423, 25878283, 25930143, 25982004, 26033864, 26085724, 26137584, 26189445, 26241305, 26293165, 26345026, 26396886, 26448746, 26500607]
the coefficients are  [-36.2 -6.93 -5.81 -5.19 -4.77 -4.44 -4.17 -3.95 -3.76 -3.59 -3.45 -3.31
 -3.19 -3.08 -2.98 -2.89 -2.8 -2.72 -2.64 -2.57 -2.5 -2.43 -2.37 -2.31
 -2.26 -2.2 -2.15 -2.1 -2.05 -2.01 -1.97 -1.92 -1.88 -1.84 -1.81 -1.77
 -1.73 -1.7 -1.67 -1.63 -1.6 -1.57 -1.54 -1.51 -1.49 -1.46 -1.43 -1.41
 -1.38 -1.36 -1.33 -1.31 -1.29 -1.26 -1.24 -1.22 -1.2 -1.18 -1.16 -1.14
 -1.12 -1.1 -1.08 -1.07 -1.05 -1.03 -1.01 -0.997 -0.98 -0.964 -0.949 -0.933
 -0.918 -0.903 -0.888 -0.874 -0.86 -0.846 -0.832 -0.819 -0.806 -0.793 -0.78
 -0.768 -0.755 -0.743 -0.731 -0.72 -0.708 -0.697 -0.686 -0.675 -0.664
 -0.653 -0.643 -0.632 -0.622 -0.612 -0.602 -0.592 -0.583 -0.573 -0.564
 -0.555 -0.546 -0.537 -0.528 -0.519 -0.511 -0.502 -0.494 -0.486 -0.478
 -0.47 -0.462 -0.454 -0.447 -0.439 -0.432 -0.424 -0.417 -0.41 -0.403 -0.396
 -0.389 -0.382 -0.375 -0.369 -0.362 -0.356 -0.349 -0.343 -0.337 -0.331
 -0.325 -0.319 -0.313 -0.307 -0.301 -0.295 -0.29 -0.284 -0.279 -0.273
 -0.268 -0.263 -0.258 -0.252 -0.247 -0.242 -0.238 -0.233 -0.228 -0.223
 -0.218 -0.214 -0.209 -0.205 -0.2 -0.196 -0.191 -0.187 -0.183 -0.179 -0.175
 -0.17 -0.167 -0.163 -0.159 -0.155 -0.151 -0.147 -0.144 -0.14 -0.136 -0.133
 -0.129 -0.126 -0.123 -0.119 -0.116 -0.113 -0.11 -0.107 -0.103 -0.1 -0.0975
 -0.0946 -0.0918 -0.0891 -0.0864 -0.0837 -0.0811 -0.0785 -0.076 -0.0736
 -0.0712 -0.0689 -0.0667 -0.0645 -0.0623 -0.0603 -0.0583 -0.0564 -0.0546
 -0.0528 -0.0511 -0.0495 -0.0479 -0.0463 -0.0447 -0.0433 -0.0419 -0.0405
 -0.0391 -0.0378 -0.0366 -0.0353 -0.0341 -0.033 -0.0318 -0.0307 -0.0296
 -0.0286 -0.0276 -0.0265 -0.0256 -0.0246 -0.0236 -0.0226 -0.0217 -0.0208
 -0.0199 -0.019 -0.0181 -0.0173 -0.0164 -0.0156 -0.0148 -0.014 -0.0132
 -0.0124 -0.0116 -0.0108 -0.01 -0.00922 -0.00844 -0.00766 -0.00685 -0.00608
 -0.00532 -0.00456 -0.00381 -0.00307 -0.00233 -0.00158 -0.000822 -7.28e-05
 0.000709 0.00149 0.00227 0.00307 0.00386 0.00461 0.00537 0.00613 0.00693
 0.00772 0.00851 0.00929 0.0101 0.0109 0.0117 0.0126 0.0134 0.0142 0.0151
 0.016 0.0169 0.0178 0.0187 0.0196 0.0206 0.0215 0.0225 0.0235 0.0245
 0.0255 0.0265 0.0275 0.0286 0.0297 0.0308 0.0321 0.0333 0.0345 0.0358
 0.0371 0.0384 0.0398 0.0413 0.0428 0.0443 0.046 0.0477 0.0494 0.0511
 0.0529 0.0547 0.0567 0.0587 0.0607 0.0628 0.065 0.0671 0.0693 0.0717
 0.0741 0.0766 0.0792 0.0819 0.0846 0.0873 0.0901 0.093 0.0959 0.0989 0.102
 0.105 0.108 0.111 0.115 0.118 0.121 0.125 0.128 0.132 0.135 0.139 0.143
 0.146 0.15 0.154 0.158 0.162 0.166 0.17 0.174 0.178 0.183 0.187 0.191
 0.195 0.2 0.204 0.209 0.214 0.218 0.223 0.228 0.233 0.238 0.243 0.248
 0.253 0.258 0.263 0.269 0.274 0.28 0.285 0.291 0.296 0.302 0.308 0.314
 0.32 0.326 0.332 0.338 0.344 0.351 0.357 0.363 0.37 0.377 0.384 0.39 0.397
 0.405 0.412 0.419 0.426 0.434 0.441 0.449 0.457 0.464 0.472 0.48 0.489
 0.497 0.505 0.514 0.522 0.531 0.54 0.549 0.558 0.567 0.577 0.586 0.596
 0.606 0.615 0.626 0.636 0.646 0.657 0.667 0.678 0.689 0.701 0.712 0.723
 0.735 0.747 0.759 0.772 0.784 0.797 0.81 0.823 0.836 0.85 0.864 0.878
 0.893 0.907 0.922 0.937 0.952 0.968 0.984 1 1.02 1.03 1.05 1.07 1.09 1.11
 1.12 1.14 1.16 1.18 1.2 1.22 1.25 1.27 1.29 1.31 1.34 1.36 1.39 1.41 1.44
 1.46 1.49 1.52 1.55 1.58 1.61 1.64 1.67 1.7 1.74 1.77 1.81 1.85 1.89 1.93
 1.97 2.01 2.06 2.1 2.15 2.2 2.26 2.31 2.37 2.43 2.5 2.57 2.64 2.71 2.8
 2.88 2.98 3.08 3.19 3.31 3.44 3.59 3.75 3.94 4.16 4.42 4.74 5.17 5.8 6.91
 36.6]
No description has been provided for this image

A sanity check with extremal values:

In [29]:
stick = np.arange(shl.n_dictionary)*nb_quant
print('Value for ones = ', rescaling(np.inf*np.ones(shl.n_dictionary), C=C))
print('Value for zeros = ', rescaling(np.zeros(shl.n_dictionary), C=C))
Value for ones =  [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
Value for zeros =  [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
In [30]:
print('Value for ones = ', quantile(P_cum, rescaling(np.inf*np.ones(shl.n_dictionary), C=C), stick))
print('Value for zeros = ', quantile(P_cum, rescaling(np.zeros(shl.n_dictionary), C=C), stick))
Value for ones =  [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
Value for zeros =  [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

a P_cum array that should not change anything to the way we chose filters:

In [31]:
P_cum_zeroeffect = np.linspace(0, 1, nb_quant, endpoint=True)[np.newaxis, :] * np.ones((n_dictionary, 1))
In [32]:
def plot_scatter_MpVsComp(sparse_vector, my_sparse_code, title='MP', xlabel='MP', ylabel='COMP'):
    fig = plt.figure(figsize=(16, 16))
    ax = fig.add_subplot(111)
    a_min = np.min((sparse_vector.min(), my_sparse_code.min()))
    a_max = np.max((sparse_vector.max(), my_sparse_code.max()))
    ax.plot(np.array([a_min, a_max]), np.array([a_min, a_max]), 'k--', lw=2)
    print(sparse_vector.shape, my_sparse_code.shape)
    ax.scatter(sparse_vector.ravel(), my_sparse_code.ravel(), alpha=0.01)
    ax.set_title(title)
    ax.set_ylabel(ylabel)
    ax.set_xlabel(xlabel)
    #ax.set_xlim(0)
    #ax.set_ylim(0)
    ax.axis('equal')
    return fig, ax

p_c = rescaling(corr, C=C, do_sym=do_sym)
fig, ax = plot_scatter_MpVsComp(p_c, quantile(P_cum_zeroeffect, p_c, stick), 
                                title='transformation with flat pcum', xlabel='input', ylabel='quantile')
(81792, 324) (81792, 324)
No description has been provided for this image

The relative mean squared error is low:

In [33]:
print('Relative difference = ', np.sum( (p_c - quantile(P_cum_zeroeffect, p_c, stick))**2) / np.sum(p_c**2))
Relative difference =  3.82686497427e-06

COMP : using modulations

let's use this new rescaling function

In [34]:
l0_sparseness = shl.l0_sparseness
def comp(data, dico, P_cum, C=C, do_sym=do_sym, verbose=0):
    if verbose!=0: t0 = time.time()
    n_samples, n_dictionary = data.shape[0], dico.shape[0]
    sparse_code = np.zeros((n_samples, n_dictionary))
    corr = (data @ dico.T)
    Xcorr = (dico @ dico.T)
    nb_quant = P_cum.shape[1]
    stick = np.arange(n_dictionary)*nb_quant
    
    for i_sample in range(n_samples):
        c = corr[i_sample, :].copy()
        if verbose!=0: ind_list=list()
        for i_l0 in range(int(l0_sparseness)):
            if P_cum is None:
                q_i = rescaling(c, C=C, do_sym=do_sym)
            else:
                q_i = quantile(P_cum, rescaling(c, C=C, do_sym=do_sym), stick, do_fast=True)
            ind  = np.argmax(q_i)
            if verbose!=0: ind_list.append(ind)

            c_ind = c[ind] / Xcorr[ind, ind]
            sparse_code[i_sample, ind] += c_ind
            c -= c_ind * Xcorr[ind, :]

        if verbose!=0 and i_sample in range(2):
            q_i = quantile(P_cum, rescaling(c, C=C, do_sym=do_sym), stick)
            print(ind_list, [q_i[i] for i in ind_list], np.median(q_i), q_i.max(), [c[i] for i in ind_list], c.min(), c.max())
    if verbose!=0:
        duration = time.time()-t0
        print('coding duration : {0}s'.format(duration))
    return sparse_code

corr = (data_test @ dico_partial_learning.dictionary.T)
C = get_rescaling(corr, nb_quant=nb_quant, verbose=True)
sparse_code = comp(data_test, dico_partial_learning.dictionary, P_cum_zeroeffect, C=C, do_sym=do_sym, verbose=1)
At indices (ranks)  [0, 51860, 103720, 155580, 207441, 259301, 311161, 363022, 414882, 466742, 518602, 570463, 622323, 674183, 726044, 777904, 829764, 881624, 933485, 985345, 1037205, 1089066, 1140926, 1192786, 1244646, 1296507, 1348367, 1400227, 1452088, 1503948, 1555808, 1607668, 1659529, 1711389, 1763249, 1815110, 1866970, 1918830, 1970690, 2022551, 2074411, 2126271, 2178132, 2229992, 2281852, 2333712, 2385573, 2437433, 2489293, 2541154, 2593014, 2644874, 2696734, 2748595, 2800455, 2852315, 2904176, 2956036, 3007896, 3059756, 3111617, 3163477, 3215337, 3267198, 3319058, 3370918, 3422778, 3474639, 3526499, 3578359, 3630220, 3682080, 3733940, 3785801, 3837661, 3889521, 3941381, 3993242, 4045102, 4096962, 4148823, 4200683, 4252543, 4304403, 4356264, 4408124, 4459984, 4511845, 4563705, 4615565, 4667425, 4719286, 4771146, 4823006, 4874867, 4926727, 4978587, 5030447, 5082308, 5134168, 5186028, 5237889, 5289749, 5341609, 5393469, 5445330, 5497190, 5549050, 5600911, 5652771, 5704631, 5756491, 5808352, 5860212, 5912072, 5963933, 6015793, 6067653, 6119513, 6171374, 6223234, 6275094, 6326955, 6378815, 6430675, 6482535, 6534396, 6586256, 6638116, 6689977, 6741837, 6793697, 6845557, 6897418, 6949278, 7001138, 7052999, 7104859, 7156719, 7208579, 7260440, 7312300, 7364160, 7416021, 7467881, 7519741, 7571602, 7623462, 7675322, 7727182, 7779043, 7830903, 7882763, 7934624, 7986484, 8038344, 8090204, 8142065, 8193925, 8245785, 8297646, 8349506, 8401366, 8453226, 8505087, 8556947, 8608807, 8660668, 8712528, 8764388, 8816248, 8868109, 8919969, 8971829, 9023690, 9075550, 9127410, 9179270, 9231131, 9282991, 9334851, 9386712, 9438572, 9490432, 9542292, 9594153, 9646013, 9697873, 9749734, 9801594, 9853454, 9905314, 9957175, 10009035, 10060895, 10112756, 10164616, 10216476, 10268336, 10320197, 10372057, 10423917, 10475778, 10527638, 10579498, 10631358, 10683219, 10735079, 10786939, 10838800, 10890660, 10942520, 10994380, 11046241, 11098101, 11149961, 11201822, 11253682, 11305542, 11357403, 11409263, 11461123, 11512983, 11564844, 11616704, 11668564, 11720425, 11772285, 11824145, 11876005, 11927866, 11979726, 12031586, 12083447, 12135307, 12187167, 12239027, 12290888, 12342748, 12394608, 12446469, 12498329, 12550189, 12602049, 12653910, 12705770, 12757630, 12809491, 12861351, 12913211, 12965071, 13016932, 13068792, 13120652, 13172513, 13224373, 13276233, 13328093, 13379954, 13431814, 13483674, 13535535, 13587395, 13639255, 13691115, 13742976, 13794836, 13846696, 13898557, 13950417, 14002277, 14054137, 14105998, 14157858, 14209718, 14261579, 14313439, 14365299, 14417159, 14469020, 14520880, 14572740, 14624601, 14676461, 14728321, 14780181, 14832042, 14883902, 14935762, 14987623, 15039483, 15091343, 15143204, 15195064, 15246924, 15298784, 15350645, 15402505, 15454365, 15506226, 15558086, 15609946, 15661806, 15713667, 15765527, 15817387, 15869248, 15921108, 15972968, 16024828, 16076689, 16128549, 16180409, 16232270, 16284130, 16335990, 16387850, 16439711, 16491571, 16543431, 16595292, 16647152, 16699012, 16750872, 16802733, 16854593, 16906453, 16958314, 17010174, 17062034, 17113894, 17165755, 17217615, 17269475, 17321336, 17373196, 17425056, 17476916, 17528777, 17580637, 17632497, 17684358, 17736218, 17788078, 17839938, 17891799, 17943659, 17995519, 18047380, 18099240, 18151100, 18202960, 18254821, 18306681, 18358541, 18410402, 18462262, 18514122, 18565982, 18617843, 18669703, 18721563, 18773424, 18825284, 18877144, 18929004, 18980865, 19032725, 19084585, 19136446, 19188306, 19240166, 19292027, 19343887, 19395747, 19447607, 19499468, 19551328, 19603188, 19655049, 19706909, 19758769, 19810629, 19862490, 19914350, 19966210, 20018071, 20069931, 20121791, 20173651, 20225512, 20277372, 20329232, 20381093, 20432953, 20484813, 20536673, 20588534, 20640394, 20692254, 20744115, 20795975, 20847835, 20899695, 20951556, 21003416, 21055276, 21107137, 21158997, 21210857, 21262717, 21314578, 21366438, 21418298, 21470159, 21522019, 21573879, 21625739, 21677600, 21729460, 21781320, 21833181, 21885041, 21936901, 21988761, 22040622, 22092482, 22144342, 22196203, 22248063, 22299923, 22351783, 22403644, 22455504, 22507364, 22559225, 22611085, 22662945, 22714806, 22766666, 22818526, 22870386, 22922247, 22974107, 23025967, 23077828, 23129688, 23181548, 23233408, 23285269, 23337129, 23388989, 23440850, 23492710, 23544570, 23596430, 23648291, 23700151, 23752011, 23803872, 23855732, 23907592, 23959452, 24011313, 24063173, 24115033, 24166894, 24218754, 24270614, 24322474, 24374335, 24426195, 24478055, 24529916, 24581776, 24633636, 24685496, 24737357, 24789217, 24841077, 24892938, 24944798, 24996658, 25048518, 25100379, 25152239, 25204099, 25255960, 25307820, 25359680, 25411540, 25463401, 25515261, 25567121, 25618982, 25670842, 25722702, 25774562, 25826423, 25878283, 25930143, 25982004, 26033864, 26085724, 26137584, 26189445, 26241305, 26293165, 26345026, 26396886, 26448746, 26500607]
the coefficients are  [-36.2 -6.93 -5.81 -5.19 -4.77 -4.44 -4.17 -3.95 -3.76 -3.59 -3.45 -3.31
 -3.19 -3.08 -2.98 -2.89 -2.8 -2.72 -2.64 -2.57 -2.5 -2.43 -2.37 -2.31
 -2.26 -2.2 -2.15 -2.1 -2.05 -2.01 -1.97 -1.92 -1.88 -1.84 -1.81 -1.77
 -1.73 -1.7 -1.67 -1.63 -1.6 -1.57 -1.54 -1.51 -1.49 -1.46 -1.43 -1.41
 -1.38 -1.36 -1.33 -1.31 -1.29 -1.26 -1.24 -1.22 -1.2 -1.18 -1.16 -1.14
 -1.12 -1.1 -1.08 -1.07 -1.05 -1.03 -1.01 -0.997 -0.98 -0.964 -0.949 -0.933
 -0.918 -0.903 -0.888 -0.874 -0.86 -0.846 -0.832 -0.819 -0.806 -0.793 -0.78
 -0.768 -0.755 -0.743 -0.731 -0.72 -0.708 -0.697 -0.686 -0.675 -0.664
 -0.653 -0.643 -0.632 -0.622 -0.612 -0.602 -0.592 -0.583 -0.573 -0.564
 -0.555 -0.546 -0.537 -0.528 -0.519 -0.511 -0.502 -0.494 -0.486 -0.478
 -0.47 -0.462 -0.454 -0.447 -0.439 -0.432 -0.424 -0.417 -0.41 -0.403 -0.396
 -0.389 -0.382 -0.375 -0.369 -0.362 -0.356 -0.349 -0.343 -0.337 -0.331
 -0.325 -0.319 -0.313 -0.307 -0.301 -0.295 -0.29 -0.284 -0.279 -0.273
 -0.268 -0.263 -0.258 -0.252 -0.247 -0.242 -0.238 -0.233 -0.228 -0.223
 -0.218 -0.214 -0.209 -0.205 -0.2 -0.196 -0.191 -0.187 -0.183 -0.179 -0.175
 -0.17 -0.167 -0.163 -0.159 -0.155 -0.151 -0.147 -0.144 -0.14 -0.136 -0.133
 -0.129 -0.126 -0.123 -0.119 -0.116 -0.113 -0.11 -0.107 -0.103 -0.1 -0.0975
 -0.0946 -0.0918 -0.0891 -0.0864 -0.0837 -0.0811 -0.0785 -0.076 -0.0736
 -0.0712 -0.0689 -0.0667 -0.0645 -0.0623 -0.0603 -0.0583 -0.0564 -0.0546
 -0.0528 -0.0511 -0.0495 -0.0479 -0.0463 -0.0447 -0.0433 -0.0419 -0.0405
 -0.0391 -0.0378 -0.0366 -0.0353 -0.0341 -0.033 -0.0318 -0.0307 -0.0296
 -0.0286 -0.0276 -0.0265 -0.0256 -0.0246 -0.0236 -0.0226 -0.0217 -0.0208
 -0.0199 -0.019 -0.0181 -0.0173 -0.0164 -0.0156 -0.0148 -0.014 -0.0132
 -0.0124 -0.0116 -0.0108 -0.01 -0.00922 -0.00844 -0.00766 -0.00685 -0.00608
 -0.00532 -0.00456 -0.00381 -0.00307 -0.00233 -0.00158 -0.000822 -7.28e-05
 0.000709 0.00149 0.00227 0.00307 0.00386 0.00461 0.00537 0.00613 0.00693
 0.00772 0.00851 0.00929 0.0101 0.0109 0.0117 0.0126 0.0134 0.0142 0.0151
 0.016 0.0169 0.0178 0.0187 0.0196 0.0206 0.0215 0.0225 0.0235 0.0245
 0.0255 0.0265 0.0275 0.0286 0.0297 0.0308 0.0321 0.0333 0.0345 0.0358
 0.0371 0.0384 0.0398 0.0413 0.0428 0.0443 0.046 0.0477 0.0494 0.0511
 0.0529 0.0547 0.0567 0.0587 0.0607 0.0628 0.065 0.0671 0.0693 0.0717
 0.0741 0.0766 0.0792 0.0819 0.0846 0.0873 0.0901 0.093 0.0959 0.0989 0.102
 0.105 0.108 0.111 0.115 0.118 0.121 0.125 0.128 0.132 0.135 0.139 0.143
 0.146 0.15 0.154 0.158 0.162 0.166 0.17 0.174 0.178 0.183 0.187 0.191
 0.195 0.2 0.204 0.209 0.214 0.218 0.223 0.228 0.233 0.238 0.243 0.248
 0.253 0.258 0.263 0.269 0.274 0.28 0.285 0.291 0.296 0.302 0.308 0.314
 0.32 0.326 0.332 0.338 0.344 0.351 0.357 0.363 0.37 0.377 0.384 0.39 0.397
 0.405 0.412 0.419 0.426 0.434 0.441 0.449 0.457 0.464 0.472 0.48 0.489
 0.497 0.505 0.514 0.522 0.531 0.54 0.549 0.558 0.567 0.577 0.586 0.596
 0.606 0.615 0.626 0.636 0.646 0.657 0.667 0.678 0.689 0.701 0.712 0.723
 0.735 0.747 0.759 0.772 0.784 0.797 0.81 0.823 0.836 0.85 0.864 0.878
 0.893 0.907 0.922 0.937 0.952 0.968 0.984 1 1.02 1.03 1.05 1.07 1.09 1.11
 1.12 1.14 1.16 1.18 1.2 1.22 1.25 1.27 1.29 1.31 1.34 1.36 1.39 1.41 1.44
 1.46 1.49 1.52 1.55 1.58 1.61 1.64 1.67 1.7 1.74 1.77 1.81 1.85 1.89 1.93
 1.97 2.01 2.06 2.1 2.15 2.2 2.26 2.31 2.37 2.43 2.5 2.57 2.64 2.71 2.8
 2.88 2.98 3.08 3.19 3.31 3.44 3.59 3.75 3.94 4.16 4.42 4.74 5.17 5.8 6.91
 36.6]
[112, 33, 278, 97, 101, 256, 30, 177, 36, 58, 59, 278, 88, 22, 278, 317] [0.66653387130965935, 0.65052856048959673, 0.59830775238067113, 0.62302444001996959, 0.64458959926104442, 0.64085691898403307, 0.61400716999153337, 0.57881440167231402, 0.0, 0.55064929939865515, 0.62791281020614631, 0.59830775238067113, 0.53729111115839245, 0.66534963023112492, 0.59830775238067113, 0.0] 0.0 0.726428524852 [0.14622114532348629, 0.11726767443482161, 0.04961471253073943, 0.075987077275044848, 0.10740013301351933, 0.10148169743045644, 0.065259558267194429, 0.034756512475519329, -0.026039702037990109, 0.019488148020597243, 0.082498085429622031, 0.04961471253073943, 0.013389142092147963, 0.14397244131880477, 0.04961471253073943, 0.0] -0.999707454888 0.287689288868
[200, 33, 73, 266, 6, 308, 278, 157, 137, 162, 194, 30, 290, 274, 212, 278] [0.52657777111988913, 0.0, 0.52231210787367488, 0.0, 0.53100666664042351, 0.50860769300083175, 0.0, 0.51572938871061569, 0.50472799517621458, 0.51739239907578116, 0.50876069676420965, 0.51956589106510176, 0.51830653729577003, 0.52029190910548584, 0.50996317355797349, 0.0] 0.0 0.550232245702 [0.0089410744039219933, -0.013120664016748315, 0.0072309348418384935, -0.00013820264737959961, 0.010729581615201201, 0.0017983779666867076, 0.0, 0.0046274909864738086, 0.00024940207037818778, 0.0052722273582512689, 0.001859008757976976, 0.006117370455199618, 0.0056274915879661741, 0.0064112040165253299, 0.0023374950692405774, 0.0] -0.07863481717 0.0192883783941
coding duration : 324.98530888557434s
In [35]:
%%timeit
corr = (data_test @ dico_partial_learning.dictionary.T)
C = get_rescaling(corr, nb_quant=nb_quant, verbose=False)
sparse_code = comp(data_test, dico_partial_learning.dictionary, P_cum_zeroeffect, C=C, do_sym=do_sym, verbose=0)
3min 24s ± 29.3 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [36]:
%%timeit
corr = (data_test @ dico_partial_learning.dictionary.T)
C = get_rescaling(corr, nb_quant=nb_quant, verbose=False)
sparse_code = comp(data_test, dico_partial_learning.dictionary, P_cum_zeroeffect, C=C, do_sym=do_sym, verbose=0)
3min 1s ± 6.55 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

testing that COMP with fixed Pcum is equivalent to MP

In [37]:
print(dico_partial_learning.P_cum)
None
In [38]:
from shl_scripts.shl_learn import get_P_cum
from shl_scripts.shl_encode import get_rescaling, rescaling
corr = (data_test @ dico_partial_learning.dictionary.T)
C = get_rescaling(corr, nb_quant=nb_quant, verbose=True)
p_c = rescaling(corr, C=C, do_sym=do_sym)
In [39]:
from shl_scripts.shl_encode import get_rescaling, mp
l0_sparseness = 50

sparse_code_mp = mp(data_test, dico_partial_learning.dictionary, P_cum=None, C=5., do_sym=do_sym, l0_sparseness=l0_sparseness, verbose=1)
sparse_code_comp = mp(data_test, dico_partial_learning.dictionary, P_cum=P_cum_zeroeffect, C=5., do_sym=do_sym, l0_sparseness=l0_sparseness, verbose=1)

#print('Difference = ', np.sum((sparse_code_mp - sparse_code_comp)**2))
#print('Relative difference = ', np.mean((sparse_code_mp - sparse_code_comp)**2))
#print('Variance = ', np.mean((sparse_code_mp)**2))
print('Relative difference = ', np.sum((sparse_code_mp - sparse_code_comp)**2)/np.sum((sparse_code_mp)**2))

fig, ax = plot_proba_histogram(sparse_code_mp)
fig.show()
fig, ax = plot_proba_histogram(sparse_code_comp)
coding duration : 90.52337574958801
coding duration : 346.8867189884186
Relative difference =  0.0
No description has been provided for this image
No description has been provided for this image
In [40]:
n_samples, nb_filter = sparse_code_mp.shape

from shl_scripts.shl_encode import get_rescaling, mp

corr = (data_test @ dico_partial_learning.dictionary.T)
C_vec = get_rescaling(corr, nb_quant=nb_quant, verbose=True)

P_cum = np.vstack((P_cum_zeroeffect, C_vec))

print('Classical MP')
sparse_code_mp = mp(data_test, dico_partial_learning.dictionary, P_cum=None, C=np.inf, do_sym=do_sym, l0_sparseness=l0_sparseness, verbose=1)
print('Homeostatic MP')
sparse_code_comp = mp(data_test, dico_partial_learning.dictionary, P_cum=P_cum, C=0., do_sym=do_sym, l0_sparseness=l0_sparseness, verbose=1)

print('Relative difference = ', np.sum((sparse_code_mp - sparse_code_comp)**2)/np.sum((sparse_code_mp)**2))
Classical MP
coding duration : 89.88441371917725
Homeostatic MP
coding duration : 536.8180718421936
Relative difference =  0.122805037453
In [41]:
fig, ax = plot_proba_histogram(sparse_code_mp)
fig.show()
fig, ax = plot_proba_histogram(sparse_code_comp)
No description has been provided for this image
No description has been provided for this image
In [42]:
fig, ax = plot_scatter_MpVsComp(sparse_code_mp, sparse_code_comp)
(81792, 324) (81792, 324)
No description has been provided for this image

gradient descent

First with one rescaling function, and varying P_cum functions.

In [43]:
from shl_scripts.shl_learn import get_P_cum
from shl_scripts.shl_encode import get_rescaling, mp

corr = (data_test @ dico_partial_learning.dictionary.T)
C_vec = get_rescaling(corr, nb_quant=nb_quant, verbose=True)

P_cum = np.linspace(0, 1, nb_quant, endpoint=True)[np.newaxis, :] * np.ones((nb_filter, 1))
print('Shape of modulation function', P_cum.shape)
P_cum = np.vstack((P_cum, C_vec))
print('Shape of modulation function', P_cum.shape)

eta_homeo = .01

for i in range(400//DEBUG_DOWNSCALE):
    
    sparse_code = mp(data_test, dico_partial_learning.dictionary, P_cum=P_cum, C=0., do_sym=do_sym, verbose=0)
    P_cum_ = get_P_cum(sparse_code, nb_quant=nb_quant, C=C_vec)
    P_cum[:-1, :] = (1-eta_homeo) * P_cum[:-1, :] + eta_homeo  * P_cum_
    
    if i % (100//DEBUG_DOWNSCALE) == 0:
        print('Learning step', i)
        fig, ax = plot_proba_histogram(sparse_code)
        plt.show()
Shape of modulation function (324, 512)
Shape of modulation function (325, 512)
Learning step 0
No description has been provided for this image
Learning step 100
No description has been provided for this image
Learning step 200
No description has been provided for this image
Learning step 300
No description has been provided for this image

now by learning the rescaling function along with the P_cum functions.

In [44]:
from shl_scripts.shl_learn import update_P_cum
from shl_scripts.shl_encode import get_rescaling, mp

corr = (data_test @ dico_partial_learning.dictionary.T)
C_vec= C_vec_init = get_rescaling(corr, nb_quant=nb_quant, verbose=True)

P_cum = np.linspace(0, 1, nb_quant, endpoint=True)[np.newaxis, :] * np.ones((nb_filter, 1))
print('Shape of modulation function', P_cum.shape)
P_cum = np.vstack((P_cum, C_vec))
print('Shape of modulation function', P_cum.shape)

eta_homeo = .01

for i in range(400//DEBUG_DOWNSCALE):
    
    sparse_code = mp(data_test, dico_partial_learning.dictionary, P_cum=P_cum, C=0., do_sym=do_sym, verbose=0)
    corr = (data_test @ dico_partial_learning.dictionary.T)
    C_vec = get_rescaling(corr, nb_quant=nb_quant, do_sym=do_sym, verbose=verbose)
    
    P_cum[-1, :]= (1 - eta_homeo) * P_cum[-1, :] + eta_homeo * C_vec
    P_cum[:-1, :] = update_P_cum(P_cum=P_cum[:-1, :], 
                                 code=sparse_code, eta_homeo=eta_homeo, 
                 C=P_cum[-1, :], nb_quant=nb_quant, do_sym=do_sym, verbose=verbose)

    if i % (100//DEBUG_DOWNSCALE) == 0:
        print('Learning step', i)
        fig, ax = plot_proba_histogram(sparse_code)
        plt.show()
Shape of modulation function (324, 512)
Shape of modulation function (325, 512)
Learning step 0
No description has been provided for this image
Learning step 100
No description has been provided for this image
Learning step 200
No description has been provided for this image
Learning step 300
No description has been provided for this image

Note how the rescaling function evolved during learning

In [45]:
fig=plt.figure(figsize=(13, 8))
ax = plt.subplot(111)
x = np.linspace(0, 30, 100)
ax.plot(x, np.ones_like(x), '--')
ax.plot(x, rescaling(x, C=C_vec_init, do_sym=do_sym), label='init')
ax.plot(x, rescaling(x, C=.5, do_sym=do_sym), label='C=5')
ax.plot(x, rescaling(x, C=C_vec, do_sym=do_sym), label='auto')
ax.set_xlabel('raw coefficient')
ax.set_ylabel('rescaled coefficient')
ax.set_xscale('log')
ax.axis('tight')
ax.legend(loc='best');
No description has been provided for this image
In [46]:
from shl_scripts.shl_tools import plot_P_cum
fig, ax = plot_P_cum(P_cum[:-1, :], verbose=False);
No description has been provided for this image
In [47]:
fig, ax = plot_P_cum(P_cum[:-1, :], verbose=False)
ax.set_ylim(0.8, 1.01);ax.set_xlim(0., 1.);
No description has been provided for this image

Similar to the previous post, one observes that the distribution gets progressively more uniform which was our goal. This implementation makes the Matching Pursuit algorithm running faster such that we can integrate it to the SHL package. This was done in this commit.

finally...

Let's try the whole learning process with the adaptive rescaling function:

In [48]:
from shl_scripts.shl_experiments import SHL
matname = tag + '_autohomeo'
shl = SHL(DEBUG_DOWNSCALE=DEBUG_DOWNSCALE,
          datapath='/tmp/database',
          data_cache='/tmp/data_cache',
          verbose=2, 
          homeo_method='HAP', homeo_params=dict(alpha_homeo=0., C=0.), 
          do_sym=do_sym, nb_quant=nb_quant)
data = shl.get_data(matname=matname)
loading the data called : /tmp/data_cache/2017-12-01_Testing_COMPs_autohomeo_data
In [49]:
test_size = data.shape[0]//2
data_training = data[:test_size, :]
data_test = data[test_size:,:]   
if DEBUG:
    test_size = data.shape[0]//20
    data_training = data[:(data.shape[0]-test_size),:].copy()
    data_test = data[:test_size, :].copy()
In [50]:
dico_partial_learning = shl.learn_dico(data=data_training, matname=matname)
loading the dico called : /tmp/data_cache/2017-12-01_Testing_COMPs_autohomeo_dico.pkl

We start off by using a short learning with no homeostasis such that we end up with a unbalanced dictionary:

In [51]:
fig, ax = shl.show_dico(dico_partial_learning, data=data, title=matname)
fig.show()
fig, ax = shl.time_plot(dico_partial_learning, variable='prob_active');
fig.show()
No description has been provided for this image
No description has been provided for this image