Source code for allensdk.internal.model.biophysical.run_passive_fit

import os
import sys
import subprocess
import numpy as np
import allensdk.internal.model.biophysical.ephys_utils as ephys_utils
from .passive_fitting import preprocess as passive_prep
import allensdk.core.json_utilities as ju
from allensdk.model.biophys_sim.config import Config
from allensdk.core.nwb_data_set import NwbDataSet
from allensdk.internal.model.biophysical.passive_fitting import neuron_passive_fit
from allensdk.internal.model.biophysical.passive_fitting import neuron_passive_fit2
from allensdk.internal.model.biophysical.passive_fitting import neuron_passive_fit_elec
from pkg_resources import resource_filename #@UnresolvedImport
import logging
import logging.config as lc


_run_passive_fit_log = logging.getLogger('allensdk.internal.model.biophysical.run_passive_fit')


def run_passive_fit(description):
    output_directory = description.manifest.get_path('WORKDIR')
    neuronal_model = ju.read(description.manifest.get_path('neuronal_model_data'))
    specimen_data = neuronal_model['specimen']
    
    is_spiny = not any(t['name'] == u'dendrite type - aspiny' for t in specimen_data['specimen_tags'])
    
    all_sweeps = specimen_data['ephys_sweeps']
    if not os.path.exists(output_directory):
        os.makedirs(output_directory)

    cap_check_sweeps, _, _ = \
        ephys_utils.get_sweeps_of_type('C1SQCAPCHK',
                                       all_sweeps)
    
    passive_fit_data = {}

    if len(cap_check_sweeps) > 0:
        data_set = NwbDataSet(description.manifest.get_path('stimulus_path'))
        d = passive_prep.get_passive_fit_data(cap_check_sweeps, data_set);

        grand_up_file = os.path.join(output_directory, 'upbase.dat')
        np.savetxt(grand_up_file, d['grand_up'])

        grand_down_file = os.path.join(output_directory, 'downbase.dat')
        np.savetxt(grand_down_file, d['grand_down'])
        
        passive_fit_data["bridge"] = d['bridge_avg']
        passive_fit_data["escape_time"] = d['escape_t']

        fit_1_file = description.manifest.get_path('fit_1_file')
        fit_1_params = subprocess.check_output([sys.executable,
                                                '-m', neuron_passive_fit.__name__, 
                                                str(d['escape_t']),
                                                os.path.realpath(description.manifest.get_path('manifest')) ])
        passive_fit_data['fit_1'] = ju.read(fit_1_file)

        fit_2_file = description.manifest.get_path('fit_2_file')

        fit_2_params = subprocess.check_output([sys.executable,
                                                '-m', neuron_passive_fit2.__name__,
                                                str(d['escape_t']),
                                                os.path.realpath(description.manifest.get_path('manifest')) ])
        passive_fit_data['fit_2'] = ju.read(fit_2_file)

        fit_3_file = description.manifest.get_path('fit_3_file')
        fit_3_params = subprocess.check_output([sys.executable,
                                                '-m', neuron_passive_fit_elec.__name__,
                                                str(d['escape_t']),
                                                str(d['bridge_avg']),
                                                str(1.0),
                                                os.path.realpath(description.manifest.get_path('manifest')) ])
        passive_fit_data['fit_3'] = ju.read(fit_3_file)
        
        # Check for potentially problematic outcomes
        cm_rel_delta = (passive_fit_data["fit_1"]["Cm"] - passive_fit_data["fit_3"]["Cm"]) / passive_fit_data["fit_1"]["Cm"]
        if passive_fit_data["fit_2"]["err"] < passive_fit_data["fit_1"]["err"]:
            _run_passive_fit_log.debug("Fixed Ri gave better results than original")
            if passive_fit_data["fit_2"]["err"] < passive_fit_data["fit_3"]["err"]:
                _run_passive_fit_log.debug("Using fixed Ri results")
                passive_fit_data["fit_for_next_step"] = passive_fit_data["fit_2"]
            else:
                _run_passive_fit_log.debug("Using electrode results")
                passive_fit_data["fit_for_next_step"] = passive_fit_data["fit_3"]
        elif abs(cm_rel_delta) > 0.1:
            _run_passive_fit_log.debug("Original and electrode fits not in sync:")
            _run_passive_fit_log.debug("original Cm: " + str(passive_fit_data["fit_1"]["Cm"]))
            _run_passive_fit_log.debug("w/ electrode Cm: " + str(passive_fit_data["fit_3"]["Cm"]))
            if passive_fit_data["fit_1"]["err"] < passive_fit_data["fit_3"]["err"]:
                _run_passive_fit_log.debug("Original has lower error")
                passive_fit_data["fit_for_next_step"] = passive_fit_data["fit_1"]
            else:
                _run_passive_fit_log.debug("Electrode has lower error")
                passive_fit_data["fit_for_next_step"] = passive_fit_data["fit_3"]
        else:
            passive_fit_data["fit_for_next_step"] = passive_fit_data["fit_1"]

        ra = passive_fit_data["fit_for_next_step"]["Ri"]
        if is_spiny:
            combo_cm = passive_fit_data["fit_for_next_step"]["Cm"]
            a1 = passive_fit_data["fit_for_next_step"]["A1"]
            a2 = passive_fit_data["fit_for_next_step"]["A2"]
            cm1 = 1.0
            cm2 = (combo_cm * (a1 + a2) - a1) / a2
        else:
            cm1 = passive_fit_data["fit_for_next_step"]["Cm"]
            cm2 = passive_fit_data["fit_for_next_step"]["Cm"]
    else:
        _run_passive_fit_log.debug("No cap check trace found")
        ra = 100.0
        cm1 = 1.0
        if is_spiny:
            cm2 = 2.0
        else:
            cm2 = 1.0

    passive_fit_data['ra'] = ra
    passive_fit_data['cm1'] = cm1
    passive_fit_data['cm2'] = cm2
    
    return passive_fit_data


def main(limit, manifest_path):
    app_config = Config()
    description = app_config.load(manifest_path)

    if 'LOG_CFG' in os.environ:
        log_config = os.environ['LOG_CFG']
    else:
        log_config = resource_filename('allensdk.model.biophysical',
                                       'logging.conf')
        os.environ['LOG_CFG'] = log_config
    lc.fileConfig(log_config)

    run_passive_fit(description)


if __name__ == "__main__":
    limit = sys.argv[-2]
    manifest_path = sys.argv[-1]
    
    main(limit, manifest_path)