compare stuff
    :param molecules:
    :param ensemble_lookup:
    :param options:
    :return:
6574
def compare(molecules, ensemble_lookup, options):
    

    print(" Analyzing differences ... ")
    print('')
    sort_order = classification.get_sort_order(molecules)

    ensemble1 = sorted(ensemble_lookup.keys())[0]
    ensemble2 = sorted(ensemble_lookup.keys())[1]

    stats = {}
    stats['header'] = [' ']
    name = os.path.basename(ensemble1).replace('.csv', '')
    stats['header'].append(name)
    name = os.path.basename(ensemble2).replace('.csv', '')
    stats['header'].append(name)
    stats['header'].append('Difference')
    stats['header'].append('95% CI')
    stats['header'].append('p-value')

    molecules1 = copy.deepcopy(molecules)
    molecules2 = copy.deepcopy(molecules)

    score_structure1 = classification.make_score_structure(molecules1, ensemble_lookup[ensemble1])
    score_structure2 = classification.make_score_structure(molecules2, ensemble_lookup[ensemble2])

    auc_structure_1 = classification.make_auc_structure(score_structure1)
    auc_structure_2 = classification.make_auc_structure(score_structure2)

    # calculate auc value differences
    auc_diff = classification.calculate_auc_diff(auc_structure_1, auc_structure_2, sort_order)

    stats['AUC'] = auc_diff

    # calculate enrichment factor differences
    fpfList = make_fpfList(options)
    for fpf in fpfList:
        fpf = float(fpf)
        ef_structure1 = classification.make_ef_structure(score_structure1, fpf, sort_order)
        ef_structure2 = classification.make_ef_structure(score_structure2, fpf, sort_order)

        if ef_structure1 and ef_structure2:
            ef_diff = classification.calculate_ef_diff(ef_structure1, ef_structure2, fpf)
            title = 'E%s' % fpf
            stats[title] = ef_diff

    # write results summary
    output.write_diff_summary(stats, options)

    # write roc curves
    if options.write_roc:
        print(" Writing ROC data ... ")
        print('')
        output.write_roc(auc_structure_1, ensemble1, options)
        output.write_roc(auc_structure_2, ensemble2, options)

    # plot
    if options.plot:
        print(" Making plots ... ")
        print('')
        plotter(molecules, ensemble_lookup, options)
712932
def create_isobaric_quant_lookup(fns, tablefn_map, featmap, pqdb,
                                 featcolnr, allquantcols, psmcolmap):
    """"""
    pqdb.store_quant_channels(
            map_psmnrcol_to_quantcol(allquantcols, psmcolmap, tablefn_map),
            psmcolmap)
    quantmap = pqdb.get_quantchannel_map()
    to_store = []
    for fn, header, pquant in tsvreader.generate_tsv_pep_protein_quants(fns):
        pqdata = get_isob_quant_data(pquant, header, tablefn_map[fn],
                                     featmap, featcolnr, quantmap)
        to_store.extend(pqdata)
        if len(to_store) > 10000:
            pqdb.store_isobaric_quants(to_store, psmcolmap)
            to_store = []
    pqdb.store_isobaric_quants(to_store, psmcolmap)
187851
def write_ensemble(ensemble, options):
    

    # set output file name
    size = len(ensemble)
    filename = '%s_%s_queries.csv' % (options.outname, size)
    file = os.path.join(os.getcwd(), filename)

    f = open(file, 'w')

    out = ', '.join(ensemble)

    f.write(out)

    f.close()
708774
def compare_ecp_pots(potential1, potential2, compare_meta=False, rel_tol=0.0):
    '''
     
    '''

    if potential1['angular_momentum'] != potential2['angular_momentum']:
        return False

    rexponents1 = potential1['r_exponents']
    rexponents2 = potential2['r_exponents']
    gexponents1 = potential1['gaussian_exponents']
    gexponents2 = potential2['gaussian_exponents']
    coefficients1 = potential1['coefficients']
    coefficients2 = potential2['coefficients']

    # integer comparison
    if rexponents1 != rexponents2:
        return False
    if not _compare_vector(gexponents1, gexponents2, rel_tol):
        return False
    if not _compare_matrix(coefficients1, coefficients2, rel_tol):
        return False
    if compare_meta:
        if potential1['ecp_type'] != potential2['ecp_type']:
            return False
        return True
    else:
        return True
874177
def _verify_configs(configs):
    
    if configs:
        scenario_names = [c.scenario.name for c in configs]
        for scenario_name, n in collections.Counter(scenario_names).items():
            if n > 1:
                msg = ("Duplicate scenario name '{}' found.  "
                       'Exiting.').format(scenario_name)
                util.sysexit_with_message(msg)

    else:
        msg = "'{}' glob failed.  Exiting.".format(MOLECULE_GLOB)
        util.sysexit_with_message(msg)
453060
def evaluate_list(molecules, ensemble_lookup, options):
    

    # create stats dictionaries to store results from each ensemble
    stats = {}  # {file name : metric_List}

    # print progress messages
    if options.write_roc:
        print(" Determining virtual screening performance and writing ROC data ... ")
        print('')
    else:
        print(" Determining virtual screening performance ...")
        print('')

    for filename in sorted(ensemble_lookup.keys()):
        metric_List = calculate_metrics(molecules, ensemble_lookup, filename, options)
        stats[filename] = metric_List

    # write results summary
    output.write_summary(stats, options, fw_type = None)

    # plot
    if options.plot:
        print(" Making plots ... ")
        print
        plotter(molecules, ensemble_lookup, options)
1330436
def rank_queries(molecules, ensemble, sort_order, options):
    

    # generate query list
    query_list = [x for x in list(molecules[0].scores.keys()) if x not in ensemble]

    results = {}

    for query in query_list:
        es = EnsembleStorage() # an ensemble storage project

        # generate test_ensemble
        test_ensemble = ensemble[0:]
        test_ensemble.append(query)
        test_ensemble = tuple(test_ensemble)
        es.set_prop('ensemble', test_ensemble)

        # calculate its performance
        score_structure = classification.make_score_structure(molecules, test_ensemble)

        # determine auc value
        auc_structure = classification.make_auc_structure(score_structure)
        auc = classification.calculate_auc(auc_structure, sort_order, 'no stats')
        es.set_prop('auc', auc)

        # if the enrichment factor was set to anything other than 1, then we're training to maximize the corresponding
        # enrichment factor
        for fpf in classification.make_fpfList(options, score_structure):
            fpf = float(fpf)
            ef_structure = classification.make_ef_structure(score_structure, fpf, sort_order)
            if ef_structure:
                ef = classification.calculate_ef(ef_structure, fpf)
                es.set_prop(ef[0], ef[1], 'ef')

        # append results to metric list
        results[test_ensemble] = es

    # peel away the best performing ensemble
    best_ensemble = screener.find_best_ensemble(results, options)

    return list(best_ensemble)
874552
def hydrophobic_interactions(atom_set_a, atom_set_b):
    
    data = namedtuple('hydroph_interaction', 'bsatom bsatom_orig_idx ligatom ligatom_orig_idx '
                                             'distance restype resnr reschain restype_l, resnr_l, reschain_l')
    pairings = []
    for a, b in itertools.product(atom_set_a, atom_set_b):
        if a.orig_idx == b.orig_idx:
            continue
        e = euclidean3d(a.atom.coords, b.atom.coords)
        if not config.MIN_DIST < e < config.HYDROPH_DIST_MAX:
            continue
        restype, resnr, reschain = whichrestype(a.atom), whichresnumber(a.atom), whichchain(a.atom)
        restype_l, resnr_l, reschain_l = whichrestype(b.orig_atom), whichresnumber(b.orig_atom), whichchain(b.orig_atom)
        contact = data(bsatom=a.atom, bsatom_orig_idx=a.orig_idx, ligatom=b.atom, ligatom_orig_idx=b.orig_idx,
                       distance=e, restype=restype, resnr=resnr,
                       reschain=reschain, restype_l=restype_l,
                       resnr_l=resnr_l, reschain_l=reschain_l)
        pairings.append(contact)
    return filter_contacts(pairings)
430688
def equalize_terminal_double_bond(mol):
    
    for i, a in mol.atoms_iter():
        if mol.neighbor_count(i) == 1:
            nb = list(mol.neighbors(i).values())[0]
            if nb.order == 2:
                nb.type = 2
318227
def similarity_lookup(self, ResponseGroup="Large", **kwargs):
        
        response = self.api.SimilarityLookup(
            ResponseGroup=ResponseGroup, **kwargs)
        root = objectify.fromstring(response)
        if root.Items.Request.IsValid == 'False':
            code = root.Items.Request.Errors.Error.Code
            msg = root.Items.Request.Errors.Error.Message
            raise SimilartyLookupException(
                "Amazon Similarty Lookup Error: '{0}', '{1}'".format(
                    code, msg))
        return [
            AmazonProduct(
                item,
                self.aws_associate_tag,
                self.api,
                region=self.region
            )
            for item in getattr(root.Items, 'Item', [])
        ]
663464
def _name_lookup(names):
    """
    mapper = dict(zip(range(len(names)), names))

    def lookup(val):
        return mapper.get(val, 'Unknown')
    return lookup
504268
def sort_basis(basis, use_copy=True):
    

    if use_copy:
        basis = copy.deepcopy(basis)

    for k, el in basis['elements'].items():
        if 'electron_shells' in el:
            el['electron_shells'] = sort_shells(el['electron_shells'], False)
        if 'ecp_potentials' in el:
            el['ecp_potentials'] = sort_potentials(el['ecp_potentials'], False)

    return sort_basis_dict(basis)
7395
def create_pep_protein_quant_lookup(fns, pqdb, poolnames, featcolnr, patterns,
                                    storefuns, isobqcolpattern=None,
                                    psmnrpattern=None):
    """"""
    tablefn_map = create_tablefn_map(fns, pqdb, poolnames)
    feat_map = pqdb.get_feature_map()
    for pattern, storefun in zip(patterns, storefuns):
        if pattern is None:
            continue
        colmap = get_colmap(fns, pattern, single_col=True)
        if colmap:
            store_single_col_data(fns, tablefn_map, feat_map,
                                  storefun, featcolnr, colmap)
    if isobqcolpattern is not None:
        isocolmap = get_colmap(fns, isobqcolpattern, antipattern=psmnrpattern)
    else:
        return
    if psmnrpattern is not None:
        psmcolmap = get_colmap(fns, psmnrpattern)
    else:
        psmcolmap = False
    create_isobaric_quant_lookup(fns, tablefn_map,
                                 feat_map, pqdb,
                                 featcolnr,
                                 isocolmap, psmcolmap)
757773
def equal(mol, query, largest_only=True, ignore_hydrogen=True):
    
    m = molutil.clone(mol)
    q = molutil.clone(query)
    if largest_only:
        m = molutil.largest_graph(m)
        q = molutil.largest_graph(q)
    if ignore_hydrogen:
        m = molutil.make_Hs_implicit(m)
        q = molutil.make_Hs_implicit(q)
    if molutil.mw(m) == molutil.mw(q):
        gm = GraphMatcher(q.graph, m.graph, node_match=atom_match)
        return gm.is_isomorphic()
    return False
5182
def basis_comparison_report(bs1, bs2, uncontract_general=False):
    '''
    
    '''

    all_bs1 = list(bs1['elements'].keys())

    if uncontract_general:
        bs1 = manip.uncontract_general(bs1)
        bs2 = manip.uncontract_general(bs2)

    not_in_bs1 = []  # Found in bs2, not in bs1
    not_in_bs2 = all_bs1.copy()  # Found in bs1, not in bs2
    no_diff = []  # Elements for which there is no difference
    some_diff = []  # Elements that are different
    big_diff = []  # Elements that are substantially different

    for k, v in bs2['elements'].items():
        if k not in all_bs1:
            not_in_bs1.append(k)
            continue

        print()
        print("-------------------------------------")
        print(" Element ", k)
        bs1_el = bs1['elements'][k]

        max_rdiff_el = 0.0
        max_rdiff_ecp = 0.0

        # Check to make sure that neither or both have ecp/electron shells
        if 'electron_shells' in v and 'electron_shells' not in bs1_el:
            print("bs2 has electron_shells, but bs1 does not")
            max_rdiff_el = float('inf')
        if 'electron_shells' in bs1_el and 'electron_shells' not in v:
            print("bs1 has electron_shells, but bs2 does not")
            max_rdiff_el = float('inf')
        if 'ecp_potentials' in v and 'ecp_potentials' not in bs1_el:
            print("bs2 has ecp_potentials, but bs1 does not")
            max_rdiff_ecp = float('inf')
        if 'ecp_potentials' in bs1_el and 'ecp_potentials' not in v:
            print("bs1 has ecp_potentials, but bs2 does not")
            max_rdiff_ecp = float('inf')

        if 'electron_shells' in v and 'electron_shells' in bs1_el:
            max_rdiff_el = max(max_rdiff_el, shells_difference(v['electron_shells'], bs1_el['electron_shells']))
        if 'ecp_potentials' in v and 'ecp_potentials' in bs1_el:
            nel1 = v['ecp_electrons']
            nel2 = bs1_el['ecp_electrons']
            if int(nel1) != int(nel2):
                print('Different number of electrons replaced by ECP ({} vs {})'.format(nel1, nel2))
                max_rdiff_ecp = float('inf')
            else:
                max_rdiff_ecp = max(max_rdiff_ecp, potentials_difference(v['ecp_potentials'],
                                                                         bs1_el['ecp_potentials']))

        max_rdiff = max(max_rdiff_el, max_rdiff_ecp)

        # Handle some differences
        if max_rdiff == float('inf'):
            big_diff.append(k)
        elif max_rdiff == 0.0:
            no_diff.append(k)
        else:
            some_diff.append(k)

        not_in_bs2.remove(k)

    print()
    print("     Not in bs1: ", _print_list(not_in_bs1))
    print("     Not in bs2: ", _print_list(not_in_bs2))
    print("  No difference: ", _print_list(no_diff))
    print("Some difference: ", _print_list(some_diff))
    print(" BIG difference: ", _print_list(big_diff))
    print()

    return (len(not_in_bs1) == 0 and len(not_in_bs2) == 0 and len(some_diff) == 0 and len(big_diff) == 0)
871126
def run(itf):
    

    if not itf:
        return 1

    # access user input
    options = SplitInput(itf)

    # check input args
    error_check(options)

    # read input files
    try:
        molecules, ensemble_lookup = ReadFiles(options)
    except:
        return 1

    if options.compare:
        compare(molecules, ensemble_lookup, options)
    else:
        evaluate_list(molecules, ensemble_lookup, options)
792718
def checkEdges(ringSet, lookup, oatoms):
    
    bondedAtoms = map( None, ringSet[:-1], ringSet[1:] )
    bondedAtoms += [ (ringSet[-1], ringSet[0]) ]

    # form a lookup for the ringSet list
    atomSet = {}
    for atomID in ringSet:
	atomSet[atomID] = 1
    results = []
    
    # for each bond in the ring, break it and find the smallest
    # rings starting on either side of the bond
    # keep the largest but rememeber to add the bond back at the
    # end
    for atom1, atom2 in bondedAtoms:
        # break a single edge in the ring
        handle1 = atom1.handle
        handle2 = atom2.handle
        oatoms1 = oatoms[handle1]
        oatoms2 = oatoms[handle2]
        index1 = oatoms1.index(atom2)
        index2 = oatoms2.index(atom1)

        # break the bond
        del oatoms1[index1]
	del oatoms2[index2]

	ring1 = getRing(atom1, atomSet, lookup, oatoms)
	ring2 = getRing(atom2, atomSet, lookup, oatoms)
	
	# keep the larger of the two rings
	if len(ring1) > len(ring2):
	    results.append((len(ring1),
                            handle1, handle2,
                            ring1))
	else:
    	    results.append((len(ring2),
                            handle2, handle1,
                            ring2))

        # retie the bond
        oatoms1.insert(index1, atom2)
        oatoms2.insert(index2, atom1)


    if not results:
	return None

    #  find the smallest ring
    size, incidentHandle, adjacentHandle, smallestRing = min(results)
    # dereference the handles
    incident, adjacent = lookup[incidentHandle], lookup[adjacentHandle]

    # break the bond between the incident and adjacent atoms
    oatomsI = oatoms[incidentHandle]
    oatomsA = oatoms[adjacentHandle]
    assert incident in oatomsA
    assert adjacent in oatomsI
    
    oatomsI.remove(adjacent)
    oatomsA.remove(incident)