"""
This module was created by Miguel Capllonch Juan on 13 November 2018.
It handles most aspects of the system's biophysics
"""

import os
import csv
import numpy as np
from neuron import h
from datetime import datetime
import planar as pl

import algebra as alg
import workspace as ws
import anatomy
import electrodes
import analysis as als
import tools as tls



class AxonGeneralModel():
	"""
	This is a class to put in all the possible properties of an 
	axon model type, including anatomical and biophysical general properties
	This class is only instantiated once per model type, so each 
	instance of a physical axon does not have to do everything 
	from scratch
	"""
	def __init__(self, model):
		self.model = model
		self.properties_ = ws.axonmodel_settings[model]

	def build_regressions(self):
		"""
		Build the necessary regression functions of the axon model from the 
		json files in the work space
		"""
		if 'regressions' in self.properties_.keys():
			regfile = os.path.join(ws.folders['models'], 
				self.properties_['regressions'])
			self.regressions = alg.parameter_regressions(regfile)


class Cable:
	"""
	This is the basic class for any cable, be this an axon or a NAELC
	"""
	def __init__(self):

		# Counter. Add one to the class instances counter and save it 
		# as an id number as well
		ws.counters['cables'] += 1
		self.cable_number = ws.counters['cables'] - 1

		# Flag indicating if this cable has recordings
		self.has_recordings = False

		# Basic geometrical properties
		if ws.settings["nerve model"] == "resistor network":
			self.x = ws.nvt.x[self.cable_number]
			self.y = ws.nvt.y[self.cable_number]
			self.r = ws.nvt.r[self.cable_number]
		elif ws.settings["nerve model"] == "simple":
			self.x = ws.nvt.x[self.cable_number + ws.nvt.nNAELC]
			self.y = ws.nvt.y[self.cable_number + ws.nvt.nNAELC]
			self.r = ws.nvt.r[self.cable_number + ws.nvt.nNAELC]

		# Create properties dictionary
		self.properties = {}

		# Find the fascicle to which this cable belongs
		self.properties["fascicle"] = ws.nvt.cables_fascicles[self.cable_number]
		
	def finish_build(self):
		"""
		Perform the last steps in building this cable model
		"""

		# Append some properties to the corresponding dictionaries
		i = self.cable_number

		# Append the cell to the list of cables in the work space
		ws.cables_hoc.append(self.cable_cell)

		# Sections and lengths
		secs = list(self.cable_cell.all)
		ws.sections[i] = secs[:] 
		ws.sectypes[i] = self.properties['sctyps']
		ws.seclens[i] = [sec.L for sec in secs]
		# print('created seclens:', i)
		
		# Segments
		for j, sec in enumerate(secs):
			ws.secsegs[(i, j)] = list(sec.allseg())[1:-1]
		segs = tls.flattenlist([list(sec.allseg())[1:-1] for sec in ws.sections[i]])
		ws.segments[i] = segs[:]
		
		# Longitudinal profile (z-profiles)
		ws.zprofs[i] = anatomy.zprofile(self.cable_cell)
		
		# Maximum number of segments
		if len(ws.zprofs[i]) > ws.maxnseg:
			ws.maxnseg = len(ws.zprofs[i])

		# Store attributes
		self.properties['zprofile'] = ws.zprofs[i]

		# Intracellular area
		try:
			self.intracellular_area = np.pi * (0.5 * self.properties["axonD"]) ** 2
		except KeyError:
			# This is a NAELC or is an axon not configured to have axonD (which will be tackled if it happens)
			self.intracellular_area = 0.

	def set_recordings(self):
		"""
		Prepare the recording vectors for the desired variables
		"""
		# ws.log("Setting recordings for cable %i"%self.cable_number)

		# Update the flag indicating if this cable has recordings
		self.has_recordings = True

		# Variables to record from the settings
		rvars = ws.settings["data records"]["variables"]

		# Dictionary with all the recordings:
		# One entry per variable. Each entry contains a list of 
		# recording identifiers (one for each variable and segment of  
		# the cable)
		self.recordings = {}

		# Which sections to record
		where_to_record = ws.settings["data records"]["record sections"]

		# If this cable is a NAELC, record only vext[0]
		if self.cabletype == 'NAELC':
			rvars = ['vext[0]']
			where_to_record = ["all"]

		# Record the variables
		timings = []
		for var in rvars:

			# Create the list
			self.recordings[var] = []

			# Create the recording vectors
			for i, seg in enumerate(ws.segments[self.cable_number]):

				record = False
				
				t0 = datetime.now()

				if where_to_record[0] == "all":
					record = True
				else:
					# for sectypename in where_to_record:
					# 	if self.properties["sctyps"][i] == sectypename:
							# record = True
					if self.properties["sctyps"][i] in where_to_record:
						record = True

				if record:
					self.recordings[var].append({
						"number": ws.counters['recordings'], 
						"segment": str(seg)
						})
					als.set_recording(seg, var)

				t1 = datetime.now()
				timings.append((t1 - t0).total_seconds())

		ws.log("\tAverage time needed by als.set_recording for cable %i: %f ms. Max: %f ms"%(self.cable_number, 1e3 * np.array(timings).mean(), 1e3 * np.array(timings).max()))


class Axon(Cable):
	"""
	This is a class to put in all the possible properties of an 
	axon model, including anatomy and biophysics
	This class inherits its properties anatomical and biophysical 
	general properties from THE INSTANCE of AxonGeneralModel which 
	describes the model type
	"""
	def __init__(self, model):

		# Inheritance
		Cable.__init__(self)
		# Type of cable
		self.cabletype = 'Axon'
		# print('***Creating axon with the model %s'%model)
		# Counter
		ws.counters['axons'] += 1
		# Specific identifier; specific for this type of cable
		self.specific_number = ws.counters['axons'] - 1
		# Properties
		self.__dict__.update(ws.ax_gral_models[model].__dict__.copy())
		# Put AxonGeneralModel.properties_ into self.properties
		self.properties.update(self.properties_)
		self.properties['fiber_radius'] = self.r
		self.properties['fiberD'] = 2. * self.properties['fiber_radius']
		# Axon model
		# self.properties['model'] = ws.nvt.models[self.cable_number]

	def build(self):
		"""
		General class that chooses between the build methods for 
		myelinated or unmyelinated axons
		"""
		build_methods = {
			'myelinated': self.build_myelinated, 
			'unmyelinated': self.build_unmyelinated
		}
		build_methods[self.properties['type']]()

		# Finish building this cable and saving its properties
		self.finish_build()


	def build_unmyelinated(self):
		"""
		Build the axon as unmyelinated
		"""
		nseg = anatomy.nseg_dlambda_rule(
			ws.length, 
			self.properties['fiber_radius'], 
			self.properties['rhoa'], 
			self.properties['cm']
			)
		self.properties['nseg'] = nseg
		self.cable_cell = h.UnMyelAxon(self.properties)
		self.properties['sctyps'] = np.zeros(nseg, dtype=np.int)

	def build_myelinated(self):
		"""
		Build the necessary properties of the axon model from the 
		json files in the work space
		"""

		# Variables that depend on a regression from some stored data
		iv = 'independent variable'
		rgs = self.regressions
		for k, func in rgs.items():
			if k != iv:
				x = self.properties[rgs[iv]]
				self.properties[k] = func(x)

		# Other properties:

		# Variable name shortenings
		deltax = self.properties['deltax']

		# 1. STIN_length
		# This is still model-specific hard-coding
		if self.model == 'MRG' or self.model == 'gaines_sensory':
			nodl = self.properties['nodelength']
			pl1 = self.properties['paralength1']
			pl2 = self.properties['paralength2']
			self.properties['STIN_length'] = deltax - nodl - 2. * (pl1 + pl2)
		elif self.model == 'MRG_STINonly':
			nodl = self.properties['nodelength']
			self.properties['STIN_length'] = deltax - nodl

		# Start point along the z-axis
		if ws.settings["nerve model"] == "resistor network":
			self.properties["start"] = deltax * ws.nvt.start_positions[self.cable_number]
		elif ws.settings["nerve model"] == "simple":
			self.properties["start"] = deltax * ws.nvt.start_positions[self.cable_number + ws.nvt.nNAELC]

		# 3. Total internodal length
		self.properties['in_l'] = deltax - self.properties['nodelength']

		# 4. Update the minimum length of one chunk of sections
		# and of the node to fiber ratio, g
		ws.deltax_min = min((ws.deltax_min, deltax))
		ws.g_min = min((ws.g_min, self.properties['g']))

		# Build the axon in hoc
		sctyps, section_counter, n_max, lengths = anatomy.all_vars(self)

		# Lengths to hoc
		lengths_hoc = {}
		for key, value in lengths.items():
			lengths_hoc[key] = h.Vector(value)

		# Number of sections of each type and section types
		nsecs = sum(section_counter.values())
		sectypes = h.List()
		for st in sctyps:
			sectypes.append(st)

		# Store variables as attributes inside the 'properties' dictionary
		self.properties['nsecs'] = nsecs
		self.properties['sectypes'] = sectypes
		self.properties['section_counter'] = section_counter
		self.properties['mx_n'] = n_max
		self.properties['lengths'] = lengths_hoc

		# Axon in hoc
		if self.model == 'MRG':
			# print(' ----------- Setting up %s'%self.model)
			self.cable_cell = h.MRGMyelAxon(self.properties)
		elif self.model == 'gaines_sensory':
			# print(' ----------- Setting up %s'%self.model)
			self.cable_cell = h.Gaines2018Sensory(self.properties)
		
		# z-profile
		zp = anatomy.zprofile(self.cable_cell)

		# If this is the smallest axon, save its longitudinal profile
		if self.properties['fiber_radius'] == ws.raxmin:
			ws.zlayers_min = zp

		# Number of segments on each section of this axon
		this_nseg = [sec.nseg for sec in self.cable_cell.all]

		# I need sctyps to reflect all the segments, not just the sections
		sctyps_str = [this_nseg[ii] * [sctyps[ii].sectype] for ii in range(len(sctyps))]
		sctyps_str = np.array(tls.flattenlist(sctyps_str))

		# Get the z-profiles of the nodes of Ranvier only
		zprofs_RN = zp[np.where(sctyps_str == 'node')[0]]
		ws.zprofs_RN[self.cable_number] = zprofs_RN

		# Store these properties
		self.properties['zprofs_RN'] = zprofs_RN
		self.properties['sctyps'] = sctyps_str


class NAELC(Cable):
	""" 
	This is a Non-Axonal Cable 
	"""
	def __init__(self, *dummyargs):
		# Inheritance
		Cable.__init__(self)
		# Type of cable
		self.cabletype = 'NAELC'
		# Counter
		ws.counters['NAELC'] += 1
		# Specific identifier; specific for this type of cable
		self.specific_number = ws.counters['NAELC'] - 1
		# Properties
		self.properties = self.properties.copy()
		self.properties.update({
					'diam': ws.cdiam, 
					'length': ws.length, 
					'extlayers': ws.NAELC_extlayers
				})

	def build(self):
		"""
		Build the cable into the hoc space
		"""
		self.properties['nseg'] = ws.NAELC_nseg
		self.properties['sctyps'] = np.zeros(ws.NAELC_nseg, dtype=np.int)
		self.cable_cell = h.Wire(self.properties)

		# Finish building this cable and saving its properties
		self.finish_build()

class ResistorNetwork():
	"""
	This is a class for the Resistor Network that describes the 
	interactions between all axons
	"""
	def __init__(self):
		pass
		
	def get_resistances(self):
		""" Calculate the values of the transverse resistances """

		self.resistors = {}
		self.res_lengs = {}
		for pair in ws.pairs:
			i, j = pair
			# Calculation of the resistor's value: 
			# shape factor of the resistor (adimensional)
			# and value in MOhm*cm
			l = ws.nvt.len_con[pair]
			if ws.EC["resistor network"]["interaxonal connections"] == "membranes":
				l -= (ws.nvt.r[i] + ws.nvt.r[j])
			elif ws.EC["resistor network"]["interaxonal connections"] == "centers":
				pass
			
			# Width
			w = ws.nvt.len_seg[pair]
			# Shape factor (adimensional (um / um))
			shape_factor = l / w
			# Large shape factors mean large resistances; not interested
			if shape_factor < 100:

				# Resistance value
				res = ws.rho['endo_trans'] * shape_factor

				# Perineurium
				resist_perineurium = 0.

				# If the pair is at opposite sides of the perineurium, 
				# add the perineurium's resistance p.u.l or 
				# resistivity (Ohm*cm)
				if ('NAELC' in ws.nvt.cables[i] and 'xon' in ws.nvt.cables[j]) \
				or ('NAELC' in ws.nvt.cables[j] and 'xon' in ws.nvt.cables[i]):
					resist_perineurium += ws.resist_perineurium

				# Also, that happens when two axons are in connected but 
				# in different fascicles
				# if 'xon' in ws.nvt.cables[i] and 'xon' in ws.nvt.cables[j]:
				if ws.cables[i].cabletype == "Axon" and ws.cables[j].cabletype == "Axon":
					# ws.log("Checking if there\'s double layer of perineurium for axons %i and %i:"%(ws.cables[i].specific_number, ws.cables[j].specific_number))
					# ws.log("%s, %s"%(ws.cables[i].properties["fascicle"], ws.cables[j].properties["fascicle"]))
					if ws.cables[i].properties["fascicle"] != ws.cables[j].properties["fascicle"]:
						resist_perineurium += 2. * ws.resist_perineurium
						# ws.log("\tAdding double perineurium between axons %i and %i"%(ws.cables[i].specific_number, ws.cables[j].specific_number))
						# ws.log("\tThese axons are placed at (%f, %f) um and (%f, %f) um, respectively"%(ws.cables[i].x, ws.cables[i].y, ws.cables[j].x, ws.cables[j].y))

				# Finally, add it up
				res += resist_perineurium / w
				# Store in the dictionary
				self.resistors[pair] = res

	def connect_resistors(self):
		""" Connect the transverse resistors with ephapses in NEURON """

		def locate_and_create(pair, positions):
			""" For a pair of cables and a set of resistor positions, 
			locate, create and connect those resistors to the cables """

			# Pair indices
			i, j = pair

			# Get resistors' lengths
			rl = get_rl(pair, positions)

			# Iterate over resistors positions
			for zr, l in zip(positions, rl):

				# For each cable, see where this falls

				# On cable i:
				i_sec, z_on_i = anatomy.locate(ws.seclens[i], zr)
				seci = ws.sections[i][i_sec]
				sgi = seci(z_on_i)

				# On cable j:
				j_sec, z_on_j = anatomy.locate(ws.seclens[j], zr)
				secj = ws.sections[j][j_sec]
				sgj = secj(z_on_j)

				# Value of the resistor (MOhm)
				ephres = res / l
				# Create the resistor
				create_eph(ephres, sgi, sgj, seci, secj, i, j, 
					zr)
				# ws.log("Created a transverse resistor at z = %f um with value %f MOhm"%(zr, ephres))

		ws.eph_counter = 0
		ws.nrn_res = {}
		for pair, res in self.resistors.items():
			i, j = pair

			# Cable
			cable_cell_i = ws.cables_hoc[i]
			cable_cell_j = ws.cables_hoc[j]

			# I have to connect the cables in a staggered manner

			# If I am only taking the nodes of Ranvier into account
			# and the connection is to be set between two axons:
			if ws.EC["resistor network"]["transverse resistor locations"] == "regular":

				# Separation between transverse resistors
				sep = float(ws.EC["resistor network"]["resistors separation"])

				# Positions of the resistors
				# res_positions = np.arange(0, ws.length, sep)
				res_positions = np.arange(0.5, ws.length, sep)

				# Create and connect the resistors
				locate_and_create(pair, res_positions)

			elif ws.EC["resistor network"]["transverse resistor locations"] == "only nodes of Ranvier" and \
			('xon' in ws.cables[i].cabletype and 'xon' in ws.cables[j].cabletype):

				# Get resistors' positions and lengths (um)

				# Merge the ws.zprofs_RN of the two axons
				# This array will contain the postions of the resistors 
				# over the z-axis
				zp_pair = np.append(ws.zprofs_RN[i], ws.zprofs_RN[j])
				zp_pair.sort()
				# Remove repeated values (can --will, actually, happen if 
				# there's aligned nodes)
				zp_pair = np.unique(zp_pair)

				# Create and connect the resistors
				locate_and_create(pair, zp_pair)

			else:
				# If I am taking all nodes into account, the compartments of only one
				# of the two cables (arbitrarily chosen; cable 'i' by default) 
				# are plugged to transverse resistors. These are then connected 
				# to the second cable wherever they must, regardless of the 
				# topology of this second cable.

				# Get resistors' lengths (um)
				# From what is explained above, in this case the z-profile of one
				# of the cables only is used as zp_pair

				rl = get_rl(pair, ws.zprofs[i])
				
				# Iterate over segments
				for ii, sgi in enumerate(ws.segments[i]):

					# Get the section and segment of cable j which sgi lies against
					zi = ws.zprofs[i][ii]
					j_sec, z_ioverj = anatomy.locate(ws.seclens[j], zi)
					seci = sgi.sec
					secj = ws.sections[j][j_sec]
					sgj = secj(z_ioverj)

					# Value of the resistor (MOhm)
					ephres = res / rl[ii]
					# Create the resistor
					create_eph(ephres, sgi, sgj, seci, secj, i, j, 
						zi)

########################################################################
# Functions

def finish_fascicles():
	"""
	Finish giving basic information to the fascicles once the cables have been set up and the resistivities have been adjusted (if they have)
	"""
	for k, fas in ws.nvt.fascicles.items():
		# Intracellular area
		ws.nvt.fascicles[k].get_intracellular_area()
		# Print properties
		fas.print_properties()
	# Now print the properties of the nerve
	ws.nvt.print_properties()

def get_rl(pair, zp_pair):
	"""
	Get the lengths (all over the z-axis) of all the transverse
	resistors for a pair of cables
	They will be used to obtain the values of each transverse resistor in MOhm.
	"""

	# Distances between nodes
	dz = zp_pair[1:] - zp_pair[:-1]
	# Get the lengths of the resistors (um)
	rl = 0.5 * (dz[1:] + dz[:-1])
	# Add one element on the left: first resistor
	# and one on the right: last resistor
	rl = np.append(np.array([zp_pair[0] + 0.5 * dz[0]]), rl)
	rl = np.append(rl, np.array([0.5 * dz[-1] + ws.length - zp_pair[-1]]))
	# Store its vale
	ws.rnet.res_lengs[pair] = rl
	return rl

# def create_eph(ephres, sgi, sgj, seci, secj, i, j, nrn_res, zr):
def create_eph(ephres, sgi, sgj, seci, secj, i, j, zr):
	"""
	Create a resistor that acts as an ephaptic connection
	"""

	# Choose the layer depending on the type of fiber or cable
	# Layer where to point to
	layeri = ws.cables[i].properties['extlayers']
	layerj = ws.cables[j].properties['extlayers']

	# If zr is not in a cuff, we will make the ephaptic connection
	# super resistive and hence disappear
	if not ws.EC["resistor network"]["outside cuffs"]:
		ephres = disconnect_under_cuff(zr, ephres)

	# If the connection is too large, do not implement it
	if ephres < 1.e2:

		eph = h.EphapDBA(2)
		eph.ex0_loc(0, sgi.x, layeri, sec=seci)
		eph.ex1_loc(1, sgj.x, layerj, sec=secj)
		eph.add_submatrix(0, 1, ephres)
		eph.install()
		# Disconnect segments from ground
		# This is redundant, but just in case
		sgi.xg[layeri - 1] = 0.
		sgj.xg[layerj - 1] = 0.

		# Add to dictionary
		ws.nrn_res = tls.append_items(ws.nrn_res, (i, j), [eph])
		ws.eph_counter += 1

def disconnect_under_cuff(zi, ephres):
	"""
	Check if a transverse resistor is not under the region of a cuff 
	electrode and if so, give it a highest value so it's not connected
	"""
	if electrodes.check_cuff_cover(zi) is None:
		ephres *= ws.rho['SuperResistive'] / ws.rho['endo_trans']
	return ephres

def build_cables():
	"""
	Build all the cables, NAELC and fibers
	"""

	ws.log("Building the cables in NEURON")
	t0 = datetime.now()

	# Dictionaries where to store information in the workspace
	ws.zprofs_RN = {}
	ws.sections = {}
	ws.sectypes = {}
	ws.secsegs = {}
	ws.seclens = {}
	ws.zprofs = {}
	ws.segments = {}
	cables = []
	ws.cables_hoc = h.List()

	# New building method
	# First, axon model properties
	# Axon general model(s)
	ws.ax_gral_models = {}
	# ws.ax_gral_models[ws.settings['axon model']] = AxonGeneralModel(ws.settings['axon model'])
	for m in ws.nvt.models_set:
		if 'NAELC' not in m:
			ws.ax_gral_models[m] = AxonGeneralModel(m)
			ws.ax_gral_models[m].build_regressions()

	# Classes of cables
	ws.cable_classes = {
		"NAELC": NAELC, 
		"Axon": Axon
	}

	# First, instantiate all the cables.
	# Then, build axons. From them, obtain the desired nseg for NAELC.
	# Finally, build the NAELC using this value

	# Instantiate everything
	nvt = ws.nvt
	for i in range(nvt.nc):

		# Create the cable(s)

		if ws.settings["nerve model"] == "resistor network":
			# Use NAELC
			# cable = ws.cable_classes[nvt.cables[i]](ws.settings['axon model'])
			cable = ws.cable_classes[nvt.cables[i]](nvt.models[i])
			# Build it only if it is an axon
			# The NAELC will be built later on
			if nvt.cables[i] == 'Axon':
				# print('building', i, nvt.cables[i])
				cable.build()
			# Append it
			cables.append(cable)

		elif ws.settings["nerve model"] == "simple":
			# Don't use NAELC
			if nvt.cables[i] == 'Axon':
				# cable = ws.cable_classes[nvt.cables[i]](ws.settings['axon model'])
				cable = ws.cable_classes[nvt.cables[i]](nvt.models[i])
				# Build the axon
				cable.build()
				# Append it
				cables.append(cable)
	
	# Obtain nseg for the NAELC if there are axons
	if ws.counters['axons'] > 0:
		# There's axons
		daxmin = 2. * ws.raxmin
		if ws.EC["resistor network"]["transverse resistor locations"] == "only nodes of Ranvier":
			nseg = int(ws.length / ws.deltax_min) + 1
		elif ws.EC["resistor network"]["transverse resistor locations"] == "regular":
			nseg = int(ws.length / float(ws.EC["resistor network"]["resistors separation"])) + 1
		else:
			# nseg for an internode
			nseg = anatomy.nseg_dlambda_rule(ws.deltax_min, ws.g_min * daxmin, 
				cables[ws.counters['NAELC']].properties['rhoa'], cm=cables[ws.counters['NAELC']].properties['cm'])
			# nseg for the whole nerve length
			nseg = int(ws.length / (ws.deltax_min / nseg)) + 1

		# Make sure it's an odd number
		# This adds 2 if it's odd, but it doesn't matter; it's even better
		nseg += nseg % 2 + 1

	else:
		# If we're here, there are no axons
		nseg = ws.nseg

	ws.NAELC_nseg = nseg
	ws.log("NAELC segment length = %f"%(ws.length / nseg))
	ws.log("nseg: %i"%nseg)
	if ws.counters['axons'] > 0:
		ws.log("Min. deltax: %f"%ws.deltax_min)

	# Now build the NAELC according to the value of nseg 
	# obtained from the axons
	for i in range(ws.counters['NAELC']):
		# print('building', i, cables[i])
		cables[i].build()

	# print(ws.nvt.cables)

	# Store the cables in the workspace
	ws.cables = cables

	# Maximum number of segments. This is a number to be found
	ws.maxnseg = 0

	t1 = datetime.now()
	ws.log("Created %i cables in %i seconds"%(ws.nvt.nc, (t1 - t0).total_seconds()))

def electrics():
	"""
	In this function, xraxial is set for every axons and all 
	connections to ground are made
	"""

	# The following loop over cables and segments and calculates:
	# - The cross-sectional area of each cable and therefore:
	#  - its value for xraxial
	#  - its connections to ground at the ends of the nerve

	for i, seglist in ws.segments.items():
		
		# Level or layer where to point to
		level = ws.cables[i].properties['extlayers'] - 1

		# Cross-sectional area and xraxial
		# Cross-sectional area assigned to the cable (cm2)
		carea = ws.nvt.free_areas[i] * 1.e-8

		# Attribute xraxial
		for j, segment in enumerate(seglist):
			# First and foremost, disconnect ALL segments from ground
			segment.xg[level] = 0.

			# If there's no EC outside the cuff-covered regions,
			# check if this point is under a cuff and connect it to ground
			if not ws.EC["resistor network"]["outside cuffs"]:
				cuff_cover = electrodes.check_cuff_cover(ws.zprofs[i][j])
				if cuff_cover is None:
					segment.xg[level] = 1.e9

			# Attribute xraxial (MOhm/cm; so I have to convert rho (in MOhm*um) 
			# to MOhm*cm)

			# Longitudinal resistivity for this cable as a weighted average between epineurium and endoneurium (the endoneurium is, in turn, a weighted average between all the fascicles intersecting the cross-sectional area of this cable)
			ctendo = ws.nvt.cables_tissues[i]['endoneurium']
			ctepi = ws.nvt.cables_tissues[i]['epineurium']
			fas_ = ws.nvt.fascicles
			resistivity = ws.rho['endo_longt'] * sum([ctendo[k] for k in fas_]) + ctepi * ws.rho['epi_longt']

			# Apply the resistivity in xraxial
			xrx =  1.e-4 * resistivity / carea
			segment.xraxial[level] = xrx

		# Connect to ground all cables at both ends of the nerve
		if ws.EC["resistor network"]["outside cuffs"]:
			
			segleft = seglist[0]
			segrght = seglist[-1]
			
			# Note: The areas of each segment need to be in cm2
			segleft.xg[level] = ws.xg_distal[i] / (segleft.area() * 1.e-8)
			segrght.xg[level] = ws.xg_proxim[i] / (segrght.area() * 1.e-8)
		else:
			# In this case, just ground the ends
			# This is regardless of the ends being cuffed or not; 
			# these are at the base of the cylinder
			seglist[0].sec(0.).xg[level] = 1.e9
			seglist[-1].sec(1.).xg[level] = 1.e9

	# Connect grounds to the nerve's walls (not the ends; these are already sorted)
	# The xg[level] in NEURON for each segment is in mho/cm2
	# All points in the contour are connected to a far ground
	for i in np.arange(ws.npc):
		# Level or layer where to point to
		level = ws.cables[i].properties['extlayers'] - 1
		for j, segment in enumerate(ws.segments[i]):

			# Convert area to cm2 (it's in um2)
			# Cable's cylindrical wall area
			area = segment.area() * 1.e-8

			# Set xg
			# This depends on whether a segment is under the cuff or not
			cuff_cover = electrodes.check_cuff_cover(ws.zprofs[i][j])
			if cuff_cover is not None:
				# xg must be in mho/cm2, so xg_rwall should be in mho; it is

				# If this segment is at one end of the nerve, it already 
				# has one connection to ground that needs to be added
				if segment in (ws.segments[i][0], ws.segments[i][-1]):
					segment.xg[level] += \
						ws.xg_rwall['with_cuff'][cuff_cover] / area
				else:
					segment.xg[level] = \
						ws.xg_rwall['with_cuff'][cuff_cover] / area
			else:
				if ws.EC["resistor network"]["outside cuffs"]:
					# If this segment is at one end of the nerve, it already 
					# has one connection to ground that needs to be added
					if segment in (ws.segments[i][0], ws.segments[i][-1]):
						segment.xg[level] += \
							ws.xg_rwall['wout_cuff'] / area
					else:
						segment.xg[level] = \
							ws.xg_rwall['wout_cuff'] / area
				else:
					# If this is not meant to be used, just ground it
					# (in a future version, I will simply crop these
					# bits of NAELC outside the cuffs)
					segment.xg[level] = 1.e9

def set_direct_extracellular_stimulation():
	""" Set up the stimulation over e_extracellular on the axons """

	# File to read the data from
	espath = os.path.join(ws.folders['data/load'], 'extstim.csv')

	# Create the time window
	# All ones for now
	# It will be cropped when the information for each pulse is fetched
	time_window = np.ones(ws.nt)

	# Time series of the extracellular field for the axons
	stim_tseries = []

	# Dictionaries with the arrays and more info
	pulses = {}

	# Open and read
	with open(espath, 'r') as f:
		reader = csv.reader(f)

		# Read line by line
		for row in reader:

			# First element of each row
			key = row[0]

			# New pulse
			if key == 'Pulse':
				k = int(row[1])
				pulses[k] = {}
				# Reset the time window
				time_window[:] = 1.

			else:

				# Parse the information in the line
				try:
					# See if this is an axon
					i = int(key)
				except ValueError:
					# If not, see if row[1] is a numerical value
					try:
						pulses[k][key] = float(row[1])
					except ValueError:
						# If not, it's a string; just save it as such 
						pulses[k][key] = row[1]
				else:
					# If it is an axon, save the array
					field = np.array([float(x) for x in row[1:]])

					# If we're here, this means we already have all the information
					# Give values to the time window (only once)
					if i == 0:
						if pulses[k]['Type'] == 'square pulse':
							time_window[np.where(ws.tarray < pulses[k]['Delay'])] = 0.
							time_window[np.where(ws.tarray >= pulses[k]['Delay'] + pulses[k]['Duration'])] = 0.

					# If this is the first pulse, create the time series for 
					# this axon
					if k == 0:
						stim_tseries.append(np.zeros((len(field), len(time_window))))

					# Now create the time series of the extracellular field 
					# for each segment of this axon and play it in NEURON
					for j, seg in enumerate(ws.segments[i]):
						stim_tseries[i][j] += field[j] * time_window
			
	# Now that each axon has its full stimulation time 
	# series, play it (fts: 'field time series')

	ws.stim_field_vectors = []

	for i in list(ws.segments.keys())[ws.counters['NAELC']:]:

		for j, seg in enumerate(ws.segments[i]):

			# Create and play the time series
			fts = h.Vector(stim_tseries[i][j])
			fts.play(seg._ref_e_extracellular, h.dt)

			# Add it to the list
			# This is necessary for this to work
			# Otherwise, fts is lost
			ws.stim_field_vectors.append(fts)