"""
Integrator interfaces.

Robert A McDougal, 2011
"""


import os
import os.path
import subprocess
from subprocess import Popen, PIPE
import platform
import time
import threading
import random
import StringIO
from errors import SimulationFailedError
from tokenizer import tokens_to_string, Variable, default_renderer, parse_string
from tokenizer import exponentiation, rpn_renderer, c_renderer, matlab_renderer
from tokenizer import simplify_expression


if platform.system() in ('Linux', 'Darwin'):
    null_device = '/dev/null'
elif platform.system() == 'Windows':
    null_device = 'NUL'
else:
    null_device = None


def _volumes(sim, renderer=default_renderer, divide_volume=True):
    """Return the volumes of each region."""
    volumes = {'': str(sim._total_volume)}
    for i in len(sim._neurons):
        volumes[str(i)] = _effective_volume(sim, i, renderer=renderer, divide_volume=divide_volume)
    return volumes
    
def _initial_value(dynamics, var, renderer=default_renderer, divide_volume=True):
    """Return the initial value of a parameter.
    
    For concentrations, multiplies by volume to get mass."""
    if var not in dynamics.concentrations:
        return _filter_initial(tokens_to_string(dynamics.algebraic[var], renderer=renderer))
    if dynamics.index is None:
        effective_volume = dynamics.simulation._total_volume
        for dyn in dynamics.simulation._neurons:
            if dyn._container is None:
                effective_volume -= float(tokens_to_string(dyn.algebraic['volume'], renderer=renderer))
    else:
        effective_volume = float(tokens_to_string(dynamics.algebraic['volume'], renderer=renderer))
        for i in reduce(lambda x, y: x + y, dynamics._compartments.values(), []):
            effective_volume -= float(tokens_to_string(dynamics.simulation._neurons[i].algebraic['volume'], renderer=renderer))
    return str(effective_volume * float(tokens_to_string(dynamics.algebraic[var], renderer=renderer)))

def _effective_volume(sim, id, renderer=default_renderer, divide_volume=True):
    """Return the effective volume of region id.
    
    Effective volume is the total volume minus the volume of any compartments.
    """
    volume_token = parse_string('volume')
    if id != '':
        id = int(id)
        neuron = sim._neurons[id]
        if 'volume' in neuron.diff_list:
            volume = ['volume' + str(id)]
        else:
            volume = [_tokens_to_string(sim._neurons[id], id, neuron.algebraic['volume'],
                      renderer=renderer, allow_external=False, divide_volume=divide_volume)]
        for i in reduce(lambda x, y: x + y, neuron._compartments.values(), []):
            volume.append('(' + _tokens_to_string(sim._neurons[i], i, volume_token, renderer=renderer, allow_external=False, divide_volume=divide_volume) + ')')
    else:
        volume = [str(sim._total_volume)]
        for i, n in enumerate(sim._neurons):
            if n._container is None:
                volume.append('(' + _tokens_to_string(n, i, volume_token, renderer=renderer, allow_external=False, divide_volume=divide_volume) + ')')
    return '(' + simplify_expression('-'.join(volume)) + ')'


def _compute_fluxes(sim, renderer=default_renderer, divide_volume=True, allow_external=True):
    """Compute all the fluxes."""
    fluxes = {'': {}}
    for i, n in enumerate(sim._neurons):
        str_i = str(i)
        fluxes[str_i] = {}
        for var in n.diff_list:
            fluxes[str_i][var] = ''
    
    for dyn in sim._shared_dynamics:
        for var in dyn.diff_list:
            fluxes[''][var] = ''
    
    for i, fs in enumerate(sim._fluxes):
        compartment = sim._neurons[i]
        container_id = str(compartment._container)
        str_i = str(i)
        if container_id == 'None':
            container_id = ''
        for f in fs:
            flux = _tokens_to_string(compartment, i, parse_string(f['rate']), renderer=renderer, divide_volume=divide_volume, allow_external=allow_external)
            fluxes[str_i][f['var']] += '-(' + flux + ')'
            fluxes[container_id][f['external_var']] += '+(' + flux + ')'

    return fluxes

def _write_dynamics_xpp(file, dynamics, fluxes, syn_inputs=[], id=''):
    """Writes a set of dynamics to file.
    
    Appends id to all of the state variables.
    """
    id = str(id)
    isyn = '+'.join([_tokens_to_string(dynamics, id, syn) for syn in syn_inputs])
    if isyn != '':
        file.write('isyn%s=%s\n' % (id, isyn))
    for var in dynamics.diff_list:
        rhs = _tokens_to_string(dynamics, id, dynamics.differential[var])
        if var == dynamics._potential_var and isyn != '':
            rhs += '-isyn' + str(id)
        if var in dynamics.concentrations and rhs != '0':
            rhs = '(' + rhs + ')*' + _effective_volume(dynamics.simulation, id)
        rhs += fluxes[id][var]
        file.write("%s%s'=%s\n" % (var, id, rhs))
        if var not in dynamics.algebraic:
            raise SimulationFailedError(0, 'Initial condition not defined for ' + var + '.')
        file.write('i %s%s=%s\n' % (var, id, _initial_value(dynamics, var)))

def _write_dynamics_c(file, dynamics, j, fluxes, syn_inputs=[], id=''):
    """Writes a set of dynamics to file.
    
    Appends id to all of the state variables.
    """
    id = str(id)
    isyn = '+'.join([_tokens_to_string(dynamics, id, syn) for syn in syn_inputs])
    for i, var in enumerate(dynamics.diff_list):
        rhs = _tokens_to_string(dynamics, id, dynamics.differential[var], c_renderer, divide_volume=False)
        if var == dynamics._potential_var and isyn != '':
            if rhs != '0':
                rhs += '-(' + isyn +')'
            else:
                rhs = '-(' + isyn + ')'
        if var in dynamics.concentrations and rhs != '0':
            rhs = '(' + rhs + ')*' + _effective_volume(dynamics.simulation, id, c_renderer, divide_volume=False)
        rhs += fluxes[id][var]
        file.write('    _deriv[%d] = %s;\n' % (j + i, rhs))

def _write_dynamics_matlab(file, dynamics, j, fluxes, syn_inputs=[], id=''):
    """Writes a set of dynamics to file.
    
    Appends id to all of the state variables.
    """
    id = str(id)
    isyn = '+'.join([_tokens_to_string(dynamics, id, syn) for syn in syn_inputs])
    for i, var in enumerate(dynamics.diff_list):
        rhs = _tokens_to_string(dynamics, id, dynamics.differential[var], matlab_renderer)
        if var == dynamics._potential_var and isyn != '':
            rhs += '-(' + isyn +')'
        if var in dynamics.concentrations and rhs != '0':
            rhs = '(' + rhs + ')*' + _effective_volume(dynamics.simulation, id, matlab_renderer)
        rhs += fluxes[id][var]
        file.write('    yprime(%d) = %s;\n' % (j + i + 1, rhs))

        
        
def _delete_file(filename, num_attempts=10, delay=3):
    """Delete the file named filename.
    
    Tries up to num_attempts times, with delay between tries measured in sec.
    Returns True on success, else false.
    """
    # sometimes os.remove seems not to work, so try just using rm
    if platform.system() in ('Linux', 'Darwin'):
        rv = os.system('rm -f ' + filename)
        if rv == 0:
            return True
        
    # try repeatedly in case file in use (Windows issue)
    for i in xrange(num_attempts):
        try:
            os.remove(filename)
        except:
            pass
        if os.path.exists(filename):
            time.sleep(delay)
        else:
            return True
    return False

def _read_data(filename, num_attempts=10, delay=1, t=None, y=None):
    """Read space separated data in from filename."""
    if t is None:
        t = []
    if y is None:
        y = []
    success = False
    for i in range(num_attempts):
        try:
            file = open(filename)
            success = True
            break
        except:
            time.sleep(delay)
    if not success:
        raise Exception('Unable to open ' + filename)
    while True:
        line = file.readline()
        if len(line) == 0:
            break
        data = [float(x) for x in line.strip().split(' ')]
        t.append(data[0])
        y.append(data[1 : ])
    file.close()
    return t, y


class _AsyncReadData(threading.Thread):
    """Asynchronous reading of data.
    
    Useful for named pipes.
    """
    def __init__(self, filename, t, y, num_attempts=10, delay=1):
        threading.Thread.__init__(self)
        self._filename = filename
        self._t = t
        self._y = y
        self._num_attempts = num_attempts
        self._delay = delay
    
    def run(self):
        _read_data(self._filename, self._num_attempts, self._delay, self._t, self._y)
    
    def shutdown(self, delay=.5):
        """Shuts down the reader if it is still waiting on input by passing it garbage data."""
        valid_data = True
        if not(self._t):
            time.sleep(delay)
            if not(self._t):
                file = open(self._filename, 'w')
                file.write('0 0\n')
                file.close()
                valid_data = False
        self.join()
        if not(valid_data):
            self._t = []
            self._y = []
        return valid_data


def _filter_initial(s):
    if s[0 : 2] == '(-' and s[-1 : ] == ')':
        return '-' + s[2 : -1]
    else:
        return s
        
def _tokens_to_string(dyn, id, tokens, renderer=default_renderer, allow_external=True, divide_volume=True):
    """
        divide_volume -- should we divide concentration variables by the volume
        
        Note:
        
        If divide_volume is False, one must do the division in advance.
    """
    print 'in _tokens_to_string'
    substitutions = dict(dyn.algebraic)
    if allow_external and 'volume' in substitutions:
        substitutions['effective_volume'] = [_effective_volume(dyn.simulation, id, renderer, divide_volume=divide_volume)]
    if id != '':
        for var in dyn.differential:
            if var in dyn.concentrations and divide_volume:
                print 'calling 268'
                substitutions[var] = ['((' + var + str(id) + ')/' + _effective_volume(dyn.simulation, id, renderer) + ')']
            else:
                if var in dyn.concentrations:
                    print 'eh?'
                substitutions[var] = [var + str(id)]
    else:
        for d in dyn.simulation._shared_dynamics:
            for var in d.diff_list:
                if var in substitutions:
                    substitutions.pop(var)
    if id != '' and allow_external:
        if dyn._container is not None:
            container_id = dyn._container
            container = dyn.simulation._neurons[container_id]
            for var in container.differential:
                external_var = var + str(container_id)
                if var in container.concentrations and divide_volume:
                    external_var = '((' + external_var + ')/' + _effective_volume(dyn.simulation, container_id, renderer) + ')'
                substitutions['external_' + var] = [external_var]
        else:
            sim = dyn.simulation
            for i, shared_dyn in enumerate(sim._shared_dynamics):
                for var in shared_dyn.differential:
                    external_var = var
                    if var in shared_dyn.concentrations and divide_volume:
                        external_var = '((' + external_var + ')/' + _effective_volume(sim, None, renderer) + ')'
                    substitutions['external_' + var] = [external_var]
    result = tokens_to_string(tokens, dyn.functions, substitutions, renderer)
    return result


def _write_noise_equations_c(file, noise, start_t, run_t, dt=10):
    if len(noise) > 1:

        t_intervals = int(run_t / dt) + 1
    
        for i in range(t_intervals):
            t = start_t + dt * i
            file.write('    if ((t >= %g) && (t < %g)) {\n' % (t, t + dt))
            for j in range(1, len(noise)):
                file.write('        noise%d = 0' % j)
                noise_data = noise[j].get_times_and_tokens(t, t + dt)
                for spike in noise_data:
                    spike_t = spike['t']
                    spike_str = tokens_to_string(spike['tokens'], {}, {}, c_renderer)
                    if spike_t < t:
                        file.write(' + ' + spike_str)
                    else:
                        file.write(' + (t > %g) * (%s)' % (spike_t, spike_str))
                file.write(';\n')
            file.write('    }\n')
 
 
 
def _get_data_xpp(ode, base_filename, end_time, num_state_var, num_attempts=3, delay=1):
    working_directory = os.path.dirname(base_filename)
    
    if not working_directory:
        working_directory = '.'

    if platform.system() in ['Darwin', 'Linux']:
        data_file = '/dev/stdout'
        xpp = Popen(['xppaut','-silent', '/dev/stdin'], cwd=working_directory, stdin=PIPE, stdout=PIPE)
        f = xpp.stdin
    else:
        data_file = base_filename + '.dat'
        f = open(base_filename + '.ode', 'w')

    f.write(ode + '\n')
    f.write('@ output=%s\n' % data_file)
    f.write('done\n')
    f.close()
    
    if platform.system() in ['Darwin', 'Linux']:
        t = []
        y = []
        while True:
            line = xpp.stdout.readline()
            if not line:
                break
            try:
                line = [float(v) for v in line.strip().split(' ')]
                if len(line) == num_state_var + 1:
                    t.append(line[0])
                    y.append(line[1 : ])
            except:
                pass
        rv = xpp.wait()
        if rv:
            raise SimulationFailedError(0, 'Could not run xppaut.')
    else:
        if null_device:
            garbage_file = null_device
        else:
            garbage_file = base_filename + '.garbage'
       
        if working_directory:
            set_dir = 'cd %s && ' % working_directory
        else:
            set_dir = ''
        
        for i in xrange(num_attempts):
            rv = os.system(set_dir + 'xppaut -silent ' + base_filename + '.ode > ' + garbage_file)
            if not rv:
                break
            time.sleep(delay)
        
        if rv:
            raise SimulationFailedError(0, 'Could not run xppaut.')

        t, y = _read_data(data_file)
        
    if not(t):
        raise SimulationFailedError(0)
    
    if max(t) < end_time:
        raise SimulationFailedError(max(t))
    
    return t, y        

def _get_data(command, data_file, end_time, num_attempts=3, delay=1):
    """Runs a command that generates a .dat file, loads the data, validates.
    
    Uses named pipes if supported by the operating system and file system,
    else uses files to transfer data.
    """
    # create a named pipe, assign a thread to receive data as it is generated
    try:
        os.mkfifo(data_file)
    except:
        use_pipes = False
    else:
        use_pipes = True
        t = []
        y = []
        reader_thread = _AsyncReadData(data_file, t, y)
        reader_thread.setDaemon(True)
        reader_thread.start()

    rv = None
    for i in range(num_attempts):
        rv = os.system(command)
        if rv == 0:
            break
        time.sleep(delay)
    
    if rv != 0:
        raise SimulationFailedError(0, 'Could not run external integrator.')

    if use_pipes:
        reader_thread.shutdown()
    else:
        t, y = _read_data(data_file)
        
    if not(t):
        raise SimulationFailedError(0)
    
    if max(t) < end_time:
        raise SimulationFailedError(max(t))
    
    return t, y

        
def _write_noise_equations_matlab(file, noise, start_t, run_t, dt=10):
    if len(noise) > 1:
        t_intervals = int(run_t / dt) + 1
        
        for i in range(t_intervals):
            t = start_t + dt * i
            file.write('    if ((t >= %g) && (t < %g))\n' % (t, t + dt))
            for j in range(1, len(noise)):
                file.write('        noise%d = 0' % j)
                noise_data = noise[j].get_times_and_tokens(t, t + dt)
                for spike in noise_data:
                    spike_t = spike['t']
                    spike_str = tokens_to_string(spike['tokens'], {}, {}, matlab_renderer)
                    if spike_t < t:
                        file.write(' + ' + spike_str)
                    else:
                        file.write(' + (t > %g) * (%s)' % (spike_t, spike_str))
                file.write(';\n')
            file.write('    end\n')



# check if we can access a ram disk
default_dir = ''
if platform.system() == 'Linux':
    if os.system('touch /dev/shm/snnet') == 0:
        default_dir = '/dev/shm/'
        _delete_file('/dev/shm/snnet')
    


class Integrator:
    """Integrator parent class. All integrators are derived from here."""
        
    def _get_base_filename(self, simulation):
        """Get a base filename for temporary files.
        
        Applies simulation suffix and any _temp_dir_handler.
        """
        base_filename = 'tmp_snnet_%d_%d%s' % (os.getpid(), random.randint(0, 10000), simulation._suffix)
        if self._temp_dir_handler is not None:
            if isinstance(self._temp_dir_handler, str):
                temp_dir = self._temp_dir_handler
            else:
                temp_dir = self._temp_dir_handler()
            if temp_dir is not None:
                base_filename = os.path.join(temp_dir, base_filename)
        return base_filename
    
    def _register_temp_files(self, base_filename, extensions):
        """Record the various temp_files being used."""
        for e in extensions:
            self._temp_files.append(base_filename + e)

    def _delete_files(self, filename, extensions):
        """Delete the temporary files"""
        for e in extensions:
            if _delete_file(filename + e):
                try:
                    self._temp_files.remove(filename + e)     
                except ValueError:
                    pass
        

class DefaultTempDirHandler(object):
    """Returns the default temp directory."""
    def __init__(self):
        self.__name__ = 'DefaultTempDirHandler()'
    def __call__(self):
        return default_dir
    def __repr__(self):
        return self.__name__
    def __str__(self):
        return repr(self)
    

default_temp_dir_handler = DefaultTempDirHandler()
        


# TODO: Modify to start integration at the correct time t.
#       This is only a problem for dynamics that explicitly depend on time.
class XppIntegrator(Integrator):
    def __init__(self, temp_dir_handler=default_temp_dir_handler, dt=.05,
                 method='qualrk', bounds=10000):
        """Initialize the xppaut integrator interface.
        
        temp_dir_handler = a function called at the beginning of each
                           simulation, returns a string with a temp directory
                           name, should also do any cleanup
       
       Integration parameters:
       
        dt = integration step size
        method = xppaut integration method
        bounds = xppaut bounds (integration stops when exceeded)
        """
        self._temp_dir_handler = temp_dir_handler
        self.dt = dt
        self.method = method
        self.bounds = bounds
    

    def __repr__(self):
        """Return a string representing the contents."""
        return "XppIntegrator(temp_dir_handler=%s, dt=%g, method='%s', bounds=%g)" % (self._temp_dir_handler.__name__, self.dt, self.method, self.bounds)
    
    def integrate(self, simulation, run_time, num_attempts=10, delay=1, temp_files=[], make_local_copy=False):
        self._temp_files = temp_files
        
        if simulation.has_synaptic_delays():
            raise SimulationFailedError(0, 'XppIntegrator does not support synaptic delays.')

        if len(simulation._noise) > 1:
            raise SimulationFailedError(0, 'XppIntegrator does not support noise.')

        base_filename = self._get_base_filename(simulation)

        extensions = ['.dat', '.ode']
        if not null_device:
            extensions.append('garbage')
        self._register_temp_files(base_filename, extensions)   

        ode, num_state_var = self._write_file(base_filename, simulation, run_time, make_local_copy)

        # run xpp
        t, y = _get_data_xpp(ode, base_filename, run_time + simulation.t, num_state_var, num_attempts, delay)
        
        # remove the temporary files
        if not make_local_copy:
            self._delete_files(base_filename, extensions)
        
        return t, y
    
    def _write_file(self, filename, simulation, run_time, make_local_copy):
        connections = simulation._connections
        
        ode = StringIO.StringIO()
        
        fluxes = _compute_fluxes(simulation)
        
        num_state_var = 0
        for i, neuron in enumerate(simulation._neurons):
            num_state_var += len(neuron.diff_list)
            
            # determine distinct synapse rules
            synapses = []
            presynaptics = []
            for syn in simulation._connections[i]:
                syn_rule = syn['rule']
                if syn_rule in synapses:
                    presynaptics[synapses.index(syn_rule)].append(syn['pre_id'])
                else:
                    synapses.append(syn_rule)
                    presynaptics.append([syn['pre_id']])
            
            # compute synaptic inputs
            syn_inputs = []
            for syn, presyn in zip(synapses, presynaptics):
                syn_inputs.append(syn.get_tokens(presyn, i, neuron._potential_var))
            
            _write_dynamics_xpp(ode, neuron, fluxes, syn_inputs, i)
        
        # write the global dynamics
        for d in simulation._shared_dynamics:
            num_state_var += len(d.diff_list)
            _write_dynamics_xpp(ode, d, fluxes)
            
        ode.write('@ total=%g, t0=%g\n' % (run_time, simulation.t))
        ode.write('@ dt=%g, bounds=%g, meth=%s, maxstor=%d\n' % (self.dt, self.bounds, self.method, int(run_time * 1.1 / self.dt) + 1))
        ode_string = ode.getvalue()
        ode.close()
        
        if make_local_copy:
            f = open(os.path.split(filename)[1] + '.ode', 'w')
            f.write(ode_string)
            f.close()
        
        return ode_string, num_state_var

        
        


class GccXppIntegrator(Integrator):
    def __init__(self, temp_dir_handler=default_temp_dir_handler, dt=.05, method='qualrk',
                 bounds=10000, gcc_flags='-m32', delay=0):
        """Initialize the gcc + xppaut integrator interface.
        
        Writes equations to a c file, compiles it, then integrates via xpp.
        
        temp_dir_handler = a function called at the beginning of each
                           simulation, returns a string with a temp directory
                           name, should also do any cleanup

        gcc_flags = flags to pass to compiler
                    (send -m32 to compile to 32 bit code)
               
        Integration parameters:
       
        dt = integration step size
        method = xppaut integration method
        bounds = xppaut bounds (integration stops when exceeded)                


        Assumes xppaut is compiled as a 32 bit application.
        """
        self._temp_dir_handler = temp_dir_handler
        self.dt = dt
        self.method = method
        self.bounds = bounds
        self.gcc_flags = gcc_flags
        self.delay = delay
    
    def __repr__(self):
        """Return a string representing the contents."""
        return "GccXppIntegrator(temp_dir_handler=%s, dt=%g, method='%s', bounds=%g, gcc_flags='%s')" % (self._temp_dir_handler.__name__, self.dt, self.method, self.bounds, self.gcc_flags)


    def integrate(self, simulation, run_time, temp_files=[], make_local_copy=False):
        self._temp_files = temp_files
        if platform.system() not in ['Darwin', 'Linux']:
            raise SimulationFailedError(0, 'GccXppIntegrator currently only works on OSX and linux. '
                                           'Figure out how to run the XPP example tstdll.ode on your system and tell me to fix.')

        base_filename = self._get_base_filename(simulation)

        extensions = ['.dat', '.c', '.so', '.ode', '.garbage']
        self._register_temp_files(base_filename, extensions)   

        ode, num_state_var = self._write_file(base_filename, simulation, run_time, make_local_copy)
        # compile the c file
        rv = None
        if platform.system() == 'Darwin':
            rv = os.system('gcc ' + self.gcc_flags +' -dynamiclib -fPIC -o ' + base_filename + '.so ' + base_filename + '.c')
        elif platform.system() == 'Linux':
            rv = os.system('gcc ' + self.gcc_flags +' -shared -fpic -o ' + base_filename + '.so ' + base_filename + '.c')
        
        if rv is not 0:
            raise SimulationFailedError(0, 'gcc failed.')
        
        # run xpp
        t, y = _get_data_xpp(ode, base_filename, run_time + simulation.t, num_state_var)
        
        # remove the temporary files
        if not make_local_copy:
            self._delete_files(base_filename, extensions)

        
        return t, y
    
    
    def _write_file(self, filename, simulation, run_time, make_local_copy):
        connections = simulation._connections
        if make_local_copy:
            filenames = [filename, os.path.split(filename)[1]]
        else:
            filenames = [filename]

        for filename in filenames:
            file = open(filename + '.c', 'w')
            file.write('#include <math.h>\n'
                       '#include <stdlib.h>\n'
                       '\n'
                       'double heav(double v) {\n'
                       '    return (v > 0);\n'
                       '}\n'
                       'double sign(double x) {\n'
                       '    if (x > 0) return(1);\n'
                       '    if (x < 0) return(-1);\n'
                       '    return(0);\n'
                       '}\n'
                       'double min(double x, double y) {\n'
                       '    return ((x < y) ? x : y);\n'
                       '}\n'
                       'double max(double x, double y) {\n'
                       '    return ((x > y) ? x : y);\n'
                       '}\n'
                       '\n'
                       'simulation(double* _a, double* _b, int _c, int _d, double* _var, double* _e) {\n'
                       '    double* _rhs = _var + 1;\n'
                       '    double t = _var[0];\n')
            # declare the variables
            j = 0;
            
            fluxes = _compute_fluxes(simulation, renderer=c_renderer, divide_volume=False)
            
            declaration_normal = []
            declaration_volume = []
            # declare neuron state variables
            for i, neuron in enumerate(simulation._neurons):
                for var in neuron.diff_list:
                    if var in neuron.concentrations:
                        declaration = '    double %s%d = _rhs[%d]/%s;\n' % (var, i, j, _effective_volume(simulation, i, c_renderer, divide_volume=False))
                    else:
                        declaration = '    double %s%d = _rhs[%d];\n' % (var, i, j)
                    j += 1
                    if var == 'volume':
                        declaration_volume.append(declaration)
                    else:
                        declaration_normal.append(declaration)
            
            # declare shared dynamics state variables
            for i, dyn in enumerate(simulation._shared_dynamics):
                for var in dyn.diff_list:
                    if var in dyn.concentrations:
                        declaration = '    double %s = _rhs[%d]/%s;\n' % (var, j, _effective_volume(simulation, '', c_renderer, divide_volume=False))
                    else:
                        declaration = '    double %s = _rhs[%d];\n' % (var, j)
                    j += 1
                    if var == 'volume':
                        declaration_volume.append(declaration)
                    else:
                        declaration_normal.append(declaration)

            for declaration in declaration_volume + declaration_normal:
                file.write(declaration)
            
            file.write('    double* _deriv = _var + %d;\n' % (j + 1))

            # declare delay variables
            if simulation.has_synaptic_delays():
                file.write('    double* _delay = _var + %d;\n' % (2 * j + 1))
            delay_count = 0
            used_ids = []
            for syns in simulation._connections:
                for syn in syns:
                    if syn['rule'].has_delay():
                        syn_id = syn['syn_id']
                        if syn_id not in used_ids:
                            file.write('    double delay%d = _delay[%d];\n' % (syn_id, delay_count))
                            delay_count += 1
                            used_ids.append(syn_id)

            
            # noise sources
            if len(simulation._noise) > 1:
                file.write('    double ' + ', '.join(['noise%d' % i for i in range(1, len(simulation._noise))]) + ';\n')
            
            # write the noise equations
            _write_noise_equations_c(file, simulation._noise, simulation.t, run_time)

            j = 0
            for i, neuron in enumerate(simulation._neurons):
                # determine distinct synapse rules
                synapses = []
                presynaptics = []
                syn_ids = []
                for syn in simulation._connections[i]:
                    syn_rule = syn['rule']
                    if syn_rule in synapses:
                        rule_index = synapses.index(syn_rule)
                        presynaptics[rule_index].append(syn['pre_id'])
                        syn_ids[rule_index].append(syn['syn_id'])
                    else:
                        synapses.append(syn_rule)
                        presynaptics.append([syn['pre_id']])
                        syn_ids.append([syn['syn_id']])
                
                # compute synaptic inputs
                syn_inputs = []
                for syn, presyn, syn_id in zip(synapses, presynaptics, syn_ids):
                    syn_inputs.append(syn.get_tokens(presyn, i, neuron._potential_var, syn_id))
                
                _write_dynamics_c(file, neuron, j, fluxes, syn_inputs, i)
                j += len(neuron.diff_list)
            
            # write the global dynamics
            for d in simulation._shared_dynamics:
                _write_dynamics_c(file, d, j, fluxes)
                j += len(d.diff_list)
            file.write('}\n')
            file.close()
            
        # write the corresponding ode file
        file = StringIO.StringIO()
        file.write("y[0..%d]'=yprime[j]\n" % (j - 1))
        file.write('yprime[0..%d]=0\n' % (j - 1))
        
        # write delay variables
        used_ids = []
        for id, syns in enumerate(simulation._connections):
            for syn in syns:
                if syn['rule'].has_delay():
                    presynaptic_id = syn['pre_id']
                    if presynaptic_id >= 0:
                        var = syn['rule']._presynaptic_var
                    else:
                        raise Exception('delay of noise not yet implemented')
                    syn_id = syn['syn_id']
                    if syn_id not in used_ids:
                        file.write('delay%d=delay(y%d, %s)\n' % (syn_id, simulation.index(var, presynaptic_id), _tokens_to_string(simulation._neurons[presynaptic_id], presynaptic_id, parse_string(syn['rule']._delay), allow_external=False)))
                        used_ids.append(syn_id)

        # initial conditions
        j = 0
        for neuron in simulation._neurons:
            for var in neuron.diff_list:
                file.write('i y%d=%s\n' % (j, _initial_value(neuron, var, renderer=c_renderer, divide_volume=False)))
                j += 1
        for dynamics in simulation._shared_dynamics:
            for var in dynamics.diff_list:
                file.write('i y%d=%s\n' % (j, _initial_value(dynamics, var, renderer=c_renderer, divide_volume=False)))
                j += 1

        # time offset from real simulation time (necessary for noise)
        file.write('p a=%g\n' % simulation.t)
        file.write('export {a}\n')
                    
        file.write('@ total=%g, t0=%g\n' % (run_time, simulation.t))
        file.write('@ dt=%g, bounds=%g, meth=%s, maxstor=%d\n' % (self.dt, self.bounds, self.method, int(run_time * 1.1 / self.dt) + 1))
        file.write('@ dll_lib=' + os.path.basename(filename) + '.so\n')
        file.write('@ dll_fun=simulation\n')
        file.write('@ nout=1\n')
        file.write('@ delay=%g\n' % self.delay)
        ode = file.getvalue()
        file.close()
        
        if make_local_copy:
            f = open(os.path.split(filename)[1] + '.ode', 'w')
            f.write(ode)
            f.close()
        return ode, j






class MatlabMexIntegrator(Integrator):
    def __init__(self, temp_dir_handler=default_temp_dir_handler, method='ode23',
                 rel_tol=1e-3, max_step=.05):
        """Initialize the matlab (with mex) integrator interface.
        
        Writes equations to a c file, compiles it, then integrates via matlab.
        
        Parameters:

            temp_dir_handler -- a function called at the beginning of each
                                simulation, returns a string with a temp directory
                                name, should also do any cleanup

            method -- MATLAB integration method (ode45, ode23, etc... could be
                      user defined)

            rel_tol -- passed to MATLAB RelTol
            
            max_step -- passed to MATLAB MaxStep
        """
        self._temp_dir_handler = temp_dir_handler
        self.method = method
        self.rel_tol = rel_tol
        self.max_step = max_step

    
    def __repr__(self):
        """Return a string representing the contents."""
        return "MatlabMexIntegrator(temp_dir_handler=%s, method='%s', rel_tol=%g, max_step=%g)" % (self._temp_dir_handler.__name__, self.method, self.rel_tol, self.max_step)
 
    
    def integrate(self, simulation, run_time, temp_files=[], make_local_copy=False):
        self._temp_files = temp_files
        if simulation.has_synaptic_delays():
            raise SimulationFailedError(0, 'MatlabMexIntegrator does not support synaptic delays.')
        
        base_filename = self._get_base_filename(simulation)
        extensions = ['.dat', '_rhs.c', '.m', '.mex*']
        if not null_device:
            extensions.append('.garbage')
        self._register_temp_files(base_filename, extensions)   
        self._write_file(base_filename, simulation, run_time, make_local_copy)
        
        # run matlab
        if null_device:
            garbage_file = null_device
        else:
            garbage_file = base_filename + '.garbage'
        base_dir = os.path.dirname(base_filename)
        if base_dir == '':
            set_correct_path = ''
        else:
            set_correct_path = 'cd ' + base_dir + '; '

        t, y = _get_data(set_correct_path + 'matlab -nodisplay < ' + base_filename + '.m > ' + garbage_file,
                         base_filename + '.dat', run_time + simulation.t)
        
        # remove the temporary files
        if not make_local_copy:
            self._delete_files(base_filename, extensions)

        
        return t, y
    
    
    def _write_file(self, filename, simulation, run_time, make_local_copy):
        connections = simulation._connections
        if make_local_copy:
            filenames = [filename, os.path.split(filename)[1]]
        else:
            filenames = [filename]

        for filename in filenames:
            file = open(filename + '_rhs.c', 'w')
            file.write('#include <math.h>\n'
                       '#include <stdlib.h>\n'
                       '#include "mex.h"\n'
                       '\n'
                       'double heav(double v) {\n'
                       '    return (v > 0);\n'
                       '}\n'
                       'double sign(double x) {\n'
                       '    if (x > 0) return(1);\n'
                       '    if (x < 0) return(-1);\n'
                       '    return(0);\n'
                       '}\n'
                       'double min(double x, double y) {\n'
                       '    return ((x < y) ? x : y);\n'
                       '}\n'
                       '\n'
                       'double max(double x, double y) {\n'
                       '    return ((x > y) ? x : y);\n'
                       '}\n'
                       '\n'
                       'void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {\n'
                       '    double* _rhs = mxGetPr(prhs[1]);\n'
                       '    double t = *mxGetPr(prhs[0]);\n')


            fluxes = _compute_fluxes(simulation, renderer=c_renderer, divide_volume=False)
            
            # declare the variables
            j = 0;
            
            declaration_normal = []
            declaration_volume = []
            # declare neuron state variables
            for i in xrange(len(simulation._neurons)):
                for var in simulation._neurons[i].diff_list:
                    if var in simulation._neurons[i].concentrations:
                        declaration = '    double %s%d = _rhs[%d]/%s;\n' % (var, i, j, _effective_volume(simulation, i, c_renderer, divide_volume=False))
                    else:
                        declaration = '    double %s%d = _rhs[%d];\n' % (var, i, j)
                    j += 1
                    if var == 'volume':
                        declaration_volume.append(declaration)
                    else:
                        declaration_normal.append(declaration)
            
            # declare shared dynamics state variables
            for i, dyn in enumerate(simulation._shared_dynamics):
                for var in dyn.diff_list:
                    if var in dyn.concentrations:
                        declaration = '    double %s = _rhs[%d]/%s;\n' % (var, j, _effective_volume(simulation, '', c_renderer, divide_volume=False))
                    else:
                        declaration = '    double %s = _rhs[%d];\n' % (var, j)
                    j += 1
                    if var == 'volume':
                        declaration_volume.append(declaration)
                    else:
                        declaration_normal.append(declaration)
                        
            for declaration in declaration_volume + declaration_normal:
                file.write(declaration)

            

            file.write('    plhs[0] = mxCreateDoubleMatrix(%d, 1, mxREAL);\n'
                       '    double* _deriv = mxGetPr(plhs[0]);\n' % (j))
            
            # noise sources
            if len(simulation._noise) > 1:
                file.write('    double ' + ', '.join(['noise%d' % i for i in range(1, len(simulation._noise))]) + ';\n')
            
            # write the noise equations
            _write_noise_equations_c(file, simulation._noise, simulation.t, run_time)

            j = 0
            for i in xrange(len(simulation._neurons)):
                neuron = simulation._neurons[i]
                
                # determine distinct synapse rules
                synapses = []
                presynaptics = []
                for syn in simulation._connections[i]:
                    syn_rule = syn['rule']
                    if syn_rule in synapses:
                        presynaptics[synapses.index(syn_rule)].append(syn['pre_id'])
                    else:
                        synapses.append(syn_rule)
                        presynaptics.append([syn['pre_id']])
                
                # compute synaptic inputs
                syn_inputs = []
                for k in xrange(len(synapses)):
                    syn_inputs.append(synapses[k].get_tokens(presynaptics[k], i, neuron._potential_var))
                
                _write_dynamics_c(file, neuron, j, fluxes, syn_inputs, i)
                j += len(neuron.diff_list)
            
            # write the global dynamics
            for d in simulation._shared_dynamics:
                _write_dynamics_c(file, d, j, fluxes)
                j += len(d.diff_list)
            file.write('}\n')
            
            # write the control code
            _write_matlab_control_code(self, filename, simulation, run_time, prepend='mex %s_rhs.c\n' % filename)



class MatlabIntegrator(Integrator):
    def __init__(self, temp_dir_handler = default_temp_dir_handler, method = 'ode23',
                 rel_tol = 1e-3, max_step = .05):
        """Initialize the MATLAB integrator interface.
        
        Writes equations to a MATLAB file, run it.
        
        Parameters:

            temp_dir_handler -- a function called at the beginning of each
                                simulation, returns a string with a temp directory
                                name, should also do any cleanup

            method -- MATLAB integration method (ode45, ode23, etc... could be
                      user defined)

            rel_tol -- passed to MATLAB RelTol
            
            max_step -- passed to MATLAB MaxStep
        """
        self._temp_dir_handler = temp_dir_handler
        self.method = method
        self.rel_tol = rel_tol
        self.max_step = max_step
    

    def __repr__(self):
        """Return a string representing the contents."""
        return "MatlabIntegrator(temp_dir_handler = %s, method = '%s', rel_tol = %g, max_step = %g)" % (self._temp_dir_handler.__name__, self.method, self.rel_tol, self.max_step)

    
    def integrate(self, simulation, run_time, temp_files, make_local_copy=False):
        self._temp_files = temp_files
        if simulation.has_synaptic_delays():
            raise SimulationFailedError(0, 'MatlabIntegrator does not support synaptic delays.')

        base_filename = self._get_base_filename(simulation)

        extensions = ['.dat', '_rhs.m', '.m']
        if not null_device:
            extensions.append('.garbage')
        self._register_temp_files(base_filename, extensions)   
        self._temp_files.append(base_filename + '.m')

        self._write_file(base_filename, simulation, run_time, make_local_copy)

        if null_device:
            garbage_file = null_device
        else:
            garbage_file = base_filename + '.garbage'
        base_dir = os.path.dirname(base_filename)
        if base_dir == '':
            set_correct_path = ''
        else:
            set_correct_path = 'cd ' + base_dir + '; '

        t, y = _get_data(set_correct_path + 'matlab -nodisplay < ' + base_filename + '.m > ' + garbage_file,
                         base_filename + '.dat', run_time + simulation.t)

        # remove the temporary files
        if not make_local_copy:
            self._delete_files(base_filename, extensions)

        
        return t, y
    
    
    
    def _write_file(self, filename, simulation, run_time, make_local_copy):
        connections = simulation._connections
        fluxes = _compute_fluxes(simulation, renderer=matlab_renderer, divide_volume=False)
        
        if make_local_copy:
            filenames = [filename, os.path.split(filename)[1]]
        else:
            filenames = [filename]

        for filename in filenames:
            # write the control code
            _write_matlab_control_code(self, filename, simulation, run_time)
            
            # write the differential equations
                       
            file = open(filename + '_rhs.m', 'w')
            file.write('function yprime = ' + os.path.basename(filename) + '_rhs(t, rhs_)\n')
            file.write('    yprime = zeros(size(rhs_));\n')
            
            # declare neuron state variables
            j = 0
            for i in range(len(simulation._neurons)):
                for var in simulation._neurons[i].diff_list:

                    j += 1
                    file.write('    %s%d = rhs_(%d);\n' % (var, i, j))
            
            # declare shared dynamics state variables
            for dyn in simulation._shared_dynamics:
                for var in dyn.diff_list:
                    j += 1
                    file.write('    %s = rhs_(%d);\n' % (var, j))
            

            # write the noise equations
            _write_noise_equations_matlab(file, simulation._noise, simulation.t, run_time)

            j = 0
            for i in range(len(simulation._neurons)):
                neuron = simulation._neurons[i]
                
                # determine distinct synapse rules
                synapses = []
                presynaptics = []
                for syn in simulation._connections[i]:
                    syn_rule = syn['rule']
                    if syn_rule in synapses:
                        presynaptics[synapses.index(syn_rule)].append(syn['pre_id'])
                    else:
                        synapses.append(syn_rule)
                        presynaptics.append([syn['pre_id']])
                
                # compute synaptic inputs
                syn_inputs = []
                for syn, presyn in zip(synapses, presynaptics):
                    syn_inputs.append(syn.get_tokens(presyn, i, neuron._potential_var))
                
                _write_dynamics_matlab(file, neuron, j, fluxes, syn_inputs, i)
                j += len(neuron.diff_list)

            # write the global dynamics
            for d in simulation._shared_dynamics:
                _write_dynamics_matlab(file, d, j, fluxes)
                j += len(d.diff_list)
            
            file.close()
        


def _write_matlab_control_code(method, filename, simulation, run_time, prepend = ''):
        file = open(filename + '.m', 'w')
        file.write(prepend)
        file.write('y0 = [];\n')
        
        # initial conditions
        for id, neuron in enumerate(simulation._neurons):
            initial = []
            for var in neuron.diff_list:
                initial.append(_initial_value(neuron, var, renderer=matlab_renderer, divide_volume=False))
            if initial == []:
                pass
            elif len(initial) == 1:
                file.write('y0 = [y0;' + initial[0] + '];\n')
            else:
                file.write('y0 = [y0;' + ';'.join(initial) + '];\n')
        for dynamics in simulation._shared_dynamics:
            initial = []
            for var in dynamics.diff_list:
                initial.append(_initial_value(dynamics, var, renderer=matlab_renderer, divide_volume=False))
            if initial == []:
                pass
            elif len(initial) == 1:
                file.write('y0 = [y0;' + initial[0] + '];\n')
            else:
                file.write('y0 = [y0;' + ';'.join(initial) + '];\n')

        # integration call
        file.write("success = false;\n"
                   "try\n"
                   "    [t, y] = %s(@%s_rhs, [%g, %g], y0, odeset('RelTol', %g, 'MaxStep', %g));\n"
                   "    success = true;\n"

                   % (method.method, os.path.basename(filename), simulation.t,
                      simulation.t + run_time, method.rel_tol, method.max_step))
        file.write("catch e\n"
                   "end\n"
                   "if not(success)\n"
                   "    f = fopen('" + filename + ".dat', 'w');\n"
                   "    fprintf(f, '0 0\\n');\n"
                   "    fclose(f);\n"
                   "end\n")
        # write dat (space separated)
        file.write("if success\n"
                   "    f = fopen('" + filename + ".dat', 'w');\n"
                   "    for i = 1 : length(t)\n"
                   "        fprintf(f, '%g', t(i));\n"
                   "        row = y(i, :);\n"
                   "        for j = 1 : length(row)\n"
                   "            fprintf(f, ' %g', row(j));\n"
                   "        end\n"
                   "        fprintf(f, '\\n');\n"
                   "    end\n"
                   "    fclose(f);\n"
                   "end\n")
                   
        file.close()
