# -* coding: utf-8 -*-
# Utility functions for moose.

from __future__ import print_function, division, absolute_import
from __future__ import print_function, division
from __future__ import absolute_import

import types
import parser
import token
import symbol
import os
import math
import warnings
from datetime import datetime
from collections import defaultdict
import re

from moose.moose_constants import *
from moose.print_utils import *
# Print and Plot utilities.
try:
    from moose.plot_utils import *
except Exception as e:
    info( "Plot utilities are not loaded due to '%s'" % e )

def create_table_path(model, graph, element, field):

    field = field[0].upper() + field[1:]

    tablePathSuffix = element.path.partition(model.path)[-1]
    if tablePathSuffix.startswith('/'):
        tablePathSuffix = tablePathSuffix[1:]

    tablePathSuffix =  tablePathSuffix.replace('/', '_') + '.' + field
    tablePathSuffix = re.sub( '.'
                            , lambda m: { '[' : '_'
                                        , ']' : '_'
                                        }.get(m.group(), m.group())
                            ,tablePathSuffix
                            )

    if tablePathSuffix.startswith("_0__"):
        tablePathSuffix = tablePathSuffix[4:]

    #tablePath = dataroot + '/' +tablePath
    tablePath = graph.path + "/" + tablePathSuffix
    return tablePath


def create_table(tablePath, element, field,tableType):
    """Create table to record `field` from element `element`

    Tables are created under `dataRoot`, the names are generally
    created by removing `/model` in the beginning of `elementPath`
    and replacing `/` with `_`. If this conflicts with an existing
    table, the id value of the target element (elementPath) is
    appended to the name.
    """
    if moose.exists(tablePath):
        table = moose.element(tablePath)
    else:
        if tableType == "Table2":
            table = moose.Table2(tablePath)
        elif tableType == "Table":
            table = moose.Table(tablePath)
        moose.connect(table, 'requestOut', element, 'get%s' % (field))
    return table

def readtable(table, filename, separator=None):
    """Reads the file specified by filename to fill the MOOSE table object.

    The file can either have one float on each line, in that case the
    table will be filled with values sequentially.
    Or, the file can have
    index value
    on each line. If the separator between index and value is anything other than
    white space, you can specify it in the separator argument."""

    in_file = open(filename)
    ii = 0
    line_no = 0
    for line in in_file:
        line_no = line_no + 1
        tokens = line.split(separator)
        if len(tokens) is 0:
            continue
        elif len(tokens) == 1:
            table[ii] = float(tokens[0])
        elif len(tokens) == 2:
            table[int(tokens[0])] = float(tokens[1])
        else:
            print("pymoose.readTable(", table, ",", filename, ",", separator
                    , ") - line#", line_no, " does not fit."
                    )

def getfields(moose_object):
    """Returns a dictionary of the fields and values in this object."""
    field_names = moose_object.getFieldNames('valueFinfo')
    fields = {}
    for name in field_names:
        fields[name] = moose_object.getField(name)
    return fields

def findAllBut(moose_wildcard, stringToExclude):
    allValidObjects = moose.wildcardFind(moose_wildcard)
    refinedList = []
    for validObject in allValidObjects:
        if validObject.path.find(stringToExclude) == -1:
            refinedList.append(validObject)

    return refinedList

def apply_to_tree(moose_wildcard, python_filter=None, value=None):
    """
    Select objects by a moose/genesis wildcard, apply a python filter on them and apply a value on them.

    moose_wildcard - this follows GENESIS convention.

    {path}/#[{condition}] returns all elements directly under {path} that satisfy condition. For example:

    '/mynetwork/mycell_0/#[TYPE=Compartment]'

    will return all Compartment objects directly under mycell_0 in mynetwork.

    '{path}/##[{condition}]' will recursively go through all the
    objects that are under {path} (i.e. children, grandchildren,
    great-grandchildren and so on up to the leaf level) and a list of
    the ones meet {condition} will be obtained.

    Thus, '/mynetwork/##[TYPE=Compartment]' will return all
    compartments under mynetwork or its children, or children thereof
    and so on.

    python_filter - if a single string, it will be taken as a
    fieldname, and value will be assigned to this field. It can also
    be a lambda function returning True or False which will be applied
    to each id in the id list returned by moose wildcard
    search. Remember, the argument to the lambda will be an Id, so it
    is up to you to wrap it into a moose object of appropriate type. An example is:

    lambda moose_id: Compartment(moose_id).diameter <  2e-6

    If your moose_wildcard selected objects of Compartment class, then
    this lambda function will select only those with diameter less
    than 2 um.

    value - can be a lambda function to apply arbitrary operations on
    the selected objects.

    If python_filter is a string it, the return
    value of applying the lambda for value() will assigned to the
    field specified by python_filter.

    But if it is value is a data object and {python_filter} is a
    string, then {value} will be assigned to the field named
    {python_filter}.


    If you want to assign Rm = 1e6 for each compartment in mycell
    whose name match 'axon_*':

    apply_to_tree('/mycell/##[Class=Compartment]',
            lambda x: 'axon_' in Neutral(x).name,
            lambda x: setattr(Compartment(x), 'Rm', 1e6))

    [you must use setattr to assign value to a field because lambda
    functions don't allow assignments].
    """
    if not isinstance(moose_wildcard, str):
        raise TypeError('moose_wildcard must be a string.')
    id_list = moose.getWildcardList(moose_wildcard, True)
    if isinstance(python_filter, types.LambdaType):
        id_list = [moose_id for moose_id in id_list if python_filter(moose_id)]
    elif isinstance(python_filter, str):
        id_list = [moose_id for moose_id in id_list if hasattr(eval('moose.%s(moose_id)' % (moose.Neutral(moose_id).className)), python_filter)]
    else:
        pass
    if isinstance(value, types.LambdaType):
        if isinstance(python_filter, str):
            for moose_id in id_list:
                moose_obj = eval('moose.%s(moose_id)' % (moose.Neutral(moose_id).className))
                setattr(moose_obj, python_filter, value(moose_id))
        else:
            for moose_id in id_list:
                value(moose_id)
    else:
        if isinstance(python_filter, str):
            for moose_id in id_list:
                moose_obj = eval('moose.%s(moose_id)' % (moose.Neutral(moose_id).className))
                setattr(moose_obj, python_filter, value)
        else:
            raise TypeError('Second argument must be a string specifying a field to assign to when third argument is a value')


def tweak_field(moose_wildcard, field, assignment_string):
    """Tweak a specified field of all objects that match the
    moose_wildcard using assignment string. All identifiers in
    assignment string must be fields of the target object.

    Example:

    tweak_field('/mycell/##[Class=Compartment]', 'Rm', '1.5 / (3.1416 * diameter * length')

    will assign Rm to every compartment in mycell such that the
    specific membrane resistance is 1.5 Ohm-m2.
    """
    if not isinstance(moose_wildcard, str):
        raise TypeError('moose_wildcard must be a string.')
    id_list = moose.getWildcardList(moose_wildcard, True)
    expression = parser.expr(assignment_string)
    expr_list = expression.tolist()
    # This is a hack: I just tried out some possible syntax trees and
    # hand coded the replacement such that any identifier is replaced
    # by moose_obj.identifier
    def replace_fields_with_value(x):
        if len(x)>1:
            if x[0] == symbol.power and x[1][0] == symbol.atom and x[1][1][0] == token.NAME:
                field = x[1][1][1]
                x[1] = [symbol.atom, [token.NAME, 'moose_obj']]
                x.append([symbol.trailer, [token.DOT, '.'], [token.NAME, field]])
            for item in x:
                if isinstance(item, list):
                    replace_fields_with_value(item)
        return x
    tmp = replace_fields_with_value(expr_list)
    new_expr = parser.sequence2st(tmp)
    code = new_expr.compile()
    for moose_id in id_list:
        moose_obj = eval('moose.%s(moose_id)' % (moose.Neutral(moose_id).className))
        value = eval(code)
        moose.setField(moose_id, field, str(value))

# 2012-01-11 19:20:39 (+0530) Subha: checked for compatibility with dh_branch
def printtree(root, vchar='|', hchar='__', vcount=1, depth=0, prefix='', is_last=False):
    """Pretty-print a MOOSE tree.

    root - the root element of the MOOSE tree, must be some derivatine of Neutral.

    vchar - the character printed to indicate vertical continuation of
    a parent child relationship.

    hchar - the character printed just before the node name

    vcount - determines how many lines will be printed between two
    successive nodes.

    depth - for internal use - should not be explicitly passed.

    prefix - for internal use - should not be explicitly passed.

    is_last - for internal use - should not be explicitly passed.

    """
    root = moose.element(root)
    # print('%s: "%s"' % (root, root.children))
    for i in range(vcount):
        print(prefix)

    if depth != 0:
        print(prefix + hchar, end=' ')
        if is_last:
            index = prefix.rfind(vchar)
            prefix = prefix[:index] + ' ' * (len(hchar) + len(vchar)) + vchar
        else:
            prefix = prefix + ' ' * len(hchar) + vchar
    else:
        prefix = prefix + vchar

    print(root.name)
    children = []
    for child_vec in root.children:
        try:
            child = moose.element(child_vec)
            children.append(child)
        except TypeError:
            pass
            # print 'TypeError:', child_vec, 'when converting to element.'
    for i in range(0, len(children) - 1):
        printtree(children[i],
                  vchar, hchar, vcount, depth + 1,
                  prefix, False)
    if len(children) > 0:
        printtree(children[-1], vchar, hchar, vcount, depth + 1, prefix, True)


def df_traverse(root, operation, *args):
    """Traverse the tree in a depth-first manner and apply the
    operation using *args. The first argument is the root object by
    default."""
    if hasattr(root, '_visited'):
        return
    operation(root, *args)
    for child in root.children:
        childNode = moose.Neutral(child)
        df_traverse(childNode, operation, *args)
    root._visited = True

def autoposition(root):
    """Automatically set the positions of the endpoints of all the
    compartments under `root`.

    This keeps x and y constant and extends the positions in
    z-direction only. This of course puts everything in a single line
    but suffices for keeping electrical properties intact.

    TODO: in future we may want to create the positions for nice
    visual layout as well. My last attempt resulted in some
    compartments overlapping in space.

    """
    compartments = moose.wildcardFind('%s/##[TYPE=Compartment]' % (root.path))
    stack = [compartment for compartment in map(moose.element, compartments)
              if len(compartment.neighbors['axial']) == 0]

    assert len(stack) == 1, 'There must be one and only one top level\
            compartment. Found %d' % len(stack)
    ret = stack[0]
    while len(stack) > 0:
        comp = stack.pop()
        parentlist = comp.neighbors['axial']
        parent = None
        if len(parentlist) > 0:
            parent = parentlist[0]
            comp.x0, comp.y0, comp.z0, = parent.x, parent.y, parent.z
        else:
            comp.x0, comp.y0, comp.z0, = (0, 0, 0)
        if comp.length > 0:
            comp.x, comp.y, comp.z, = comp.x0, comp.y0, comp.z0 + comp.length
        else:
            # for spherical compartments x0, y0, z0 are centre
            # position nad x,y,z are on the surface
            comp.x, comp.y, comp.z, = comp.x0, comp.y0, comp.z0 + comp.diameter/2.0
        # We take z == 0 as an indicator that this compartment has not
        # been processed before - saves against inadvertent loops.
        stack.extend([childcomp for childcomp in map(moose.element, comp.neighbors['raxial']) if childcomp.z == 0])
    return ret

def loadModel(filename, target,method='ee'):
    moose.loadModel(filename,target)
    moose.mooseAddChemSolver(target,method)
    if moose.exists(target+'/kinetics/info'):
        moose.element(target+'/kinetics/info').solver = method

def readcell_scrambled(filename, target, method='ee'):
    """A special version for handling cases where a .p file has a line
    with specified parent yet to be defined.

    It creates a temporary file with a sorted version based on
    connectivity, so that parent is always defined before child."""
    pfile = open(filename, "r")
    tmpfilename = filename + ".tmp"
    graph = defaultdict(list)
    data = {}
    error = None
    root = None
    ccomment_started = False
    current_compt_params = []
    for line in pfile:
        tmpline = line.strip()
        if not tmpline or tmpline.startswith("//"):
            continue
        elif tmpline.startswith("/*"):
            ccomment_started = True
        if tmpline.endswith('*/'):
            ccomment_started = False
        if ccomment_started:
            continue
        if tmpline.startswith('*set_compt_param'):
            current_compt_params.append(tmpline)
            continue
        node, parent, rest, = tmpline.partition(' ')
        print('22222222', node, parent)
        if (parent == "none"):
            if (root is None):
                root = node
            else:
                raise ValueError("Duplicate root elements: ", root, node, "> Cannot process any further.")
                break
        graph[parent].append(node)
        data[node] = '\n'.join(current_compt_params)

    tmpfile = open(tmpfilename, "w")
    stack = [root]
    while stack:
        current = stack.pop()
        children = graph[current]
        stack.extend(children)
        print('#########"', current, '": ', data[current])
        tmpfile.write(data[current])
    tmpfile.close()
    ret = moose.loadModel(tmpfilename, target, method)
    return ret

## Subha: In many scenarios resetSim is too rigid and focussed on
## neuronal simulation.  The setDefaultDt and
## assignTicks/assignDefaultTicks keep the process of assigning dt to
## ticks and assigning ticks to objects separate. reinit() should be
## explicitly called by user just before running a simulation, and not
## when setting it up.
def updateTicks(tickDtMap):
    """Try to assign dt values to ticks.

    Parameters
    ----------
    tickDtMap: dict
    map from tick-no. to dt value. if it is empty, then default dt
    values are assigned to the ticks.

    """
    for tickNo, dt in list(tickDtMap.items()):
        if tickNo >= 0 and dt > 0.0:
            moose.setClock(tickNo, dt)
    if all([(v == 0) for v in list(tickDtMap.values())]):
        setDefaultDt()

def assignTicks(tickTargetMap):
    """
    Assign ticks to target elements.

    Parameters
    ----------
    tickTargetMap:
    Map from tick no. to target path and method. The path can be wildcard expression also.
    """
    if len(tickTargetMap) == 0:
        assignDefaultTicks()
    for tickNo, target in list(tickTargetMap.items()):
        if not isinstance(target, str):
            if len(target) == 1:
                moose.useClock(tickNo, target[0], 'process')
            elif len(target) == 2:
                moose.useClock(tickNo, target[0], target[1])
        else:
            moose.useClock(tickNo, target, 'process')

    # # This is a hack, we need saner way of scheduling
    # ticks = moose.vec('/clock/tick')
    # valid = []
    # for ii in range(ticks[0].localNumField):
    #     if ticks[ii].dt > 0:
    #         valid.append(ii)
    # if len(valid) == 0:
    #     assignDefaultTicks()

def setDefaultDt(elecdt=1e-5, chemdt=0.01, tabdt=1e-5, plotdt1=1.0, plotdt2=0.25e-3):
    """Setup the ticks with dt values.

    Parameters
    ----------

    elecdt: dt for ticks used in computing electrical biophysics, like
    neuronal compartments, ion channels, synapses, etc.

    chemdt: dt for chemical computations like enzymatic reactions.

    tabdt: dt for lookup tables

    plotdt1: dt for chemical kinetics plotting

    plotdt2: dt for electrical simulations

    """
    moose.setClock(0, elecdt)
    moose.setClock(1, elecdt)
    moose.setClock(2, elecdt)
    moose.setClock(3, elecdt)
    moose.setClock(4, chemdt)
    moose.setClock(5, chemdt)
    moose.setClock(6, tabdt)
    moose.setClock(7, tabdt)
    moose.setClock(8, plotdt1) # kinetics sim
    moose.setClock(9, plotdt2) # electrical sim

def assignDefaultTicks(modelRoot='/model', dataRoot='/data', solver='hsolve'):
    if isinstance(modelRoot, moose.melement) or isinstance(modelRoot, moose.vec):
        modelRoot = modelRoot.path
    if isinstance(dataRoot, moose.melement) or isinstance(dataRoot, moose.vec):
        dataRoot = dataRoot.path
    if solver != 'hsolve' or len(moose.wildcardFind('%s/##[ISA=HSolve]' % (modelRoot))) == 0:
        moose.useClock(0, '%s/##[ISA=Compartment]' % (modelRoot), 'init')
        moose.useClock(1, '%s/##[ISA=Compartment]'  % (modelRoot), 'process')
        moose.useClock(2, '%s/##[ISA=HHChannel]'  % (modelRoot), 'process')
        # moose.useClock(2, '%s/##[ISA=ChanBase]'  % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=IzhikevichNrn]' % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=GapJunction]' % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=HSolve]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=LeakyIaF]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=IntFire]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=SpikeGen]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=PulseGen]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=StimulusTable]'  % (modelRoot), 'process')
    moose.useClock(1, '%s/##[ISA=TimeTable]'  % (modelRoot), 'process')
    moose.useClock(2, '%s/##[ISA=HHChannel2D]'  % (modelRoot), 'process')
    moose.useClock(2, '%s/##[ISA=SynChan]'  % (modelRoot), 'process')
    moose.useClock(2, '%s/##[ISA=MgBlock]'  % (modelRoot), 'process')
    moose.useClock(3, '%s/##[ISA=CaConc]'  % (modelRoot), 'process')
    moose.useClock(3, '%s/##[ISA=Func]' % (modelRoot), 'process')
    # The voltage clamp circuit depends critically on the dt used for
    # computing soma Vm and need to be on a clock with dt=elecdt.
    moose.useClock(0, '%s/##[ISA=DiffAmp]'  % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=VClamp]' % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=PIDController]' % (modelRoot), 'process')
    moose.useClock(0, '%s/##[ISA=RC]' % (modelRoot), 'process')
    # Special case for kinetics models
    kinetics = moose.wildcardFind('%s/##[FIELD(name)=kinetics]' % modelRoot)
    if len(kinetics) > 0:
        # Do nothing for kinetics models - until multiple scheduling issue is fixed.
        moose.useClock(4, '%s/##[ISA!=PoolBase]' % (kinetics[0].path), 'process')
        moose.useClock(5, '%s/##[ISA==PoolBase]' % (kinetics[0].path), 'process')
        moose.useClock(18, '%s/##[ISA=Table2]' % (dataRoot), 'process')
    else:
        # input() function is called in Table. process() which gets
        # called at each timestep. When a message is connected
        # explicitly to input() dest field, it is driven by the sender
        # and process() adds garbage value to the vector. Hence not to
        # be scheduled.
        for tab in moose.wildcardFind('%s/##[ISA=Table]' % (dataRoot)):
            if len(tab.neighbors['input']) == 0:
                moose.useClock(9, tab.path, 'process')

def stepRun(simtime, steptime, verbose=True, logger=None):
    """Run the simulation in steps of `steptime` for `simtime`."""
    clock = moose.element('/clock')
    if verbose:
        msg = 'Starting simulation for %g' % (simtime)
        if logger is None:
            print(msg)
        else:
            logger.info(msg)
    ts = datetime.now()
    while clock.currentTime < simtime - steptime:
        ts1 = datetime.now()
        moose.start(steptime)
        te = datetime.now()
        td = te - ts1
        if verbose:
            msg = 'Simulated till %g. Left: %g. %g of simulation took: %g s' % (clock.currentTime, simtime - clock.currentTime, steptime, td.days * 86400 + td.seconds + 1e-6 * td.microseconds)
            if logger is None:
                print(msg)
            else:
                logger.info(msg)

    remaining = simtime - clock.currentTime
    if remaining > 0:
        if verbose:
            msg = 'Running the remaining %g.' % (remaining)
            if logger is None:
                print(msg)
            else:
                logger.info(msg)
        moose.start(remaining)
    te = datetime.now()
    td = te - ts
    dt = min([t for t in moose.element('/clock').dts if t > 0.0])
    if verbose:
        msg = 'Finished simulation of %g with minimum dt=%g in %g s' % (simtime, dt, td.days * 86400 + td.seconds + 1e-6 * td.microseconds)
        if logger is None:
            print(msg)
        else:
            logger.info(msg)



############# added by Aditya Gilra -- begin ################

def resetSim(simpaths, simdt, plotdt, simmethod='hsolve'):
    """ For each of the MOOSE paths in simpaths, this sets the clocks and finally resets MOOSE.
    If simmethod=='hsolve', it sets up hsolve-s for each Neuron under simpaths, and clocks for hsolve-s too. """
    print('Solver:', simmethod)
    moose.setClock(INITCLOCK, simdt)
    moose.setClock(ELECCLOCK, simdt) # The hsolve and ee methods use clock 1
    moose.setClock(CHANCLOCK, simdt) # hsolve uses clock 2 for mg_block, nmdachan and others.
    moose.setClock(POOLCLOCK, simdt) # Ca/ion pools & funcs use clock 3
    moose.setClock(STIMCLOCK, simdt) # Ca/ion pools & funcs use clock 3
    moose.setClock(PLOTCLOCK, plotdt) # for tables
    for simpath in simpaths:
        ## User can connect [qty]Out of an element to input of Table or
        ## requestOut of Table to get[qty] of the element.
        ## Scheduling the Table to a clock tick, will call process() of the Table
        ## which will send a requestOut and overwrite any value set by input(),
        ## thus adding garbage value to the vector. Hence schedule only if
        ## input message is not connected to the Table.
        for table in moose.wildcardFind(simpath+'/##[TYPE=Table]'):
            if len(table.neighbors['input']) == 0:
                moose.useClock(PLOTCLOCK, table.path, 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=PulseGen]', 'process')
        moose.useClock(STIMCLOCK, simpath+'/##[TYPE=DiffAmp]', 'process')
        moose.useClock(STIMCLOCK, simpath+'/##[TYPE=VClamp]', 'process')
        moose.useClock(STIMCLOCK, simpath+'/##[TYPE=PIDController]', 'process')
        moose.useClock(STIMCLOCK, simpath+'/##[TYPE=RC]', 'process')
        moose.useClock(STIMCLOCK, simpath+'/##[TYPE=TimeTable]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=LeakyIaF]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=IntFire]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=IzhikevichNrn]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=SpikeGen]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=Interpol]', 'process')
        moose.useClock(ELECCLOCK, simpath+'/##[TYPE=Interpol2D]', 'process')
        moose.useClock(CHANCLOCK, simpath+'/##[TYPE=HHChannel2D]', 'process')
        moose.useClock(CHANCLOCK, simpath+'/##[TYPE=SynChan]', 'process')
        ## If simmethod is not hsolve, set clocks for the biophysics,
        ## else just put a clock on the hsolve:
        ## hsolve takes care of the clocks for the biophysics
        if 'hsolve' not in simmethod.lower():
            print('Using exp euler')
            moose.useClock(INITCLOCK, simpath+'/##[TYPE=Compartment]', 'init')
            moose.useClock(ELECCLOCK, simpath+'/##[TYPE=Compartment]', 'process')
            moose.useClock(CHANCLOCK, simpath+'/##[TYPE=HHChannel]', 'process')
            moose.useClock(POOLCLOCK, simpath+'/##[TYPE=CaConc]', 'process')
            moose.useClock(POOLCLOCK, simpath+'/##[TYPE=Func]', 'process')
        else: # use hsolve, one hsolve for each Neuron
            print('Using hsolve')
            element = moose.Neutral(simpath)
            for childid in element.children:
                childobj = moose.Neutral(childid)
                classname = childobj.className
                if classname in ['Neuron']:
                    neuronpath = childobj.path
                    h = moose.HSolve( neuronpath+'/solve' )
                    h.dt = simdt
                    h.target = neuronpath
                    moose.useClock(INITCLOCK, h.path, 'process')
    moose.reinit()

def setupTable(name, obj, qtyname, tables_path=None, threshold=None, spikegen=None):
    """ Sets up a table with 'name' which stores 'qtyname' field from 'obj'.
    The table is created under tables_path if not None, else under obj.path . """
    if tables_path is None:
        tables_path = obj.path+'/data'
    ## in case tables_path does not exist, below wrapper will create it
    tables_path_obj = moose.Neutral(tables_path)
    qtyTable = moose.Table(tables_path_obj.path+'/'+name)
    ## stepMode no longer supported, connect to 'input'/'spike' message dest to record Vm/spiktimes
    # qtyTable.stepMode = TAB_BUF
    if spikegen is None:
        if threshold is None:
            ## below is wrong! reads qty twice every clock tick!
            #moose.connect( obj, qtyname+'Out', qtyTable, "input")
            ## this is the correct method
            moose.connect( qtyTable, "requestOut", obj, 'get'+qtyname)
        else:
            ## create new spikegen
            spikegen = moose.SpikeGen(tables_path_obj.path+'/'+name+'_spikegen')
            ## connect the compartment Vm to the spikegen
            moose.connect(obj,"VmOut",spikegen,"Vm")
            ## spikegens for different synapse_types can have different thresholds
            spikegen.threshold = threshold
            spikegen.edgeTriggered = 1 # This ensures that spike is generated only on leading edge.
    else:
        moose.connect(spikegen,'spikeOut',qtyTable,'input') ## spikeGen gives spiketimes
    return qtyTable

def connectSynapse(compartment, synname, gbar_factor):
    """
    Creates a synname synapse under compartment, sets Gbar*gbar_factor, and attaches to compartment.
    synname must be a synapse in /library of MOOSE.
    """
    synapseid = moose.copy(moose.SynChan('/library/'+synname),compartment,synname)
    synapse = moose.SynChan(synapseid)
    synapse.Gbar = synapse.Gbar*gbar_factor
    synapse_mgblock = moose.Mstring(synapse.path+'/mgblockStr')
    if synapse_mgblock.value=='True': # If NMDA synapse based on mgblock, connect to mgblock
        mgblock = moose.Mg_block(synapse.path+'/mgblock')
        compartment.connect("channel", mgblock, "channel")
    else:
        compartment.connect("channel", synapse, "channel")
    return synapse

def printNetTree():
    """ Prints all the cells under /, and recursive prints the cell tree for each cell. """
    root = moose.Neutral('/')
    for id in root.children: # all subelements of 'root'
        if moose.Neutral(id).className == 'Cell':
            cell = moose.Cell(id)
            print("-------------------- CELL : ",cell.name," ---------------------------")
            printCellTree(cell)

def printCellTree(cell):
    """
    Prints the tree under MOOSE object 'cell'.
    Assumes cells have all their compartments one level below,
    also there should be nothing other than compartments on level below.
    Apart from compartment properties and messages,
    it displays the same for subelements of compartments only one level below the compartments.
    Thus NMDA synapses' mgblock-s will be left out.

    FIXME: no lenght cound on compartment.
    """
    for compartmentid in cell.children: # compartments
        comp = moose.Compartment(compartmentid)
        print("  |-",comp.path, 'l=',comp.length, 'd=',comp.diameter, 'Rm=',comp.Rm, 'Ra=',comp.Ra, 'Cm=',comp.Cm, 'EM=',comp.Em)
        #for inmsg in comp.inMessages():
        #    print "    |---", inmsg
        #for outmsg in comp.outMessages():
        #    print "    |---", outmsg
        printRecursiveTree(compartmentid, level=2) # for channels and synapses and recursively lower levels

## Use printCellTree which calls this
def printRecursiveTree(elementid, level):
    """ Recursive helper function for printCellTree,
    specify depth/'level' to recurse and print subelements under MOOSE 'elementid'. """
    spacefill = '  '*level
    element = moose.Neutral(elementid)
    for childid in element.children:
        childobj = moose.Neutral(childid)
        classname = childobj.className
        if classname in ['SynChan','KinSynChan']:
            childobj = moose.SynChan(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'Gbar=',childobj.Gbar, 'numSynapses=', childobj.numSynapses)
            return # Have yet to figure out the children of SynChan, currently not going deeper
        elif classname in ['HHChannel', 'HHChannel2D']:
            childobj = moose.HHChannel(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'Gbar=',childobj.Gbar, 'Ek=',childobj.Ek)
        elif classname in ['CaConc']:
            childobj = moose.CaConc(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'thick=',childobj.thick, 'B=',childobj.B)
        elif classname in ['Mg_block']:
            childobj = moose.Mg_block(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'CMg',childobj.CMg, 'KMg_A',childobj.KMg_A, 'KMg_B',childobj.KMg_B)
        elif classname in ['SpikeGen']:
            childobj = moose.SpikeGen(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'threshold',childobj.threshold)
        elif classname in ['Func']:
            childobj = moose.Func(childid)
            print(spacefill+"|--", childobj.name, childobj.className, 'expr',childobj.expr)
        elif classname in ['Table']: # Table gives segfault if printRecursiveTree is called on it
            return # so go no deeper
        #for inmsg in childobj.inMessages():
        #    print spacefill+"  |---", inmsg
        #for outmsg in childobj.outMessages():
        #    print spacefill+"  |---", outmsg
        if len(childobj.children)>0:
            printRecursiveTree(childid, level+1)

def setup_vclamp(compartment, name, delay1, width1, level1, gain=0.5e-5):
    """
    Sets up a voltage clamp with 'name' on MOOSE 'compartment' object:
    adapted from squid.g in DEMOS (moose/genesis)
    Specify the 'delay1', 'width1' and 'level1' of the voltage to be applied to the compartment.
    Typically you need to adjust the PID 'gain'
    For perhaps the Davison 4-compartment mitral or the Davison granule:
    0.5e-5 optimal gain - too high 0.5e-4 drives it to oscillate at high frequency,
    too low 0.5e-6 makes it have an initial overshoot (due to Na channels?)
    Returns a MOOSE table with the PID output.
    """
    ## If /elec doesn't exists it creates /elec and returns a reference to it.
    ## If it does, it just returns its reference.
    moose.Neutral('/elec')
    pulsegen = moose.PulseGen('/elec/pulsegen'+name)
    vclamp = moose.DiffAmp('/elec/vclamp'+name)
    vclamp.saturation = 999.0
    vclamp.gain = 1.0
    lowpass = moose.RC('/elec/lowpass'+name)
    lowpass.R = 1.0
    lowpass.C = 50e-6 # 50 microseconds tau
    PID = moose.PIDController('/elec/PID'+name)
    PID.gain = gain
    PID.tau_i = 20e-6
    PID.tau_d = 5e-6
    PID.saturation = 999.0
    # All connections should be written as source.connect('',destination,'')
    pulsegen.connect('outputSrc',lowpass,'injectMsg')
    lowpass.connect('outputSrc',vclamp,'plusDest')
    vclamp.connect('outputSrc',PID,'commandDest')
    PID.connect('outputSrc',compartment,'injectMsg')
    compartment.connect('VmSrc',PID,'sensedDest')

    pulsegen.trigMode = 0 # free run
    pulsegen.baseLevel = -70e-3
    pulsegen.firstDelay = delay1
    pulsegen.firstWidth = width1
    pulsegen.firstLevel = level1
    pulsegen.secondDelay = 1e6
    pulsegen.secondLevel = -70e-3
    pulsegen.secondWidth = 0.0

    vclamp_I = moose.Table("/elec/vClampITable"+name)
    vclamp_I.stepMode = TAB_BUF #TAB_BUF: table acts as a buffer.
    vclamp_I.connect("inputRequest", PID, "output")
    vclamp_I.useClock(PLOTCLOCK)

    return vclamp_I

def setup_iclamp(compartment, name, delay1, width1, level1):
    """
    Sets up a current clamp with 'name' on MOOSE 'compartment' object:
    Specify the 'delay1', 'width1' and 'level1' of the current pulse to be applied to the compartment.
    Returns the MOOSE pulsegen that sends the current pulse.
    """
    ## If /elec doesn't exists it creates /elec and returns a reference to it.
    ## If it does, it just returns its reference.
    moose.Neutral('/elec')
    pulsegen = moose.PulseGen('/elec/pulsegen'+name)
    iclamp = moose.DiffAmp('/elec/iclamp'+name)
    iclamp.saturation = 1e6
    iclamp.gain = 1.0
    pulsegen.trigMode = 0 # free run
    pulsegen.baseLevel = 0.0
    pulsegen.firstDelay = delay1
    pulsegen.firstWidth = width1
    pulsegen.firstLevel = level1
    pulsegen.secondDelay = 1e6 # to avoid repeat
    pulsegen.secondLevel = 0.0
    pulsegen.secondWidth = 0.0
    pulsegen.connect('output',iclamp,'plusIn')
    iclamp.connect('output',compartment,'injectMsg')
    return pulsegen

def get_matching_children(parent, names):
    """ Returns non-recursive children of 'parent' MOOSE object
    with their names containing any of the strings in list 'names'. """
    matchlist = []
    for childID in parent.children:
        child = moose.Neutral(childID)
        for name in names:
            if name in child.name:
                matchlist.append(childID)
    return matchlist

def underscorize(path):
    """ Returns: / replaced by underscores in 'path'.
    But async13 branch has indices in the path like [0],
    so just replacing / by _ is not enough,
    should replace [ and ] also by _ """
    return path.replace('/','_').replace('[','-').replace(']','-')

def blockChannels(cell, channel_list):
    """
    Sets gmax to zero for channels of the 'cell' specified in 'channel_list'
    Substring matches in channel_list are allowed
    e.g. 'K' should block all K channels (ensure that you don't use capital K elsewhere in your channel name!)
    """
    for compartmentid in cell.children: # compartments
        comp = moose.Compartment(compartmentid)
        for childid in comp.children:
            child = moose.Neutral(childid)
            if child.className in ['HHChannel', 'HHChannel2D']:
                chan = moose.HHChannel(childid)
                for channame in channel_list:
                    if channame in chan.name:
                        chan.Gbar = 0.0

def get_child_Mstring(mooseobject,mstring):
    for child in mooseobject.children:
        if child.className=='Mstring' and child.name==mstring:
            child = moose.Mstring(child)
            return child
    return None

def connect_CaConc(compartment_list, temperature=None):
    """ Connect the Ca pools and channels within each of the compartments in compartment_list
     Ca channels should have a child Mstring named 'ion' with value set in MOOSE.
     Ca dependent channels like KCa should have a child Mstring called 'ionDependency' with value set in MOOSE.
     Call this only after instantiating cell so that all channels and pools have been created. """
    for compartment in compartment_list:
        caconc = None
        for child in compartment.children:
            neutralwrap = moose.Neutral(child)
            if neutralwrap.className == 'CaConc':
                caconc = moose.CaConc(child)
                break
        if caconc is not None:
            child = get_child_Mstring(caconc,'phi')
            if child is not None:
                caconc.B = float(child.value) # B = phi by definition -- see neuroml 1.8.1 defn
            else:
                ## B has to be set for caconc based on thickness of Ca shell and compartment l and dia,
                ## OR based on the Mstring phi under CaConc path.
                ## I am using a translation from Neuron for mitral cell, hence this method.
                ## In Genesis, gmax / (surfacearea*thick) is set as value of B!
                caconc.B = 1 / (2*FARADAY) / \
                    (math.pi*compartment.diameter*compartment.length * caconc.thick)
            for child in compartment.children:
                neutralwrap = moose.Neutral(child)
                if neutralwrap.className == 'HHChannel':
                    channel = moose.HHChannel(child)
                    ## If child Mstring 'ion' is present and is Ca, connect channel current to caconc
                    for childid in channel.children:
                        # in async13, gates which have not been created still 'exist'
                        # i.e. show up as a child, but cannot be wrapped.
                        try:
                            child = moose.element(childid)
                            if child.className=='Mstring':
                                child = moose.Mstring(child)
                                if child.name=='ion':
                                    if child.value in ['Ca','ca']:
                                        moose.connect(channel,'IkOut',caconc,'current')
                                        #print 'Connected IkOut of',channel.path,'to current of',caconc.path
                                ## temperature is used only by Nernst part here...
                                if child.name=='nernst_str':
                                    nernst = moose.Nernst(channel.path+'/nernst')
                                    nernst_params = child.value.split(',')
                                    nernst.Cout = float(nernst_params[0])
                                    nernst.valence = float(nernst_params[1])
                                    nernst.Temperature = temperature
                                    moose.connect(nernst,'Eout',channel,'setEk')
                                    moose.connect(caconc,'concOut',nernst,'ci')
                                    #print 'Connected Nernst',nernst.path
                        except TypeError:
                            pass

                if neutralwrap.className == 'HHChannel2D':
                    channel = moose.HHChannel2D(child)
                    ## If child Mstring 'ionDependency' is present, connect caconc Ca conc to channel
                    for childid in channel.children:
                        # in async13, gates which have not been created still 'exist'
                        # i.e. show up as a child, but cannot be wrapped.
                        try:
                            child = moose.element(childid)
                            if child.className=='Mstring' and child.name=='ionDependency':
                                child = moose.Mstring(child)
                                if child.value in ['Ca','ca']:
                                    moose.connect(caconc,'concOut',channel,'concen')
                                    #print 'Connected concOut of',caconc.path,'to concen of',channel.path
                        except TypeError:
                            pass
