This function calculates the  STA of LFP, extracts kernels and recontructs the LFP from kernels.
    
    Arguments
    ::
      transient : the time in milliseconds, after which the analysis should begin
                so as to avoid any starting transients
      X : id of presynaptic trigger population

def fig_kernel_lfp_EITN_II(savefolders, params, transient=200, T=[800., 1000.], X='L5E',
                   lags=[20, 20], channels=[0,3,7,11,13]):

    '''
    This function calculates the  STA of LFP, extracts kernels and recontructs the LFP from kernels.
    
    Arguments
    ::
      transient : the time in milliseconds, after which the analysis should begin
                so as to avoid any starting transients
      X : id of presynaptic trigger population
       
    '''

    # Electrode geometry
    zvec = np.r_[params.electrodeParams['z']]

    alphabet = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'



    ana_params.set_PLOS_2column_fig_style(ratio=0.5)
    # Start the figure
    fig = plt.figure()
    fig.subplots_adjust(left=0.06, right=0.95, bottom=0.08, top=0.9, hspace=0.23, wspace=0.55) 

    # create grid_spec
    gs = gridspec.GridSpec(len(channels), 7)


    ###########################################################################
    # spikegen "network" activity
    ############################################################################
    
    
    # path to simulation files
    savefolder = 'simulation_output_spikegen'
    params.savefolder = os.path.join(os.path.split(params.savefolder)[0],
                                         savefolder)
    params.figures_path = os.path.join(params.savefolder, 'figures')
    params.spike_output_path = os.path.join(params.savefolder,
                                                       'processed_nest_output')
    params.networkSimParams['spike_output_path'] = params.spike_output_path
    
    
    # Get the spikegen LFP:
    f = h5py.File(os.path.join('simulation_output_spikegen', 'LFPsum.h5'))
    srate = f['srate'].value
    tvec = np.arange(f['data'].shape[1]) * 1000. / srate
    
    # slice
    inds = (tvec < params.tstop) & (tvec >= transient)
    
    data_sg_raw = f['data'].value.astype(float)
    data_sg = data_sg_raw[:, inds]
    f.close()
    
    # kernel width
    kwidth = 20
    
    # create some dummy spike times
    activationtimes = np.array([x*100 for x in range(3,11)] + [200])
    networkSimSpikegen = CachedNetwork(**params.networkSimParams)

    x, y = networkSimSpikegen.get_xy([transient, params.tstop])

    
    ############################################################################
    ## Part A: spatiotemporal kernels, all presynaptic populations
    #############################################################################
    #
    #titles = ['TC',
    #          'L23E/I',
    #          'LFP kernels \n L4E/I',
    #          'L5E/I',
    #          'L6E/I',
    #          ]
    #
    #COUNTER = 0 
    #for i, X__ in enumerate(([['TC']]) + zip(params.X[1::2], params.X[2::2])):        
    #    ax = fig.add_subplot(gs[:len(channels), i])
    #    if i == 0:
    #        phlp.annotate_subplot(ax, ncols=7, nrows=4, letter=alphabet[0], linear_offset=0.02)
    #
    #    for j, X_ in enumerate(X__):
    #        # create spikegen histogram for population Y
    #        cinds = np.arange(activationtimes[np.arange(-1, 8)][COUNTER]-kwidth,
    #                          activationtimes[np.arange(-1, 8)][COUNTER]+kwidth+2)
    #        x0_sg = np.histogram(x[X_], bins=cinds)[0].astype(float)
    #        
    #        if X_ == ('TC'):
    #            color='r'
    #        else:
    #            color=('r', 'b')[j]
    #        
    #        
    #        # plot kernel as correlation of spikegen LFP signal with delta spike train
    #        xcorr, vlimround = plotting_correlation(params,
    #                                                x0_sg/x0_sg.sum()**2, 
    #                                                data_sg_raw[:, cinds[:-1]]*1E3,
    #                                                ax, normalize=False,
    #                                                lag=kwidth,
    #                                                color=color,
    #                                                scalebar=False)
    #        if i > 0:
    #            ax.set_yticklabels([])
    #        
    #        ## Create scale bar
    #        ax.plot([kwidth, kwidth],
    #            [-1500 + j*3*100, -1400 + j*3*100], lw=2, color=color,
    #            clip_on=False)
    #        ax.text(kwidth*1.08, -1450 + j*3*100, '%.1f $\mu$V' % vlimround,
    #                    rotation='vertical', va='center')
    #
    #        ax.set_xlim((-5, kwidth))
    #        ax.set_xticks([-20, 0, 20])
    #        ax.set_xticklabels([-20, 0, 20])
    #        
    #        COUNTER += 1
    #        
    #    ax.set_title(titles[i])
    

    ################################################
    # Iterate over savefolders
    ################################################

    for i, (savefolder, lag) in enumerate(zip(savefolders, lags)):
        
        # path to simulation files
        params.savefolder = os.path.join(os.path.split(params.savefolder)[0],
                                         savefolder)
        params.figures_path = os.path.join(params.savefolder, 'figures')
        params.spike_output_path = os.path.join(params.savefolder,
                                                'processed_nest_output')
        params.networkSimParams['spike_output_path'] = params.spike_output_path
        
        #load spike as database inside function to avoid buggy behaviour
        networkSim = CachedNetwork(**params.networkSimParams)
        
        
        
        # Get the Compound LFP: LFPsum : data[nchannels, timepoints ]
        f = h5py.File(os.path.join(params.savefolder, 'LFPsum.h5'))
        data_raw = f['data'].value
        srate = f['srate'].value
        tvec = np.arange(data_raw.shape[1]) * 1000. / srate
        # slice
        inds = (tvec < params.tstop) & (tvec >= transient)
        data = data_raw[:, inds]
        # subtract mean
        dataT = data.T - data.mean(axis=1)
        data = dataT.T
        f.close()
        
        # Get the spikegen LFP:
        f = h5py.File(os.path.join('simulation_output_spikegen', 'LFPsum.h5'))
        data_sg_raw = f['data'].value
        
        f.close()
        #
        #
        #
        #
        #########################################################################
        ## Part B: STA LFP
        #########################################################################
        #
        #titles = ['staLFP(%s)\n(spont.)' % X, 'staLFP(%s)\n(AC. mod.)' % X]
        #ax = fig.add_subplot(gs[:len(channels), 5 + i])
        #if i == 0:
        #    phlp.annotate_subplot(ax, ncols=15, nrows=4, letter=alphabet[i+1],
        #                          linear_offset=0.02)
        #  
        #collect the spikes x is the times, y is the id of the cell.
        x, y = networkSim.get_xy([0,params.tstop])
        #
        ## Get the spikes for the population of interest given as 'Y'
        bins = np.arange(0, params.tstop+2) + 0.5
        x0_raw = np.histogram(x[X], bins=bins)[0]
        x0 = x0_raw[inds].astype(float)
        #
        ## correlation between firing rate and LFP deviation
        ## from mean normalized by the number of spikes  
        #xcorr, vlimround = plotting_correlation(params,
        #                                        x0/x0.sum(), 
        #                                        data*1E3,
        #                                        ax, normalize=False, 
        #                                        #unit='%.3f mV',
        #                                        lag=lag,
        #                                        scalebar=False,
        #                                        color='k',
        #                                        title=titles[i],
        #                                        )
        #
        ## Create scale bar
        #ax.plot([lag, lag],
        #    [-1500, -1400], lw=2, color='k',
        #    clip_on=False)
        #ax.text(lag*1.08, -1450, '%.1f $\mu$V' % vlimround,
        #            rotation='vertical', va='center')
        #
        #
        #[Xind] = np.where(np.array(networkSim.X) == X)[0]        
        #        
        ## create spikegen histogram for population Y
        #x0_sg = np.zeros(x0.shape, dtype=float)
        #x0_sg[activationtimes[Xind]] += params.N_X[Xind]
        #
        #
        #ax.set_yticklabels([])
        #ax.set_xticks([-lag, 0, lag])
        #ax.set_xticklabels([-lag, 0, lag])

        
        ###########################################################################
        # Part C, F: LFP and reconstructed LFP
        ############################################################################

        # create grid_spec
        gsb = gridspec.GridSpec(len(channels), 8)

      
        ax = fig.add_subplot(gsb[:, (i*4):(i*4+2)])
        phlp.annotate_subplot(ax, ncols=8/2., nrows=4, letter=alphabet[i*3+2],
                              linear_offset=0.02)
          
        # extract kernels, force negative lags to be zero
        kernels = np.zeros((len(params.N_X), 16, kwidth*2))
        for j in range(len(params.X)):
            kernels[j, :, kwidth:] = data_sg_raw[:, (j+2)*100:kwidth+(j+2)*100]/params.N_X[j]
    
        LFP_reconst_raw = np.zeros(data_raw.shape)
    
        for j, pop in enumerate(params.X):
            x0_raw = np.histogram(x[pop], bins=bins)[0].astype(float)
            for ch in range(kernels.shape[1]):
                LFP_reconst_raw[ch] += np.convolve(x0_raw, kernels[j, ch],
                                                   'same')
        
        # slice
        LFP_reconst = LFP_reconst_raw[:, inds]
        # subtract mean
        LFP_reconstT = LFP_reconst.T - LFP_reconst.mean(axis=1)
        LFP_reconst = LFP_reconstT.T
        vlimround = plot_signal_sum(ax, params,
                                    fname=os.path.join(params.savefolder,
                                                           'LFPsum.h5'),
                    unit='mV', scalebar=True,
                        T=T, ylim=[-1550, 50],
                        color='k', label='$real$',
                        rasterized=False)
    
        plot_signal_sum(ax, params, fname=LFP_reconst_raw,
                    unit='mV', scaling_factor= 1., scalebar=False,
                        vlimround=vlimround,
                        T=T, ylim=[-1550, 50],
                        color='r', label='$reconstr$',
                        rasterized=False)
        ax.set_title('LFP & population \n rate predictor')
        if i > 0:
            ax.set_yticklabels([])


        
        ###########################################################################
        # Part D,G: Correlation coefficient
        ############################################################################
        
        ax = fig.add_subplot(gsb[:, i*4+2:i*4+3])
        phlp.remove_axis_junk(ax)
        phlp.annotate_subplot(ax, ncols=8./1, nrows=4, letter=alphabet[i*3+3],
                              linear_offset=0.02)
            
        cc = np.zeros(len(zvec))
        for ch in np.arange(len(zvec)):
            cc[ch] = np.corrcoef(data[ch], LFP_reconst[ch])[1, 0]
        
        ax.barh(zvec, cc, height=90, align='center', color='1', linewidth=0.5)
        ax.set_ylim([-1550, 50])
        ax.set_yticklabels([])
        ax.set_yticks(zvec)
        ax.set_xlim([0.0, 1.])
        ax.set_xticks([0.0, 0.5, 1])
        ax.yaxis.tick_left()
        ax.set_xlabel('$cc$ (-)', labelpad=0.1)
        ax.set_title('corr. \n coef.')
    
        print 'correlation coefficients:'
        print cc


        ###########################################################################
        # Part E,H: Power spectra
        ############################################################################
        

        #compute PSDs ratio between ground truth and estimate
        freqs, PSD_data = calc_signal_power(params, fname=data,
                                            transient=transient, Df=None, mlab=True,
                                            NFFT=256, noverlap=128,
                                            window=plt.mlab.window_hanning)
        freqs, PSD_LFP_reconst = calc_signal_power(params, fname=LFP_reconst,
                                            transient=transient, Df=None, mlab=True,
                                            NFFT=256, noverlap=128,
                                            window=plt.mlab.window_hanning)

        zv = np.r_[params.electrodeParams['z']]
        zv = np.r_[zv, zv[-1] + np.diff(zv)[-1]]
        inds = freqs >= 1  # frequencies greater than 1 Hz  

        for j, ch in enumerate(channels):
        
            ax = fig.add_subplot(gsb[j, (i*4+3):(i*4+4)])
            if j == 0:
                phlp.annotate_subplot(ax, ncols=8./1, nrows=4.5*len(channels),
                                      letter=alphabet[i*3+4], linear_offset=0.02)
                ax.set_title('PSD')
                

            phlp.remove_axis_junk(ax)
            ax.loglog(freqs[inds], PSD_data[ch, inds], 'k', label='LFP', clip_on=True)
            ax.loglog(freqs[inds], PSD_LFP_reconst[ch, inds], 'r', label='predictor', clip_on=True)
            ax.set_xlim([4E0,4E2])
            ax.set_ylim([1E-8, 1E-4])
            ax.tick_params(axis='y', which='major', pad=0)
            ax.set_yticks([1E-8,1E-6,1E-4])
            ax.yaxis.set_minor_locator(plt.NullLocator())
            ax.text(0.8, 0.9, 'ch. %i' % (ch+1),
                    horizontalalignment='left',
                    verticalalignment='center',
                    fontsize=6, 
                    transform=ax.transAxes)
            
            if j == 0:
                ax.set_ylabel('(mV$^2$/Hz)', labelpad=0.)
            if j > 0:
                ax.set_yticklabels([])
            if j == len(channels)-1:
                ax.set_xlabel(r'$f$ (Hz)', labelpad=0.)
            else:
                ax.set_xticklabels([])



    return fig, PSD_LFP_reconst, PSD_data