ObjSim/python/objsimpy/som_results.py

# -*- coding: utf-8 -*-
from . import eval_sim_data
import os
from .sim_data import SimInfo, read_config
from objsimpy.stim.stim_sequence import read_input_file, InputTrace
from . import spike_file
import numpy as np
import numpy.testing as npt
import pylab as plt
#import pymeas.plot3d
from mpl_toolkits.axes_grid import AxesGrid
from . import colormaps
from .config import SIMDATA_BASE_DIR, MOVIE_DIR

def main():
    for i in range(26):
        simname = 'taurec_9__2'
        data_dir = os.path.join(SIMDATA_BASE_DIR, 'csim', 'simresults', simname)
        evaluate_single_sim(data_dir=data_dir,
                        test_nr=i, 
                        show=False,
                        reevaluate=False)

def evaluate_single_sim(data_dir=os.path.join(SIMDATA_BASE_DIR, 'csim', 'som02', 'asymgauss__1'), 
                        stim_dir=MOVIE_DIR,
                        test_nr=30,
                        reevaluate=False,
                        evaluate_spikes=False,
                        show=False,
                        fig_format='.png',                        
                        ):

    sim_info = SimInfo(os.path.join(data_dir, "Sim.SimInfo"))
    layer_name = "RepresentationLayer"
    input_name = "MovieScanInput0"
    layer_info = [l for l in sim_info.get_layers() if l.name==layer_name][0]
    
    plt.close("all")
    
    if evaluate_spikes:
        spike_file_name = layer_name + "spikes.datTEST"+str(test_nr)+".bin"      
        input_file_name = input_name + ".inp.datTEST"+str(test_nr)  
        spikes = spike_file.read(os.path.join(data_dir, spike_file_name))
        input = InputTrace(os.path.join(data_dir, input_file_name))
        
        print(len(spikes[0]))
        #max_time = max(spikes[0]) + 1
        #n_neurons = max(spikes[1]) + 1
        n_neurons = layer_info.LayerDimensions['N']
        #spikes = create_dummy_spikes(n_neurons=n_neurons, max_time_ms=max_time)
        #return
        response_strength_pickle_file = os.path.join(data_dir, spike_file_name + ".response.pickle")
        response_strength = eval_sim_data.calc_response_strength2_PICKLED(response_strength_pickle_file, 
                                                                          spikes, input, 430,540, n_neurons=n_neurons, reevaluate=reevaluate)
        response_selectivity = calc_response_selectivity(response_strength)
        
        print("response_strength: ")
        print(response_strength.shape)
        plt.plot(response_strength[:,0])
    
        stim_onsets = input.get_stim_onset()
        stim_numbers = [onset[1] for onset in stim_onsets]
    
        if False:
            fig = plt.figure()
            plt.hist(stim_numbers, bins=list(range(-2,410)))
        
        nx = sim_info.get_inputs()[0].NXPara
        ny = sim_info.get_inputs()[0].NYPara
        pickle_file = os.path.join(data_dir, spike_file_name +".pref_xy_sel.pickle")
        preferred_stimuli_x, preferred_stimuli_y, selectivity = calc_preferred_stimuli_from_response_PICKLED(pickle_file, response_strength, nx, ny, reevaluate=reevaluate)
        fig = show_preferred_stimulus_maps(preferred_stimuli_x, preferred_stimuli_y, selectivity, nx, ny)
    
    
    weight_file = os.path.join(data_dir, "conForward0weights.dat" + str(test_nr))
    stim_file = os.path.join(stim_dir, "gauss20x20.idlmov")
    #fig = plt.figure()
    stim_nx = 20
    stim_ny = 20
    pickle_file = os.path.join(data_dir, weight_file + ".pref_stim_sel_inc_map" + ".pickle")
    fname_stimcorr_pickle = os.path.join(data_dir, weight_file + ".stimcorr" + ".pickle")
    pref_x, pref_y, selectivity, inc_sum = eval_sim_data.calc_pref_stim_maps_from_weights_PICKLED(
                 pickle_file,
                 weight_file,
                 stim_file,
                 stim_nx,
                 stim_ny,
                 reevaluate=reevaluate,
                 fname_stimcorr_pickle=fname_stimcorr_pickle,
                 stimcorr_reeval=reevaluate)

    fig = show_preferred_stimulus_maps(pref_x, pref_y, selectivity, inc_sum, stim_nx, stim_ny)
    
    # patch_sizes_x, patch_sizes_y = eval_sim_data.calc_patch_sizes(weight_file, stim_file, stim_nx, stim_ny)

    if False:
        show_neuron_vs_stimuli(response_strength)
        for x in [10]:
            for y in [10]:
                pass
                show_stimuli_vs_neuron(response_strength, x_0=x, y_0=y)
    
    if False:
        show_total_psth(spikes, input)
        total_stim_response = response_strength.sum(axis=0)
        plt.matshow(total_stim_response.reshape(20,20), cmap=plt.cm.gray)
    
    if False:
        fig = plt.figure()
        response_to_all_stimuli = response_strength.sum(axis=1)
        plt.imshow(response_to_all_stimuli.reshape(100,100), origin='lower', interpolation='nearest', cmap=plt.cm.gray)
    
    if show:
        plt.show()

    save_prefmap(pref_x, stim_nx, 'preferred stimuli x', os.path.join(data_dir, 'prefmap_x_' + str(test_nr) + fig_format))
    save_prefmap(pref_y, stim_ny, 'preferred stimuli y', os.path.join(data_dir, 'prefmap_y_' + str(test_nr) + fig_format))
    
def create_dummy_spikes(n_neurons=10, max_time_ms=1000, spike_rate_herz=10):
    n_spikes = n_neurons * max_time_ms / (1000*spike_rate_herz)
    spike_times = round(np.linspace(0, max_time_ms-1, n_spikes))
    

def calc_response_selectivity(response_strength):
    total_mean_response = response_strength.mean()
    no_response = response_strength.max(axis=1) < total_mean_response
    print("no response n=", no_response.sum())

def show_preferred_stimulus_maps(preferred_stimuli_x, 
                                 preferred_stimuli_y, 
                                 selectivity=None, 
                                 inc_sum=None, 
                                 nx=20, 
                                 ny=20,
                                 figtext=None):

    fig = plt.figure(figsize=(10,11))

    if figtext is not None:
        plt.figtext(0.5, 0.965, figtext, ha='center', color='black', weight='bold', size='large')
    
    show_preferred_stimuli_x_y(preferred_stimuli_x, preferred_stimuli_y, nx, ny, fig, 421, 422)
    if selectivity is not None:
        show_selectivity_plot(selectivity, fig, 423)
    if inc_sum is not None:
        show_incomming_weights_sum_plot(inc_sum, fig, 424)
    show_preferred_stimuli_histograms(preferred_stimuli_x, preferred_stimuli_y, nx, ny, fig, 413, 414)
    return fig

def save_prefmap(pref, n_stim, title, fname):
    plt.figure()
    plot_prefmap(pref, n_stim, title)
    plt.savefig(fname)
    
def plot_prefmap(pref, n_stim, title):
    circular_cm = colormaps.cmap_color_circle()
    norm = plt.matplotlib.colors.Normalize(0, n_stim)
    im = plt.imshow(pref, origin='lower', interpolation='nearest', cmap=circular_cm)
    im.set_norm(norm)
    plt.title(title)
    plt.colorbar()

def print_y_hist(preferred_stimuli_y):
    y_hist = {}
    for y in preferred_stimuli_y.flat:
        if not y in y_hist:
            y_hist[y] = 1
        else:
            y_hist[y] += 1
    for y in y_hist:
        print("y = " + str(y) + " h=" + str(y_hist[y]))
    

def show_preferred_stimuli_x_y(preferred_stimuli_x, preferred_stimuli_y, nx, ny, fig, pos_x, pos_y):
    fig.canvas.set_window_title("preferred stimuli") 
    
    fig.add_subplot(pos_x)
    plot_prefmap(preferred_stimuli_x, nx, "preferred stimuli x")

    fig.add_subplot(pos_y)
    plot_prefmap(preferred_stimuli_y, ny, "preferred stimuli y")

def show_preferred_stimuli_histograms(preferred_stimuli_x, preferred_stimuli_y, nx, ny, fig, pos_x, pos_y):    
    sub = fig.add_subplot(pos_x)
    plt.hist(preferred_stimuli_x.flat, bins=np.arange(-2,20)+0.5, facecolor='green', alpha=0.75)
    plt.title("x histogram")
    plt.xlim([-1,nx])

    sub = fig.add_subplot(pos_y)
    plt.hist(preferred_stimuli_y.flat, bins=np.arange(-2,20)+0.5, facecolor='blue', alpha=0.75)
    plt.title("y histogram")
    plt.xlim([-1,ny])


def show_selectivity_plot(selectivity, fig, pos):
    sub = fig.add_subplot(pos)
    im_selectivity = plt.imshow(selectivity, origin='lower', interpolation='nearest', cmap=plt.cm.gray)
    plt.title("selectivity")
    norm_0_1 = plt.matplotlib.colors.Normalize(0, 1)
    im_selectivity.set_norm(norm_0_1)
    plt.colorbar()
    
def show_incomming_weights_sum_plot(inc_sum, fig, pos):
    sub = fig.add_subplot(pos)
    im_inc_sum = plt.imshow(inc_sum, origin='lower', interpolation='nearest', cmap=plt.cm.gray)
    plt.title("incomming weight sum")
    plt.colorbar()
    

def show_neuron_vs_stimuli(response_strength, n_x=5, n_y=5, x_0=0, y_0=0):
    columns = list(range(n_x))
    rows = list(range(n_y))
    fig = plt.figure(figsize=(n_x, n_y))
    grid = AxesGrid(fig, 111, nrows_ncols = (n_x, n_y), axes_pad=0.1)
    title = "neuron_vs_stimuli x0=" + str(x_0) + " y0=" + str(y_0)
    fig.canvas.set_window_title(title)        
    for x in columns:
        for y in rows:
            neuron_nr = (x_0 + x) + 100*(y_0 + y)
            sgl_response_map = response_strength[neuron_nr,:].reshape(20,20) 
            grid[x+n_x*y].matshow(sgl_response_map, cmap=plt.cm.gray)

def show_stimuli_vs_neuron(response_strength, n_x=5, n_y=5, x_0=0, y_0=0):
    columns = list(range(n_x))
    rows = list(range(n_y))
    show_neuron_vs_stimuli(response_strength)   
    fig = plt.figure()
    grid = AxesGrid(fig, 111, nrows_ncols = (n_x, n_y), axes_pad=0.1)
    title = "x0="+str(x_0)+", y0="+str(y_0)
    fig.canvas.set_window_title(title)               
    for x in columns:
        for y in rows:
            stim_nr = x_0 + x + 20 * (y_0 + y)
            sgl_stim_response = response_strength[:,stim_nr].reshape(100,100) 
            grid[x+n_x*y].matshow(sgl_stim_response, cmap=plt.cm.gray)                

def show_total_psth(spikes, input):
    total_psth = eval_sim_data.calc_total_psth(spikes, input, 1, 1000)
    fig = plt.figure()
    plt.plot(total_psth, 'o-')

def show_pref_stim_maps_from_weights(weight_file, stim_file, fig, pos1, pos2):
    stim_nx = 20
    stim_ny = 20
    pref_x, pref_y, selectivity = eval_sim_data.calc_pref_stim_maps_from_weights(weight_file, stim_file, stim_nx, stim_ny)
    show_preferred_stimuli_x_y(pref_x, pref_y, stim_nx, stim_ny, fig, pos1, pos2)




if __name__=="__main__":
    main()