text: from __future__ import division
  
  import numpy as np
  
  from pycog import tasktools
  import matplotlib.pyplot as plt         # Alfred
  from matplotlib import cm as cm         # Alfred
  import seaborn as sb
  import shutil
  import os
  import cPickle as pickle
  import random
  
  import scipy.linalg as la
  from scipy.signal import hilbert, chirp
  from sklearn.decomposition import PCA
  
  from scipy.io import savemat
  
  
  #-----------------------------------------------------------------------------------------
  # Network structure
  #-----------------------------------------------------------------------------------------
  
  Nin = 1
  N = 100
  Nout1 = 8
  Nout2 = 8
  # Nout = Nout1 + Nout2
  Nout = max(Nout1, Nout2)
  
  # E/I
  ei, EXC, INH = tasktools.generate_ei(N)
  
  # Time scaling
  ts = 1.0 # default is 1.0 by Alfred
  
  # Time constant
  tau = 50#*np.ones(N)        # modified by Alfred
  #tau[:80] = 50        # added by Alfred
  
  # Time step duration
  dt = 10.0/ts # default is 10 by Alfred
  
  # hidden_activation = 'softplus'	# added by Alfred
  
  lambda1_out = 2
  learning_rate = 1e-1
  
  # Setting the fontsize of Python plots by Alfred
  # plt.rcParams.update({'font.size': 14})  # for phase plots
  plt.rcParams.update({'font.size': 18})  # for wij plots
  
  #-----------------------------------------------------------------------------------------
  # Noise
  #-----------------------------------------------------------------------------------------
  
  var_rec = 0.01**2  # commented by Alfred
  
  # The following lines were added by Alfred for biological correctness
  Cout = np.zeros((Nout,N))
  for i in range(Nout):
      for j in range(N):
          if(i==j-0*Nout):
              Cout[i,j] = 1
  
  # Crec = np.ones((N,N))
  # for i in range(int(0.8*N)):
  #    for j in range(int(0.8*N)):
  #        if((random.uniform(0,1)>0.5) or i==j):
  #            Crec[i,j] = 0
  
  def rectify(x):
      return x*(x > 0)   # commented by Alfred for extracting x
  
  
  def generate_trial(rng, dt, params):
      T = 10000/ts    # changed to 3000 by Alfred
      TT = 800
      signal_time = int(rng.uniform(100, 110)/ts)#+signal_shift # (default=int(rng.uniform(100, 110)) or 0 for chirp inputs)
      delay = TT/ts#T-(100/ts) # changed from 800 to 750 by Alfred (default=900)
      width1 = (10.0/Nout1)*20/ts #+ np.random.uniform(-2,2,Nout) # (default=20)"" # noise added for some jitter # Alfred
      width2 = (10.0/Nout2)*20/ts #+ np.random.uniform(-2,2,Nout) # (default=20)"" # noise added for some jitter # Alfred
      magnitude = 3.0  #(default=3)
      slope = 0.5
      add_T = 3000/ts # added by Alfred (default=0)
  
      epochs = {}
      epochs['T'] = T
      t, e  = tasktools.get_epochs_idx(dt, epochs) # Time, task epochs in discrete time
      trial = {'t': t, 'epochs': epochs}           # Trial
  
      trial['info'] = {}
  
      signal_time /= dt
      delay /= dt
      width1 /= dt    # added newly by Alfred for two seqences case
      width2 /= dt    # added newly by Alfred for two seqences case
      add_T /= dt     # added newly by Alfred for 2nd delayed input
  
      LT = len(t)
      X = np.zeros((LT, Nin))
      inp = np.zeros((LT, Nin))
      rnd_freq = int(rng.uniform(40, 45)/ts) + 0# random frequency (Alfred) # (default=int(rng.uniform(40, 45)))
      rnd_bin = 0#rng.binomial(1,0.5)   # generate a random binomial term to decide increasing/descreasing sequence - added by Alfred
  
      for tt in range(int(LT*TT/T + signal_time)):
          if tt >= signal_time: #and tt<=1.1*signal_time:
              # inp[tt][0] = (1-slope)*(LT*TT/T-tt+signal_time)/6 + ((((1-slope)*(LT*TT/T-tt+signal_time)/6)*100)%100)/80	# Linear sawtooth input
              inp[tt][0] = 7#32*np.exp((LT*TT/T-10*tt+signal_time)/20) #(1-slope)*(LT*TT/T-tt+signal_time)/6 + ((((1-slope)*(LT*TT/T-tt+signal_time)/6)*100)%100)/80	# Linear sawtooth input
  
  
      if(rnd_bin==0): # decreasing case (default)
          for tt in range(int(LT*TT/T + signal_time)):
              if tt >= signal_time:
                  X[tt][0] = inp[tt][0]	# Linear sawtooth input
  
      else:   # increasing case
          for tt in range(int(LT*TT/T + signal_time)):
              if tt >= signal_time:
                  X[tt][0] = inp[int(np.floor(LT*TT/T + 2*int(signal_time) - tt))][0]	# Linear sawtooth input
  
  # For chirp as 2nd trigger
  #     for tt1 in range(int(LT*TT/T)):
  #         tt = int(tt1+add_T)
  #         if tt >= add_T:
  # #            X[tt][0] += (1.0+np.sin((tt - add_T)*(tt - add_T)/125))*chirp_mag	# Chirp input
  #             # X[tt][0] = inp[tt1][0]*1.7#np.sin(tt1*tt1/70)*chirp_mag + (1-slope)*(LT*500/T-tt1)/6	# Chirp input with slop downward
  #             X[tt][0] = inp[int(np.floor(LT*TT/T + 2*int(signal_time) - tt1))][0]*2#np.sin(tt1*tt1/70)*chirp_mag + (1-slope)*(LT*500/T-tt1)/6    # Chirp input with slop downward
  # #                X[tt][0] = np.sin((LT-tt1)*(LT-tt1)/70)*chirp_mag	+ tt1*slope/1# Chirp input with slope upward
  # #                X[tt][0] = (np.sin((LT - tt + add_T)*(LT - tt + add_T)/70))*chirp_mag	# Chirp input
  
  
      trial['inputs'] = X
  
      #---------------------------------------------------------------------------------
      # Target output
      #---------------------------------------------------------------------------------
  
      if params.get('target_output', False):
          Y = np.zeros((len(t), Nout)) # Output matrix
          M = np.zeros((len(t), Nout)) # Mask matrix
          out = np.zeros((len(t), Nout)) # Temp Output matrix
  
          for i in range(Nout1):
              for tt in range(int(LT*TT/T + signal_time)):
                  if tt >= signal_time:
                      out[tt][i] = np.exp( -(tt - (signal_time + delay / (Nout1+1) * (i+1)))**2 / (2 * width1**2)) * magnitude  # forward output: 1 to Nout1  for decay-rate-based outputs
  
          if(rnd_bin==0): # decreasing case (default)
              for i in range(Nout1):
                  for tt in range(int(LT*TT/T + signal_time)):
                      if tt >= signal_time:
                          Y[tt][i] = out[tt][i]#np.exp( -(tt - (signal_time + delay / (Nout1+1) * (i+1)))**2 / (2 * width1**2)) * magnitude  # forward output: 1 to Nout1  for decay-rate-based outputs
                          M[tt][i] = 1
          else:   # increasing case
              for j in range(Nout2):
                  for tt in range(int(LT*TT/T + signal_time)):
                      if tt >= signal_time:
                          # Y[tt][j] = out[int(np.floor(LT*TT/T + 2*int(signal_time) - tt))][j]#np.exp( -(tt - (signal_time + delay / (Nout2+1) * (Nout2-i)))**2 / (2 * width2**2)) * magnitude      # reverse output: Nout2 to 1
                          # M[tt][j] = 1
                          Y[tt][Nout-j-1] = out[tt][j]
                          M[tt][Nout-j-1] = 1
  
          trial['outputs'] = Y
  
  
      # heat_map = sb.heatmap(Y.T)
      # plt.title('Heat map of target sequential activation of neurons')
      # plt.ylabel('Target output neural nodes')
      # plt.xlabel('Time')
      # plt.show()
  
  
      # plt.plot(t/(20*tau), X[:,0])
      # legend = ['Input']
      # for j in range(Nout):
      #     plt.plot(t/(20*tau), Y[:,j])
      #     legend.append('Output {}'.format(j+1))
      # plt.title('Target sequential activation of neurons')
      # plt.ylabel('Target output neural nodes')
      # plt.xlabel('Time')
      # plt.legend(legend)
      # plt.show()
  
      return trial
  
  
  
  
  min_error = 0.1    # commented by Alfred # (default=0.1)
  n_gradient = 40  # uncommented by Alfred # (default=20)
  mode         = 'batch' # (default mode)
  
  
  if __name__ == '__main__':
      from pycog          import RNN
      from pycog.figtools import Figure
  
      rng = np.random.RandomState(1234)   # Added by Alfred
  
      # savefile = 'examples/work/data/cluster/two_sequences_sawtooth/two_sequences_sawtoothNout8_lr_01_seed60_inpswap.pkl'
      savefile = 'examples/work/data/cluster/two_sequences_sawtooth/two_sequences_sawtoothNout8_lr_01_seed30_2.pkl' # used for dc input/outputs
      # savefile = 'examples/work/data/cluster/two_sequences_sawtooth2/two_sequences_sawtooth2Nout8_lr_01_seed80_v2.pkl' # used for 1to8&9to16 cases in paper
      # savefile = 'examples/work/data/cluster/two_sequences_sawtooth/two_sequences_sawtoothNout8_lr_01_seed80_2.pkl'
      # savefile = 'examples/work/data/two_sequences_sawtooth2/two_sequences_sawtooth2Nout_1to8f_9to16r_v2.pkl'
  
      dt = 10.0/ts # added by Alfred
      factorwrec = 1.0#0.7 + np.arange(14)/20.0 # added by Alfred
      factorwout = 1.0 # added by Alfred
      factor = 1.0 # added by Alfred
      rnn  = RNN(savefile, {'dt': dt, 'var_rec': 0.01**2})  # added by Alfred
      # rnn  = RNN(savefile, {'dt': dt, 'var_rec': 0.0})  # added by Alfred
      trial_args = {}
  
      info1 = rnn.run(inputs=(generate_trial, trial_args), seed=200)
      Z0 = rnn.z
  
      node_drop_errors = np.zeros([1,N])
      node_drop_sums = np.zeros([1,N])
  
      rnn_zs = np.zeros([N,Nout,len(rnn.z[0])])
  
      r1 = random.randint(0,8)
      r2 = random.randint(0,8)
      while(r1!=r2):
          r2 = random.randint(0,8)
      temp1 = rnn.Wrec[r1,r2]
      temp2 = rnn.Wrec[r2,r1]
  
      for i in range(1):
          rnn  = RNN(savefile, {'dt': dt, 'var_rec': 0.01**2})  # added by Alfred
          # rnn  = RNN(savefile, {'dt': dt, 'var_rec': 0.0})  # added by Alfred
          trial_args = {}
          # wrec_old = rnn.Wrec
          # rnn.Wrec = factorwrec*rnn.Wrec #- 1.0*eyeN.astype(np.float32)
          # rnn.Wrec[:int(0.8*N),:int(0.8*N)] = rnn.Wrec[:int(0.8*N),:int(0.8*N)]*0.75#*1.01005#*factor#*factorwrec  # scaling  E to E weights
          # rnn.Wrec[int(0.8*N):,int(0.8*N):] = rnn.Wrec[int(0.8*N):,int(0.8*N):]*2.0#*factor#*factorwrec  # scaling  I to I weights
          # rnn.Wrec[:int(0.8*N),int(0.8*N):] = rnn.Wrec[:int(0.8*N),int(0.8*N):]*1.2#*factor#*factorwrec  # scaling  I to E weights
          # rnn.Wrec[int(0.8*N):,:int(0.8*N)] = rnn.Wrec[int(0.8*N):,:int(0.8*N)]*3.0#*factor#*factorwrec  # scaling  E to I weights
          for ii in range(80):
              for jj in range(80):
                  # rnn.Wrec[ii,jj] = 0
                  temp = rng.uniform(0,1)
                  if(temp>1.0):
                      rnn.Wrec[ii,jj] = 0
          # rnn.Wrec[r2,r1] = -1
      # rnn.Wrec[r1,r2] = temp1
      # rnn.Wrec[r2,r1] = temp2
  
  
      info1 = rnn.run(inputs=(generate_trial, trial_args), seed=200)
      rnn_Wrec_r = np.dot(rnn.Wrec,rnn.r)
      for j in range(Nout):
          # rnn_zs[i,j,:] = rnn.z[j]/np.max(rnn.z[j])
          rnn_zs[i,j,:] = rnn.r[j]/1
  
  
  
  
      # print "z is: ",rnn.z
      # ix = 0
      # # # np.savetxt("x1.csv", rnn_zs[0,:,:], delimiter=",")
      # # print "x is: ",rnn_zs[0,:,:]
      # # xdic = {"u{}".format(itrial): rnn.u,"x{}".format(itrial): rnn_zs[0,:,:]}
      xdic = {"u": rnn.u,"x": rnn_zs[0,:,:]}
      savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/x_dc7.mat", xdic)
      # xdic = {"Wrec": rnn.Wrec}
      # savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/two_seq2/Wrec_1to8_9to16.mat", xdic)
      # savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/two_seq2/Wrec_dc.mat", xdic)
      # # zdic = {"u": rnn.u,"z": rnn_zs[0,:,:]}
      # # savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/z_delta.mat", zdic)
  
  
      # inp = np.linspace(-7,7,15)
      # L = len(inp)
  
      # out = np.zeros([8,L])
      # out_mean = inp*0#np.zeros([1,L])
      # out_std = inp*0#np.zeros([1,L])
      # data = []
  
      # for i in range(L//2+1):
      #     temp = []
      #     # for j in range(8-abs(7-i)):
      #     print "i is: ",i
      #     for j in range(i+1):
      #         # print "p is: ",j,"; q is: ",j+7-i
      #         temp.append(rnn.Wrec[j,j+7-i])
      #     print "temp is ",temp
      #     data.append(temp)
      #     out_mean[i] = np.mean(temp)
      #     out_std[i] = np.std(temp)
  
      # for i in range(L//2):
      #     temp = []
      #     # for j in range(8-abs(7-i)):
      #     print "i is: ",i
      #     for j in range(abs(7-i)):
      #         # print "p is: ",j+i+1,"; q is: ",j
      #         temp.append(rnn.Wrec[j+i+1,j])
      #     print "temp is ",temp
      #     data.append(temp)
      #     out_mean[L//2+1+i] = np.mean(temp)
      #     out_std[L//2+1+i] = np.std(temp)
      # # for i in range(8):
      # #     for j in range(L):
      # #         out[i][j] = rnn.Wrec[8-j,]
  
      # # plt.errorbar(inp, out_mean, yerr=np.std(out,axis=0), fmt='.k',ecolor='lightgray',elinewidth=3);
      # plt.errorbar(inp, out_mean, yerr=out_std, fmt='.k',ecolor='lightgray',elinewidth=3)
      # plt.plot(inp, out_mean, '--')
      # plt.xlabel("(i-j) of $W^{ij}_{rec}$'s indices")
      # plt.ylabel("$W_{ij}$")
      # plt.show()
  
  
      # # np.savetxt("x1.csv", rnn_zs[0,:,:], delimiter=",")
      # print "x is: ",rnn_zs[0,:,:]
      # xdic = {"i_minus_j":inp,"wij_mean": out_mean,"wij_std": out_std, "wij_data":data}
      # savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/two_seq2/wij_stat_twoseq.mat", xdic)
      # zdic = {"u": rnn.u,"z": rnn_zs[0,:,:]}
      # savemat("exp4/Best_model_experiments/3tobase_experiments/23a_best_model/z_delta.mat", zdic)
  
      # print rnn.Wrec
      # heat_map = sb.heatmap(rnn.Wrec)
      # plt.title('Heat map of $W_{rec}$ weights matrix')
      # plt.ylabel('Rows')
      # plt.xlabel('Columns')
      # plt.show()
  
      # for i in range(len(rnn.Wrec[0])):
      #     print rnn.Wrec[i,:]
      # heat_map = sb.heatmap(rnn.Wrec)
      # plt.title('Heat map of $W_{rec}$ weights matrix')
      # plt.ylabel('Rows')
      # plt.xlabel('Columns')
      # plt.show()
  
      # heat_map = sb.heatmap(rnn.Wout)
      # plt.title('Heat map of $W_{out}$ weights matrix')
      # plt.ylabel('Rows')
      # plt.xlabel('Columns')
      # plt.show()
  
      # for i in range(len(rnn_zs[0,1,:100])):
      #     # print len(rnn_zs[0,1,:])
      #     print rnn_zs[0,1,i]
      tstart = 0
      tend = 200
  
  
      for i in range(1):
          num_ticks = 3#Nout+1
          xticks = np.linspace(0, len(rnn.t[tstart:-1])-1, num_ticks, dtype=np.int)
          yticks = np.linspace(1, Nout, num_ticks, dtype=np.int)
          xticklabels = [int(idx*dt/1000.0 + 0.1) for idx in xticks]
          yticklabels = [idy for idy in yticks]
          heat_map = sb.heatmap(rnn_zs[i,:,tstart:], xticklabels=xticklabels, yticklabels=yticklabels)
          heat_map.set_xticks(xticks)
          heat_map.set_yticks(yticks)
          plt.ylabel('Nodes')
          plt.xlabel('Time (sec)')
          plt.show()
  
  
      #plt.rc('legend', fontsize=7)    # legend fontsize
      lw = 2
      plt.plot(rnn.t[tstart:-1]/(20*tau), rnn.u[0,tstart:-1], '--', linewidth=lw)
      legend = ['$u(t)$']
      for j in range(Nout-1):
          plt.plot(rnn.t[tstart:-1]/(20*tau), rnn_zs[0,j,tstart:-1], linewidth=lw)
          # plt.plot(rnn.t[tstart:-1]/(20*tau), rnn.u[0,tstart:-1]*j, '--', linewidth=lw)
          legend.append('$z_{}(t)$'.format(j+1))
      plt.plot(rnn.t[tstart:-1]/(20*tau), rnn_zs[0,Nout-1,tstart:-1], linewidth=lw)
      legend.append('$z_{8}(t)$')
      plt.ylabel('$u(t)$ or $z(t)$')
      # plt.ylabel('DC Input, $u(t)$')
      plt.xlabel('Time, t (sec)')
  #    plt.legend(legend)#, loc="upper left")
      plt.show()
  
  
  
          # Limit cycle figures saved
     # results_dir = 'exp4/Best_model_experiments/3tobase_experiments/23a_best_model'
     # if not os.path.isdir(results_dir):
     #     os.makedirs(results_dir)
  
      tb = 315
      for i in range(1):
          figure = plt.figure()
  #        plt.quiver(rnn_zs[0,i,113:-1], rnn_zs[0,i+1,113:-1], rnn_zs[0,i,114:]-rnn_zs[0,i,113:-1], rnn_zs[0,i+1,114:]-rnn_zs[0,i+1,113:-1], scale_units='xy', angles='xy', scale=1)
          plt.quiver(rnn_zs[0,i,tb:-1], rnn_zs[0,i+1,tb:-1], rnn_zs[0,i,tb+1:]-rnn_zs[0,i,tb:-1], rnn_zs[0,i+1,tb+1:]-rnn_zs[0,i+1,tb:-1], scale_units='xy', angles='xy', scale=1)
          plt.title('Phase diagram of $x(t)$')
          plt.ylabel('$x_{}(t)$'.format(i+2))
          plt.xlabel('$x_{}(t)$'.format(i+1))
          # plt.xlim((-7.5,-1))
          # plt.ylim((-7.5,-1))
  #        figure.savefig(results_dir+'/phase_x_{}_x{}_2nd_trigger_4thloc.png'.format(i+1,i+2),dpi=200)
  #        figure.savefig(results_dir+'/phase_x_{}_x{}.png'.format(i+1,i+2),dpi=200)
          plt.show()
  