def get_arguments():
    info_txt = 'Usage: python plot_activity_as_colormap.py FOLDER_NAME CELL_TYPE [PATTERN_NR]'
    try:
        folder = sys.argv[1]
        params_fn = os.path.abspath(folder) + '/Parameters/simulation_parameters.json'
        param_tool = simulation_parameters.parameter_storage(params_fn=params_fn)
    except:
        print info_txt
        print 'Taking the parameters currently in simulation_parameters.py'
        param_tool = simulation_parameters.parameter_storage()

    params = param_tool.params
    print 'debug n_cells', params['n_cells']

    try:
        cell_type = sys.argv[2]
    except:
        print 'Missing cell_type argument'
        print info_txt
        exit(1)

    try:
        pn = int(sys.argv[3])
    except:
        print info_txt
        print 'Plotting pattern 0'
        pn = 0

    return params, cell_type, pn

def return_plot(cell_gids=[], subplot_code=111, fig=None, input_fn_base=None, motion_params=None):
    network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
    params = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
    if input_fn_base == None:
        input_fn_base = params['input_rate_fn_base']

    fn = params['tuning_prop_means_fn']
    if motion_params == None:
        mp = params['motion_params']
    else:
        mp = motion_params
    tp = np.loadtxt(fn)

    if len(cell_gids) == 0:
        n_cells =15 
        cell_gids = np.random.randint(0, params['n_exc'], n_cells)
    n_cells = len(cell_gids)
    ms = 5# markersize for scatterplots
    bg_color = 'w'
    pylab.rcParams['lines.markeredgewidth'] = 0

#    input_sum = np.zeros(n_cells)
#    for i, gid in enumerate(cell_gids):
#        input_fn = input_fn_base + str(gid) + '.dat'
#        rate = np.loadtxt(input_fn)
#        input_sum[i] = rate.sum()
#    input_max = input_sum.max()

    if fig == None:
        fig = pylab.figure(facecolor=bg_color)
    ax = fig.add_subplot(subplot_code)
    colors = ['b', 'g']
    for i, gid in enumerate(cell_gids):
        x, y, u, v = tp[gid, :]
#        print 'tp[%d]:' % (gid), tp[gid, :]

#        h = 240.
#        l = 1. - 0.5 * input_sum[i] / input_max
#        s = 1. # saturation
#        assert (0 <= h and h < 360)
#        assert (0 <= l and l <= 1)
#        assert (0 <= s and s <= 1)
#        (r, g, b) = utils.convert_hsl_to_rgb(h, s, l)
#        ax.plot(x, y, 'o', c=(r,g,b), markersize=ms)

        ax.plot(x, y, 'o', c=colors[i%len(colors)], markersize=ms)
        ax.quiver(x, y, u, v, angles='xy', scale_units='xy', scale=1, width=0.02)#, headwidth=6)


    # plot stimulus
    stim_color = 'y'
    ax.quiver(mp[0], mp[1], mp[2], mp[3], angles='xy', scale_units='xy', scale=1, color=stim_color, headwidth=6, width=0.02)
    ax.annotate('Stimulus', (mp[0]+.5*mp[2], mp[1]+0.1), fontsize=12, color=stim_color)

    ax.set_xlim((0, 1))
    ax.set_ylim((0, 1))
    #output_fn_fig = 'delme_test.png'
    #print "Saving figure: ", output_fn_fig
    #pylab.savefig(output_fn_fig)#, facecolor=bg_color)
    return ax

    def __init__(self, params=None, **kwargs):

        print 'BasicPlotter'
        if params == None:
            self.network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
            self.params = self.network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        else:
            self.params = params

        self.subfig_cnt = 1
        self.n_fig_x = kwargs.get('n_fig_x', 2)
        self.n_fig_y = kwargs.get('n_fig_y', 2)
        self.tuning_prop = np.loadtxt(self.params['tuning_prop_means_fn'])
        assert (self.tuning_prop[:, 0].size == self.params['n_exc']), 'Number of cells does not match in %s and simulation_parameters!\n Wrong tuning_prop file?' % self.params['tuning_prop_means_fn']

        # figure details
        fig_width_pt = 800.0  # Get this from LaTeX using \showthe\columnwidth
        inches_per_pt = 1.0/72.27               # Convert pt to inch
        golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
        fig_width = fig_width_pt*inches_per_pt  # width in inches
        fig_height = fig_width*golden_mean      # height in inches
#        fig_size =  [fig_width,fig_height]
        fig_size =  [fig_height,fig_width]
        params = {#'backend': 'png',
                  'titel.fontsize': 16,
                  'axes.labelsize' : 12,
                  'text.fontsize': 12,
                  'figure.figsize': fig_size}
        pylab.rcParams.update(params)

def plot_conductances_vs_time(params=None, comm=None, data_fn=None, inh_spikes = None):

    t_start = time.time()
    if params== None:
        network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
#        P = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        params = network_params.params

    sim_cnt = 0
    if data_fn == None:
        data_fn = params['exc_spiketimes_fn_merged'] + '%d.ras' % (sim_cnt)

#    if inh_spikes == None:
#        inh_spikes = params['inh_spiketimes_fn_merged'] + '%d.ras' % (sim_cnt)

    plotter = P.PlotConductances(params, comm, data_fn)
    plotter.load_spiketimes()

    if plotter.no_spikes:
        return
    output_fn_base = '%s%s_wsigmaX_%.2f_wsigmaV%.2f_wthresh%.1e' % (params['grouped_actitivty_fig_fn_base'], params['connectivity_code'], \
            params['w_sigma_x'], params['w_sigma_v'], params['w_thresh_connection'])

    # fig 1
    # neuronal level
    plotter.create_fig()  # create an empty figure
    plotter.plot_rasterplot('exc', 1)               # 1 
    plotter.plot_rasterplot('inh', 2)               # 2 
    plotter.plot_group_spikes_vs_time(3)            # 3
    output_fn = output_fn_base + '_0.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn)

#    output_fn = '%sgoodcell_connections_%s_wsigmaX_%.2f_wsigmaV%.2f_wthresh%.1e.png' % (params['figures_folder'], params['connectivity_code'], \
#            params['w_sigma_x'], params['w_sigma_v'], params['w_thresh_connection'])
#    plotter.create_fig()  # create an empty figure
#    plotter.plot_good_cell_connections(1) # subplot 1 + 2
#    print 'Saving figure to:', output_fn
#    pylab.savefig(output_fn)

    # fig 2
    plotter.create_fig()  # create an empty figure
    plotter.plot_input_cond(1) # subplot 1 + 2
    plotter.plot_conductances()
    output_fn = output_fn_base + '_1.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn)


    plotter.create_fig()  # create an empty figure
    plotter.plot_input_cond()

    t_stop = time.time()
    t_run = t_stop - t_start
    print "PlotConductance duration: %d sec or %.1f min for %d cells (%d exc, %d inh)" % (t_run, (t_run)/60., \
            params['n_cells'], params['n_exc'], params['n_inh'])

    def __init__(self, params=None, comm=None, data_fn=None, sim_cnt=0):

        if params == None:
            self.network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
            self.params = self.network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        else:
            self.params = params
        self.no_spikes = False
        self.comm = comm

        self.n_fig_x = 2
        self.n_fig_y = 2
        self.tuning_params = np.loadtxt(self.params['tuning_prop_means_fn'])

        # define parameters
        self.n_cells = self.params['n_exc']
        self.time_binsize = 20 # [ms]
        self.n_bins = int((self.params['t_sim'] / self.time_binsize) )
        self.time_bins = [self.time_binsize * i for i in xrange(self.n_bins)]
        self.t_axis = np.arange(0, self.n_bins * self.time_binsize, self.time_binsize)

        self.n_good = self.params['n_exc'] * .10 # fraction of 'interesting' cells
        print 'Number of cells with \'good\' tuning_properties = ', self.n_good

        # create data structures
        self.nspikes = np.zeros(self.n_cells)                                   # summed activity
        self.nspikes_binned = np.zeros((self.n_cells, self.n_bins))             # binned activity over time
        self.spiketrains = [[] for i in xrange(self.n_cells)]

        self.tuning_prop = np.loadtxt(self.params['tuning_prop_means_fn'])
        # sort the cells by their proximity to the stimulus into 'good_gids' and the 'rest'
        # cell in 'good_gids' should have the highest response to the stimulus
        print 'utils.sort_gids_by_distance_to_stimulus'
        all_gids, all_distances = utils.sort_gids_by_distance_to_stimulus(self.tuning_prop, self.params['motion_params'], self.params) 
        self.good_gids, self.good_distances = all_gids[0:self.n_good], all_distances[0:self.n_good]
        print 'Saving gids to record to', self.params['gids_to_record_fn']
        np.savetxt(self.params['gids_to_record_fn'], np.array(self.good_gids), fmt='%d')
        self.rest_gids = range(self.n_cells)
        for gid in self.good_gids:
            self.rest_gids.remove(gid)

        fig_width_pt = 800.0  # Get this from LaTeX using \showthe\columnwidth
        inches_per_pt = 1.0/72.27               # Convert pt to inch
        golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
        fig_width = fig_width_pt*inches_per_pt  # width in inches
        fig_height = fig_width*golden_mean      # height in inches
        fig_size =  [fig_width,fig_height]
        params = {#'backend': 'png',
#                  'axes.labelsize': 10,
#                  'text.fontsize': 10,
                  'legend.fontsize': 10,
#                  'xtick.labelsize': 8,
#                  'ytick.labelsize': 8,
#                  'text.usetex': True,
                  'figure.figsize': fig_size}
        pylab.rcParams.update(params)

def sort_cells_by_distance_to_stimulus(n_cells, verbose=True):
    import simulation_parameters
    network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
    params = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
    tp = np.loadtxt(params['tuning_prop_means_fn'])
    mp = params['motion_params']
    indices, distances = sort_gids_by_distance_to_stimulus(tp , mp, params) # cells in indices should have the highest response to the stimulus
    print 'Motion parameters', mp
    print 'GID\tdist_to_stim\tx\ty\tu\tv\t\t'
    if verbose:
        for i in xrange(n_cells):
            gid = indices[i]
            print gid, '\t', distances[i], tp[gid, :]
    return indices, distances

    def __init__(self, param_fn=None, spiketimes_fn=None):
        """
        params : dictionary or NeuroTools.parameters ParameterSet
        """
        if params == None:
            self.network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
            self.params = self.network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        else:
            self.params = params
        self.no_spikes = False

        self.n_fig_x = 2
        self.n_fig_y = 2

        self.n_cells = self.params['n_exc']
        self.nspikes = np.zeros(self.n_cells)                                   # summed activity
        self.spiketrains = [[] for i in xrange(self.n_cells)]
        self.load_spiketimes(data_fn)

def plot_input_colormap(params=None, data_fn=None, inh_spikes = None):

    if params== None:
        network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
#        P = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        params = network_params.params

    if data_fn == None:
        if params.has_key('merged_input_spiketrains_fn'):
            output_fn = params['merged_input_spiketrains_fn']
        else:
            params['merged_input_spiketrains_fn'] = "%sinput_spiketrain_merged.dat" % (params['input_folder'])
            data_fn = params['merged_input_spiketrains_fn']
    if not os.path.exists(data_fn):
        merge_input_spiketrains(params)
    

    plotter = P.PlotPrediction(params, data_fn)
    pylab.rcParams['axes.labelsize'] = 14
    pylab.rcParams['axes.titlesize'] = 16
    if plotter.no_spikes:
        return

    plotter.compute_v_estimates()
#    plotter.compute_position_estimates()
#    plotter.compute_theta_estimates()

    # fig 1
    # neuronal level
    output_fn_base = params['figures_folder'] + 'input_colormap.png'

    plotter.create_fig()  # create an empty figure
    pylab.subplots_adjust(left=0.07, bottom=0.07, right=0.97, top=0.93, wspace=0.3, hspace=.2)
    plotter.n_fig_x = 2
    plotter.n_fig_y = 2
#    plotter.plot_rasterplot('exc', 1)               # 1 
#    plotter.plot_rasterplot('inh', 2)               # 2 
    plotter.plot_vx_grid_vs_time(1)              # 3 
    plotter.plot_vy_grid_vs_time(2)              # 4 
    plotter.plot_x_grid_vs_time(3, ylabel='x-position of stimulus')
    plotter.plot_y_grid_vs_time(4, ylabel='y-position of stimulus')
    output_fn = output_fn_base + '_0.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn)

    def __init__(self, params=None):
        if params == None:
            network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
    #        P = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
            self.params = network_params.params
        else:
            self.params = params

        self.n_exc = self.params['n_exc']
        self.output = []
        self.g_in_histograms = []
        self.output_fig = self.params['conductances_fig_fn_base']
        self.n_good = self.params['n_exc'] * .05 # fraction of 'good' (well-tuned) cells
        print 'Number of \'good\' (well-tuned) cells:', self.n_good

        self.no_spikes = False

        self.load_nspikes()

        self.conn_dict = {}
        for conn_type in self.params['conn_types']:
            print 'Calling utils.get_conn_dict(..., %s)' % conn_fn
            conn_fn = self.params['conn_list_%s_fn' % conn_type]
            self.conn_dict[conn_type] = utils.get_conn_dict(self.params, conn_fn)

        fig_width_pt = 800.0  # Get this from LaTeX using \showthe\columnwidth
        inches_per_pt = 1.0/72.27               # Convert pt to inch
        golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
        fig_width = fig_width_pt*inches_per_pt  # width in inches
        fig_height = fig_width*golden_mean      # height in inches
        fig_size =  [fig_width,fig_height]
        params = {#'backend': 'png',
                  'axes.labelsize': 12,
#                  'text.fontsize': 14,
#                  'legend.fontsize': 10,
#                  'xtick.labelsize': 8,
#                  'ytick.labelsize': 8,
#                  'text.usetex': True,
                  'figure.figsize': fig_size}
        pylab.rcParams.update(params)
        pylab.subplots_adjust(bottom=0.30)

    def __init__(self, argv):
        if len(argv) > 1:
            if argv[1].isdigit():
                gid = int(argv[1])
            else:
                param_fn = argv[1]
                if os.path.isdir(param_fn):
                    param_fn += '/Parameters/simulation_parameters.json'
                print '\nLoading parameters from %s\n' % (param_fn)
                f = file(param_fn, 'r')
                params = json.load(f)
        else:
            print '\nLoading the default paremters...\n'
            import simulation_parameters
            ps = simulation_parameters.parameter_storage()
            params = ps.params

        self.params = params
        self.conn_list_loaded = [False, False, False, False]
        self.conn_mat_loaded = [False, False, False, False]
        self.conn_lists = {}

    def __init__(self, params=None, comm=None):

        if params == None:
            network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
            params = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
            print 'Merging connlists ...'
            os.system('python merge_connlists.py')
        else:
            self.params = params
            print 'Assuming that \n\tpython merge_connlists.py \nhas been called before in the directory %s' % params['folder_name']

        self.comm = comm
        if comm != None:
            self.pc_id, self.n_proc = comm.rank, comm.size
        self.conn_lists = {}
        self.n_fig_x = 1
        self.n_fig_y = 1

        # cell markers
        self.markersize_cell = 10
        self.markersize_min = 3
        self.markersize_max = 12
        self.shaft_width = 0.005
        self.conn_type_dict = {'e' : 'excitatory', 'i' : 'inhibitory'}

    def __init__(self, param_fn=None, spiketimes_fn=None):
        """
        params : dictionary or NeuroTools.parameters ParameterSet
        """

        print "debug", type(param_fn)
        if param_fn == None:
            print "Loading default parameters stored in simulation_parameters.py"
            self.network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
            self.params = self.network_params.load_params()                       # params stores cell numbers, etc as a dictionary
            self.params_fn = self.params['params_fn']
        else:
            if type(param_fn) != type(""): raise TypeError("File name expected for param_fn")
            self.params_fn = param_fn
            self.params = ntp.ParameterSet(param_fn)
            self.params
    
        self.spiketimes_fn = spiketimes_fn
        print os.path.abspath(self.params_fn)
        print os.path.abspath(self.spiketimes_fn)

        fig_width_pt = 800.0  # Get this from LaTeX using \showthe\columnwidth
        inches_per_pt = 1.0/72.27               # Convert pt to inch
        golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
        fig_width = fig_width_pt*inches_per_pt  # width in inches
        fig_height = fig_width*golden_mean      # height in inches
        fig_size =  [fig_width,fig_height]
        params = {#'backend': 'png',
#                  'axes.labelsize': 10,
#                  'text.fontsize': 10,
#                  'legend.fontsize': 10,
#                  'xtick.labelsize': 8,
#                  'ytick.labelsize': 8,
#                  'text.usetex': True,
                  'figure.figsize': fig_size}
        pylab.rcParams.update(params)

    # EPTH -> OB connections are not affected by the pattern
    print "Creating connections: orn -> mit"
    ConnectionClass = CreateObConnections.CreateObConnections(params)
    ConnectionClass.connect_orn_mit()
    ConnectionClass.connect_orn_pg()
    ConnectionClass.connect_pg_mit_serial()
    ConnectionClass.connect_pg_mit_reciprocal()
    ConnectionClass.connect_mt_gran_local()
    ConnectionClass.connect_mt_gran_global()


if __name__ == '__main__':
    print info_txt
    # ------------ I N I T -----------------------------
    # The simulation_parameters module defines a class for simulation parameter storage
    param_tool = simulation_parameters.parameter_storage()
    # params is the dictionary with all parameters
    params = param_tool.params

    if not (params['test_pattern_rivalry']):
        print '\n\n\tThis scipt is not intended to be used with the train_pattern_rivalry flag!\nWill now quit'
        exit(1)

    print 'New folder:', params['folder_name']
    ok = raw_input('\nContinue to create this folder structure? Parameters therein will be overwritten\n\ty / Y / blank = OK; anything else --> exit\n')
    if not ((ok == '') or (ok.capitalize() == 'Y')):
        print 'quit'
        exit(1)

    OrnParamClass = CreateOrnParameters.CreateOrnParameters(params) # patterns for ORN activation must be recreated to add noise

def plot_prediction(params=None, data_fn=None, inh_spikes = None):

    if params== None:
        network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
#        P = network_params.load_params()                       # params stores cell numbers, etc as a dictionary
        params = network_params.params

    if data_fn == None:
        data_fn = params['exc_spiketimes_fn_merged'] + '.ras'

#    if inh_spikes == None:
#        inh_spikes = params['inh_spiketimes_fn_merged'] + '.ras'

#    params['t_sim'] = 1200

    plotter = P.PlotPrediction(params, data_fn)
    pylab.rcParams['axes.labelsize'] = 14
    pylab.rcParams['axes.titlesize'] = 16
    if plotter.no_spikes:
        return

    plotter.compute_v_estimates()
#    plotter.compute_position_estimates() --> happening in compute_v_estimates
    plotter.compute_theta_estimates()
    plotter.compute_orientation_estimates()

    # fig 1
    # neuronal level
    output_fn_base = '%s%s_wsigmaX_%.2f_wsigmaV%.2f_delayScale%d_connRadius%.2f_wee%.2f' % (params['prediction_fig_fn_base'], params['connectivity_code'], \
            params['w_sigma_x'], params['w_sigma_v'], params['delay_scale'], params['connectivity_radius'], params['w_tgt_in_per_cell_ee'])


    plotter.create_fig()  # create an empty figure
    pylab.subplots_adjust(left=0.07, bottom=0.07, right=0.97, top=0.93, wspace=0.3, hspace=.2)
    plotter.n_fig_x = 2
    plotter.n_fig_y = 3
    plotter.plot_rasterplot('exc', 1)               # 1 
    plotter.plot_rasterplot('inh', 2)               # 2 
    plotter.plot_vx_grid_vs_time(3)              # 3 
    plotter.plot_vy_grid_vs_time(4)              # 4 
    plotter.plot_x_grid_vs_time(5)
    plotter.plot_y_grid_vs_time(6)
    output_fn = output_fn_base + '_0.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn, dpi=200)
#    output_fn = output_fn_base + '_0.pdf'
#    print 'Saving figure to:', output_fn
#    pylab.savefig(output_fn, dpi=200)

    # orientation figure
    plotter.create_fig()  # create an empty figure
    plotter.n_fig_x = 3
    plotter.n_fig_y = 2 # to have the ~same figure sizes as for the others
    plotter.plot_orientation_grid_vs_time(1)
    plotter.plot_orientation_estimates(2)
    plotter.plot_orientation_diff(3)
    output_fn = output_fn_base + '_orientation.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn, dpi=200)

#    output_fn = output_fn_base + '_0.eps'
#    print 'Saving figure to:', output_fn
#    pylab.savefig(output_fn, dpi=200)



    # poplation level, short time-scale
    plotter.n_fig_x = 3
    plotter.n_fig_y = 2
    plotter.create_fig()
    pylab.rcParams['legend.fontsize'] = 12
    pylab.subplots_adjust(left=0.07, bottom=0.07, right=0.97, top=0.93, wspace=0.3, hspace=.3)
    plotter.plot_vx_estimates(1)
    plotter.plot_vy_estimates(2)
    plotter.plot_vdiff(3)
    plotter.plot_x_estimates(4)
    plotter.plot_y_estimates(5) 
    plotter.plot_xdiff(6)
    output_fn = output_fn_base + '_1.png'
    print 'Saving figure to:', output_fn
    pylab.savefig(output_fn, dpi=200)
#    output_fn = output_fn_base + '_1.pdf'
#    print 'Saving figure to:', output_fn
#    pylab.savefig(output_fn, dpi=200)
#    output_fn = output_fn_base + '_1.eps'
#    print 'Saving figure to:', output_fn
#    pylab.savefig(output_fn, dpi=200)



#    plotter.plot_theta_estimates(5)


    # fig 3
    # population level, long time-scale

#    plotter.n_fig_x = 1
#    plotter.n_fig_y = 4
#    pylab.rcParams['legend.fontsize'] = 10
#    pylab.subplots_adjust(hspace=0.5)
#.........这里部分代码省略.........

    Usage:
        python plot_response_curve.py [FOLDER] [CELLTYPE] 
        or
        python plot_response_curve.py [FOLDER] [CELLTYPE] [PATTERN_NUMBER]
    """
    assert (len(sys.argv) > 2), 'ERROR: folder and cell_type not given\n' + info_txt
    folder = sys.argv[1]
    cell_type = sys.argv[2]
    try:
        pn = int(sys.argv[3])
    except:
        print 'WARNING: Using the default pattern number 0'
        pn = 0

    params_fn = os.path.abspath(folder) + '/Parameters/simulation_parameters.json'
    param_tool = simulation_parameters.parameter_storage(params_fn=params_fn)
    params = param_tool.params

    if cell_type == 'all':
#        cell_types = params['cell_types']
        cell_types = ['mit', 'pg', 'gran']
    else:
        cell_types = [cell_type]

    for cell_type in cell_types:
        print 'Plotting raster for:', cell_type
        plot_raster_for_celltype(params, cell_type, title='%s spikes pattern %d' % (cell_type.upper(), pn))
#        plot_raster_for_celltype(params, cell_type, title='%s spikes pattern %d' % (cell_type.upper(), pn))

    pylab.show()

import sys


if len(sys.argv) > 1:
    param_fn = sys.argv[1]
    if os.path.isdir(param_fn):
        param_fn += '/Parameters/simulation_parameters.json'
    import json
    f = file(param_fn, 'r')
    print 'Loading parameters from', param_fn
    params = json.load(f)

else:
    print '\n NOT successfull\nLoading the parameters currently in simulation_parameters.py\n'
    import simulation_parameters
    network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
    params = network_params.load_params()                       # params stores cell numbers, etc as a dictionary

# E -> E 
tmp_fn = 'delme_tmp_%d' % (np.random.randint(0, 1e8))
cat_cmd = 'cat %s* > %s' % (params['conn_list_ee_fn_base'], tmp_fn)
sort_cmd = 'sort -gk 1 -gk 2 %s > %s' % (tmp_fn, params['merged_conn_list_ee'])
rm_cmd = 'rm %s' % (tmp_fn)

print cat_cmd
os.system(cat_cmd)
print sort_cmd
os.system(sort_cmd)
print rm_cmd
os.system(rm_cmd)

def return_plot(iteration=None, fig=None):
    network_params = simulation_parameters.parameter_storage()  # network_params class containing the simulation parameters
    params = network_params.load_params()                       # params stores cell numbers, etc as a dictionary

    fn = params['tuning_prop_means_fn']
    tp = np.loadtxt(fn)
    if iteration == None:
        mp = params['motion_params']
    else:
        motion_params_fn = params['parameters_folder'] + 'input_params.txt'
        all_mp = np.loadtxt(motion_params_fn)
        mp = all_mp[iteration, :]

    n_cells = tp[:, 0].size
    #n_cells = 10
    ms = 10 # markersize for scatterplots
    bg_color = 'w'
    pylab.rcParams['lines.markeredgewidth'] = 0

    input_sum = np.zeros(n_cells)
    if iteration == None:
        input_fn_base = params['input_rate_fn_base'] 
        fn_ending = '.npy'
    else:
        input_fn_base = params['folder_name'] + 'TrainingInput_%d/abstract_input_' % iteration
        fn_ending = '.dat'

    for i in xrange(n_cells):
        input_fn = input_fn_base + str(i) + fn_ending
        if iteration == None:
            rate = np.load(input_fn)
        else:
            rate = np.loadtxt(input_fn)
        input_sum[i] = rate.sum()

    input_max = input_sum.max()
    print 'input_max', input_max
    idx = input_sum.argsort()
    n_mac = int(round(params['n_exc'] * 0.05))
    mac = idx[-n_mac:]
    print 'motion stimulus', mp
    print 'most activated cells:'
    for i in mac:
        dist = tp[i, :] - mp
#        print i, tp[i, :], np.sqrt(np.dot(dist, dist)), input_sum[i], input_max
    if fig == None:
        fig = pylab.figure(facecolor=bg_color)
    ax = fig.add_subplot(211)
#    h = 240.
#    s = 1. # saturation

    o_min = 0.
    o_max = input_max
    norm = matplotlib.mpl.colors.Normalize(vmin=o_min, vmax=o_max)
    m = matplotlib.cm.ScalarMappable(norm=norm, cmap=cm.binary)#PuBu)#jet)
    m.set_array(np.arange(o_min, o_max, (o_max-o_min)/1000.))

    for i in idx:
        if i in mac:
            print 'debug', i, input_sum[i], input_max, input_sum[i] / input_max
        c = m.to_rgba(input_sum[i])
#        l = 1. - 0.5 * input_sum[i] / input_max
#        assert (0 <= h and h < 360)
#        assert (0 <= l and l <= 1)
#        assert (0 <= s and s <= 1)
#        (r, g, b) = utils.convert_hsl_to_rgb(h, s, l)
        x, y, u, v = tp[i, :]
        ax.plot(x, y, 'o', c=c, markersize=ms)
#        ax.plot(x, y, 'o', c=(r,g,b), markersize=ms)
        
#        if l < .7:
        if i in mac:
            ax.annotate('%d' % i, (x + np.random.rand() * 0.02, y + np.random.rand() * 0.02), fontsize=10)

    stim_color = 'k'
    ax.quiver(mp[0], mp[1], mp[2], mp[3], angles='xy', scale_units='xy', scale=1, color=stim_color, headwidth=4, pivot='tail')
    ax.annotate('Stimulus', (mp[0]+.5*mp[2], mp[1]+0.1), fontsize=12, color=stim_color)

    fig.colorbar(m)

    ax.set_xlim((-.05, 1.05))
    ax.set_ylim((-.05, 1.05))
    if iteration == None:
        iteration = 0
    ax.set_title('Abstract activtation iteration %d' % iteration)

    ax = fig.add_subplot(212)
    ax.bar(range(params['n_exc']), input_sum, width=1)

    output_fn = params['figures_folder'] + 'abstract_activation_%d.png' % (iteration)
    print "Saving figure: ", output_fn
    pylab.savefig(output_fn)#, facecolor=bg_color)
    return ax

times = []
times.append(time.time())

full_system = False
if len(sys.argv) > 1:
    if sys.argv[1] == 'full':
        full_system = True

tau_dict = {'tau_zi' : 50.,    'tau_zj' : 5., 
            'tau_ei' : 50.,   'tau_ej' : 50., 'tau_eij' : 50.,
            'tau_pi' : 500.,  'tau_pj' : 500., 'tau_pij' : 500.,
            }

PS = simulation_parameters.parameter_storage()
params = PS.params
PS.create_folders()
PS.write_parameters_to_file()

n_cells = params['n_exc']
my_units = utils.distribute_n(n_cells, n_proc, pc_id)

mp = params['motion_params']

# P R E P A R E     T U N I N G    P R O P E R T I E S
tuning_prop = utils.set_tuning_prop(params, mode='hexgrid', v_max=params['v_max'])
np.savetxt(params['tuning_prop_means_fn'],tuning_prop)
#exit(1)

# load

import pylab
import numpy as np
import simulation_parameters as sp
import matplotlib
from matplotlib import cm

import time

PS = sp.parameter_storage()
params = PS.load_params()
tp = np.loadtxt(params['tuning_prop_means_fn'])
input_params = np.loadtxt(params['parameters_folder'] + 'input_params.txt')

#network_activity_fn = 'AndersWij/activity_test.dat'
#output_folder_fig  = 'AndersWij/test_activity/'

#network_activity_fn = 'Abstract/ANNActivity/ann_activity_40iterations_no_rec.dat'
network_activity_fn = 'Abstract/ANNActivity/ann_activity_40iterations.dat'
output_folder_fig  = params['figures_folder']

print 'Loading ', network_activity_fn
network_activity = np.loadtxt(network_activity_fn)
network_activity = np.exp(network_activity)

scale = 3
#n_time_steps = d[:, 0].size
n_time_steps_per_iteration = 300
n_cells = params['n_exc']

n_iteration = 40