text:
from conn import DB
import numpy as np
from query import Type
from sklearn.metrics import f1_score

class QueryBag:
    """intermediate structure when formalizing query matrix representation"""
    def __init__(self):
        self.id_set = set()
        self.query_list = []

    def add_q(self, new_q):
        if type(new_q) is list:
            self.query_list.extend(new_q)
        else:
            self.query_list.append(new_q)

    def add_ids(self, ids):
        self.id_set.update(set(ids))

class Query:
    """the query class"""

    query_type = Type.QueryType.WCQ

    def __init__(self):
        self.index = -1
        self.cond_list = []

        self.query_type = -1

        self.query_bags = []  # intermediate holder for each condition/predicate
        self.query_matrix = []  # W
        self.domain_hist = []  # x
        self.true_answer = []  # the true counting of queries using matrix mechanism
        self.noisy_answer = []  # noisy counting
        self.lap_noise = []  # a list of laplace noise

        self.cardinality = -1  # cardinality of the table, will assign when get_cond_list

        # WCQ related
        self.answer_diff = []

        # ICQ related
        self.icq_c = -1
        self.selected_cond_index = []
        self.count_poking = -1
        self.real_count_poking = []  # a list of true poking times per predicate

        # TCQ related
        self.tcq_k = -1
        self.selected_arg_list = []

        # data set related
        self.table_name = ""
        self.data_size = 1
        self.db = ""

        self.m_type = ""
        self.query_key = -1

        self.Wx_cache_W = dict()
        self.Wx_cache_x = dict()
        self.Wx_cache_N = dict()
        self.eps_cache = dict()

    # set the query index
    def set_index(self, i):
        self.index = i

    # set query predicates
    def set_cond_list(self, p):
        self.cond_list = p

    # set the table and db connection
    def set_data_and_size(self, ds, s):
        if ds == 'census':
            self.table_name = 'census.income'
        elif ds == 'location':
            self.table_name = 'location.trip'
        else:
            print("ERROR: not supported data set ", ds)
        # set the db connection and size
        self.db = DB.DB(ds)
        self.data_size = s

    def set_query_type(self, qt, icq_c=-1, mpm_p=-1, tcq_k=-1):
        self.query_type = qt
        if qt == Type.QueryType.ICQ:
            assert icq_c != -1
            self.icq_c = icq_c
            self.count_poking = mpm_p
            self.real_count_poking = 0
        elif qt == Type.QueryType.TCQ:
            assert tcq_k != -1
            self.tcq_k = tcq_k

    def set_mechanism_type(self, m):
        self.m_type = m

    def set_cache(self, w, x, n, e):
        self.Wx_cache_W = w
        self.Wx_cache_x = x
        self.Wx_cache_N = n
        self.eps_cache = e

    # collect the ids for each condition to build histogram
    def get_cond_ids(self):
        counter = 0
        table_name = self.table_name + str(self.data_size)

        for cond in self.cond_list:
            # initialize an object to hold the ids and its belonging queries
            crnt_query_bag = QueryBag()
            crnt_query_bag.add_q(counter)

            # prepare the query to query ids
            crnt_query = 'select id from ' + table_name + ' where ' + cond
            crnt_query_rs = self.db.run(crnt_query)

            # store the row ids
            crnt_query_bag.add_ids([i[0] for i in crnt_query_rs])
            self.query_bags.append(crnt_query_bag)
            # print(len(crnt_query_rs))
            counter = counter + 1

        # get the total count from the table
        crnt_query = 'select count(*) from ' + table_name
        self.cardinality = self.db.run(crnt_query)[0][0]
        self.db.close_conn()

    def get_cond_counts(self):
        table_name = self.table_name + str(self.data_size)
        counts_list = []

        for cond in self.cond_list:
            # prepare the query to query ids
            crnt_query = 'select count(*) from ' + table_name + ' where ' + cond
            counts_list.append(self.db.run(crnt_query)[0][0])

        # get the total count from the table
        crnt_query = 'select count(*) from ' + table_name
        self.cardinality = self.db.run(crnt_query)[0][0]
        self.db.close_conn()

        return counts_list

    # construct the matrix representation of the query
    def to_matrix(self):

        # check if the results was previously cached
        # key qw1_100_1
        self.query_key = self.index.__name__ + str(len(self.cond_list)) + str(self.data_size)

        if self.query_key in self.Wx_cache_W:
            self.domain_hist = self.Wx_cache_W.get(self.query_key)
            self.query_matrix = self.Wx_cache_x.get(self.query_key)
            self.cardinality = self.Wx_cache_N.get(self.query_key)
            print("DEBUG: reuse cache query_key=", self.query_key)
            return
        else:
            print("DEBUG: no cache query_key=", self.query_key)

        # if the query is qt_2 or qt_4, otherwise just count
        if self.index.__name__ in ['qt_2', 'qt_4']:  # HD
            # first of all, collect the ids for each condition
            self.get_cond_ids()

            disjoint_query_bags = list()
            print("len(self.query_bags)= ", len(self.query_bags))
            disjoint_query_bags.append(self.query_bags[0])

            for i in range(1, len(self.cond_list)):
                # use the current query bag to iterate the existing disjoint bags
                crnt_query_bag = self.query_bags[i]
                disjoint_query_bags_size = len(disjoint_query_bags)

                for j in range(0, disjoint_query_bags_size):
                    existing_query_bag = disjoint_query_bags[j]
                    # find the intersection
                    set_intersect = existing_query_bag.id_set & crnt_query_bag.id_set
                    if len(set_intersect) > 0:
                        # intersection is not empty
                        intersect_query_bag = QueryBag()

                        # add the intersected ids into new query bag
                        intersect_query_bag.add_ids(set_intersect)

                        # set the queries from both current and previous query bags
                        intersect_query_bag.add_q(existing_query_bag.query_list)
                        intersect_query_bag.add_q(crnt_query_bag.query_list)

                        # extract the intersection from previous two bags
                        existing_query_bag.id_set = existing_query_bag.id_set - set_intersect
                        crnt_query_bag.id_set = crnt_query_bag.id_set - set_intersect

                        # append the intersected query bag into disjoint bag list
                        disjoint_query_bags.append(intersect_query_bag)

                # add the remaining crnt_query_bag into disjoint bag list
                if len(crnt_query_bag.id_set) > 0:
                    disjoint_query_bags.append(crnt_query_bag)

            # clean the disjoin query bags
            caught_count = 0
            for crnt_bag in disjoint_query_bags:
                if len(crnt_bag.id_set) == 0:
                    disjoint_query_bags.remove(crnt_bag)
                caught_count += len(crnt_bag.id_set)

            # add the uncaught ids into the last position
            if caught_count < self.cardinality:
                last_query_bag = QueryBag()
                last_query_bag.add_ids(range(0, self.cardinality - caught_count))
                disjoint_query_bags.append(last_query_bag)

            disjoint_query_bags_size = len(disjoint_query_bags)

            # generate the domain histogram
            self.domain_hist = [len(crnt_bag.id_set) for crnt_bag in disjoint_query_bags]
            print("domain_hist: \n", self.domain_hist)
            print("domain_hist_len: \n", len(self.domain_hist))

            # test disjoint query bags
            # for i in range(0, len(disjoint_query_bags)):
            #     crnt_bag = disjoint_query_bags[i]
            #     print("query_list= ", crnt_bag.query_list)
            #     print("id_set= ", crnt_bag.id_set)

            # initialize a matrix
            self.query_matrix = np.zeros((len(self.cond_list), disjoint_query_bags_size))
            for i in range(0, disjoint_query_bags_size):
                crnt_query_list = disjoint_query_bags[i].query_list
                for j in range(0, len(crnt_query_list)):
                    self.query_matrix[crnt_query_list[j]][i] = 1

            print("query matrix: \n", self.query_matrix)

        else:
            # for 1D and 2D histogram, just count(*)

            # get all the counts for each predicate
            self.domain_hist = self.get_cond_counts()

            # for i in range(0, len(self.domain_hist)):
            #     print(i, "\t", self.domain_hist[i])

            if self.index.__name__ in ['qw_1', 'qi_3', 'qt_1', 'qw_4', 'qi_2', 'qt_3']:  # 1D/2D histogram
                is_hist = True
            else:
                is_hist = False
                assert self.index.__name__ in ['qw_2', 'qw_3', 'qi_1', 'qi_4']  # 1D prefix

            if is_hist:
                # for histogram query
                count_row_by_predicates = sum(self.domain_hist)
            else:
                # for prefix query
                count_list = list(self.domain_hist)
                count_row_by_predicates = self.domain_hist[-1]
                for i in range(1, len(self.domain_hist)):
                    self.domain_hist[i] = count_list[i] - count_list[i - 1]
                    # print("after minus:\t", i, "\t", self.domain_hist[i])
                    assert self.domain_hist[i] >= 0

            # now the domain_hist represents non-overlapping sets
            count_remaining = self.cardinality - count_row_by_predicates
            if count_remaining > 0:
                self.domain_hist.append(count_remaining)
            self.query_matrix = np.zeros((len(self.cond_list), len(self.domain_hist)))

            if is_hist:
                for i in range(0, len(self.cond_list)):
                    self.query_matrix[i][i] = 1
            else:
                for i in range(0, len(self.cond_list)):
                    for j in range(0, i+1):
                        self.query_matrix[i][j] = 1

        # cache the results
        self.Wx_cache_W[self.query_key] = self.domain_hist
        self.Wx_cache_x[self.query_key] = self.query_matrix
        self.Wx_cache_N[self.query_key] = self.cardinality

    # return the sensitivity of matrix
    def get_sensitivity(self):
        return max([sum(a) for a in self.query_matrix])

    # to count f1 for ICQ and TCQ, max error for WCQ
    def get_accuracy(self):

        if self.query_type == Type.QueryType.WCQ:
            # print("DEBUG: max_error", max([abs(crnt_answer_diff) for crnt_answer_diff in self.answer_diff]))
            max_errpr = 1.0 - max([abs(crnt_answer_diff) for crnt_answer_diff in self.answer_diff]) / self.cardinality
            return [max_errpr, 'N/A']

        true_list = [0] * len(self.cond_list)
        pred_list = [0] * len(self.cond_list)

        c = 0
        if self.query_type == Type.QueryType.ICQ:
            c = self.icq_c
            for i in self.selected_cond_index:
                pred_list[i] = 1
            for i in range(0, len(self.cond_list)):
                if self.true_answer[i] > self.icq_c:
                    true_list[i] = 1

        elif self.query_type == Type.QueryType.TCQ:
            for i in self.selected_arg_list:
                pred_list[i] = 1

            cp_true = list(self.true_answer)
            cp_true.sort(reverse=True)
            c = cp_true[self.tcq_k - 1]

            for i in range(0, len(self.cond_list)):
                if self.true_answer[i] >= c:
                    true_list[i] = 1

        # return sum(true_list), sum(pred_list), f1_score(true_list, pred_list)
        f1 = f1_score(true_list, pred_list)

        alpha_hat = 0
        for i in range(0, len(self.cond_list)):
            if true_list[i] != pred_list[i]:
                alpha_hat = max(alpha_hat, abs(c - self.true_answer[i]))

        error = alpha_hat / self.cardinality

        return [error, f1]

max_stars_repo_name: cgebest/APEx

max_issues_repo_name: cgebest/APEx

max_forks_repo_name: cgebest/APEx


text:
from __future__ import print_function
from icecube import dataclasses

from icecube.load_pybindings import load_pybindings
load_pybindings(__name__,__path__)

#helper function to build the tau_dnde vector needed /fhl
def tau_dnde_builder(it):
    tmp = []
    for i in it:
        if len(i) != 2:
            raise ValueError("All elements must have exactly two parameters (tau and dN/dE)!")
        if i[0] <= 0:
            raise ValueError("tau must be greater than 0")
        tmp.append(dataclasses.make_pair(i[0],i[1]))
    return dataclasses.I3VectorDoubleDouble(tmp)

def MakeCLSimPropagator(DetectorParams, UseCPUs=False, UseGPUs=True,
                        CopyConfig=False):
    """
    Configure ppc to for use with clsim.

    :param DetectorParams: configuration dict returned by clsim.traysegments.common.setupDetector()
    :param UseCPUs: use CPUs
    :param UseGPUs: use GPUs
    :param CopyConfig: copy the cfg.txt file and update first 2 rows instead of
    generating it from scratch
    """
    from icecube.icetray import I3Units
    from icecube import clsim
    from os.path import expandvars, join, isfile
    from os import environ
    import tempfile, shutil
    import numpy

    propagator = I3CLSimStepToPhotonConverterPPC.instance
    propagator.SetGeometry(DetectorParams['Geometry'])
    propagator.SetWlenBias(DetectorParams['WavelengthGenerationBias'])

    tmpdir = tempfile.mkdtemp()

    if CopyConfig:
        shutil.copy(join(DetectorParams['IceModelLocation'], 'cfg.txt'), tmpdir)
        with open(join(tmpdir, 'cfg.txt'), 'r') as f:
            cfg = f.readlines()
        if cfg[0][0] == '#':
            first = 1
        else:
            first = 0
        cfg[first] = str(DetectorParams['DOMOversizeFactor']) + '\n'
        cfg[first + 1] = str(DetectorParams['UnshadowedFraction']*DetectorParams['MediumProperties'].efficiency) + '\n'
        with open(join(tmpdir, 'cfg.txt'), 'w') as f:
            f.writelines(cfg)
    else:
        # create config file from scratch
        with open(join(tmpdir, 'cfg.txt'), 'w') as f:
            write = lambda v: f.write('{}\n'.format(v))
            write(DetectorParams['DOMOversizeFactor'])
            write(DetectorParams['UnshadowedFraction']*DetectorParams['MediumProperties'].efficiency)

            # FIXME: extract scattering parameters from ice model
            icemodel = DetectorParams['MediumProperties']
            icemodel.GetScatteringCosAngleDistribution()
            write(0.35) # 0=HG; 1=SAM
            write(0.9) # g=<cos(theta)>

    shutil.copy(expandvars('$I3_BUILD/ppc/resources/ice/rnd.txt'), tmpdir)

    # configure unit angular acceptance, disabling angular downsampling in PPC
    with open(join(tmpdir, 'as.dat'), 'w') as f:
        for _ in range(2):
            f.write('{}\n'.format(1.0))

    # convert generation bias to cumulative distribution of photon wavelengths
    with open(join(tmpdir, 'wv.dat'), 'w') as f:
        # integrate Cherenkov spectrum weighted with given bias
        bias = DetectorParams['WavelengthGenerationBias']
        nPhase = DetectorParams['MediumProperties'].GetPhaseRefractiveIndex(0)
        def cherenkovDist(wlen, beta=1):
            # dN/dxdwlen
            return (2.*numpy.pi/(137.*(wlen*wlen)))*(1. - 1./ ( (beta*nPhase.GetValue(wlen))**2 ) )
        wavelength = numpy.array(list(map(bias.GetEntryWavelength, range(bias.GetNumEntries()))))
        value = numpy.array(list(map(cherenkovDist, wavelength)))*numpy.array(list(map(bias.GetEntryValue, range(bias.GetNumEntries()))))
        f.write('{:.8f} {:.0f}\n'.format(0, (wavelength[0] - bias.GetWavelengthStepping()/2)/I3Units.nanometer))
        for wvl, v in zip(wavelength + bias.GetWavelengthStepping()/2, numpy.cumsum(value)/sum(value)):
            f.write('{:.8f} {:.0f}\n'.format(v, wvl/I3Units.nanometer))

    # copy the ice model itself over blindly
    shutil.copy(join(DetectorParams['IceModelLocation'], 'icemodel.par'), tmpdir)
    shutil.copy(join(DetectorParams['IceModelLocation'], 'icemodel.dat'), tmpdir)

    for ext in 'par', 'dat':
        fname = join(DetectorParams['IceModelLocation'], 'tilt.'+ext)
        if isfile(fname):
            shutil.copy(fname, tmpdir)

    environ['PPCTABLESDIR'] = tmpdir
    if UseCPUs and not UseGPUs:
        environ['OCPU'] = '1'
    elif UseGPUs and not UseCPUs:
        environ['OGPU'] = '1'
    elif not UseGPUs and not UseCPUs:
        raise ValueError("Need to use at least one device type")

    try:
        propagator.Initialize()
    finally:
        shutil.rmtree(tmpdir)
    return propagator

ext: py

reason: import numpy


text:
import numpy
import pytest
from matchms import Scores

class DummySimilarityFunction:
    def __init__(self):
        """constructor"""

    def __call__(self, reference, query):
        """call method"""
        s = reference + query
        return s, len(s)

class DummySimilarityFunctionParallel:
    def __init__(self):
        """constructor"""

    def __call__(self, references, queries):
        """call method"""
        shape = len(references), len(queries)
        s = numpy.empty(shape, dtype="object")
        for index_reference, reference in enumerate(references):
            for index_query, query in enumerate(queries):
                rq = reference + query
                s[index_reference, index_query] = rq, len(rq)
        return s

def test_scores_calculate():
    dummy_similarity_function = DummySimilarityFunction()
    scores = Scores(references=["r0", "r1", "r2"],
                    queries=["q0", "q1"],
                    similarity_function=dummy_similarity_function)
    scores.calculate()
    actual = list(scores)
    expected = [
        ("r0", "q0", "r0q0", 4),
        ("r0", "q1", "r0q1", 4),
        ("r1", "q0", "r1q0", 4),
        ("r1", "q1", "r1q1", 4),
        ("r2", "q0", "r2q0", 4),
        ("r2", "q1", "r2q1", 4)
    ]
    assert actual == expected

def test_scores_calculate_parallel():
    dummy_similarity_function = DummySimilarityFunctionParallel()
    scores = Scores(references=["r0", "r1", "r2"],
                    queries=["q0", "q1"],
                    similarity_function=dummy_similarity_function)
    scores.calculate_parallel()
    actual = list(scores)
    expected = [
        ("r0", "q0", "r0q0", 4),
        ("r0", "q1", "r0q1", 4),
        ("r1", "q0", "r1q0", 4),
        ("r1", "q1", "r1q1", 4),
        ("r2", "q0", "r2q0", 4),
        ("r2", "q1", "r2q1", 4)
    ]
    assert actual == expected

def test_scores_init_with_list():

    dummy_similarity_function = DummySimilarityFunction()
    scores = Scores(references=["r0", "r1", "r2"],
                    queries=["q0", "q1"],
                    similarity_function=dummy_similarity_function)
    assert scores.scores.shape == (3, 2)

def test_scores_init_with_numpy_array():

    dummy_similarity_function = DummySimilarityFunction()
    scores = Scores(references=numpy.asarray(["r0", "r1", "r2"]),
                    queries=numpy.asarray(["q0", "q1"]),
                    similarity_function=dummy_similarity_function)
    assert scores.scores.shape == (3, 2)

def test_scores_init_with_queries_dict():

    dummy_similarity_function = DummySimilarityFunction()
    with pytest.raises(AssertionError) as msg:
        _ = Scores(references=["r0", "r1", "r2"],
                   queries=dict(k0="q0", k1="q1"),
                   similarity_function=dummy_similarity_function)

    assert str(msg.value) == "Expected input argument 'queries' to be list or tuple or numpy.ndarray."

def test_scores_init_with_references_dict():

    dummy_similarity_function = DummySimilarityFunction()
    with pytest.raises(AssertionError) as msg:
        _ = Scores(references=dict(k0="r0", k1="r1", k2="r2"),
                   queries=["q0", "q1"],
                   similarity_function=dummy_similarity_function)

    assert str(msg.value) == "Expected input argument 'references' to be list or tuple or numpy.ndarray."

def test_scores_init_with_tuple():

    dummy_similarity_function = DummySimilarityFunction()
    scores = Scores(references=("r0", "r1", "r2"),
                    queries=("q0", "q1"),
                    similarity_function=dummy_similarity_function)
    assert scores.scores.shape == (3, 2)

def test_scores_next():

    dummy_similarity_function = DummySimilarityFunction()
    scores = Scores(references=["r", "rr", "rrr"],
                    queries=["q", "qq"],
                    similarity_function=dummy_similarity_function).calculate()

    actual = list(scores)
    expected = [
        ("r", "q", "rq", 2),
        ("r", "qq", "rqq", 3),
        ("rr", "q", "rrq", 3),
        ("rr", "qq", "rrqq", 4),
        ("rrr", "q", "rrrq", 4),
        ("rrr", "qq", "rrrqq", 5)
    ]
    assert actual == expected


text:
"""
Provides the main Move class which allows definition of moves
which alter the positions of subsets of atoms in a context during a BLUES
simulation, in order to increase sampling.
Also provides functionality for CombinationMove definitions which consist of
a combination of other pre-defined moves such as via instances of Move.

Authors: Samuel C. Gill

Contributors: Nathan M. Lim, Kalistyn Burley, David L. Mobley
"""

import copy
import math
import random
import sys
import traceback

import mdtraj
import numpy
import parmed
from simtk import unit

try:
    import openeye.oechem as oechem
    if not oechem.OEChemIsLicensed():
        print('ImportError: Need License for OEChem! SideChainMove class will be unavailable.')
    try:
        import oeommtools.utils as oeommtools
    except ImportError:
        print('ImportError: Could not import oeommtools. SideChainMove class will be unavailable.')
except ImportError:
    print('ImportError: Could not import openeye-toolkits. SideChainMove class will be unavailable.')

class Move(object):
    """This is the base Move class. Move provides methods for calculating properties
    and applying the move on the set of atoms being perturbed in the NCMC simulation.
    """

    def __init__(self):
        """Initialize the Move object
        Currently empy.
        """

    def initializeSystem(self, system, integrator):
        """If the system or integrator needs to be modified to perform the move
        (ex. adding a force) this method is called during the start
        of the simulation to change the system or integrator to accomodate that.

        Parameters
        ----------
        system : openmm.System
            System to be modified.
        integrator : openmm.Integrator
            Integrator to be modified.

        Returns
        -------
        system : openmm.System
            The modified System object.
        integrator : openmm.Integrator
            The modified Integrator object.

        """
        new_sys = system
        new_int = integrator
        return new_sys, new_int

    def beforeMove(self, context):
        """This method is called at the start of the NCMC portion if the
        context needs to be checked or modified before performing the move
        at the halfway point.

        Parameters
        ----------
        context : openmm.Context
            Context containing the positions to be moved.

        Returns
        -------
        context : openmm.Context
            The same input context, but whose context were changed by this function.

        """
        return context

    def afterMove(self, context):
        """This method is called at the end of the NCMC portion if the
        context needs to be checked or modified before performing the move
        at the halfway point.

        Parameters
        ----------
        context : openmm.Context
            Context containing the positions to be moved.

        Returns
        -------
        context : openmm.Context
            The same input context, but whose context were changed by this function.

        """

        return context

    def _error(self, context):
        """This method is called if running during NCMC portion results
        in an error. This allows portions of the context, such as the
        context parameters that would not be fixed by just reverting the
        positions/velocities of the context.

        Parameters
        ----------
        context : openmm.Context
            Context containing the positions to be moved.

        Returns
        -------
        context : openmm.Context
            The same input context, but whose context were changed by this function.

        """

        return context

    def move(self, context):
        """This method is called at the end of the NCMC portion if the
        context needs to be checked or modified before performing the move
        at the halfway point.

        Parameters
        ----------
        context : openmm.Context
            Context containing the positions to be moved.

        Returns
        -------
        context : openmm.Context
            The same input context, but whose context were changed by this function.
        """
        return context

class RandomLigandRotationMove(Move):
    """RandomLightRotationMove that provides methods for calculating properties on the
    object 'model' (i.e ligand) being perturbed in the NCMC simulation.
    Current methods calculate the object's atomic masses and center of masss.
    Calculating the object's center of mass will get the positions
    and total mass.

    Parameters
    ----------
    resname : str
        String specifying the residue name of the ligand.
    structure: parmed.Structure
        ParmEd Structure object of the relevant system to be moved.
    random_state : integer or numpy.RandomState, optional
        The generator used for random numbers. If an integer is given, it fixes the seed. Defaults to the global numpy random number generator.

    Attributes
    ----------
    structure : parmed.Structure
        The structure of the ligand or selected atoms to be rotated.
    resname : str, default='LIG'
        The residue name of the ligand or selected atoms to be rotated.
    topology : openmm.Topology
        The topology of the ligand or selected atoms to be rotated.
    atom_indices : list
        Atom indicies of the ligand.
    masses : list
        Particle masses of the ligand with units.
    totalmass : int
        Total mass of the ligand.
    center_of_mass : numpy.array
        Calculated center of mass of the ligand in XYZ coordinates. This should
        be updated every iteration.
    positions : numpy.array
        Ligands positions in XYZ coordinates. This should be updated
        every iteration.

    Examples
    --------
    >>> from blues.move import RandomLigandRotationMove
    >>> ligand = RandomLigandRotationMove(structure, 'LIG')
    >>> ligand.resname
        'LIG'
    """

    def __init__(self, structure, resname='LIG', random_state=None):
        self.structure = structure
        self.resname = resname
        self.random_state = random_state
        self.atom_indices = self.getAtomIndices(structure, self.resname)
        self.topology = structure[self.atom_indices].topology
        self.totalmass = 0
        self.masses = []

        self.center_of_mass = None
        self.positions = structure[self.atom_indices].positions
        self._calculateProperties()

    def getAtomIndices(self, structure, resname):
        """
        Get atom indices of a ligand from ParmEd Structure.

        Parameters
        ----------
        resname : str
            String specifying the residue name of the ligand.
        structure: parmed.Structure
            ParmEd Structure object of the atoms to be moved.

        Returns
        -------
        atom_indices : list of ints
            list of atoms in the coordinate file matching lig_resname
        """
        # TODO: Add option for resnum to better select residue names
        atom_indices = []
        topology = structure.topology
        for atom in topology.atoms():
            if str(resname) in atom.residue.name:
                atom_indices.append(atom.index)
        return atom_indices

    def getMasses(self, topology):
        """
        Returns a list of masses of the specified ligand atoms.

        Parameters
        ----------
        topology: parmed.Topology
            ParmEd topology object containing atoms of the system.

        Returns
        -------
        masses: 1xn numpy.array * simtk.unit.dalton
            array of masses of len(self.atom_indices), denoting
            the masses of the atoms in self.atom_indices
        totalmass: float * simtk.unit.dalton
            The sum of the mass found in masses
        """
        masses = unit.Quantity(numpy.zeros([int(topology.getNumAtoms()), 1], numpy.float32), unit.dalton)
        for idx, atom in enumerate(topology.atoms()):
            masses[idx] = atom.element._mass
        totalmass = masses.sum()
        return masses, totalmass

    def getCenterOfMass(self, positions, masses):
        """Returns the calculated center of mass of the ligand as a numpy.array

        Parameters
        ----------
        positions: nx3 numpy array * simtk.unit compatible with simtk.unit.nanometers
            ParmEd positions of the atoms to be moved.
        masses : numpy.array
            numpy.array of particle masses

        Returns
        -------
        center_of_mass: numpy array * simtk.unit compatible with simtk.unit.nanometers
            1x3 numpy.array of the center of mass of the given positions
        """
        coordinates = numpy.asarray(positions._value, numpy.float32)
        center_of_mass = parmed.geometry.center_of_mass(coordinates, masses) * positions.unit
        return center_of_mass

    def _calculateProperties(self):
        """Calculate the masses and center of mass for the object. This function
        is called upon initailization of the class."""
        self.masses, self.totalmass = self.getMasses(self.topology)
        self.center_of_mass = self.getCenterOfMass(self.positions, self.masses)

    def move(self, context):
        """Function that performs a random rotation about the
        center of mass of the ligand.

        Parameters
        ----------
        context: simtk.openmm.Context object
            Context containing the positions to be moved.

        Returns
        -------
        context: simtk.openmm.Context object
            The same input context, but whose positions were changed by this function.
        """
        positions = context.getState(getPositions=True).getPositions(asNumpy=True)

        self.positions = positions[self.atom_indices]
        self.center_of_mass = self.getCenterOfMass(self.positions, self.masses)
        reduced_pos = self.positions - self.center_of_mass

        # Define random rotational move on the ligand
        rand_quat = mdtraj.utils.uniform_quaternion(size=None, random_state=self.random_state)
        rand_rotation_matrix = mdtraj.utils.rotation_matrix_from_quaternion(rand_quat)
        #multiply lig coordinates by rot matrix and add back COM translation from origin
        rot_move = numpy.dot(reduced_pos, rand_rotation_matrix) * positions.unit + self.center_of_mass

        # Update ligand positions in nc_sim
        for index, atomidx in enumerate(self.atom_indices):
            positions[atomidx] = rot_move[index]
        context.setPositions(positions)
        positions = context.getState(getPositions=True).getPositions(asNumpy=True)
        self.positions = positions[self.atom_indices]
        return context

class MoveEngine(object):
    """MoveEngine provides perturbation functions for the context during the NCMC
    simulation.

    Parameters
    ----------
    moves : blues.move object or list of N blues.move objects
        Specifies the possible moves to be performed.

    probabilities: list of floats, optional, default=None
        A list of N probabilities, where probabilities[i] corresponds to the
        probaility of moves[i] being selected to perform its associated
        move() method. If None, uniform probabilities are assigned.

    Attributes
    ----------
    moves : blues.move or list of N blues.move objects
        Possible moves to be performed.
    probabilities : list of floats
        Normalized probabilities for each move.
    selected_move : blues.move
        Selected move to be performed.
    move_name : str
        Name of the selected move to be performed

    Examples
    --------
    Load a parmed.Structure, list of moves with probabilities, initialize
    the MoveEngine class, and select a move from our list.

    >>> import parmed
    >>> from blues.moves import *
    >>> structure = parmed.load_file('tests/data/eqToluene.prmtop', xyz='tests/data/eqToluene.inpcrd')
    >>> probabilities = [0.25, 0.75]
    >>> moves = [SideChainMove(structure, [111]),RandomLigandRotationMove(structure, 'LIG')]
    >>> mover = MoveEngine(moves, probabilities)
    >>> mover.moves
    [<blues.moves.SideChainMove at 0x7f2eaa168470>,
     <blues.moves.RandomLigandRotationMove at 0x7f2eaaaa51d0>]
    >>> mover.selectMove()
    >>> mover.selected_move
    <blues.moves.RandomLigandRotationMove at 0x7f2eaaaa51d0>
    """

    def __init__(self, moves, probabilities=None):
        #make a list from moves if not a list
        if isinstance(moves, list):
            self.moves = moves
        else:
            self.moves = [moves]
        #normalize probabilities
        if probabilities is None:
            single_prob = 1. / len(self.moves)
            self.probabilities = [single_prob for x in (self.moves)]
        else:
            prob_sum = float(sum(probabilities))
            self.probabilities = [x / prob_sum for x in probabilities]
        #if move and probabilitiy lists are different lengths throw error
        if len(self.moves) != len(self.probabilities):
            print('moves and probability list lengths need to match')
            raise IndexError
        #use index in selecting move
        self.selected_move = None

    def selectMove(self):
        """Chooses the move which will be selected for a given NCMC
        iteration
        """
        rand_num = numpy.random.choice(len(self.probabilities), p=self.probabilities)
        self.selected_move = self.moves[rand_num]
        self.move_name = self.selected_move.__class__.__name__

    def runEngine(self, context):
        """Selects a random Move object based on its
        assigned probability and and performs its move() function
        on a context.

        Parameters
        ----------
        context : openmm.Context
            OpenMM context whose positions should be moved.

        Returns
        -------
        context : openmm.Context
            OpenMM context whose positions have been moved.
        """
        try:
            new_context = self.selected_move.move(context)
        except Exception as e:
            #In case the move isn't properly implemented, print out useful info
            print('Error: move not implemented correctly, printing traceback:')
            ex_type, ex, tb = sys.exc_info()
            traceback.print_tb(tb)
            print(e)
            raise SystemExit

        return new_context

########################
## UNDER DEVELOPMENT ###
########################

class SideChainMove(Move):
    """**NOTE:** Usage of this class requires a valid OpenEye license.

    SideChainMove provides methods for calculating properties needed to
    rotate a sidechain residue given a parmed.Structure. Calculated properties
    include: backbone atom indicies, atom pointers and indicies of the residue
    sidechain, bond pointers and indices for rotatable heavy bonds in
    the sidechain, and atom indices upstream of selected bond.

    The class contains functions to randomly select a bond and angle to be rotated
    and applies a rotation matrix to the target atoms to update their coordinates on the
    object 'model' (i.e sidechain) being perturbed in the NCMC simulation.

    Parameters
    ----------
    structure : parmed.Structure
        The structure of the entire system to be simulated.
    residue_list : list of int
        List of the residue numbers of the sidechains to be rotated.
    verbose : bool, default=False
        Enable verbosity to print out detailed information of the rotation.
    write_move : bool, default=False
        If True, writes a PDB of the system after rotation.

    Attributes
    ----------
    structure : parmed.Structure
        The structure of the entire system to be simulated.
    molecule : oechem.OEMolecule
        The OEMolecule containing the sidechain(s) to be rotated.
    residue_list : list of int
        List containing the residue numbers of the sidechains to be rotated.
    all_atoms : list of int
        List containing the atom indicies of the sidechains to be rotated.
    rot_atoms : dict
        Dictionary of residues, bonds and atoms to be rotated
    rot_bonds : dict of oechem.OEBondBase
        Dictionary containing the bond pointers of the rotatable bonds.
    qry_atoms : dict of oechem.OEAtomBase
        Dictionary containing all the atom pointers (as OpenEye objects) that
        make up the given residues.

    Examples
    --------
    >>> from blues.move import SideChainMove
    >>> sidechain = SideChainMove(structure, [1])

    """

    def __init__(self, structure, residue_list, verbose=False, write_move=False):
        self.structure = structure
        self.molecule = self._pmdStructureToOEMol()
        self.residue_list = residue_list
        self.all_atoms = [atom.index for atom in self.structure.topology.atoms()]
        self.rot_atoms, self.rot_bonds, self.qry_atoms = self.getRotBondAtoms()
        self.atom_indices = self.rot_atoms
        self.verbose = verbose
        self.write_move = write_move

    def _pmdStructureToOEMol(self):
        """Helper function for converting the parmed structure into an OEMolecule."""
        top = self.structure.topology
        pos = self.structure.positions
        molecule = oeommtools.openmmTop_to_oemol(top, pos, verbose=False)
        oechem.OEPerceiveResidues(molecule)
        oechem.OEFindRingAtomsAndBonds(molecule)

        return molecule

    def getBackboneAtoms(self, molecule):
        """Takes an OpenEye Molecule, finds the backbone atoms and
        returns the indicies of the backbone atoms.

        Parameters
        ----------
        molecule : oechem.OEMolecule
            The OEmolecule of the simulated system.

        Returns
        -------
        backbone_atoms : list of int
            List containing the atom indices of the backbone atoms.

        """

        backbone_atoms = []
        pred = oechem.OEIsBackboneAtom()
        for atom in molecule.GetAtoms(pred):
            bb_atom_idx = atom.GetIdx()
            backbone_atoms.append(bb_atom_idx)

        return backbone_atoms

    def getTargetAtoms(self, molecule, backbone_atoms, residue_list):
        """Takes an OpenEye molecule and a list of residue numbers then
        generates a dictionary containing all the atom pointers and indicies for the
        non-backbone, atoms of those target residues, as well as a list of backbone atoms.
        Note: The atom indicies start at 0 and are thus -1 from the PDB file indicies

        Parameters
        ----------
        molecule : oechem.OEMolecule
            The OEmolecule of the simulated system.
        backbone_atoms : list of int
            List containing the atom indices of the backbone atoms.
        residue_list : list of int
            List containing the residue numbers of the sidechains to be rotated.

        Returns
        -------
        backbone_atoms : list of int
            List containing the atom indices of the backbone atoms to be rotated.
        qry_atoms : dict of oechem.OEAtomBase
            Dictionary containing all the atom pointers (as OpenEye objects) that
            make up the given residues.

        """

        # create and clear dictionary to store atoms that make up residue list
        qry_atoms = {}
        qry_atoms.clear()

        reslib = []

        #print('Searching residue list for atoms...')
        # loop through all the atoms in the PDB OEGraphMol structure
        for atom in molecule.GetAtoms():
            # check if the atom is in backbone
            if atom.GetIdx() not in backbone_atoms:
                # if heavy, find what residue it is associated with
                myres = oechem.OEAtomGetResidue(atom)
                # check if the residue number is amongst the list of residues
                if myres.GetResidueNumber() in residue_list and myres.GetName() != "HOH":
                    # store the atom location in a query atom dict keyed by its atom index
                    qry_atoms.update({atom: atom.GetIdx()})
                    #print('Found atom %s in residue number %i %s'%(atom,myres.GetResidueNumber(),myres.GetName()))
                    if myres not in reslib:
                        reslib.append(myres)

        return qry_atoms, backbone_atoms

    def findHeavyRotBonds(self, pdb_OEMol, qry_atoms):
        """Takes in an OpenEye molecule as well as a dictionary of atom locations (keys)
        and atom indicies.  It loops over the query atoms and identifies any heavy bonds associated with each atom.
        It stores and returns the bond indicies (keys) and the two atom indicies for each bond in a dictionary
        Note: atom indicies start at 0, so are offset by 1 compared to pdb)

        Parameters
        ----------
        pdb_OEMol : oechem.OEMolecule
            The OEmolecule of the simulated system generated from a PDB file.
        qry_atoms : dict of oechem.OEAtomBase
            Dictionary containing all the atom pointers (as OpenEye objects) that
            make up the given residues.

        Returns
        -------
        rot_bonds : dict of oechem.OEBondBase
            Dictionary containing the bond pointers of the rotatable bonds.

        """
        # create and clear dictionary to store bond and atom indicies that are rotatable + heavy
        rot_bonds = {}
        rot_bonds.clear()

        for atom in qry_atoms.keys():
            myres = oechem.OEAtomGetResidue(atom)
            for bond in atom.GetBonds():
                # retrieve the begnning and ending atoms
                begatom = bond.GetBgn()
                endatom = bond.GetEnd()
                # if begnnning and ending atoms are not Hydrogen, and the bond is rotatable
                if endatom.GetAtomicNum() > 1 and begatom.GetAtomicNum() > 1 and bond.IsRotor():
                    # if the bond has not been added to dictionary already..
                    # (as would happen if one of the atom pairs was previously looped over)
                    if bond not in rot_bonds:
                        #print('Bond number',bond, 'is rotatable, non-terminal, and contains only heavy atoms')
                        # store bond pointer (key) and atom indicies in dictionary if not already there
                        #rot_bonds.update({bond : {'AtomIdx_1' : bond.GetBgnIdx(), 'AtomIdx_2': bond.GetEndIdx()}})
                        rot_bonds.update({bond: myres.GetResidueNumber()})

        return rot_bonds

    def getRotAtoms(self, rotbonds, molecule, backbone_atoms):
        """Function identifies and stores neighboring, upstream atoms for a given sidechain bond.

        Parameters
        ----------
        rot_bonds : dict of oechem.OEBondBase
            Dictionary containing the bond pointers of the rotatable bonds.
        molecule : oechem.OEMolecule
            The OEmolecule of the simulated system.
        backbone_atoms : list of int
            List containing the atom indices of the backbone atoms.

        Returns
        -------
        rot_atom_dict : dict of oechem.OEAtomBase
            Dictionary containing the atom pointers for a given sidechain bond.

        """
        backbone = backbone_atoms
        query_list = []
        idx_list = []
        rot_atom_dict = {}
        rot_atom_dict.clear()

        for bond in rotbonds.keys():
            idx_list.clear()
            query_list.clear()
            resnum = (rotbonds[bond])
            thisbond = bond
            ax1 = bond.GetBgn()
            ax2 = bond.GetEnd()

            if resnum in rot_atom_dict.keys():
                rot_atom_dict[resnum].update({thisbond: []})
            else:
                rot_atom_dict.update({resnum: {thisbond: []}})

            idx_list.append(ax1.GetIdx())
            idx_list.append(ax2.GetIdx())

            if ax1 not in query_list and ax1.GetIdx() not in backbone_atoms:
                query_list.append(ax1)
            if ax2 not in query_list and ax2.GetIdx() not in backbone_atoms:
                query_list.append(ax2)

            for atom in query_list:
                checklist = atom.GetAtoms()
                for candidate in checklist:
                    if candidate not in query_list and candidate.GetIdx() not in backbone and candidate != ax2:
                        query_list.append(candidate)
                        if candidate.GetAtomicNum() > 1:
                            can_nbors = candidate.GetAtoms()
                            for can_nbor in can_nbors:
                                if can_nbor not in query_list and candidate.GetIdx(
                                ) not in backbone and candidate != ax2:
                                    query_list.append(can_nbor)

            for atm in query_list:
                y = atm.GetIdx()
                if y not in idx_list:
                    idx_list.append(y)

            rot_atom_dict[resnum].update({thisbond: list(idx_list)})
            #print("Moving these atoms:", idx_list)

        return rot_atom_dict

    def getRotBondAtoms(self):
        """This function is called on class initialization.

        Takes in a PDB filename (as a string) and list of residue numbers.  Returns
        a nested dictionary of rotatable bonds (containing only heavy atoms), that are keyed by residue number,
        then keyed by bond pointer, containing values of atom indicies [axis1, axis2, atoms to be rotated]
        Note: The atom indicies start at 0, and are offset by -1 from the PDB file indicies

        Returns
        -------
        rot_atoms : dict
            Dictionary of residues, bonds and atoms to be rotated
        rot_bonds : dict of oechem.OEBondBase
            Dictionary containing the bond pointers of the rotatable bonds.
        qry_atoms : dict of oechem.OEAtomBase
            Dictionary containing all the atom pointers (as OpenEye objects) that
            make up the given residues.

        """
        backbone_atoms = self.getBackboneAtoms(self.molecule)

        # Generate dictionary containing locations and indicies of heavy residue atoms
        #print('Dictionary of all query atoms generated from residue list\n')
        qry_atoms, backbone_atoms = self.getTargetAtoms(self.molecule, backbone_atoms, self.residue_list)

        # Identify bonds containing query atoms and return dictionary of indicies
        rot_bonds = self.findHeavyRotBonds(self.molecule, qry_atoms)

        # Generate dictionary of residues, bonds and atoms to be rotated
        rot_atoms = self.getRotAtoms(rot_bonds, self.molecule, backbone_atoms)
        return rot_atoms, rot_bonds, qry_atoms

    def chooseBondandTheta(self):
        """This function is called on class initialization.

        Takes a dictionary containing nested dictionary, keyed by res#,
        then keyed by bond_ptrs, containing a list of atoms to move, randomly selects a bond,
        and generates a random angle (radians).  It returns the atoms associated with the
        the selected bond, the pointer for the selected bond and the randomly generated angle

        Returns
        -------
        theta_ran :

        targetatoms :

        res_choice :

        bond_choice :

        """

        res_choice = random.choice(list(self.rot_atoms.keys()))
        bond_choice = random.choice(list(self.rot_atoms[res_choice].keys()))
        targetatoms = self.rot_atoms[res_choice][bond_choice]
        theta_ran = random.random() * 2 * math.pi

        return theta_ran, targetatoms, res_choice, bond_choice

    def rotation_matrix(self, axis, theta):
        """Function returns the rotation matrix associated with counterclockwise rotation
        about the given axis by theta radians.

        Parameters
        ----------
        axis :

        theta : float
            The angle of rotation in radians.
        """
        axis = numpy.asarray(axis)
        axis = axis / math.sqrt(numpy.dot(axis, axis))
        a = math.cos(theta / 2.0)
        b, c, d = -axis * math.sin(theta / 2.0)
        aa, bb, cc, dd = a * a, b * b, c * c, d * d
        bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
        return numpy.array([[aa + bb - cc - dd, 2 * (bc + ad),
                             2 * (bd - ac)], [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                            [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

    def move(self, context, verbose=False):
        """Rotates the target atoms around a selected bond by angle theta and updates
        the atom coordinates in the parmed structure as well as the ncmc context object

        Parameters
        ----------
        context: simtk.openmm.Context object
            Context containing the positions to be moved.
        verbose : bool, default=False
            Enable verbosity to print out detailed information of the rotation.

        Returns
        -------
        context: simtk.openmm.Context object
            The same input context, but whose positions were changed by this function.

        """

        # determine the axis, theta, residue, and bond + atoms to be rotated
        theta, target_atoms, res, bond = self.chooseBondandTheta()
        print('Rotating bond: %s in resnum: %s by %.2f radians' % (bond, res, theta))

        #retrieve the current positions
        initial_positions = context.getState(getPositions=True).getPositions(asNumpy=True)
        nc_positions = copy.deepcopy(initial_positions)

        model = copy.copy(self.structure)

        # set the parmed model to the same coordinates as the context
        for idx, atom in enumerate(self.all_atoms):
            if self.verbose:
                print('Before:')
                print(atom, idx)
                print(nc_positions[atom], model.positions[atom])

            model.atoms[atom].xx = nc_positions[atom][0].value_in_unit(unit.angstroms)
            model.atoms[atom].xy = nc_positions[atom][1].value_in_unit(unit.angstroms)
            model.atoms[atom].xz = nc_positions[atom][2].value_in_unit(unit.angstroms)

            if self.verbose:
                print('After:')
                print(nc_positions[atom], model.positions[atom])

        positions = model.positions

        # find the rotation axis using the updated positions
        axis1 = target_atoms[0]
        axis2 = target_atoms[1]
        rot_axis = (positions[axis1] - positions[axis2]) / positions.unit

        #calculate the rotation matrix
        rot_matrix = self.rotation_matrix(rot_axis, theta)

        # apply the rotation matrix to the target atoms
        for idx, atom in enumerate(target_atoms):

            my_position = positions[atom]

            if self.verbose:
                print('The current position for %i is: %s' % (atom, my_position))

            # find the reduced position (substract out axis)
            red_position = (my_position - model.positions[axis2])._value
            # find the new positions by multiplying by rot matrix
            new_position = numpy.dot(rot_matrix, red_position) * positions.unit + positions[axis2]

            if self.verbose: print("The new position should be:", new_position)

            positions[atom] = new_position
            # Update the parmed model with the new positions
            model.atoms[atom].xx = new_position[0] / positions.unit
            model.atoms[atom].xy = new_position[1] / positions.unit
            model.atoms[atom].xz = new_position[2] / positions.unit

            #update the copied ncmc context array with the new positions
            nc_positions[atom][0] = model.atoms[atom].xx * nc_positions.unit / 10
            nc_positions[atom][1] = model.atoms[atom].xy * nc_positions.unit / 10
            nc_positions[atom][2] = model.atoms[atom].xz * nc_positions.unit / 10

            if self.verbose:
                print('The updated position for this atom is:', model.positions[atom])

        # update the actual ncmc context object with the new positions
        context.setPositions(nc_positions)

        # update the class structure positions
        self.structure.positions = model.positions

        if self.write_move:
            filename = 'sc_move_%s_%s_%s.pdb' % (res, axis1, axis2)
            mod_prot = model.save(filename, overwrite=True)
        return context

class SmartDartMove(RandomLigandRotationMove):
    """**WARNING:** This class has not been completely tested. Use at your own risk.

    Move object that allows center of mass smart darting moves to be performed on a ligand,
    allowing translations of a ligand between pre-defined regions in space. The
    `SmartDartMove.move()` method translates the ligand to the locations of the ligand
    found in the coord_files. These locations are defined in terms of the basis_particles.
    These locations are picked with a uniform probability. Based on Smart Darting Monte Carlo [smart-dart]_

    Parameters
    ----------
    structure: parmed.Structure
        ParmEd Structure object of the relevant system to be moved.
    basis_particles: list of 3 ints
        Specifies the 3 indices of the protein whose coordinates will be used
        to define a new set of basis vectors.
    coord_files: list of str
        List containing paths to coordinate files of the whole system for smart darting.
    topology: str, optional, default=None
        A path specifying a topology file matching the files in coord_files. Not
        necessary if the coord_files already contain topologies (ex. PDBs).
    dart_radius: simtk.unit float object compatible with simtk.unit.nanometers unit,
        optional, default=0.2*simtk.unit.nanometers
        The radius of the darting region around each dart.
    self_dart: boolean, optional, default='False'
        When performing the center of mass darting in `SmartDartMove.move()`,this
        specifies whether or not to include the darting region where the center
        of mass currently resides as an option to dart to.
    resname : str, optional, default='LIG'
        String specifying the residue name of the ligand.

    References
    ----------
    .. [smart-dart] I. Andricioaei, J. E. Straub, and A. F. Voter, J. Chem. Phys. 114, 6994 (2001).
        https://doi.org/10.1063/1.1358861

    """

    def __init__(self,
                 structure,
                 basis_particles,
                 coord_files,
                 topology=None,
                 dart_radius=0.2 * unit.nanometers,
                 self_dart=False,
                 resname='LIG'):

        super(SmartDartMove, self).__init__(structure, resname=resname)

        if len(coord_files) < 2:
            raise ValueError('You should include at least two files in coord_files ' +
                             'in order to benefit from smart darting')
        self.dartboard = []
        self.n_dartboard = []
        self.particle_pairs = []
        self.particle_weights = []
        self.basis_particles = basis_particles
        self.dart_radius = dart_radius
        self.calculateProperties()
        self.self_dart = self_dart
        self.dartsFromParmEd(coord_files, topology)

    def dartsFromParmEd(self, coord_files, topology=None):
        """
        Used to setup darts from a generic coordinate file, through MDtraj using the basis_particles to define
        new basis vectors, which allows dart centers to remain consistant through a simulation.
        This adds to the self.n_dartboard, which defines the centers used for smart darting.

        Parameters
        ----------
        coord_files: list of str
            List containing coordinate files of the whole system for smart darting.
        topology: str, optional, default=None
            A path specifying a topology file matching the files in coord_files. Not
            necessary if the coord_files already contain topologies.

        """

        n_dartboard = []
        dartboard = []
        #loop over specified files and generate parmed structures from each
        #then the center of masses of the ligand in each structureare found
        #finally those center of masses are added to the `self.dartboard`s to
        #be used in the actual smart darting move to define darting regions
        for coord_file in coord_files:
            if topology == None:
                #if coord_file contains topology info, just load coord file
                temp_md = parmed.load_file(coord_file)
            else:
                #otherwise load file specified in topology
                temp_md = parmed.load_file(topology, xyz=coord_file)
            #get position values in terms of nanometers
            context_pos = temp_md.positions.in_units_of(unit.nanometers)
            lig_pos = numpy.asarray(context_pos._value)[self.atom_indices] * unit.nanometers
            particle_pos = numpy.asarray(context_pos._value)[self.basis_particles] * unit.nanometers
            #calculate center of mass of ligand
            self.calculateProperties()
            center_of_mass = self.getCenterOfMass(lig_pos, self.masses)
            #get particle positions
            new_coord = self._findNewCoord(particle_pos[0], particle_pos[1], particle_pos[2], center_of_mass)
            #old_coord should be equal to com
            old_coord = self._findOldCoord(particle_pos[0], particle_pos[1], particle_pos[2], new_coord)
            numpy.testing.assert_almost_equal(old_coord._value, center_of_mass._value, decimal=1)
            #add the center of mass in euclidian and new basis set (defined by the basis_particles)
            n_dartboard.append(new_coord)
            dartboard.append(old_coord)
        self.n_dartboard = n_dartboard
        self.dartboard = dartboard

    def move(self, context):
        """
        Function for performing smart darting move with darts that
        depend on particle positions in the system.

        Parameters
        ----------
        context: simtk.openmm.Context object
            Context containing the positions to be moved.

        Returns
        -------
        context: simtk.openmm.Context object
            The same input context, but whose positions were changed by this function.

        """

        atom_indices = self.atom_indices
        if len(self.n_dartboard) == 0:
            raise ValueError('No darts are specified. Make sure you use ' +
                             'SmartDartMove.dartsFromParmed() before using the move() function')

        #get state info from context
        stateinfo = context.getState(True, True, False, True, True, False)
        oldDartPos = stateinfo.getPositions(asNumpy=True)
        #get the ligand positions
        lig_pos = numpy.asarray(oldDartPos._value)[self.atom_indices] * unit.nanometers
        #updates the darting regions based on the current position of the basis particles
        self._findDart(context)
        #find the ligand's current center of mass position
        center = self.getCenterOfMass(lig_pos, self.masses)
        #calculate the distance of the center of mass to the center of each darting region
        selected_dart, changevec = self._calc_from_center(com=center)
        #selected_dart is the selected darting region

        #if the center of mass was within one darting region, move the ligand to another region
        if selected_dart != None:
            newDartPos = numpy.copy(oldDartPos)
            #find the center of mass in the new darting region
            dart_switch = self._reDart(selected_dart, changevec)
            #find the vector that will translate the ligand to the new darting region
            vecMove = dart_switch - center
            #apply that vector to the ligand to actually translate the coordinates
            for atom in atom_indices:
                newDartPos[atom] = newDartPos[atom] + vecMove._value
            #set the positions after darting
            context.setPositions(newDartPos)

            return context

    def _calc_from_center(self, com):
        """
        Helper function that finds the distance of the current center of
        mass to each dart center in self.dartboard

        Parameters
        --------
        com: 1x3 numpy.array*simtk.unit.nanometers
            Current center of mass coordinates of the ligand.

        Returns
        -------
        selected_dart: simtk.unit.nanometers, or None
            The distance of a dart to a center. Returns
            None if the distance is greater than the darting region.
        changevec: 1x3 numpy.array*simtk.unit.nanometers,
            The vector from the ligand center of mass
            to the center of a darting region.

        """

        distList = []
        diffList = []
        indexList = []
        #Find the distances of the COM to each dart, appending
        #the results to distList
        for dart in self.dartboard:
            diff = com - dart
            dist = numpy.sqrt(numpy.sum((diff) * (diff))) * unit.nanometers
            distList.append(dist)
            diffList.append(diff)
        selected_dart = []
        #Find the dart(s) less than self.dart_radius
        for index, entry in enumerate(distList):
            if entry <= self.dart_radius:
                selected_dart.append(index)
                diff = diffList[index]
                indexList.append(index)
        #Dart error checking
        #to ensure reversibility the COM should only be
        #within self.dart_radius of one dart
        if len(selected_dart) == 1:
            return selected_dart[0], diffList[indexList[0]]
        elif len(selected_dart) == 0:
            return None, diff
        elif len(selected_dart) >= 2:
            #COM should never be within two different darts
            raise ValueError(' The spheres defining two darting regions have overlapped, ' +
                             'which results in potential problems with detailed balance. ' +
                             'We are terminating the simulation. Please check the size and ' +
                             'identity of your darting regions defined by dart_radius.')
            #TODO can treat cases using appropriate probablility correction
            #see https://doi.org/10.1016/j.patcog.2011.02.006

    def _findDart(self, context):
        """
        Helper function to dynamically update dart positions based on the current positions
        of the basis particles.

        Parameters
        ---------
        context: Context object from simtk.openmm
            Context from the ncmc simulation.

        Returns
        -------
        dart_list list of 1x3 numpy.arrays in units.nm
            new dart positions calculated from the particle_pairs
            and particle_weights.

        """

        basis_particles = self.basis_particles
        #make sure there's an equal number of particle pair lists
        #and particle weight lists
        dart_list = []
        state_info = context.getState(True, True, False, True, True, False)
        temp_pos = state_info.getPositions(asNumpy=True)
        part1 = temp_pos[basis_particles[0]]
        part2 = temp_pos[basis_particles[1]]
        part3 = temp_pos[basis_particles[2]]
        for dart in self.n_dartboard:
            old_center = self._findOldCoord(part1, part2, part3, dart)
            dart_list.append(old_center)
        self.dartboard = dart_list[:]
        return dart_list

    def _reDart(self, selected_dart, changevec):
        """
        Helper function to choose a random dart and determine the vector
        that would translate the COM to that dart center + changevec.
        This is called reDart in the sense that it helps to switch
        the ligand to another darting region.

        Parameters
        ---------
        selected_dart :
        changevec: 1x3 numpy.array * simtk.unit.nanometers
            The vector difference of the ligand center of mass
            to the closest dart center (if within the dart region).

        Returns
        -------
        dart_switch: 1x3 numpy.array * simtk.unit.nanometers

        """
        dartindex = list(range(len(self.dartboard)))
        if self.self_dart == False:
            dartindex.pop(selected_dart)
        dartindex = numpy.random.choice(dartindex)
        dvector = self.dartboard[dartindex]
        dart_switch = dvector + changevec
        return dart_switch

    def _changeBasis(self, a, b):
        """
        Changes positions of a particle (b) in the regular basis set to
        another basis set (a). Used to recalculate the center of mass
        in terms of the local coordinates defined by self.basis_particles.
        Used to change between the basis sets defined from the basis_particles
        and the normal euclidian basis set.

        Parameters
        ----------
        a: 3x3 numpy.array
            Defines vectors that will create the new basis.
        b: 1x3 numpy.array
            Defines position of particle to be transformed into
            new basis set.

        Returns
        -------
        changed_coord: 1x3 numpy.array
            Coordinates of b in new basis.

        """

        ainv = numpy.linalg.inv(a.T)
        changed_coord = numpy.dot(ainv, b.T) * unit.nanometers
        return changed_coord

    def _undoBasis(self, a, b):
        """
        Transforms positions in a transformed basis (b) to the regular
        basis set. Used to transform the dart positions in the local
        coordinate basis set to the cartesian basis set.

        Parameters
        ----------
        a: 3x3 numpy.array
            Defines vectors that defined the new basis.
        b: 1x3 numpy.array
            Defines position of particle to be transformed into
            regular basis set.

        Returns
        -------
        changed_coord: 1x3 numpy.array
            Coordinates of b in new basis.
        """

        a = a.T
        changed_coord = numpy.dot(a, b.T) * unit.nanometers
        return changed_coord

    def _normalize(self, vector):
        """Normalize a given vector

        Parameters
        ----------
        vector: 1xn numpy.array
            Vector to be normalized.

        Returns
        -------
        unit_vec: 1xn numpy.array
            Normalized vector.

        """

        magnitude = numpy.sqrt(numpy.sum(vector * vector))
        unit_vec = vector / magnitude
        return unit_vec

    def _localCoord(self, particle1, particle2, particle3):
        """
        Defines a new coordinate system using 3 particles
        returning the new basis set vectors

        Parameters
        ----------
        particle1, particle2, particle3: 1x3 numpy.array
            numpy.array corresponding to a given particle's positions

        Returns
        -------
        vec1, vec2, vec3: 1x3 numpy.array
            Basis vectors of the coordinate system defined
            by particles1-3.

        """

        part2 = particle2 - particle1
        part3 = particle3 - particle1
        vec1 = part2
        vec2 = part3
        vec3 = numpy.cross(vec1, vec2) * unit.nanometers
        return vec1, vec2, vec3

    def _findNewCoord(self, particle1, particle2, particle3, center):
        """
        Finds the coordinates of a given center in the standard basis
            in terms of a new basis defined by particles1-3

        Parameters
        ----------
        particle1, particle2, particle3: 1x3 numpy.array
            numpy.array corresponding to a given particle's positions
        center: 1x3 numpy.array * simtk.unit compatible with simtk.unit.nanometers
            Coordinate of the center of mass in the standard basis set.

        Returns
        -------
        new_coord : numpy.array
            Updated coordinates in terms of new basis.
        """

        #calculate new basis set
        vec1, vec2, vec3 = self._localCoord(particle1, particle2, particle3)
        basis_set = numpy.zeros((3, 3)) * unit.nanometers
        basis_set[0] = vec1
        basis_set[1] = vec2
        basis_set[2] = vec3
        #since the origin is centered at particle1 by convention
        #subtract to account for this
        recenter = center - particle1
        #find coordinate in new coordinate system
        new_coord = self._changeBasis(basis_set, recenter)
        return new_coord

    def _findOldCoord(self, particle1, particle2, particle3, center):
        """
        Finds the coordinates of a given center (defined by a different basis
        given by particles1-3) back in the euclidian coordinates

        Parameters
        ----------
        particle1, particle2, particle3: 1x3 numpy.array
            numpy.array corresponding to a given particle's positions
        center: 1x3 numpy.array * simtk.unit compatible with simtk.unit.nanometers
            Coordinate of the center of mass in the non-standard basis set.

        Returns
        -------
        adjusted_center : numpy.array
            Corrected coordinates of new center in euclidian coordinates.

        """

        vec1, vec2, vec3 = self._localCoord(particle1, particle2, particle3)
        basis_set = numpy.zeros((3, 3)) * unit.nanometers
        basis_set[0] = vec1
        basis_set[1] = vec2
        basis_set[2] = vec3
        #since the origin is centered at particle1 by convention
        #subtract to account for this
        old_coord = self._undoBasis(basis_set, center)
        adjusted_center = old_coord + particle1
        return adjusted_center

class CombinationMove(Move):
    """**WARNING:** This class has not been completely tested. Use at your own risk.

    Move object that allows Move object moves to be performed according to
    the order in move_list. To ensure detailed balance, the moves have an equal
    chance to be performed in listed or reverse order.

    Parameters
    ----------
    moves : list of blues.move.Move

    """

    def __init__(self, moves):
        self.moves = moves

    def move(self, context):
        """Performs the move() functions of the Moves in move_list on
        a context.

        Parameters
        ----------
        context: simtk.openmm.Context object
            Context containing the positions to be moved.

        Returns
        -------
        context: simtk.openmm.Context object
            The same input context, but whose positions were changed by this function.

        """
        rand = numpy.random.random()
        #to maintain detailed balance this executes both
        #the forward and reverse order moves with equal probability
        if rand > 0.5:
            for single_move in self.move_list:
                single_move.move(context)
        else:
            for single_move in reverse(self.move_list):
                single_move.move(context)


text:
def extract_balise(dataset):
    """
    Extract balise gene and their associated value,
    take the mean of the column
    """

    ## importation
    import pandas as pd

    ## parameters
    target_genes = [
        "NUP210L",
        "SPIRE1",
        "GATAD1",
        "HVCN1",
        "ENO1",
        "FLNA"
    ]

    ## load dataset
    df = pd.read_csv(dataset)

    ## try to extract target genes
    try:
        df = df[target_genes]
    except:
        print("[!] ERROR, missing genes in "+str(dataset))
        return 0

    ## compute mean of target genes
    gene_to_mean = {}
    for gene in target_genes:
        gene_to_mean[gene] = df[gene].mean()

    ## return
    return gene_to_mean

def save_balise(ref_dataset, target_dataset):
    """
    Extract gene balise for ref_dataset and target_dataset,
    save results in pickle stuff
    """

    ## importation
    import pickle

    ## extract gene to mean to ref
    gene_to_mean_ref = extract_balise(ref_dataset)

    ## extract gene to mean to target
    gene_to_mean_target = extract_balise(target_dataset)

    ## save dict in picle format
    with open('ressources/gene_to_mean_ref.pickle', 'wb') as handle:
        pickle.dump(gene_to_mean_ref, handle, protocol=pickle.HIGHEST_PROTOCOL)
    with open('ressources/gene_to_mean_target.pickle', 'wb') as handle:
        pickle.dump(gene_to_mean_target, handle, protocol=pickle.HIGHEST_PROTOCOL)

def interpolate(x):
    """
    Intertpolate value x from target to projection domain
    """

    ## importation
    import pickle
    import numpy as np

    ## load dictionnary
    with open('ressources/gene_to_mean_ref.pickle', 'rb') as handle:
        gene_to_mean_ref = pickle.load(handle)
    with open('ressources/gene_to_mean_target.pickle', 'rb') as handle:
        gene_to_mean_target = pickle.load(handle)

    ## locate x interval
    a_gene = "NA"
    b_gene = "NA"
    cmpt = 0
    gene_list = list(gene_to_mean_target.keys())
    for gene in gene_list:
        if(cmpt+1 < len(gene_list)):
            a_value = gene_to_mean_target[gene]
            b_value = gene_to_mean_target[gene_list[cmpt+1]]
            if(x >= a_value and x <= b_value):
                a_gene = gene
                b_gene = gene_list[cmpt+1]
        cmpt +=1

    ## compute interpolation
    if(a_gene != "NA" and b_gene != "NA"):
        y = gene_to_mean_ref[a_gene] + (gene_to_mean_ref[b_gene] - gene_to_mean_ref[a_gene])*((x - gene_to_mean_target[a_gene])/(gene_to_mean_target[b_gene]-gene_to_mean_target[a_gene]))

    elif(x <= gene_to_mean_target[gene_list[0]]):
        y = gene_to_mean_ref[gene_list[0]]
    elif(x >= gene_to_mean_target[gene_list[-1]]):
        y =  gene_to_mean_ref[gene_list[-1]]
    elif(np.isnan(x)):
        y = np.nan

    ## return interpolated value
    return y

def run_interpolation(target_dataset, output_filename, verbose_mode):
    """
    Run interpolation on target dataset

        - target_dataset is a string, the name of the dataset to interpolate
        - output_filename is a string, the name of the interpolated dataset
        - verbose_mode is a boolean, usually set to True

    -> extract selected genes
    -> load target dataset
    -> apply interpolate function on each entry of the target dataset
    -> save interpolated dataset
    """

    ## importation
    import pandas as pd
    import numpy as np

    ## parameters
    feature_file = "ressources/features.csv"

    ## display information if needed
    if(verbose_mode):
        print("[+] Run interpolation")

    ## extract selected genes
    selected_genes = pd.read_csv(feature_file)
    selected_genes_list = list(selected_genes['FEATURES'])

    ## load dataset
    df = pd.read_csv(target_dataset)

    ## select variables
    df_selected = pd.DataFrame()
    for var in selected_genes_list:
        try:
            df_selected[var] = df[var]
        except:
            df_selected[var] = np.nan
    df_selected['ID'] = df['ID']
    df_selected = df_selected.set_index('ID')

    ## apply interpolate function on each entry of the dataset
    df_selected = df_selected.applymap(interpolate)

    ## save dataset
    df_selected.to_csv(output_filename)

    ## display information if needed
    if(verbose_mode):
        print("[+] Interpolation done")

max_stars_repo_name: LBAI-InfoLab/SjTree

max_issues_repo_name: LBAI-InfoLab/SjTree

max_forks_repo_name: LBAI-InfoLab/SjTree


text:
import numpy as np
import tempfile
import torch.onnx

from contextlib import contextmanager
from typing import Any, Dict, List, Optional

from dnnv.nn import operations
from dnnv.nn.visitors import EnsureSupportVisitor
from dnnv.verifiers.common.base import Parameter, Verifier
from dnnv.verifiers.common.results import SAT, UNSAT, UNKNOWN

from .errors import VerinetError, VerinetTranslatorError

class VeriNet(Verifier):
    EXE = "verinet.py"
    translator_error = VerinetTranslatorError
    verifier_error = VerinetError
    parameters = {
        "max_proc": Parameter(int, help="Maximum number of processes to use."),
        "no_split": Parameter(bool, help="Whether or not to do splitting."),
    }

    def build_inputs(self, prop):
        if prop.input_constraint.num_variables > 1:
            raise self.translator_error(
                "Unsupported network: More than 1 input variable"
            )

        with tempfile.NamedTemporaryFile(
            mode="w+", suffix=".onnx", delete=False
        ) as onnx_model_file:
            op_graph = prop.suffixed_op_graph()
            supported_operations = [
                operations.Concat,
                operations.Conv,
                operations.Flatten,
                operations.Gather,
                operations.Gemm,
                operations.Input,
                operations.Relu,
                operations.Reshape,
                operations.Shape,
                operations.Sigmoid,
                operations.Tanh,
                operations.Transpose,
                operations.Unsqueeze,
            ]
            op_graph.walk(
                EnsureSupportVisitor(supported_operations, self.translator_error)
            )
            op_graph.output_operations[0].b = np.hstack(
                [
                    op_graph.output_operations[0].b,
                    np.zeros(
                        op_graph.output_operations[0].b.shape,
                        dtype=op_graph.output_operations[0].b.dtype,
                    ),
                ]
            )

            torch.onnx.export(
                op_graph.as_pytorch(),
                tuple(
                    torch.ones(tuple((d if d > 0 else 1) for d in shape))
                    for shape in op_graph.input_shape
                ),
                onnx_model_file.name,
            )

        lb = prop.input_constraint.lower_bounds[0].flatten().copy()
        ub = prop.input_constraint.upper_bounds[0].flatten().copy()

        input_bounds = np.array(list(zip(lb, ub))).reshape(
            tuple(prop.input_constraint.lower_bounds[0].shape) + (2,)
        )
        with tempfile.NamedTemporaryFile(
            mode="w+", suffix=".npy", delete=False
        ) as input_bounds_file:
            np.save(input_bounds_file.name, input_bounds)

        with tempfile.NamedTemporaryFile(
            mode="w+", suffix=".npy", delete=False
        ) as output_file:
            self._tmp_output_file = output_file
        args = (
            "verinet.py",
            onnx_model_file.name,
            input_bounds_file.name,
            "-o",
            self._tmp_output_file.name,
        )
        if "max_proc" in self.parameters and self.parameters["max_proc"] is not None:
            value = self.parameters["max_proc"]
            args += (f"--max_procs={value}",)
        if "no_split" in self.parameters and self.parameters["no_split"]:
            args += ("--no_split",)
        return args

    def parse_results(self, prop, results):
        try:
            status, cex = np.load(self._tmp_output_file.name, allow_pickle=True)
            if status == "Safe":
                return UNSAT, None
            elif status == "Unsafe":
                input_shape, input_dtype = prop.op_graph.input_details[0]
                return SAT, cex.reshape(input_shape)
            elif status == "Unknown":
                return UNKNOWN, None
            elif status == "Undecided":
                raise self.verifier_error("Undecided")
            raise self.translator_error(f"Unknown verification result: {result_str}")
        finally:
            del self._tmp_output_file


text:
# Copyright (c) Microsoft Corporation. All rights reserved.
# Licensed under the MIT License.
#

import os
import io
import shutil
import argparse

import numpy as np

from utils import get_songs

parser = argparse.ArgumentParser(description="histogram evaluation for melody.")

parser.add_argument('--lyric-file', type=str,
                    help="The original lyric file")
parser.add_argument('--melody-file', type=str,
                    help="The original melody file")
parser.add_argument('--song-id-file', type=str,
                    help="The original song id file")
parser.add_argument('--generated-melody-file', type=str,
                    help="The generated melody file by")
parser.add_argument('--metric', choices=['pitch', 'duration'], default='pitch',
                    help="Calculate the pitch/duration distribution similarity")

def get_pitch_count(x):
    cnt = [0] * 129
    for i in x:
        cnt[int(i[0])] += 1
    return np.array(cnt)

def get_duration_count(x):
    cnt = [0 for _ in range(25, 3325, 25)]
    for i in x:
        cnt[int(i[1]) - 129] += 1
    return np.array(cnt)

def measure_pitch_simiarlity(targets, hypos):
    song_num = len(hypos)

    similarity = 0

    def get_pitch_histo(x):
        x = get_pitch_count(x)
        x_sum = np.sum(x) if np.sum(x) > 0 else 1
        x_dist = x.astype(np.float32) / x_sum
        return x_dist

    for i in range(song_num):
        hypo_histo = get_pitch_histo(hypos[i])
        target_histo = get_pitch_histo(targets[i])

        pitch_diff = np.abs(hypo_histo - target_histo)
        pitch_overlap_i = (hypo_histo + target_histo - pitch_diff) / 2

        similarity += np.sum(pitch_overlap_i)

    return similarity / song_num

def measure_duration_similarity(targets, hypos):
    song_num = len(hypos)

    similarity = 0

    def get_duration_histo(x):
        x = get_duration_count(x)
        x_sum = np.sum(x) if np.sum(x) > 0 else 1
        x_dist = x.astype(np.float32) / x_sum
        return x_dist

    for i in range(song_num):
        hypo_histo = get_duration_histo(hypos[i])
        target_histo = get_duration_histo(targets[i])

        duration_diff = np.abs(hypo_histo - target_histo)
        duration_overlap_i = (hypo_histo + target_histo - duration_diff) / 2

        similarity += np.sum(duration_overlap_i)

    return similarity / song_num

def main(args):
    hypos = get_songs(
        args.generated_melody_file,
        args.lyric_file,
        args.song_id_file,
        is_generated=True,
        get_last=True,
        find_structure=True,
        cut_exceed_sent=True,
    )

    targets = get_songs(
        args.melody_file,
        args.lyric_file,
        args.song_id_file,
    )

    metric_func = {
        'pitch': measure_pitch_simiarlity,
        'duration': measure_duration_similarity,
    }

    similarity = metric_func[args.metric](targets, hypos)
    print('The {} distribution similarity is {}'.format(args.metric, similarity))

if __name__ == '__main__':
    args = parser.parse_args()
    main(args)

max_stars_repo_name: hongwen-sun/muzic

max_issues_repo_name: hongwen-sun/muzic

max_forks_repo_name: hongwen-sun/muzic


text:
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

data = pd.read_csv("sensor_data_700.txt", delimiter=" ", header=None, names=("data","time","ir","lidar"))

# 12時から16時までのデータを抽出
d = data[ (data["time"] < 160000) & (data["time"] >= 120000) ]
d = d.loc[:, ["ir", "lidar"]]

sns.jointplot(d["ir"], d["lidar"], d, kind="kde")
plt.show()

# In[]:
print("光センサの計測値の分散:", d.ir.var())
print("LiDARの計測値の分散", d.lidar.var())

diff_ir = d.ir - d.ir.mean()
diff_lidar = d.lidar - d.lidar.mean()

a = diff_ir * diff_lidar
print("共分散:", sum(a)/(len(d-1)))

d.mean()

# In[]:
d.cov()

# In[]:
from scipy.stats import multivariate_normal

# 多次元ガウス分布のオブジェクト生成
irlidar = multivariate_normal(mean=d.mean().values.T, cov=d.cov().values)

# In[]:
import numpy as np

# 2次元平面に均等にx座標，y座標を作成
x, y = np.mgrid[0:40, 710:750]
# 40×40×2のリストを作成（xは40×40のリスト，さらにもう一次元を追加）
pos = np.empty(x.shape + (2,))
# 加えた3次元目にx,yを代入する
pos[:, :, 0] = x
pos[:, :, 1] = y
# x, y座標とそれに対応する密度を算出
cont = plt.contour(x, y, irlidar.pdf(pos))
# 等高線に値を書き込む
cont.clabel(fmt='%1.1e')

plt.show()

# In[]:
print("x座標:", x)
print("y座標:", y)

# In[]:
c = d.cov().values + np.array([[0, 20], [20,0]])
tmp = multivariate_normal(mean=d.mean().values.T, cov=c)
count = plt.contour(x, y, tmp.pdf(pos))
plt.show()

reason: import numpy,from scipy


text:
from sympy import *

r = symbols("r")

f = exp(-r/5)
#f = exp(-r*r/5)

s = integrate(f, (r, 0.0, oo))
print("s    = ", s)
print("N(s) = ", N(s))

reason: from sympy


text:
# Copyright (c) 2019 Alibaba Inc. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================

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

import tensorflow as tf
import numpy as np
import gym

from easy_rl.agents import agents
from easy_rl.models import EvolutionStrategy
from easy_rl.utils.window_stat import WindowStat

MODEL_CONFIG = dict(
    # specific
    type="ES",

    # common
    init_lr=0.01,
    lr_strategy_spec={
        'type': 'exponential_decay',
        'decay_steps': 50,
        'decay_rate': 0.9
    },
    global_norm_clip=40)

AGENT_CONFIG = dict(
    type="Agent",
    # how many trials needed for one optimization
    sample_batch_size=100,
)
np.random.seed(0)

class MyESmodel(EvolutionStrategy):
    def _encode_obs(self, input_obs, scope="encode_obs"):
        with tf.variable_scope(name_or_scope=scope):
            h1 = tf.layers.dense(
                input_obs,
                units=64,
                activation=tf.nn.relu,
                kernel_initializer=tf.random_normal_initializer(
                    mean=0.0, stddev=0.1, seed=0))
            h2 = tf.layers.dense(
                h1,
                units=64,
                activation=tf.nn.relu,
                kernel_initializer=tf.random_normal_initializer(
                    mean=0.0, stddev=0.1, seed=0))
            logits = tf.layers.dense(
                h2,
                units=2,
                activation=None,
                kernel_initializer=tf.random_normal_initializer(
                    mean=0.0, stddev=0.1, seed=0))
            return logits

def main():
    env = gym.make("CartPole-v0")
    env.seed(0)

    agent_class = agents[AGENT_CONFIG["type"]]
    agent = agent_class(
        env.observation_space,
        env.action_space,
        AGENT_CONFIG,
        MODEL_CONFIG,
        distributed_spec={},
        checkpoint_dir="ckpt_dir",
        custom_model=MyESmodel)

    reward_window = WindowStat("reward", 50)
    length_window = WindowStat("length", 50)

    init_perturbation_scale = 0.1

    seeds, rewards, perturbation_scales = list(), list(), list()
    is_positive_direction = list()

    # how many episodes needed for one trial
    episode_per_perturbation = 1

    returns = list()

    for i in range(4000):
        ob = env.reset()
        done = False
        episode_reward = .0
        episode_len = 0

        if i % episode_per_perturbation == 0:
            # perturb parameters every `episode_per_seed` episodes
            is_positive = True if len(
                is_positive_direction
            ) == 0 else is_positive_direction[-1] != True

            # each seed twice
            seed = np.random.randint(1000000) if is_positive else seeds[-1]
            perturbation_scale = max(
                init_perturbation_scale * (1 - i / 2000.0), 0.02)

            feed = agent.model.perturbation_feed
            fetch = [agent.model.reset_perturbation_op]

            agent.executor.run(
                fetches=fetch,
                feed_dict={
                    feed['perturbation_seeds']: [seed],
                    feed['perturbation_scales']: [perturbation_scale],
                    feed['positive_perturbation']: is_positive
                })

            if is_positive:
                seeds.append(seed)
                perturbation_scales.append(perturbation_scale)
            is_positive_direction.append(is_positive)

        while not done:
            action, result = agent.act([ob], True, use_perturbed_action=True)

            next_ob, reward, done, info = env.step(action[0])

            ob = next_ob
            episode_reward += reward
            episode_len += 1

        rewards.append(episode_reward)
        reward_window.push(episode_reward)
        length_window.push(episode_len)
        if len(rewards) == episode_per_perturbation:
            returns.append(np.mean(rewards))
            rewards = []
            if len(returns) == 2 * agent.config.get('sample_batch_size', 100):
                print(reward_window)
                assert len(seeds) == (len(returns) / 2)
                assert len(perturbation_scales) == (len(returns) / 2)
                agent.learn(
                    batch_data=dict(
                        perturbation_seeds=seeds,
                        perturbation_scales=perturbation_scales,
                        returns=np.reshape(returns, [-1, 2])))
                seeds = []
                perturbation_scales = []
                returns = []
                is_positive_direction = []

                # evaluation 20 episodes
                test_rewards = list()
                for j in range(20):
                    done = False
                    ob = env.reset()
                    episode_reward = 0
                    episode_len = 0
                    while not done:
                        action, result = agent.act(
                            [ob], True, use_perturbed_action=False)

                        next_ob, reward, done, info = env.step(action[0])

                        ob = next_ob
                        episode_reward += reward
                        episode_len += 1
                    test_rewards.append(episode_reward)
                print("[evaluation] average reward of 20 episodes:",
                      np.mean(test_rewards))

                print('train at ', i)

    agent.export_saved_model(export_dir="dump_dir")
    print("Done.")

if __name__ == "__main__":
    main()

max_stars_repo_name: simonoso/EasyRL

max_issues_repo_name: simonoso/EasyRL

max_forks_repo_name: simonoso/EasyRL


text:
import numpy as np

# Copyright (c) 2018-present, Facebook, Inc.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.
#

import numpy as np
#import h5py

mpii_metadata = {
    'layout_name': 'mpii',
    'num_joints': 16,
    'keypoints_symmetry': [
        [3, 4, 5, 13, 14, 15],
        [0, 1, 2, 10, 11, 12],
    ]
}

coco_metadata = {
    'layout_name': 'coco',
    'num_joints': 17,
    'keypoints_symmetry': [
        [1, 3, 5, 7, 9, 11, 13, 15],
        [2, 4, 6, 8, 10, 12, 14, 16],
    ]
}

h36m_metadata = {
    'layout_name': 'h36m',
    'num_joints': 17,
    'keypoints_symmetry': [
        [4, 5, 6, 11, 12, 13],
        [1, 2, 3, 14, 15, 16],
    ]
}

humaneva15_metadata = {
    'layout_name': 'humaneva15',
    'num_joints': 15,
    'keypoints_symmetry': [
        [2, 3, 4, 8, 9, 10],
        [5, 6, 7, 11, 12, 13]
    ]
}

humaneva20_metadata = {
    'layout_name': 'humaneva20',
    'num_joints': 20,
    'keypoints_symmetry': [
        [3, 4, 5, 6, 11, 12, 13, 14],
        [7, 8, 9, 10, 15, 16, 17, 18]
    ]
}

def suggest_metadata(name):
    names = []
    for metadata in [mpii_metadata, coco_metadata, h36m_metadata, humaneva15_metadata, humaneva20_metadata]:
        if metadata['layout_name'] in name:
            return metadata
        names.append(metadata['layout_name'])
    raise KeyError('Cannot infer keypoint layout from name "{}". Tried {}.'.format(name, names))

def import_detectron_poses(path):
    # Latin1 encoding because Detectron runs on Python 2.7
    data = np.load(path, encoding='latin1')
    kp = data['keypoints']
    bb = data['boxes']
    results = []
    for i in range(len(bb)):
        if len(bb[i][1]) == 0:
            assert i > 0
            # Use last pose in case of detection failure
            results.append(results[-1])
            continue
        best_match = np.argmax(bb[i][1][:, 4])
        keypoints = kp[i][1][best_match].T.copy()
        results.append(keypoints)
    results = np.array(results)
    return results[:, :, 4:6]  # Soft-argmax
    # return results[:, :, [0, 1, 3]] # Argmax + score

def import_cpn_poses(path):
    data = np.load(path)
    kp = data['keypoints']
    return kp[:, :, :2]

def import_sh_poses(path):
    with h5py.File(path) as hf:
        positions = hf['poses'].value
    return positions.astype('float32')

def suggest_pose_importer(name):
    if 'detectron' in name:
        return import_detectron_poses
    if 'cpn' in name:
        return import_cpn_poses
    if 'sh' in name:
        return import_sh_poses
    raise KeyError('Cannot infer keypoint format from name "{}". Tried detectron, cpn, sh.'.format(name))

def fetch(subjects, keypoints, dataset, downsample, action_filter=None, cam_filter=None, subset=1, parse_3d_poses=True):
    out_poses_3d = []
    out_poses_2d = []
    out_camera_params = []
    for subject in subjects:
        for action in keypoints[subject].keys():
            action_split = action.split(' ')[0]
            if action_filter is not None:
                found = False
                # distinguish the actions:'Sitting' and 'SittingDown'
                for act in action_filter:
                    act = act.split(' ')[0]
                    if action_split == act:
                        found = True
                        break
                if not found:
                    continue
            poses_2d = keypoints[subject][action]
            index = np.random.randint(0,4)
            if cam_filter==[5]: #random camera index
                out_poses_2d.append(poses_2d[index])
                print('choose a camera index for each action:', index)

            elif cam_filter:
                for j in cam_filter: # Select by some camera viewpoints
                    out_poses_2d.append(poses_2d[j])

            else:
                for i in range(len(poses_2d)):  # Iterate across cameras
                    out_poses_2d.append(poses_2d[i])

            if subject in dataset.cameras():
                cams = dataset.cameras()[subject]
                assert len(cams) == len(poses_2d), 'Camera count mismatch'
                if cam_filter==[5]:
                    cam = cams[index]
                    if 'intrinsic' in cam:
                        use_params = {}
                        use_params['intrinsic'] = cam['intrinsic']
                        if 'normalization_params' in dataset[subject][action]:
                            use_params['normalization_params'] = \
                                dataset[subject][action]['normalization_params'][index]
                        out_camera_params.append(use_params)
                elif cam_filter:
                    for j in cam_filter:
                        for i, cam in enumerate(cams):
                            if j == i:
                                if 'intrinsic' in cam:
                                    use_params = {}
                                    use_params['intrinsic'] = cam['intrinsic']
                                    if 'normalization_params' in dataset[subject][action]:
                                        use_params['normalization_params'] = \
                                        dataset[subject][action]['normalization_params'][i]
                                    out_camera_params.append(use_params)
                else:
                    for i, cam in enumerate(cams):
                        if 'intrinsic' in cam:
                            use_params = {}
                            use_params['intrinsic'] = cam['intrinsic']
                            if 'normalization_params' in dataset[subject][action]:
                                use_params['normalization_params'] = dataset[subject][action]['normalization_params'][i]
                            out_camera_params.append(use_params)

            if parse_3d_poses and 'positions_3d' in dataset[subject][action]:
                poses_3d = dataset[subject][action]['positions_3d']
                assert len(poses_3d) == len(poses_2d), 'Camera count mismatch'
                if cam_filter==[5]:
                    out_poses_3d.append(poses_3d[index])
                elif cam_filter :
                    for j in cam_filter:
                        out_poses_3d.append(poses_3d[j])
                else:
                    for i in range(len(poses_3d)):  # Iterate across cameras
                        out_poses_3d.append(poses_3d[i])

    if len(out_camera_params) == 0:
        out_camera_params = None
    if len(out_poses_3d) == 0:
        out_poses_3d = None

    stride = downsample
    if subset < 1:
        for i in range(len(out_poses_2d)):
            n_frames = int(round(len(out_poses_2d[i]) // stride * subset) * stride)
            start = deterministic_random(0, len(out_poses_2d[i]) - n_frames + 1, str(len(out_poses_2d[i])))
            out_poses_2d[i] = out_poses_2d[i][start:start + n_frames:stride]
            if out_poses_3d is not None:
                out_poses_3d[i] = out_poses_3d[i][start:start + n_frames:stride]
    elif stride > 1:
        # Downsample as requested
        for i in range(len(out_poses_2d)):
            out_poses_2d[i] = out_poses_2d[i][::stride]
            if out_poses_3d is not None:
                out_poses_3d[i] = out_poses_3d[i][::stride]

    return out_camera_params, out_poses_3d, out_poses_2d

max_stars_repo_name: ailingzengzzz/Split-and-Recombine-Net

max_issues_repo_name: ailingzengzzz/Split-and-Recombine-Net

max_forks_repo_name: ailingzengzzz/Split-and-Recombine-Net