<commit_before>import numpy as np
import sys
import os
sys.path.append('/Users/palmer/Documents/python_codebase/')


def get_variables(json_filename):
    import json
    config = json.loads(open(json_filename).read())
    # maintain compatibility with previous versions
    if 'clean_im' not in config['image_generation']:
        config['image_generation']['clean_im']=True
        config['image_generation']['interp']=''
    return config


### We simulate a mass spectrum for each sum formula/adduct combination. This generates a set of isotope patterns (see http://www.mi.fu-berlin.de/wiki/pub/ABI/QuantProtP4/isotope-distribution.pdf) which can provide additional informaiton on the molecule detected. This gives us a list of m/z centres for the molecule
def calculate_isotope_patterns(sum_formulae, adduct='', isocalc_sig=0.01, isocalc_resolution=200000.,
                                   isocalc_do_centroid=True, charge='1'):
    from pyMS.pyisocalc import pyisocalc
    ### Generate a mz list of peak centroids for each sum formula with the given adduct
    mz_list = {}
    for n, sum_formula in enumerate(sum_formulae):
        sf = pyisocalc.complex_to_simple(sum_formula+adduct)
        if sf == None: #negative atoms as a result of simplification
            print 'negative adduct for {} : {}'.format(sum_formula,adduct)
            continue
        isotope_ms = pyisocalc.isodist(sf, plot=False, sigma=isocalc_sig, charges=charge,
                                       resolution=isocalc_resolution, do_centroid=isocalc_do_centroid)
        if not sum_formula in mz_list:
            mz_list[sum_formula] = {}
        mz_list[sum_formula][adduct] = isotope_ms.get_spectrum(source='centroids')
    return mz_list


def generate_isotope_patterns(config):
    from pySpatialMetabolomics.parse_databases import parse_databases

    import pickle
    # Extract variables from config dict
    db_filename = config['file_inputs']['database_file']
    db_dump_folder = config['file_inputs']['database_load_folder']
    isocalc_sig = float(config['isotope_generation']['isocalc_sig'])
    isocalc_resolution = float(config['isotope_generation']['isocalc_resolution'])
    if len(config['isotope_generation']['charge']) > 1:
        print 'Warning: only first charge state currently accepted'
    charge = int('{}{}'.format(config['isotope_generation']['charge'][0]['polarity'],
                               config['isotope_generation']['charge'][0][
                                   'n_charges']))  # currently only supports first charge!!
    adducts = [a['adduct'] for a in config['isotope_generation']['adducts']]

    # Read in molecules
    sum_formulae = parse_databases.read_generic_csv(db_filename)
    if '' in sum_formulae:
        print 'empty sf removed from list'
        del sum_formulae['']
    # Check if already genrated and load if possible, otherwise calculate fresh   
    db_name = os.path.splitext(os.path.basename(db_filename))[0]
    mz_list = {}
    for adduct in adducts:
        load_file = '{}/{}_{}_{}_{}.dbasedump'.format(db_dump_folder, db_name, adduct, isocalc_sig, isocalc_resolution)
        if os.path.isfile(load_file):
            print "{} -> loading".format(load_file)
            mz_list_tmp = pickle.load(open(load_file, 'r'))
        else:
            print "{} -> generating".format(load_file)
            mz_list_tmp = calculate_isotope_patterns(sum_formulae, adduct=adduct, isocalc_sig=isocalc_sig,
                                                     isocalc_resolution=isocalc_resolution, charge=charge)
            if db_dump_folder != "":
                pickle.dump(mz_list_tmp, open(load_file, 'w'))
        # add patterns to total list
        for sum_formula in sum_formulae:
            if sum_formula not in mz_list_tmp:# could be missing if [M-a] would have negative atoms
                continue
            if sum_formula not in mz_list:
                mz_list[sum_formula]={}
            ## this limit of 4 is hardcoded to reduce the number of calculations
            n = np.min([4,len(mz_list_tmp[sum_formula][adduct][0])])
            mz_list[sum_formula][adduct] = [mz_list_tmp[sum_formula][adduct][0][0:n],mz_list_tmp[sum_formula][adduct][1][0:n]]
    print 'all isotope patterns generated and loaded'
    return sum_formulae, adducts, mz_list


def hot_spot_removal(xics, q):
    for xic in xics:
        xic_q = np.percentile(xic, q)
        xic[xic > xic_q] = xic_q
    return xics


def run_search(config, IMS_dataset, sum_formulae, adducts, mz_list):
    from pyIMS.image_measures import level_sets_measure, isotope_image_correlation, isotope_pattern_match
    import time
    ### Runs the main pipeline
    # Get sum formula and predicted m/z peaks for molecules in database
    ppm = config['image_generation']['ppm']  # parts per million -  a measure of how accuracte the mass spectrometer is
    nlevels = config['image_generation']['nlevels']  # parameter for measure of chaos
    q = config['image_generation']['q']
    clean_im = config['image_generation']['clean_im']
    interp = config['image_generation']['interp']
    measure_value_score = {}
    iso_correlation_score = {}
    iso_ratio_score = {}
    t0 = time.time()
    t_el = 0
    for adduct in adducts:
        print 'searching -> {}'.format(adduct)
        for ii,sum_formula in enumerate(sum_formulae):
            if adduct not in mz_list[sum_formula]:#adduct may not be present if it would make an impossible formula, is there a better way to handle this?
                # print '{} adduct not found for {}'.format(adduct, mz_list[sum_formula])
                continue
            if time.time() - t_el > 10.:
                t_el = time.time()
                print '{:3.2f} done in {:3.0f} seconds'.format(float(ii)/len(sum_formulae),time.time()-t0)
            # Allocate dicts if required
            if not sum_formula in measure_value_score:
                    measure_value_score[sum_formula] = {}
            if not sum_formula in iso_correlation_score:
                    iso_correlation_score[sum_formula] = {}
            if not sum_formula in iso_ratio_score:
                    iso_ratio_score[sum_formula] = {}
            try:
                # 1. Generate ion images
                ion_datacube = IMS_dataset.get_ion_image(mz_list[sum_formula][adduct][0],
                                                         ppm)  # for each spectrum, sum the intensity of all peaks within tol of mz_list
                ion_datacube.xic = hot_spot_removal(ion_datacube.xic, q)
                # 2. Spatial Chaos
                measure_value_score[sum_formula][adduct] = level_sets_measure.measure_of_chaos(
                    ion_datacube.xic_to_image(0), nlevels, interp=interp, clean_im=clean_im)[0]
                if measure_value_score[sum_formula][adduct] == 1:
                    measure_value_score[sum_formula][adduct] = 0
                # 3. Score correlation with monoiso
                if len(mz_list[sum_formula][adduct][1]) > 1:
                    iso_correlation_score[sum_formula][adduct] = isotope_image_correlation.isotope_image_correlation(
                        ion_datacube.xic, weights=mz_list[sum_formula][adduct][1][1:])
                else:  # only one isotope peak, so correlation doesn't make sense
                    iso_correlation_score[sum_formula][adduct] = 1
                # 4. Score isotope ratio
                iso_ratio_score[sum_formula][adduct] = isotope_pattern_match.isotope_pattern_match(ion_datacube.xic,
                                                                                                   mz_list[sum_formula][
                                                                                                       adduct][1])
            except KeyError:<commit_msg>image processing exists as a separate function - currently on quantile scaling implemented (smoothing still in level_sets) updated config handling of default parameters for compatability with older config files<commit_after>import numpy as np
import sys
import os
sys.path.append('/Users/palmer/Documents/python_codebase/')


def d_update(d, u):
    import collections
    #http://stackoverflow.com/questions/3232943/update-value-of-a-nested-dictionary-of-varying-depth
    for k, v in u.iteritems():
        if isinstance(v, collections.Mapping):
            r = d_update(d.get(k, {}), v)
            d[k] = r
        else:
            d[k] = u[k]
    return d


def get_variables(json_filename):
    import json
    config = json.loads(open(json_filename).read())
    # maintain compatibility with previous versions
    # defaults are be how everything *should* be -> update makes sure that whatever is loaded conforms to this
    compatability_defaults = {'image_generation':{'smooth': True, 'smooth_params':()} }
    config = d_update(compatability_defaults,config)
    return compatability_defaults


### We simulate a mass spectrum for each sum formula/adduct combination. This generates a set of isotope patterns (see http://www.mi.fu-berlin.de/wiki/pub/ABI/QuantProtP4/isotope-distribution.pdf) which can provide additional informaiton on the molecule detected. This gives us a list of m/z centres for the molecule
def calculate_isotope_patterns(sum_formulae, adduct='', isocalc_sig=0.01, isocalc_resolution=200000.,
                                   isocalc_do_centroid=True, charge='1'):
    from pyMS.pyisocalc import pyisocalc
    ### Generate a mz list of peak centroids for each sum formula with the given adduct
    mz_list = {}
    for n, sum_formula in enumerate(sum_formulae):
        sf = pyisocalc.complex_to_simple(sum_formula+adduct)
        if sf == None: #negative atoms as a result of simplification
            print 'negative adduct for {} : {}'.format(sum_formula,adduct)
            continue
        isotope_ms = pyisocalc.isodist(sf, plot=False, sigma=isocalc_sig, charges=charge,
                                       resolution=isocalc_resolution, do_centroid=isocalc_do_centroid)
        if not sum_formula in mz_list:
            mz_list[sum_formula] = {}
        mz_list[sum_formula][adduct] = isotope_ms.get_spectrum(source='centroids')
    return mz_list


def generate_isotope_patterns(config):
    from pySpatialMetabolomics.parse_databases import parse_databases

    import pickle
    # Extract variables from config dict
    db_filename = config['file_inputs']['database_file']
    db_dump_folder = config['file_inputs']['database_load_folder']
    isocalc_sig = float(config['isotope_generation']['isocalc_sig'])
    isocalc_resolution = float(config['isotope_generation']['isocalc_resolution'])
    if len(config['isotope_generation']['charge']) > 1:
        print 'Warning: only first charge state currently accepted'
    charge = int('{}{}'.format(config['isotope_generation']['charge'][0]['polarity'],
                               config['isotope_generation']['charge'][0][
                                   'n_charges']))  # currently only supports first charge!!
    adducts = [a['adduct'] for a in config['isotope_generation']['adducts']]

    # Read in molecules
    sum_formulae = parse_databases.read_generic_csv(db_filename)
    if '' in sum_formulae:
        print 'empty sf removed from list'
        del sum_formulae['']
    # Check if already genrated and load if possible, otherwise calculate fresh   
    db_name = os.path.splitext(os.path.basename(db_filename))[0]
    mz_list = {}
    for adduct in adducts:
        load_file = '{}/{}_{}_{}_{}.dbasedump'.format(db_dump_folder, db_name, adduct, isocalc_sig, isocalc_resolution)
        if os.path.isfile(load_file):
            print "{} -> loading".format(load_file)
            mz_list_tmp = pickle.load(open(load_file, 'r'))
        else:
            print "{} -> generating".format(load_file)
            mz_list_tmp = calculate_isotope_patterns(sum_formulae, adduct=adduct, isocalc_sig=isocalc_sig,
                                                     isocalc_resolution=isocalc_resolution, charge=charge)
            if db_dump_folder != "":
                pickle.dump(mz_list_tmp, open(load_file, 'w'))
        # add patterns to total list
        for sum_formula in sum_formulae:
            if sum_formula not in mz_list_tmp:# could be missing if [M-a] would have negative atoms
                continue
            if sum_formula not in mz_list:
                mz_list[sum_formula]={}
            ## this limit of 4 is hardcoded to reduce the number of calculations
            n = np.min([4,len(mz_list_tmp[sum_formula][adduct][0])])
            mz_list[sum_formula][adduct] = [mz_list_tmp[sum_formula][adduct][0][0:n],mz_list_tmp[sum_formula][adduct][1][0:n]]
    print 'all isotope patterns generated and loaded'
    return sum_formulae, adducts, mz_list


def hot_spot_removal(xics, q):
    for xic in xics:
        xic_q = np.percentile(xic, q)
        xic[xic > xic_q] = xic_q
    return xics

def apply_image_processing(config, ion_datacube):
    """
    Function to apply pre-defined image processing methods to ion_datacube
    #todo: expose parameters in config
    :param ion_datacube:
    :return:
    """
    q = config['image_generation']['q']
    if q > 0:
        ion_datacube.xic = hot_spot_removal(ion_datacube.xic, q)
    smooth_method = config['image_generation']['smooth']
    smooth_params = config['image_generation']['smooth_params']
    if smooth_method != '':
        raise NotImplementedError
        from pyIMS import smoothing
        for ii in range(n_im):
            im = ion_datacube.xic_to_image(ii)
            im = smoothing.smooth(smooth_method,*smooth_params)
        #todo put xic back into vecotr


def run_search(config, IMS_dataset, sum_formulae, adducts, mz_list):
    from pyIMS.image_measures import level_sets_measure, isotope_image_correlation, isotope_pattern_match
    import time
    ### Runs the main pipeline
    # Get sum formula and predicted m/z peaks for molecules in database
    ppm = config['image_generation']['ppm']  # parts per million -  a measure of how accuracte the mass spectrometer is
    nlevels = config['image_generation']['nlevels']  # parameter for measure of chaos
    do_preprocessing = config['image_generation']['do_preprocessing']
    interp = config['image_generation']['smooth']
    measure_value_score = {}
    iso_correlation_score = {}
    iso_ratio_score = {}
    t0 = time.time()
    t_el = 0
    for adduct in adducts:
        print 'searching -> {}'.format(adduct)
        for ii,sum_formula in enumerate(sum_formulae):
            if adduct not in mz_list[sum_formula]:#adduct may not be present if it would make an impossible formula, is there a better way to handle this?
                # print '{} adduct not found for {}'.format(adduct, mz_list[sum_formula])
                continue
            if time.time() - t_el > 10.:
                t_el = time.time()
                print '{:3.2f} done in {:3.0f} seconds'.format(float(ii)/len(sum_formulae),time.time()-t0)
            # Allocate dicts if required
            if not sum_formula in measure_value_score:
                    measure_value_score[sum_formula] = {}
            if not sum_formula in iso_correlation_score:
                    iso_correlation_score[sum_formula] = {}
            if not sum_formula in iso_ratio_score:
                    iso_ratio_score[sum_formula] = {}
            try:
                # 1. Generate ion images
                ion_datacube = IMS_dataset.get_ion_image(mz_list[sum_formula][adduct][0],
                                                         ppm)  # for each spectrum, sum the intensity of all peaks within tol of mz_list
                if do_preprocessing:
                    apply_image_processing(config,ion_datacube) #currently just supports hot-spot removal
                # 2. Spatial Chaos
                measure_value_score[sum_formula][adduct] = level_sets_measure.measure_of_chaos(
                    ion_datacube.xic_to_image(0), nlevels, interp=interp, clean_im=False)[0]
                if measure_value_score[sum_formula][adduct] == 1:
                    measure_value_score[sum_formula][adduct] = 0
                # 3. Score correlation with monoiso
                if len(mz_list[sum_formula][adduct][1]) > 1:
                    iso_correlation_score[sum_formula][adduct] = isotope_image_correlation.isotope_image_correlation(
                        ion_datacube.xic, weights=mz_list[sum_formula][adduct][1][1:])
                else:  # only one isotope peak, so correlation doesn't make sense
                    iso_correlation_score[sum_formula][adduct] = 1
                # 4. Score isotope ratio
                iso_ratio_score[sum_formula][adduct] = isotope_pattern_match.isotope_pattern_match(ion_datacube.xic,
                                                                                                   mz_list[sum_formula][
                                                                                                       adduct][1])
            except KeyError as e:
                print str(e)
{'fraction_non_alphanumeric': 0.09328336358845231, 'fraction_numerical': 0.005855829478245593, 'mean_word_length': 2.554954204829309, 'pattern_counts': {'":': 2, '<': 3, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4704027', 'n_tokens_mistral': 5262, 'n_tokens_neox': 4860, 'n_words': 1181}