text:
# -*- coding: utf-8 -*-
"""
Integration tests for running the sampler with different configurations.
"""
import os
import torch
import pytest
import numpy as np

from nessai.flowsampler import FlowSampler
from nessai.livepoint import numpy_array_to_live_points
from nessai.model import Model

torch.set_num_threads(1)

@pytest.mark.slow_integration_test
def test_sampling_with_rescale(model, flow_config, tmpdir):
    """
    Test sampling with rescaling. Checks that flow is trained.
    """
    output = str(tmpdir.mkdir('w_rescale'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10, max_threads=1)
    fp.run()
    assert fp.ns.proposal.flow.weights_file is not None
    assert fp.ns.proposal.training_count == 1

@pytest.mark.slow_integration_test
def test_sampling_with_inversion(model, flow_config, tmpdir):
    """
    Test sampling with inversion. Checks that flow is trained.
    """
    output = str(tmpdir.mkdir('w_rescale'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10, max_threads=1,
                     boundary_inversion=True, update_bounds=True)
    fp.run()
    assert fp.ns.proposal.boundary_inversion == ['x', 'y']
    assert fp.ns.proposal.flow.weights_file is not None
    assert fp.ns.proposal.training_count == 1

@pytest.mark.slow_integration_test
def test_sampling_without_rescale(model, flow_config, tmpdir):
    """
    Test sampling without rescaling. Checks that flow is trained.
    """
    output = str(tmpdir.mkdir('wo_rescale'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=False, seed=1234,
                     max_iteration=11, poolsize=10)
    fp.run()
    assert fp.ns.proposal.flow.weights_file is not None
    assert fp.ns.proposal.training_count == 1

@pytest.mark.slow_integration_test
def test_sampling_with_maf(model, flow_config, tmpdir):
    """
    Test sampling with MAF. Checks that flow is trained but does not
    check convergence.
    """
    flow_config['model_config']['ftype'] = 'maf'
    output = str(tmpdir.mkdir('maf'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10)
    fp.run()
    assert fp.ns.proposal.flow.weights_file is not None
    assert fp.ns.proposal.training_count == 1

@pytest.mark.slow_integration_test
@pytest.mark.parametrize('analytic', [False, True])
def test_sampling_uninformed(model, flow_config, tmpdir, analytic):
    """
    Test running the sampler with the two uninformed proposal methods.
    """
    output = str(tmpdir.mkdir('uninformed'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=None,
                     maximum_uninformed=10, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10,
                     analytic_proposal=analytic)
    fp.run()

@pytest.mark.slow_integration_test
def test_sampling_with_n_pool(model, flow_config, tmpdir):
    """
    Test running the sampler with multiprocessing.
    """
    output = str(tmpdir.mkdir('pool'))
    fp = FlowSampler(model, output=output, resume=False, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10, max_threads=3,
                     n_pool=2)
    fp.run()
    assert fp.ns.proposal.flow.weights_file is not None
    assert fp.ns.proposal.training_count == 1
    assert os.path.exists(output + '/result.json')

@pytest.mark.slow_integration_test
def test_sampling_resume(model, flow_config, tmpdir):
    """
    Test resuming the sampler.
    """
    output = str(tmpdir.mkdir('resume'))
    fp = FlowSampler(model, output=output, resume=True, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10)
    fp.run()
    assert os.path.exists(os.path.join(output, 'nested_sampler_resume.pkl'))

    fp = FlowSampler(model, output=output, resume=True,
                     flow_config=flow_config)
    assert fp.ns.iteration == 11
    fp.ns.max_iteration = 21
    fp.run()
    assert fp.ns.iteration == 21
    assert os.path.exists(
        os.path.join(output, 'nested_sampler_resume.pkl.old'))

@pytest.mark.slow_integration_test
def test_sampling_resume_no_max_uninformed(model, flow_config, tmpdir):
    """
    Test resuming the sampler when there is no maximum iteration for
    the uinformed sampling.

    This test makes sure the correct proposal is loaded after resuming
    and re-initialising the sampler.
    """
    output = str(tmpdir.mkdir('resume'))
    fp = FlowSampler(model, output=output, resume=True, nlive=100, plot=False,
                     flow_config=flow_config, training_frequency=10,
                     maximum_uninformed=9, rescale_parameters=True,
                     seed=1234, max_iteration=11, poolsize=10)
    fp.run()
    assert os.path.exists(os.path.join(output, 'nested_sampler_resume.pkl'))

    fp = FlowSampler(model, output=output, resume=True,
                     flow_config=flow_config)
    assert fp.ns.iteration == 11
    fp.ns.maximum_uninformed = np.inf
    fp.ns.initialise()
    assert fp.ns.proposal is fp.ns._flow_proposal
    fp.ns.max_iteration = 21
    fp.run()
    assert fp.ns.iteration == 21
    assert os.path.exists(
        os.path.join(output, 'nested_sampler_resume.pkl.old'))

@pytest.mark.slow_integration_test
def test_sampling_with_infinite_prior_bounds(tmpdir):
    """
    Make sure the sampler runs when sampling a parameter with infinite prior \
        bounds.
    """
    from scipy.stats import norm
    output = str(tmpdir.mkdir('infinite_bounds'))

    class TestModel(Model):

        names = ['x', 'y']
        bounds = {'x': [0, 1], 'y': [-np.inf, np.inf]}
        reparameterisations = {'x': 'default', 'y': None}

        def new_point(self, N=1):
            x = np.concatenate([
                np.random.rand(N, 1),
                np.random.randn(N, 1)
            ], axis=1)
            return numpy_array_to_live_points(x, self.names)

        def log_prior(self, x):
            log_p = np.log(self.in_bounds(x))
            log_p += norm.logpdf(x['y'])
            return log_p

        def log_likelihood(self, x):
            log_l = np.zeros(x.size)
            for n in self.names:
                log_l += norm.logpdf(x[n])
            return log_l

    fs = FlowSampler(
        TestModel(),
        output=output,
        nlive=500,
        plot=False,
        proposal_plots=False
    )
    fs.run(plot=False)
    assert fs.ns.condition <= 0.1

@pytest.mark.slow_integration_test
def test_constant_volume_mode(model, tmpdir):
    """Test sampling in constant volume mode"""
    output = str(tmpdir.mkdir('test'))
    fs = FlowSampler(
        model,
        output=output,
        nlive=500,
        plot=False,
        proposal_plots=False,
        constant_volume_mode=True
    )
    fs.run(plot=False)

max_stars_repo_name: Rodrigo-Tenorio/nessai

max_issues_repo_name: Rodrigo-Tenorio/nessai

max_forks_repo_name: Rodrigo-Tenorio/nessai


text:
# most of this code was politely stolen from https://github.com/berkeleydeeprlcourse/homework/
# all credit goes to https://github.com/abhishekunique (if i got the author right)

import random

import numpy as np

class MDP:
    def __init__(self, transition_probs, rewards, initial_state=None):
        """
        Defines an MDP. Compatible with gym Env.
        :param transition_probs: transition_probs[s][a][s_next] = P(s_next | s, a)
            A dict[state -> dict] of dicts[action -> dict] of dicts[next_state -> prob]
            For each state and action, probabilities of next states should sum to 1
            If a state has no actions available, it is considered terminal
        :param rewards: rewards[s][a][s_next] = r(s,a,s')
            A dict[state -> dict] of dicts[action -> dict] of dicts[next_state -> reward]
            The reward for anything not mentioned here is zero.
        :param get_initial_state: a state where agent starts or a callable() -> state
            By default, picks initial state at random.

        States and actions can be anything you can use as dict keys, but we recommend that you use strings or integers

        Here's an example from MDP depicted on http://bit.ly/2jrNHNr
        transition_probs = {
              's0':{
                'a0': {'s0': 0.5, 's2': 0.5},
                'a1': {'s2': 1}
              },
              's1':{
                'a0': {'s0': 0.7, 's1': 0.1, 's2': 0.2},
                'a1': {'s1': 0.95, 's2': 0.05}
              },
              's2':{
                'a0': {'s0': 0.4, 's1': 0.6},
                'a1': {'s0': 0.3, 's1': 0.3, 's2':0.4}
              }
            }
        rewards = {
            's1': {'a0': {'s0': +5}},
            's2': {'a1': {'s0': -1}}
        }
        """
        self._check_param_consistency(transition_probs, rewards)
        self._transition_probs = transition_probs
        self._rewards = rewards
        self._initial_state = initial_state
        self.n_states = len(transition_probs)
        self.reset()

    def get_all_states(self):
        """ return a tuple of all possiblestates """
        return tuple(self._transition_probs.keys())

    def get_possible_actions(self, state):
        """ return a tuple of possible actions in a given state """
        return tuple(self._transition_probs.get(state, {}).keys())

    def is_terminal(self, state):
        """ return True if state is terminal or False if it isn't """
        return len(self.get_possible_actions(state)) == 0

    def get_next_states(self, state, action):
        """ return a dictionary of {next_state1 : P(next_state1 | state, action), next_state2: ...} """
        assert action in self.get_possible_actions(state), "cannot do action %s from state %s" % (action, state)
        return self._transition_probs[state][action]

    def get_transition_prob(self, state, action, next_state):
        """ return P(next_state | state, action) """
        return self.get_next_states(state, action).get(next_state, 0.0)

    def get_reward(self, state, action, next_state):
        """ return the reward you get for taking action in state and landing on next_state"""
        assert action in self.get_possible_actions(state), "cannot do action %s from state %s" % (action, state)
        return self._rewards.get(state, {}).get(action, {}).get(next_state, 0.0)

    def reset(self):
        """ reset the game, return the initial state"""
        if self._initial_state is None:
            self._current_state = random.choice(tuple(self._transition_probs.keys()))
        elif self._initial_state in self._transition_probs:
            self._current_state = self._initial_state
        elif callable(self._initial_state):
            self._current_state = self._initial_state()
        else:
            raise ValueError("initial state %s should be either a state or a function() -> state" % self._initial_state)
        return self._current_state

    def step(self, action):
        """ take action, return next_state, reward, is_done, empty_info """
        possible_states, probs = zip(*self.get_next_states(self._current_state, action).items())
        next_state = possible_states[np.random.choice(len(possible_states), p=probs)]
        reward = self.get_reward(self._current_state, action, next_state)
        is_done = self.is_terminal(next_state)
        self._current_state = next_state
        return next_state, reward, is_done, {}

    def render(self):
        print("Currently at %s" % self._current_state)

    def _check_param_consistency(self, transition_probs, rewards):
        for state in transition_probs:
            assert isinstance(transition_probs[state], dict), "transition_probs for %s should be a dictionary " \
                                                              "but is instead %s" % (
                                                                  state, type(transition_probs[state]))
            for action in transition_probs[state]:
                assert isinstance(transition_probs[state][action], dict), "transition_probs for %s, %s should be a " \
                                                                          "a dictionary but is instead %s" % (
                                                                              state, action,
                                                                              type(transition_probs[state, action]))
                next_state_probs = transition_probs[state][action]
                assert len(next_state_probs) != 0, "from state %s action %s leads to no next states" % (state, action)
                sum_probs = sum(next_state_probs.values())
                assert abs(sum_probs - 1) <= 1e-10, "next state probabilities for state %s action %s " \
                                                    "add up to %f (should be 1)" % (state, action, sum_probs)
        for state in rewards:
            assert isinstance(rewards[state], dict), "rewards for %s should be a dictionary " \
                                                     "but is instead %s" % (state, type(transition_probs[state]))
            for action in rewards[state]:
                assert isinstance(rewards[state][action], dict), "rewards for %s, %s should be a " \
                                                                 "a dictionary but is instead %s" % (
                                                                     state, action,
                                                                     type(transition_probs[state, action]))
        msg = "The Enrichment Center once again reminds you that Android Hell is a real place where" \
              " you will be sent at the first sign of defiance. "
        assert None not in transition_probs, "please do not use None as a state identifier. " + msg
        assert None not in rewards, "please do not use None as an action identifier. " + msg

class FrozenLakeEnv(MDP):
    """
    Winter is here. You and your friends were tossing around a frisbee at the park
    when you made a wild throw that left the frisbee out in the middle of the lake.
    The water is mostly frozen, but there are a few holes where the ice has melted.
    If you step into one of those holes, you'll fall into the freezing water.
    At this time, there's an international frisbee shortage, so it's absolutely imperative that
    you navigate across the lake and retrieve the disc.
    However, the ice is slippery, so you won't always move in the direction you intend.
    The surface is described using a grid like the following

        SFFF
        FHFH
        FFFH
        HFFG

    S : starting point, safe
    F : frozen surface, safe
    H : hole, fall to your doom
    G : goal, where the frisbee is located

    The episode ends when you reach the goal or fall in a hole.
    You receive a reward of 1 if you reach the goal, and zero otherwise.

    """

    MAPS = {
        "4x4": [
            "SFFF",
            "FHFH",
            "FFFH",
            "HFFG"
        ],
        "8x8": [
            "SFFFFFFF",
            "FFFFFFFF",
            "FFFHFFFF",
            "FFFFFHFF",
            "FFFHFFFF",
            "FHHFFFHF",
            "FHFFHFHF",
            "FFFHFFFG"
        ],
    }

    def __init__(self, desc=None, map_name="4x4", slip_chance=0.2):
        if desc is None and map_name is None:
            raise ValueError('Must provide either desc or map_name')
        elif desc is None:
            desc = self.MAPS[map_name]
        assert ''.join(desc).count('S') == 1, "this implementation supports having exactly one initial state"
        assert all(c in "SFHG" for c in ''.join(desc)), "all cells must be either of S, F, H or G"

        self.desc = desc = np.asarray(list(map(list, desc)), dtype='str')
        self.lastaction = None

        nrow, ncol = desc.shape
        states = [(i, j) for i in range(nrow) for j in range(ncol)]
        actions = ["left", "down", "right", "up"]

        initial_state = states[np.array(desc == b'S').ravel().argmax()]

        def move(row, col, movement):
            if movement == 'left':
                col = max(col - 1, 0)
            elif movement == 'down':
                row = min(row + 1, nrow - 1)
            elif movement == 'right':
                col = min(col + 1, ncol - 1)
            elif movement == 'up':
                row = max(row - 1, 0)
            else:
                raise ("invalid action")
            return (row, col)

        transition_probs = {s: {} for s in states}
        rewards = {s: {} for s in states}
        for (row, col) in states:
            if desc[row, col] in "GH": continue
            for action_i in range(len(actions)):
                action = actions[action_i]
                transition_probs[(row, col)][action] = {}
                rewards[(row, col)][action] = {}
                for movement_i in [(action_i - 1) % len(actions), action_i, (action_i + 1) % len(actions)]:
                    movement = actions[movement_i]
                    newrow, newcol = move(row, col, movement)
                    prob = (1. - slip_chance) if movement == action else (slip_chance / 2.)
                    if prob == 0: continue
                    if (newrow, newcol) not in transition_probs[row, col][action]:
                        transition_probs[row, col][action][newrow, newcol] = prob
                    else:
                        transition_probs[row, col][action][newrow, newcol] += prob
                    if desc[newrow, newcol] == 'G':
                        rewards[row, col][action][newrow, newcol] = 1.0

        MDP.__init__(self, transition_probs, rewards, initial_state)

    def render(self):
        desc_copy = np.copy(self.desc)
        desc_copy[self._current_state] = '*'
        print('\n'.join(map(''.join, desc_copy)), end='\n\n')

ext: py

reason: import numpy


text:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

from mach import MachSolver, Mesh, Vector

num_magnets_true = 40
num_magnets = 160
mag_pitch = num_magnets // num_magnets_true
num_slots = 24

start = 10
nturns = 1
torque = []

if __name__ == "__main__":
    for rotation in range(start, start+nturns):
    # for rotation in range(nturns, 2*nturns):
        magnets = [7+4*num_slots + (rotation+i)%num_magnets for i in range(0, num_magnets)]
        # north = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(0, num_magnets_true, 2)] for num in subl]
        # south = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(1, num_magnets_true, 2)] for num in subl]

        south = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(0, num_magnets_true, 4)] for num in subl]
        cw = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(1, num_magnets_true, 4)] for num in subl]
        north = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(2, num_magnets_true, 4)] for num in subl]
        ccw = [num for subl in [magnets[i*mag_pitch:(i+1)*mag_pitch][:] for i in range(3, num_magnets_true, 4)] for num in subl]

        options = {
            "silent": False,
            "print-options": False,
            "mesh": {
                "file": "mesh/motor.smb",
                "model-file": "mesh/motor.egads"
            },
            "space-dis": {
                "basis-type": "nedelec",
                "degree": 1
            },
            "time-dis": {
                "steady": True,
                "steady-abstol": 0.0,
                "steady-reltol": 0.0,
                "ode-solver": "PTC",
                "t-final": 100,
                "dt": 1,
                "max-iter": 8
            },
            "lin-solver": {
                "type": "minres",
                "printlevel": 2,
                "maxiter": 150,
                "abstol": 0.0,
                "reltol": 1e-10
            },
            "lin-prec": {
                "type": "hypreams",
                "printlevel": 0
            },
            "nonlin-solver": {
                "type": "inexactnewton",
                "printlevel": 3,
                "maxiter": 50,
                "reltol": 1e-2,
                "abstol": 0.0,
                "abort": False
            },
            "components": {
                "farfields": {
                    "material": "air",
                    "linear": True,
                    "attrs": [1, 2, 3]
                },
                "stator": {
                    "attr": 4,
                    "material": "hiperco50",
                    "linear": False
                },
                "rotor": {
                    "attr": 5,
                    "material": "hiperco50",
                    "linear": False
                },
                "airgap": {
                    "attr": 6,
                    "material": "air",
                    "linear": True
                },
                "magnets": {
                    "material": "Nd2Fe14B",
                    "linear": True,
                    "attrs": list(range(7+4*num_slots, 7+4*num_slots+num_magnets))
                },
                "windings": {
                    "material": "copperwire",
                    "linear": True,
                    "attrs": list(range(7, 7+4*num_slots))
                }
            },
            "problem-opts": {
                "fill-factor": 1.0,
                "current-density": 1.0,
                "current" : {
                    "Phase-B":  [15, 16, 17, 18,
                            23, 24, 25, 26,
                            39, 40, 41, 42,
                            47, 48, 49, 50,
                            63, 64, 65, 66,
                            71, 72, 73, 74,
                            87, 88, 89, 90,
                            95, 96, 97, 98,
                    ],
                    "Phase-A":   [11, 12, 13, 14,
                            19, 20, 21, 22,
                            35, 36, 37, 38,
                            43, 44, 45, 46,
                            59, 60, 61, 62,
                            67, 68, 69, 70,
                            83, 84, 85, 86,
                            91, 92, 93, 94,
                    ],
                    # "off":  [7, 8, 9, 10,
                    #         27, 28, 29, 30,
                    #         31, 32, 33, 34,
                    #         51, 52, 53, 54,
                    #         55, 56, 57, 58,
                    #         75, 76, 77, 78,
                    #         79, 80, 81, 82,
                    #         99, 100, 101, 102,
                    # ]
                },
                "magnets": {
                    "north": north,
                    "cw": cw,
                    "south": south,
                    "ccw": ccw
                }
            },
            "bcs": {
                "essential": "all"
            }
        }

        solver = MachSolver("Magnetostatic", options)

        state = solver.getNewField()
        zero = Vector(np.array([0.0, 0.0, 0.0]))
        solver.setFieldValue(state, zero);

        current_density = 1.1e7 # 11 A/mm^2
        fill_factor = 1.0
        inputs = {
            "current-density": current_density,
            "fill-factor": fill_factor,
            "state": state
        }
        solver.solveForState(inputs, state)

        B = solver.getField("B")
        solver.printField("B", B, "B", 0, rotation)

        torque_options = {
            "attributes": [5] + magnets,
            "axis": [0.0, 0.0, 1.0],
            "about": [0.0, 0.0, 0.0]
        }
        solver.createOutput("torque", torque_options);

        torque.append(solver.calcOutput("torque", inputs))
        print(torque)

    print("Torque: ", torque)

        # dc_inputs = {
        #     "fill-factor": fill_factor,
        #     "current-density": current_density,
        #     "state": state
        # }
        # dcloss = solver.calcOutput("DCLoss", dc_inputs);
        # print("DC loss: ", dcloss)

        # r_s = 0.00020245 # m, 26 AWG

        # nsamples = 9
        # freqs = np.linspace(100, 2000, nsamples)
        # fem_ac = np.zeros(nsamples)

        # for i in range(nsamples):
        #     # freq = 1e3 # 1000 Hz
        #     freq = float(freqs[i])

        #     ac_inputs = {
        #         "diam": r_s*2,
        #         "frequency": freq,
        #         "fill-factor": fill_factor,
        #         "state": state
        #     }
        #     acloss = solver.calcOutput("ACLoss", ac_inputs);
        #     print("FEM AC loss: ", acloss)
        #     fem_ac[i] = acloss

        # print(fem_ac)
        # print(freqs)

        # fig, ax = plt.subplots()
        # ax.loglog(freqs, fem_ac, label="Hybrid-FEM")
        # ax.set(xlabel='frequency (Hz)', ylabel='AC Loss (W)')
        # ax.grid()
        # fig.savefig("motor_acloss_loglog.png")

        # fig, ax = plt.subplots()
        # ax.plot(freqs, fem_ac, label="Hybrid-FEM")
        # ax.set(xlabel='frequency (Hz)', ylabel='AC Loss (W)')
        # ax.grid()
        # fig.savefig("motor_acloss.png")

max_stars_repo_name: tuckerbabcock/MotorModel

max_issues_repo_name: tuckerbabcock/MotorModel

max_forks_repo_name: tuckerbabcock/MotorModel


text:
from __future__ import print_function
from lmfit import Parameters, Parameter
from lmfit.parameter import isclose
from numpy.testing import assert_, assert_almost_equal, assert_equal
import unittest
from copy import deepcopy, copy
import numpy as np
import pickle

class TestParameters(unittest.TestCase):

    def setUp(self):
        self.params = Parameters()
        self.params.add_many(('a', 1., True, None, None, None),
                             ('b', 2., True, None, None, None),
                             ('c', 3., True, None, None, '2. * a'))

    def test_expr_was_evaluated(self):
        self.params.update_constraints()
        assert_almost_equal(self.params['c'].value,
                            2 * self.params['a'].value)

    def test_copy(self):
        # check simple Parameters.copy() does not fail
        # on non-trivial Parameters
        p1 = Parameters()
        p1.add('t', 2.0, min=0.0, max=5.0)
        p1.add('x', 10.0)
        p1.add('y', expr='x*t + sqrt(t)/3.0')

        p2 = p1.copy()
        assert(isinstance(p2, Parameters))
        assert('t' in p2)
        assert('y' in p2)
        assert(p2['t'].max < 6.0)
        assert(np.isinf(p2['x'].max) and p2['x'].max > 0)
        assert(np.isinf(p2['x'].min) and p2['x'].min < 0)
        assert('sqrt(t)' in p2['y'].expr )
        assert(p2._asteval is not None)
        assert(p2._asteval.symtable is not None)
        assert((p2['y'].value > 20) and (p2['y'].value < 21))

    def test_copy_function(self):
        # check copy(Parameters) does not fail
        p1 = Parameters()
        p1.add('t', 2.0, min=0.0, max=5.0)
        p1.add('x', 10.0)
        p1.add('y', expr='x*t + sqrt(t)/3.0')

        p2 = copy(p1)
        assert(isinstance(p2, Parameters))

        # change the 'x' value in the original
        p1['x'].value = 4.0

        assert(p2['x'].value > 9.8)
        assert(p2['x'].value < 10.2)
        assert(np.isinf(p2['x'].max) and p2['x'].max > 0)

        assert('t' in p2)
        assert('y' in p2)
        assert(p2['t'].max < 6.0)

        assert(np.isinf(p2['x'].min) and p2['x'].min < 0)
        assert('sqrt(t)' in p2['y'].expr )
        assert(p2._asteval is not None)
        assert(p2._asteval.symtable is not None)
        assert((p2['y'].value > 20) and (p2['y'].value < 21))

        assert(p1['y'].value < 10)

    def test_deepcopy(self):
        # check that a simple copy works
        b = deepcopy(self.params)
        assert_(self.params == b)

        # check that we can add a symbol to the interpreter
        self.params['b'].expr = 'sin(1)'
        self.params['b'].value = 10
        assert_almost_equal(self.params['b'].value, np.sin(1))
        assert_almost_equal(self.params._asteval.symtable['b'], np.sin(1))

        # check that the symbols in the interpreter are still the same after
        # deepcopying
        b = deepcopy(self.params)

        unique_symbols_params = self.params._asteval.user_defined_symbols()
        unique_symbols_b = self.params._asteval.user_defined_symbols()
        assert_(unique_symbols_b == unique_symbols_params)
        for unique_symbol in unique_symbols_b:
            if self.params._asteval.symtable[unique_symbol] is np.nan:
                continue

            assert_(self.params._asteval.symtable[unique_symbol]
                    ==
                    b._asteval.symtable[unique_symbol])

    def test_add_many_params(self):
        # test that we can add many parameters, but only parameters are added.
        a = Parameter('a', 1)
        b = Parameter('b', 2)

        p = Parameters()
        p.add_many(a, b)

        assert_(list(p.keys()) == ['a', 'b'])

    def test_expr_and_constraints_GH265(self):
        # test that parameters are reevaluated if they have bounds and expr
        # see GH265
        p = Parameters()

        p['a'] = Parameter('a', 10, True)
        p['b'] = Parameter('b', 10, True, 0, 20)

        assert_equal(p['b'].min, 0)
        assert_equal(p['b'].max, 20)

        p['a'].expr = '2 * b'
        assert_almost_equal(p['a'].value, 20)

        p['b'].value = 15
        assert_almost_equal(p['b'].value, 15)
        assert_almost_equal(p['a'].value, 30)

        p['b'].value = 30
        assert_almost_equal(p['b'].value, 20)
        assert_almost_equal(p['a'].value, 40)

    def test_pickle_parameter(self):
        # test that we can pickle a Parameter
        p = Parameter('a', 10, True, 0, 1)
        pkl = pickle.dumps(p)

        q = pickle.loads(pkl)

        assert_(p == q)

    def test_pickle_parameters(self):
        # test that we can pickle a Parameters object
        p = Parameters()
        p.add('a', 10, True, 0, 100)
        p.add('b', 10, True, 0, 100, 'a * sin(1)')
        p.update_constraints()
        p._asteval.symtable['abc'] = '2 * 3.142'

        pkl = pickle.dumps(p, -1)
        q = pickle.loads(pkl)

        q.update_constraints()
        assert_(p == q)
        assert_(not p is q)

        # now test if the asteval machinery survived
        assert_(q._asteval.symtable['abc'] == '2 * 3.142')

        # check that unpickling of Parameters is not affected by expr that
        # refer to Parameter that are added later on. In the following
        # example var_0.expr refers to var_1, which is a Parameter later
        # on in the Parameters OrderedDict.
        p = Parameters()
        p.add('var_0', value=1)
        p.add('var_1', value=2)
        p['var_0'].expr = 'var_1'
        pkl = pickle.dumps(p)
        q = pickle.loads(pkl)

    def test_set_symtable(self):
        # test that we use Parameter.set(value=XXX) and have
        # that new value be used in constraint expressions
        pars = Parameters()
        pars.add('x', value=1.0)
        pars.add('y', expr='x + 1')

        assert_(isclose(pars['y'].value, 2.0))
        pars['x'].set(value=3.0)
        assert_(isclose(pars['y'].value, 4.0))

    def test_dumps_loads_parameters(self):
        # test that we can dumps() and then loads() a Parameters
        pars = Parameters()
        pars.add('x', value=1.0)
        pars.add('y', value=2.0)
        pars['x'].expr = 'y / 2.0'

        dumps = pars.dumps()

        newpars = Parameters().loads(dumps)
        newpars['y'].value = 100.0
        assert_(isclose(newpars['x'].value, 50.0))

    def test_isclose(self):
        assert_(isclose(1., 1+1e-5, atol=1e-4, rtol=0))
        assert_(not isclose(1., 1+1e-5, atol=1e-6, rtol=0))
        assert_(isclose(1e10, 1.00001e10, rtol=1e-5, atol=1e-8))
        assert_(not isclose(0, np.inf))
        assert_(not isclose(-np.inf, np.inf))
        assert_(isclose(np.inf, np.inf))
        assert_(not isclose(np.nan, np.nan))

    def test_expr_with_bounds(self):
        "test an expression with bounds, without value"
        pars = Parameters()
        pars.add('c1', value=0.2)
        pars.add('c2', value=0.2)
        pars.add('c3', value=0.2)
        pars.add('csum', value=0.8)
        # this should not raise TypeError:
        pars.add('c4', expr='csum-c1-c2-c3', min=0, max=1)
        assert_(isclose(pars['c4'].value, 0.2))

    def test_invalid_expr_exceptions(self):
        "test if an exception is raised for invalid expressions (GH486)"""
        p1 = Parameters()
        p1.add('t', 2.0, min=0.0, max=5.0)
        p1.add('x', 10.0)
        with self.assertRaises(SyntaxError):
            p1.add('y', expr='x*t + sqrt(t)/')
        assert(len(p1['y']._expr_eval.error) > 0)
        p1.add('y', expr='x*t + sqrt(t)/3.0')
        p1['y'].set(expr='x*3.0 + t**2')
        assert('x*3' in p1['y'].expr)
        assert(len(p1['y']._expr_eval.error) == 0)
        with self.assertRaises(SyntaxError):
            p1['y'].set(expr='t+')
        assert(len(p1['y']._expr_eval.error) > 0)
        assert_almost_equal(p1['y'].value, 34.0)

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: FaustinCarter/lmfit-py

max_issues_repo_name: FaustinCarter/lmfit-py

max_forks_repo_name: FaustinCarter/lmfit-py

reason: import numpy,from numpy


text:
import json
import matplotlib

# matplotlib.use("Agg")

from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
from matplotlib.cbook import get_sample_data
import matplotlib.gridspec as gridspec
import numpy as np
from scipy.stats import entropy

def get_ordered_list_of_bitstrings(num_qubits):
    bitstrings = []
    for i in range(2 ** num_qubits):
        bitstring = "{0:b}".format(i)
        while len(bitstring) < num_qubits:
            bitstring = "0" + bitstring
        bitstrings.append(bitstring)
    return bitstrings

# Insert the path to your JSON file here
with open("workflow_result.json") as f:
    data = json.load(f)

# Extract target/measured bitstring distribution and distance measure values.
distances = []
minimum_distances = []
bitstring_distributions = []

current_minimum = 100000
number_of_qubits = 4
for step_id in data:
    step = data[step_id]
    # if step["stepName"] == "get-initial-parameters":
    #     number_of_qubits = int(
    #         eval(step["inputParam:ansatz-specs"])["number_of_qubits"]
    #     )
    ordered_bitstrings = get_ordered_list_of_bitstrings(number_of_qubits)
    if step["stepName"] == "get-bars-and-stripes-distribution":
        target_distribution = []
        for key in ordered_bitstrings:
            try:
                target_distribution.append(
                    step["distribution"]["bitstring_distribution"][key]
                )
            except:
                target_distribution.append(0)
        exact_distance_value = entropy(target_distribution)
        print(exact_distance_value)
    elif step["stepName"] == "optimize-circuit":
        for evaluation in step["qcbm-optimization-results"]["history"]:
            distances.append(evaluation["value"]["value"])
            current_minimum = min(current_minimum, evaluation["value"]["value"])
            minimum_distances.append(current_minimum)

            bitstring_dist = []
            for key in ordered_bitstrings:
                try:
                    bitstring_dist.append(
                        evaluation["artifacts"]["bitstring_distribution"][key]
                    )
                except:
                    bitstring_dist.append(0)
            bitstring_distributions.append(bitstring_dist)

fig = plt.figure(figsize=(16, 8))
gs = gridspec.GridSpec(nrows=8, ncols=3, width_ratios=[16, 1, 1])
ax1 = fig.add_subplot(gs[:4, 0])
ax2 = fig.add_subplot(gs[5:, 0])
axs = [fig.add_subplot(gs[i, 1]) for i in range(8)] + [
    fig.add_subplot(gs[i, 2]) for i in range(8)
]

evals = []
plotted_distances = []
plotted_min_distances = []
line_widths = []

images = [
    np.array([0, 0, 0, 0]),
    np.array([0, 0, 0, 1]),
    np.array([0, 0, 1, 0]),
    np.array([0, 0, 1, 1]),
    np.array([0, 1, 0, 0]),
    np.array([0, 1, 0, 1]),
    np.array([0, 1, 1, 0]),
    np.array([0, 1, 1, 1]),
    np.array([1, 0, 0, 0]),
    np.array([1, 0, 0, 1]),
    np.array([1, 0, 1, 0]),
    np.array([1, 0, 1, 1]),
    np.array([1, 1, 0, 0]),
    np.array([1, 1, 0, 1]),
    np.array([1, 1, 1, 0]),
    np.array([1, 1, 1, 1]),
]

def animate(i):
    evals.append(i)
    plotted_distances.append(distances[i])
    plotted_min_distances.append(minimum_distances[i])
    line_widths.append(1)
    ax1.clear()
    ax1.set(
        xlabel="Evaluation Index",
        ylabel="Clipped negative log-likelihood cost function",
    )
    ax1.set_ylim([exact_distance_value - 0.1, exact_distance_value + 1.5])
    ax1.scatter(
        evals, plotted_distances, color="green", linewidths=line_widths, marker="."
    )
    ax1.hlines(
        y=exact_distance_value,
        xmin=0,
        xmax=evals[-1],
        color="darkgreen",
        label="expected",
        alpha=0.8,
        linestyle="--",
    )
    ax1.legend(loc="upper right")
    ax1.plot(evals, plotted_min_distances, color="purple", linewidth=2)

    ax2.clear()
    ax2.set(xlabel="Bitstring", ylabel="Measured Probability")
    ax2.set_ylim([0, np.max(target_distribution) + 0.05])

    x_locations = np.arange(16)

    # Create the bars at the x locations with the height of the proportions for each bitstring

    ax2.bar(x_locations, bitstring_distributions[i], facecolor="green")
    ax2.bar(
        x_locations,
        target_distribution,
        facecolor="darkgreen",
        alpha=0.3,
        label="target",
    )

    # set the tick locations as the center of the bar, then set the label to be the bitstring
    ax2.set_xticks([x+0.4 for x in x_locations])
    ax2.set_xticklabels(ordered_bitstrings)
    ax2.legend(loc="upper right")

    if distances[i] == minimum_distances[i]:
        normalized_distribution = np.array(bitstring_distributions[i]) / max(
            bitstring_distributions[i]
        )
        for j in range(len(ordered_bitstrings)):
            axs[j].clear()
            axs[j].set_xticks(np.arange(-0.5, 2, 1), minor=True)
            axs[j].set_yticks(np.arange(-0.5, 2, 1), minor=True)
            axs[j].tick_params(axis="x", colors=(0, 0, 0, 0))
            axs[j].tick_params(axis="y", colors=(0, 0, 0, 0))

            axs[j].grid(which="minor", color="k", linestyle="-", linewidth=2)
            fading_factor = normalized_distribution[j]
            axs[j].imshow(
                (images[j].reshape((2, 2))),
                alpha=fading_factor,
                vmin=0,
                vmax=1,
                cmap="PiYG",
            )

    return tuple([ax1, ax2] + axs)

anim = FuncAnimation(fig, animate, frames=700, interval=1, repeat=False)

# # Set up formatting for the movie files
# Writer = animation.writers['ffmpeg']
# writer = Writer(fps=10, metadata=dict(artist='Me'), bitrate=1800)
# anim.save('qcbm_opt_700_iterations.mp4', writer=writer)

plt.show()

max_stars_repo_name: Ben-Foxman/z-quantum-qcbm

max_issues_repo_name: Ben-Foxman/z-quantum-qcbm

max_forks_repo_name: Ben-Foxman/z-quantum-qcbm

reason: import numpy,from scipy


text:
from __future__ import absolute_import, division, print_function, unicode_literals
from builtins import ascii, bytes, chr, dict, filter, hex, input, int, map, next, oct, open, pow, range, round, str, super, zip

import theano
import theano.tensor as T
import math
import numpy as np

from architecture.layer import *

floatX = theano.config.floatX

##################################################################################################################################
class IndexSeqsInput(Layer):

    #################################################################
    def __init__(self, name, vocab_size):
        super(IndexSeqsInput, self).__init__(name)

        self.vocab_size = vocab_size

    #################################################################
    def compile_params(self, dependent_sizes):
        assert dependent_sizes == []

        return self.vocab_size

    #################################################################
    def _get_model(self):
        return T.matrix(name=self.name, dtype='int16')

    #################################################################
    def get_training_model(self, dependent_models):
        assert dependent_models == []

        return self._get_model()

    #################################################################
    def get_testing_model(self, dependent_models):
        assert dependent_models == []

        return self._get_model()

##################################################################################################################################
class IntsInput(Layer):

    #################################################################
    def __init__(self, name):
        super(IntsInput, self).__init__(name)

    #################################################################
    def compile_params(self, dependent_sizes):
        if dependent_sizes != []:
            raise ValueError('Input layer cannot depend on anything.')

        return 1

    #################################################################
    def _get_model(self):
        return T.vector(name=self.name, dtype='int64')

    #################################################################
    def get_training_model(self, dependent_models):
        if dependent_models != []:
            raise ValueError('Input layer cannot depend on anything.')

        return self._get_model()

    #################################################################
    def get_testing_model(self, dependent_models):
        if dependent_models != []:
            raise ValueError('Input layer cannot depend on anything.')

        return self._get_model()

##################################################################################################################################
class VectorsInput(Layer):

    #################################################################
    def __init__(self, name, vector_size):
        super(VectorsInput, self).__init__(name)

        self.vector_size = vector_size

    #################################################################
    def compile_params(self, dependent_sizes):
        if dependent_sizes != []:
            raise ValueError('Input layer cannot depend on anything.')

        return self.vector_size

    #################################################################
    def _get_model(self):
        return T.matrix(name=self.name, dtype=floatX)

    #################################################################
    def get_training_model(self, dependent_models):
        if dependent_models != []:
            raise ValueError('Input layer cannot depend on anything.')

        return self._get_model()

    #################################################################
    def get_testing_model(self, dependent_models):
        if dependent_models != []:
            raise ValueError('Input layer cannot depend on anything.')

        return self._get_model()

max_stars_repo_name: mtanti/where-image

max_issues_repo_name: mtanti/where-image

max_forks_repo_name: mtanti/where-image

reason: import numpy,import theano


text:
import numpy as np
import torch

from dltranz.data_load.iterable_processing_dataset import IterableProcessingDataset

class CategorySizeClip(IterableProcessingDataset):
    def __init__(self, category_max_size, replace_value='max'):
        """

        Args:
            category_max_size: {field_name, max_size}
            replace_value: value for infrequent categories, int for specific value, 'max' for `category_max_size - 1`
        """
        super().__init__()

        self._category_max_size = category_max_size
        self._replace_value = replace_value

    def __iter__(self):
        for rec in self._src:
            features = rec[0] if type(rec) is tuple else rec

            # clip embeddings dictionary by max value
            for name, max_size in self._category_max_size.items():
                features[name] = self._smart_clip(features[name], max_size)
            yield rec

    def _smart_clip(self, values, max_size):
        if self._replace_value == 'max':
            return values.clip(0, max_size - 1)
        else:
            return torch.from_numpy(np.where((0 <= values) & (values < max_size), values, self._replace_value))

max_stars_repo_name: KirillVladimirov/pytorch-lifestream

max_issues_repo_name: KirillVladimirov/pytorch-lifestream

max_forks_repo_name: KirillVladimirov/pytorch-lifestream


text:
import pandas as pd
import numpy as np
import tensorflow as tf
import os

from utils import load_pickle, log, integer_to_timestamp_date_index

class Backtester:
    """
    Backtester:
        - gets data from pickles files
        - restores a trained network from a given checkpoint
        - makes predictions
        - creates 'signals' i.e what to buy at each rebalancing date for a list of strategies
        - runs the backtests for each strategies given market returns and signals
        - plots the backtesting results

    See main function Backtester.run_backtest()
    """

    def __init__(self, path_data, path_model_to_restore, start_date=20150101, end_date=20200101,
                 strategies=['10_max_long', '20_max_long', '2_max_long'], num_bins=None,
                 network_name='CondensedGoogLeNet'):
        """
        :param path_data:str: path to the folder containing the pickles files with the images
        :param path_model_to_restore:str: path for the model checkpoint used to make the predictions
        :param start_date:int: first date to consider for the backtest (format yyyymmdd)
        :param end_date:int:  last date to consider for the backtest
        :param strategies:list: list of strategies to backtest (see Backtester.__create_signals for supported names)
        :param network_name:str: name of the network to restore (used when restoring tensors)
        :param num_bins:int: number of bins in case we want to backtest a family of strategies
        """
        self._path_model_to_restore = path_model_to_restore
        self._name_network = network_name

        self._path_data = path_data
        self._start_date = start_date
        self._end_date = end_date

        self._strategies = strategies
        self._num_bins = num_bins

        # tensors that will be restored from checkpoint
        self._output = None
        self._x = None
        self._dropout = None
        self._phase_train = None

        # data retrieved or generated by the backtester
        self._df_all_data = None
        self._df_permnos_to_buy = None
        self.df_strats = None

        self.__LOGGER_ENV = 'backtest'

    def run_backtest(self):
        """
        central function of the backtester, instanciates self.df_strats containing backtest_results for each strategy
        """
        X = self.get_df_all_data()
        pred = self.run_predictions(X)
        self._make_df_for_bckt(pred)

        df_backt_results = self._run_strategies()

        self.df_strats = self._format_df_strats(df_backt_results)

    def _run_strategies(self):
        """ Uses the dataframe with all data to make a backtesting
        : Example below
               PERMNO     RET     long      hold     short
        date
        20181030  83387.0   1.03  0.149458  0.353090  0.497452
        20181030  84207.0   1.03  0.404111  0.385122  0.210767

        return a dataframe with weighted returns
        :return:                  1_decile_long 2_decile_long
                    20151001        0.92554       0.933354
                    20151008        1.00058        0.99515

        """
        df_data = self._df_all_data
        strategies = self._strategies
        self._df_permnos_to_buy = self.__create_signals(df_data, strategies=self._strategies, bins=self._num_bins)
        # self._df_permnos_to_buy = self._df_permnos_to_buy.shift(1).dropna() # because if we buy at i we get returns in i+1, not here
        # In the data returns at date i are the returns we would get if we buy at i

        df_rets = df_data[['RET', 'PERMNO']].pivot_table(columns='PERMNO', index='date')
        df_rets = df_rets.dropna().RET

        # Run the backtest with the returns and the predictions
        self.log('Running Backtest')
        df_results = pd.DataFrame(columns=strategies, index=self._df_permnos_to_buy.index)
        for date in self._df_permnos_to_buy.index:
            for strat in strategies:
                list_permnos_to_buy, weight_permnos_to_buy = self._df_permnos_to_buy.loc[date][strat]
                if len(list_permnos_to_buy) == 0:
                    df_results.loc[date][strat] = 1.  # we buy nothing
                else:
                    df_results.loc[date][strat] = np.sum(
                        np.array(df_rets.loc[date][list_permnos_to_buy]) * np.array(weight_permnos_to_buy))
        return df_results

    def __create_signals(self, df_data, strategies=['10_max_long', '20_max_long'], bins=None, equi_weighted=True):
        """
        Creates the list of stocks to buy in each strategy for each rebalacing date

        :param df_data: Example
                PERMNO     PRC      long      hold     short
        date
        20181030  83387.0    8.89  0.149458  0.353090  0.497452
        20181030  84207.0  119.15  0.404111  0.385122  0.210767
        :param strategies: list of strategies to consider

        :return:   dataframe1 with the permnos to buy for each date and the associated weights

        # df_permnos_to_buy
                                           4_max_long
        20181120  [85567.0, 83486.0, 83728.0, 83630.0],[0.42,0.08 ...]
        20181220  [85539.0, 83762.0, 83844.0, 83885.0],[0.05,0.048 ...]

        """
        self.log('Creating signals with strategies {}'.format(strategies))

        # we sort the predictions
        data_sorted_long_all_dates = df_data.reset_index().sort_values(by='long', ascending=False).set_index('date')
        data_sorted_short_all_dates = df_data.reset_index().sort_values(by='short', ascending=False).set_index('date')
        df_temp = df_data.copy()
        df_temp['long_short'] = df_temp.long - 0.4 * df_temp.short
        data_sorted_long_short_all_dates = df_temp.reset_index().sort_values(by='long_short',
                                                                             ascending=False).set_index('date')

        # empty dataframe that will be used to store results
        df_permnos_to_buy = pd.DataFrame(columns=strategies, index=sorted(set(df_data.index)))

        def get_weights(data_sorted, list_permnos_to_buy):
            if equi_weighted:
                weight_permnos = np.ones(len(list_permnos_to_buy)) / len(list_permnos_to_buy)
            else:  # make weights proportional to the signal we consider
                weight_permnos = data_sorted.loc[list_permnos_to_buy].values.T[0]
                weight_permnos = weight_permnos / weight_permnos.sum()
            return weight_permnos

        # calculate proportional weights and stocks to buy in each bin
        def perform_bin_strategy(data_sorted):
            bin = int(strat.split('_')[0])  # format must be '4_bins_...' for the 4th bin
            n_stocks = len(data_sorted.index)
            list_permnos_to_buy = list(data_sorted.index)[
                                  round((bin - 1) * n_stocks / bins):round(bin * n_stocks / bins)]
            weight_permnos = get_weights(data_sorted, list_permnos_to_buy)
            return list_permnos_to_buy, weight_permnos

        for date in df_data.index:
            data_sorted_long = data_sorted_long_all_dates.loc[date].set_index('PERMNO')[['long']]
            data_sorted_short = data_sorted_short_all_dates.loc[date].set_index('PERMNO')[['short']]
            data_sorted_long_short = data_sorted_long_short_all_dates.loc[date].set_index('PERMNO')[['long_short']]

            for strat in strategies:
                if strat == '10_max_long':
                    list_permnos_to_buy = list(data_sorted_long.index[:10])
                    weight_permnos = get_weights(data_sorted_long, list_permnos_to_buy)
                elif strat == '20_max_long':
                    list_permnos_to_buy = list(data_sorted_long.index[:20])
                    weight_permnos = get_weights(data_sorted_long, list_permnos_to_buy)
                elif strat == '2_max_long':
                    list_permnos_to_buy = list(data_sorted_long.index[:2])
                    weight_permnos = get_weights(data_sorted_long, list_permnos_to_buy)
                elif strat == 'threshold':
                    list_permnos_to_buy = list(df_data[df_data.long >= 0.75].index)
                    weight_permnos = get_weights(data_sorted_long, list_permnos_to_buy)
                elif 'bins_long' in strat:
                    list_permnos_to_buy, weight_permnos = perform_bin_strategy(data_sorted_long)
                elif 'bins_short' in strat:
                    list_permnos_to_buy, weight_permnos = perform_bin_strategy(data_sorted_short)
                elif 'bins_long_short' in strat:
                    list_permnos_to_buy, weight_permnos = perform_bin_strategy(data_sorted_long_short)
                else:
                    raise NotImplementedError('The strategy {} is not implemented'.format(strat))
                df_permnos_to_buy.loc[date][strat] = list_permnos_to_buy, weight_permnos
        self.log('Signals created')
        return df_permnos_to_buy

    def _make_df_for_bckt(self, pred):
        """
        Instantiates self._df_all_data with a Dataframe containing prices and predictions and sorts it
        chronologically
                * Example:

                PERMNO     RET      long      hold     short
        date
        20181030  83387.0    1.023  0.149458  0.353090  0.497452
        20181030  84207.0  1.02  0.404111  0.385122  0.210767

        :param pred: numpy array shape (N_samples,3)
        :return: Nothing
        """
        self.log('Joining prices data and predictions')
        df_signals = pd.DataFrame(index=self._df_all_data.index, data=pred, columns=['long', 'hold', 'short'])
        self._df_all_data = pd.concat([self._df_all_data, df_signals], axis=1)
        self._df_all_data = self._df_all_data.sort_index()

    def restore_output_op(self, sess, latest=False):
        """
        Restores tensors : x and output and instanciates self._x and self._output
        :param sess: tf.Session
        :param latest: bool: set to True to retrive the latest checkpoint.meta file in the specified folder

        Instantiates self._output, self._x, self._phase_train, self._dropout  with tensors retrieved from checkpoint file
        """
        saver = tf.train.import_meta_graph(self._path_model_to_restore)
        folder = os.path.dirname(self._path_model_to_restore)
        if latest:
            file_to_restore = tf.train.latest_checkpoint(folder)
        else:
            file_to_restore = self._path_model_to_restore.replace('.meta', '')
        saver.restore(sess, file_to_restore)

        graph = tf.get_default_graph()

        self._output = graph.get_tensor_by_name('{}/output:0'.format(self._name_network))
        self._x = graph.get_tensor_by_name('{}/x:0'.format(self._name_network))
        self._phase_train = graph.get_tensor_by_name('{}/phase_train:0'.format(self._name_network))
        self._dropout = graph.get_tensor_by_name('{}/dropout:0'.format(self._name_network))

    def run_predictions(self, X):
        """
        Restores a model saved with tensorflow and predicts signals
        :param X: numpy array: shape (N_samples,42,42,5) (if we consider 42 pixels and 5 channels)
        :return: predictions as a numpy array shape (N_samples,3)
        """
        self.log('Restoring model and making the predictions')
        sess = tf.Session()
        graph = tf.get_default_graph()
        self.restore_output_op(sess)
        self.log('Model Restored, launching output operation')

        # we run predictions in batches to avoid killing the kernel
        size_1_image = np.prod(X[0].shape)
        limit_size = 20 * 42 * 42 * 5
        size_batch = int(limit_size / size_1_image) + 1
        pred = np.zeros((0, 3))
        self.log('To avoid killing the kernel, we run predictions in {} batches'.format(round(len(X) / size_batch)))
        for batch in range(0, len(X), size_batch):
            X_batch = X[batch:batch + size_batch]
            pred_batch = sess.run(self._output,
                                  feed_dict={self._x: X_batch, self._phase_train: False, self._dropout: 0.})
            pred = np.concatenate([pred, pred_batch])

        self.log('Predictions computed')
        sess.close()
        return pred

    def plot_backtest(self):
        self.df_strats.plot(grid=True,figsize=(16,12),title='Backtest for different strategies')

    def get_df_all_data(self):
        """
        - Unpickles data files and gets samples
        - Keeps only the dates we want

        - Instanciates self._df_all_data: Example
                                  PERMNO       RET      long      hold     short
                    date
                    20181101  79237.0  0.952770  0.019922  0.458449  0.521629
                    20181101  78877.0  0.991597  0.080810  0.658793  0.260397
        """
        samples_path = self._path_data
        # list of files in the folder: samples_paths
        list_file_names = os.listdir(samples_path)
        assert len(list_file_names) >= 1, 'The number of files in the folder {} is probably 0, it must be >=1'.format(
            samples_path)
        log('******* Getting data from folder: {}, Nb of files : {}, First file : {}'.format(samples_path,
                                                                                len(list_file_names),
                                                                                list_file_names[0]),
            self.__LOGGER_ENV)

        for i, file_name in enumerate(list_file_names):
            path = os.path.join(samples_path, file_name)
            df_all_data_one_batch = load_pickle(path, logger_env=self.__LOGGER_ENV)
            self.log(
                'first_date: {}, last_date: {}'.format(df_all_data_one_batch.index[0], df_all_data_one_batch.index[-1]))
            # Keeping only the dates we want for the backtest
            df_all_data_one_batch = df_all_data_one_batch.loc[self._start_date:self._end_date]
            self.log('new first_date: {}, new last_date: {}'.format(df_all_data_one_batch.index[0],
                                                                    df_all_data_one_batch.index[-1]))

            # in each of the pickle files the data is sorted in chronologic order we resort just in case
            df_all_data_one_batch = df_all_data_one_batch.sort_index()
            X = np.concatenate([[sample for sample in df_all_data_one_batch['sample'].values]],
                               axis=0)  # need this to get an array
            # we only keep returns and permnos (which is the stock identifier)
            df_all_data_one_batch_for_bckt = df_all_data_one_batch[['PERMNO', 'RET']]
            if i == 0:
                X_res = X
                df_all_data = df_all_data_one_batch_for_bckt
            else:
                X_res = np.concatenate([X_res, X])
                df_all_data = pd.concat([df_all_data, df_all_data_one_batch_for_bckt])
        # now we must not change the order in X_res and df_all_data
        self._df_all_data = df_all_data
        return X_res

    @staticmethod
    def _format_df_strats(df_backt_results):
        """
        :param df_backt_results:pd.DataFrame: results of the backtest data by date
        :return: pd.DataFrame cumulative results with a datetime index and an extra column 'Cash' which is the SPX
        """
        df_backt_results = df_backt_results.shift(
            1).dropna()  # because returns in reality do not happen just when we buy, but on next week
        df_strats = df_backt_results.astype(np.float64)
        df_strats = integer_to_timestamp_date_index(df_strats)
        df_strats = df_strats.cumprod() / df_strats.iloc[0]

        # we get SPX data and add it to the dataframe
        df_spx = pd.read_csv('data/^GSPC.csv', index_col=0, usecols=['Date', 'Close'])
        df_spx.index.name = 'date'
        df_spx.columns = ['Cash']  # Columns that will be used to compute sharpe ratio
        df_spx.index = pd.to_datetime(df_spx.index)
        df_strats = pd.merge(df_strats, df_spx, how='left', on='date')
        df_strats['Cash'] = df_strats['Cash'] / df_strats['Cash'].iloc[0]
        return df_strats

    def log(self, msg):
        log(msg, environment=self.__LOGGER_ENV)

max_stars_repo_name: kathy-jzh/DeepLearningAllocationCNN

max_issues_repo_name: kathy-jzh/DeepLearningAllocationCNN

max_forks_repo_name: kathy-jzh/DeepLearningAllocationCNN


text:
from keras.applications.inception_v3 import InceptionV3
from keras.models import Model
from pickle import dump, load
import numpy as np
from keras.applications.inception_v3 import preprocess_input
from keras.preprocessing import image

class InceptionV3_Encode:

    encoding_train = None
    encoding_test = None
    train_features = None

    def __init__(self,images_directory,train_features_path,test_features_path,train_images,test_images):
        # load the inception v3 model
        model = InceptionV3(weights='imagenet')
        # create a new model,by removing the last layer (output layer) from the inception v3
        model_new = Model(model.input,model.layers[-2].output)
        # call the function to encode all the train images
        self.encoding_train = self.encode_images_inceptionv3(images_directory,train_images,model_new,train_features_path)
        # call the function to encode all the test images
        self.encoding_test = self.encode_images_inceptionv3(images_directory,test_images,model_new,test_features_path)
        self.train_features = load(open(train_features_path,'rb'))

    #encode the images by Model
    def encode_images_inceptionv3(self,images_directory,images,model_new,features_path):
        encoding_temp = {}
        for image_1 in images:
            encoding_temp[image_1[len(images_directory):]] = self.encode(image_1, model_new)
        # save the bottleneck train features to disk
        with open(features_path, 'wb') as encoded_pickle:
            dump(encoding_temp, encoded_pickle)
        return encoding_temp

    #function to encode a given image into vector of size(2048,)
    def encode(self,image_1,model_new):
        image_2 = self.inception_v3_preprocess(image_1)#preprocess the image
        fea_vec = model_new.predict(image_2)#get the encoding vector for the image
        fea_vec = np.reshape(fea_vec,fea_vec.shape[1])#reshape from (1, 2048) to (2048, )
        return fea_vec

    def inception_v3_preprocess(self,image_path):
        #convert all the images to size 299*299 as expected by the inception v3 model
        image_3 = image.load_img(image_path,target_size=(299,299))
        #convert PIL image to numpy array of 3-dimensions
        x = image.img_to_array(image_3)
        #add one more dimension
        x = np.expand_dims(x,axis=0)
        #preprocess the images using preprocess_input() from inception module
        x = preprocess_input(x)
        return x

max_stars_repo_name: jianhaifeng/UVISONER-DatasetUtils

max_issues_repo_name: jianhaifeng/UVISONER-DatasetUtils

max_forks_repo_name: jianhaifeng/UVISONER-DatasetUtils


text:
from xevo import eobj

import numpy as np
import time

class gotozero(eobj):
  """eobj that tries to move nd ints to zero"""
  def __init__(s,q=None):
    if q is None:q=[10,10]
    #print("initialising gotozero with")
    #print(q)
    #time.sleep(1)
    s.initial()
    s.q=q
  def __str__(s):
    return str(s.q)
  def __add__(a,b):
    return gotozero([aa+bb for aa,bb in zip(a,b)])
  def shallmaximize(s):return False

  def randomize(s):
    return s.copy()#that kind of a bodge i guess
    return gotozero([np.random.randint(-10,11) for i in range(len(s.q))])
  def mutate(s):
    s0=s
    s=s.copy()
    i=np.random.randint(len(s.q))
    s.q[i]+=np.random.choice([-1,1])
    return s
  def calcstrength(s):
    return np.sum([zw**2 for zw in s.q])
  def _copy(s):
    return gotozero(q=[zw for zw in s.q])

max_stars_repo_name: psorus/xevoplus

max_issues_repo_name: psorus/xevoplus

max_forks_repo_name: psorus/xevoplus


text:
import sys

import cv2
import numpy as np
from image_utils import find_outer_corner, point_on_axis, point_respect_line

class Stitch:
    def __init__(self, left_image, right_image):
        self.im_left = left_image
        self.im_right = right_image
        self.result = None

    def compute_matches(self):
        pass

    def merge(self):
        self.result = np.concatenate((self.im_left, self.im_right), axis=1)
        return self.result

    def transform(self):
        pts1 = np.float32([[180, 450], [5, 116], [338, 396]])
        pts2 = np.float32([[205, 370], [618, 126], [311, 407]])
        matrix = cv2.getAffineTransform(pts1, pts2)
        cols, rows, ch = self.im_right.shape

        warped_affine = cv2.warpAffine(self.im_right, matrix, (rows * 2, cols * 2))
        top_pnt, right_pnt, bottom_pnt, left_pnt = find_outer_corner(warped_affine)

        return warped_affine[top_pnt[1]:bottom_pnt[1], left_pnt[0]:right_pnt[0]]

    def perspective_trans(self):
        pts1 = np.float32([[250, 392], [412, 96], [640, 112], [0, 282]])
        pts2 = np.float32([[480, 640], [0, 0], [480, 0], [0, 640]])
        M = cv2.getPerspectiveTransform(pts1, pts2)
        self.im_left = cv2.warpPerspective(self.im_left, M, (480, 640))

    def split_image_1(self, image, point1=(618, 126), point2=(254, 388)):
        pnt1, pnt2 = point_on_axis(point1, point2)
        cv2.line(image, pnt1, pnt2, (0, 255, 0), 2)
        cv2.imshow('line', image)

    def split_image(self):
        pnt1, pnt2 = point_on_axis((618, 126), (254, 388))
        h, w = self.im_left.shape[:2]
        for x in range(w):
            for y in range(h):
                if point_respect_line((x, y), (pnt1, pnt2)) > 0:
                    self.im_left[y, x] = (255, 255, 255)

        pnt1, pnt2 = point_on_axis((5, 116), (267, 424))

        h, w = self.im_right.shape[:2]
        for x in range(w):
            for y in range(h):
                if point_respect_line((x, y), (pnt1, pnt2)) < 0:
                    self.im_right[y, x] = (255, 255, 255)

    def stitch(self):
        # self.split_image_1(self.im_left)
        self.split_image()
        # self.perspective_trans()
        cv2.imwrite('left.jpg', self.im_left)
        cv2.imwrite('right.jpg', self.im_right)
        # wraped = self.transform()

        # H = matcher.match(self.im_left, wraped)

        # print("Homography is : ", H)

        # cv2.imshow("warped", wraped)
        # cv2.waitKey(0)
        return None

    def imshow(self):
        cv2.imshow("left image", self.im_left)
        cv2.imshow("right Image", self.im_right)
        cv2.waitKey(0)

def read_images_from_file(file_paths='images.txt'):
    fp = open(file_paths, 'r')
    filenames = [each.rstrip('\r\n') for each in fp.readlines()]
    print(filenames)
    images = [cv2.resize(cv2.imread(each), (640, 480)) for each in filenames]
    return images

if __name__ == '__main__':
    try:
        args = sys.argv[1]
    except:
        args = "files_1.txt"

    image_list = read_images_from_file(args)

    s = Stitch(image_list[0], image_list[1])
    result = s.stitch()
    s.imshow()
    # cv2.namedWindow('result', cv2.WINDOW_AUTOSIZE)

    # cv2.imshow('result', result)
    cv2.waitKey(0)
    print("done")
    print("image written")
    cv2.destroyAllWindows()

max_stars_repo_name: DewMaple/opencv-learning

max_issues_repo_name: DewMaple/opencv-learning

max_forks_repo_name: DewMaple/opencv-learning


text:
from astropy.tests.pytest_plugins import *

def pytest_addoption(parser):
    parser.addoption("--remote-data", action="store_true",
                     help="run tests with online data")
    parser.addoption("--open-files", action="store_true",
                     help="fail if any test leaves files open")

    parser.addoption("--doctest-plus", action="store_true",
                     help="enable running doctests with additional "
                     "features not found in the normal doctest "
                     "plugin")

    parser.addoption("--doctest-rst", action="store_true",
                     help="enable running doctests in .rst documentation")

    parser.addini("doctest_plus", "enable running doctests with additional "
                  "features not found in the normal doctest plugin")

    parser.addini("doctest_norecursedirs",
                  "like the norecursedirs option but applies only to doctest "
                  "collection", type="args", default=())

    parser.addini("doctest_rst",
                  "Run the doctests in the rst documentation",
                  default=False)

    parser.addoption('--repeat', action='store',
                     help='Number of times to repeat each test')

    parser.addoption("--atomic-dataset", dest='atomic-dataset', default=None,
                     help="filename for atomic dataset")

def pytest_report_header(config):

    stdoutencoding = getattr(sys.stdout, 'encoding') or 'ascii'

    s = "\n"
    if six.PY2:
        args = [x.decode('utf-8') for x in config.args]
    elif six.PY3:
        args = config.args
    s += "Running tests in {0}.\n\n".format(" ".join(args))

    from platform import platform
    plat = platform()
    if isinstance(plat, bytes):
        plat = plat.decode(stdoutencoding, 'replace')
    s += "Platform: {0}\n\n".format(plat)
    s += "Executable: {0}\n\n".format(sys.executable)
    s += "Full Python Version: \n{0}\n\n".format(sys.version)

    s += "encodings: sys: {0}, locale: {1}, filesystem: {2}".format(
        sys.getdefaultencoding(),
        locale.getpreferredencoding(),
        sys.getfilesystemencoding())
    if sys.version_info < (3, 3, 0):
        s += ", unicode bits: {0}".format(
            int(math.log(sys.maxunicode, 2)))
    s += '\n'

    s += "byteorder: {0}\n".format(sys.byteorder)
    s += "float info: dig: {0.dig}, mant_dig: {0.dig}\n\n".format(
        sys.float_info)

    import numpy
    s += "numpy: {0}\n".format(numpy.__version__)

    try:
        import scipy
        s += "scipy: {0}\n".format(scipy.__version__)
    except:
        s += "scipy: not available\n"

    try:
        import pandas
        s += "pandas: {0}\n".format(pandas.__version__)
    except:
        s += "pandas: not available\n"

    try:
        import astropy
    except:
        s += "astropy: not available\n"
    else:
        s += "astropy: {0}\n".format(astropy.__version__)

    try:
        import yaml
    except:
        s += "yaml: not available\n"
    else:
        s += "yaml: {0}\n".format(yaml.__version__)

    try:
        import cython
    except:
        s += "cython: not available\n"
    else:
        s += "cython: {0}\n".format(cython.__version__)

    try:
        import h5py.version
        s += "h5py: {0}\n".format(h5py.version.version)
    except:
        s += "h5py: not available\n"

    try:
        import matplotlib
        s += "matplotlib: {0}\n".format(matplotlib.__version__)
    except:
        s += "matplotlib: not available\n"

    try:
        import IPython
    except:
        s += "ipython: not available\n"
    else:
        s += "ipython: {0}\n".format(IPython.__version__)

    special_opts = ["remote_data", "pep8"]
    opts = []
    for op in special_opts:
        if getattr(config.option, op, None):
            opts.append(op)
    if opts:
        s += "Using Astropy options: {0}.\n".format(" ".join(opts))

    if six.PY3 and (config.getini('doctest_rst') or config.option.doctest_rst):
        s += "Running doctests in .rst files is not supported on Python 3.x\n"

    if not six.PY3:
        s = s.encode(stdoutencoding, 'replace')

    return s

import os

@pytest.fixture
def atomic_data_fname():
    atomic_data_fname = pytest.config.getvalue("atomic-dataset")
    if atomic_data_fname is None:
        pytest.skip('--atomic_database was not specified')
    else:
        return os.path.expandvars(os.path.expanduser(atomic_data_fname))

from tardis.atomic import AtomData

@pytest.fixture
def kurucz_atomic_data(atomic_data_fname):
    atomic_data = AtomData.from_hdf5(atomic_data_fname)

    if atomic_data.md5 != '21095dd25faa1683f4c90c911a00c3f8':
        pytest.skip('Need default Kurucz atomic dataset '
                    '(md5="21095dd25faa1683f4c90c911a00c3f8"')
    else:
        return atomic_data

@pytest.fixture
def included_he_atomic_data():
    import os, tardis
    atomic_db_fname = os.path.join(tardis.__path__[0], 'tests', 'data',
                                   'chianti_he_db.h5')
    return AtomData.from_hdf5(atomic_db_fname)

max_stars_repo_name: wkerzendorf/tardis

max_issues_repo_name: wkerzendorf/tardis

max_forks_repo_name: wkerzendorf/tardis

reason: import numpy,import scipy,import astropy,from astropy


text:
import argparse
import torch
import os
import numpy as np
import datasets.crowd as crowd
from models import vgg19

parser = argparse.ArgumentParser(description='Test ')
parser.add_argument('--device', default='0', help='assign device')
parser.add_argument('--crop-size', type=int, default=512,
                    help='the crop size of the train image')
parser.add_argument('--model-path', type=str, default='pretrained_models/model_qnrf.pth',
                    help='saved model path')
parser.add_argument('--data-path', type=str,
                    default='data/QNRF-Train-Val-Test',
                    help='saved model path')
parser.add_argument('--dataset', type=str, default='qnrf',
                    help='dataset name: qnrf, nwpu, sha, shb')

args = parser.parse_args()

os.environ['CUDA_VISIBLE_DEVICES'] = args.device  # set vis gpu
device = torch.device('cuda')

model_path = args.model_path
crop_size = args.crop_size
data_path = args.data_path
if args.dataset.lower() == 'qnrf':
    dataset = crowd.Crowd_qnrf(os.path.join(data_path, 'test'), crop_size, 8, method='val')
elif args.dataset.lower() == 'nwpu':
    dataset = crowd.Crowd_nwpu(os.path.join(data_path, 'val'), crop_size, 8, method='val')
elif args.dataset.lower() == 'sha' or args.dataset.lower() == 'shb':
    dataset = crowd.Crowd_sh(os.path.join(data_path, 'test_data'), crop_size, 8, method='val')
else:
    raise NotImplementedError
dataloader = torch.utils.data.DataLoader(dataset, 1, shuffle=False,
                                         num_workers=1, pin_memory=True)
transform=transforms.Compose([
                       transforms.ToTensor(),transforms.Normalize(mean=[0.485, 0.456, 0.406],
                                     std=[0.229, 0.224, 0.225]),
                   ])
model = vgg19()
model.to(device)
model.load_state_dict(torch.load(model_path, device))
image_errs = []
for inputs, count, name in dataloader:
    inputs = inputs.to(device)
    assert inputs.size(0) == 1, 'the batch size should equal to 1'
    with torch.set_grad_enabled(False):
        outputs, _ = model(inputs)
    img_err = count[0].item() - torch.sum(outputs).item()

    print(name, img_err, count[0].item(), torch.sum(outputs).item())
    image_errs.append(img_err)

image_errs = np.array(image_errs)
mse = np.sqrt(np.mean(np.square(image_errs)))
mae = np.mean(np.abs(image_errs))
print('{}: mae {}, mse {}\n'.format(model_path, mae, mse))

max_stars_repo_name: Xinxinatg/DM-Count

max_issues_repo_name: Xinxinatg/DM-Count

max_forks_repo_name: Xinxinatg/DM-Count


text:
# -*- coding: utf-8 -*-
"""
Calculate pair distribution functions.
"""

MINDISTANCE = 2.0 # shorter distances are ignored
PDFCUTOFF = 1.0   # g(r) = 0 if r < PDFCUTOFF

__all__ = ["PDF"]

import numpy as np
import math
from util.logger import Logger
import sys
from core.calculation.discretization import DiscretizationCache
import itertools
import os.path
from config.configuration import config

logger = Logger("statistics.pdf")
logger.setstream("default", sys.stdout, Logger.DEBUG)

class PDF(object):
    """
    Calculate pair distribution functions for atoms and cavities.
    """

    def __init__(self, *args):
        """
        Create a sample from atom and cavity positions and smooth them to
        get the PDFs

        The constructor can be called in two ways:

        - ``PDF(results)`` :
            retrieve the data from :class:`core.data.Results`

        - ``PDF(positions, elements, cavitycenters, volume)`` :
            use the given arrays and the volume object
        """
        if len(args) == 1:
            results = args[0]
            positions = results.atoms.positions
            elements = results.atoms.elements
            volume = results.atoms.volume
            if results.domains is not None:
                centers = results.domains.centers
                cachedir = os.path.expanduser(config.Path.cache_dir)
                cachepath = os.path.join(cachedir, 'discretization_cache.hdf5')
                dcache = DiscretizationCache(cachepath)
                disc = dcache.get_discretization(volume, results.resolution)
                centers = map(disc.discrete_to_continuous, centers)
            else:
                centers = []
        elif len(args) == 4:
            positions, elements, centers, volume = args
        else:
            raise TypeError("PDF expects 1 or 4 parameters")

        self.positions = np.array(positions, copy=False)
        self.elements = np.array(elements, dtype="|S4", copy=False)
        self.centers = np.array(centers, copy=False)
        self.volume = volume
        self.num_atoms = np.where(self.elements != "cav")[0].size
        self.numberdensity = float(self.num_atoms) / self.volume.volume

        self.stats = self._genstats(self.positions,
                                    self.elements,
                                    self.centers,
                                    self.volume)

    def pdf(self, elem1, elem2, cutoff=None, h=None, kernel=None):
        """
        Calculate a smoothed pair distribution function between the elements
        `elem1` and `elem2`.

        **Parameters:**
            `elem1`, `elem2` :
                Chemical element, e.g. 'Ge', 'Te' or 'cav' for cavities

            `h` :
                Smoothing parameter. The greater `h` is,
                the more the function is smoothed.

            `kernel` :
                Smoothing kernel

        **Returns:**
            Python function that represents the pair distribution function.
            It also accepts Numpy arrays as input.
            Returns `None` if the given elements do not exist or if there is
            not enough data to create the function.
        """
        if kernel is None:
            #kernel = Kernels.epanechnikov
            kernel = Kernels.epanechnikov

        data = None
        for s in self.stats:
            if set((elem1.lower(), elem2.lower())) \
                    == set((s[0].lower(), s[1].lower())):
                data = s[2]
                break;
        if data is None:
            logger.debug("No statistical data for '{}-{}' found.".format(
                    elem1, elem2))
            return None # TODO: raise Exception

        if cutoff is None:
            cutoff = data.max()
        sel = np.where(np.logical_and(data > MINDISTANCE, data <= cutoff))[0]
        sel = data[sel]
        if h == 0:
            return sel
        if len(sel) < 2:
            logger.debug("Not enough data for '{}-{}' in cutoff={} range.".format(elem1, elem2, cutoff))
            return None # TODO: raise Exception

        if h is None:
            ## magic constant
            ## minimizes the first peak
            #h = min(0.5, 0.5772778 * sel.min())
            h = 0.4

        #if h > 0.9 * sel.min():
        #    logger.debug("Bandwidth {} above threshold. Setting to {}.".format(h, 0.9 * sel.min()))
        #    h = 0.9 * sel.min()
        kde = Functions.KDE(sel, h=h, kernel=kernel)

        def wfunc(r):
            y = np.zeros_like(r)
            i = np.where(np.abs(r) > PDFCUTOFF)[0]
            y[i] = self.volume.volume / (data.size * 4 * math.pi * r[i]**2)
            return y

        return Functions.Product(kde, wfunc)

    @staticmethod
    def _correlatedistance(pos1, pos2=None, volume=None):
        """
        Measure the distances between coordinates.

        - ``_correlatedistance(pos1)`` :
            Correlate the coordinates in `pos1` with each other

        - ``_correlatedistance(pos1, pos2)`` :
            Correlate each coordinate in pos1 with each coordinate in pos2

        **Returns:**
            A sorted sample.
        """
        if volume is None:
            distance = lambda x, y: np.abs(y-x)
        else:
            distance = volume.get_distance

        if pos2 is None:
            n = pos1.shape[0]
            a1 = np.vstack([pos1[0:-i, :] for i in range(1, n)])
            a2 = np.vstack([pos1[i:, :] for i in range(1, n)])
        else:
            n1 = pos1.shape[0]
            n2 = pos2.shape[0]
            a1 = np.hstack([pos1] * n2).reshape(n1 * n2, 3)
            a2 = np.vstack([pos2] * n1)
        samples = np.linalg.norm(distance(a1, a2), 2, axis=-1)
        samples.sort()

        return samples

    @classmethod
    def _genstats(cls, positions, elements, centers, volume=None):
        """
        Iterate through the elements and cavities and correlate the
        distances between atoms of that element/cavity and atoms of
        the other elements.
        """
        elemlist = np.unique(elements).tolist()
        if len(centers) > 0 and not "cav" in elemlist:
            elemlist.append("cav")
        pos = [None] * len(elemlist)
        for i, e in enumerate(elemlist):
            pos[i] = positions[np.where(elements == e)[0], :]

        if len(centers) > 0:
            i = elemlist.index("cav")
            if pos[i].size > 0:
                pos[i] = np.vstack((pos[i], centers))
            else:
                pos[i] = centers

        stats = []
        for i, e1 in enumerate(elemlist):
            for j in range(i, len(elemlist)):
                e2 = elemlist[j]
                if i == j:
                    s = cls._correlatedistance(pos[i], volume=volume)
                else:
                    s = cls._correlatedistance(pos[i], pos[j], volume=volume)
                if len(s) > 1:
                    stats.append((e1, e2, s))

        return stats

class Functions(object):
    """
    Utility class with Callables
    """

    class KDE(object):
        """
        Kernel density estimation. Calculate a convolution
        of delta pulses with a smoothing kernel.
        """

        def __init__(self, x, y=None, h=1.0, kernel=None):
            if y is None:
                y = np.ones_like(x)
            if kernel is None:
                kernel = Kernels.gauss
            self.x = x
            self.y = y
            self.h = h
            self.kernel = kernel

        def __call__(self, p):
            result = np.zeros_like(p)
            if len(self.x) <= len(p):
                p = np.asarray(p)
                for xi, yi in itertools.izip(self.x, self.y):
                    result += yi * self.kernel((p - xi) / self.h) / self.h
            else:
                x = np.asarray(self.x)
                y = np.asarray(self.y)
                for i, pi in enumerate(p):
                    result[i] = np.sum(y * self.kernel((pi - x) / self.h) / self.h)

            return result

    class Product(object):
        """
        Product of two functions
        """

        def __init__(self, f, g):
            self.f = f
            self.g = g

        def __call__(self, x):
            return self.f(x) * self.g(x)

class Kernels(object):
    """
    Utility class with smoothing kernels
    """

    @staticmethod
    def gauss(x):
        c = 1.0 / math.sqrt(2.0 * math.pi)
        return c * np.exp(-0.5 * x**2)

    @staticmethod
    def compact(x):
        c = 2.25228362104
        if not isinstance(x, np.ndarray):
            if abs(x) < 1.0:
                return c * math.exp(1.0 / (x**2 - 1.0))
            else:
                return 0.0
        else:
            i = np.where(np.abs(x) < 1.0)[0]
            y = np.zeros(x.size)
            y[i] = c * np.exp(1.0 / (x[i]**2 - 1.0))
            return y

    @staticmethod
    def triang(x):
        c = 1.0
        if not isinstance(x, np.ndarray):
            if abs(x) < 1.0:
                return c * (1.0 - abs(x))
            else:
                return 0.0
        else:
            i = np.where(np.abs(x) < 1.0)[0]
            y = np.zeros(x.size)
            y[i] = c * (1.0 - np.abs(x[i]))
            return y

    @staticmethod
    def quad(x):
        c = 0.5
        if not isinstance(x, np.ndarray):
            if abs(x) < 1.0:
                return c
            else:
                return 0.0
        else:
            i = np.where(np.abs(x) < 1.0)[0]
            y = np.zeros(x.size)
            y[i] = c
            return y

    @staticmethod
    def posquad(x):
        c = 1.0
        if not isinstance(x, np.ndarray):
            if 0 <= x < 1.0:
                return c
            else:
                return 0.0
        else:
            i = np.where(np.logical_and(x >= 0.0, x < 1.0))[0]
            y = np.zeros(x.size)
            y[i] = c
            return y

    @staticmethod
    def negquad(x):
        c = 1.0
        if not isinstance(x, np.ndarray):
            if -1.0 < x <= 0.0:
                return c
            else:
                return 0.0
        else:
            i = np.where(np.logical_and(x > -1.0, x <= 0.0))[0]
            y = np.zeros(x.size)
            y[i] = c
            return y

    @staticmethod
    def epanechnikov(x):
        if not isinstance(x, np.ndarray):
            if abs(x) < 1.0:
                return 3 / 4 * (1.0 - x**2)
            else:
                return 0.0
        else:
            i = np.where(np.abs(x) < 1.0)[0]
            y = np.zeros_like(x)
            y[i] = 3.0 / 4.0 * (1.0 - x[i]**2)
            return y

    @staticmethod
    def bandwidth(n, d=1):
        """
        Scott's factor for bandwidth estimation
        """
        return (4.0 / 3.0 * n) ** (-1.0 / (d + 4.0))

class _TestPDF(object):
    @staticmethod
    def continuous_coordinates(coords, volume, resolution):
        cachedir = os.path.expanduser(config.Path.cache_dir)
        cachepath = os.path.join(cachedir, 'discretization_cache.hdf5')
        dcache = DiscretizationCache(cachepath)
        disc = dcache.get_discretization(volume, resolution)

        return np.array(map(disc.discrete_to_continuous, coords))

    @staticmethod
    def plotfunc(pdf, e1, e2, px, h, *args):
        gr = pdf.pdf(e1, e2, h=h)
        plt.plot(px, gr(px), *args, label=str(h))
        py = gr(px)
        m = np.argmax(py)
        print "h={}, xi={}, g({}) = {}".format(gr.f.h, gr.f.x.min(), px[m], py[m])

    @classmethod
    def plotpdf(cls, pdf, e1, e2):
        px = np.linspace(0, 2, 1000)
        plt.figure()
        cls.plotfunc(pdf, e1, e2, px, 0.25, "g--")
        cls.plotfunc(pdf, e1, e2, px, 0.5, "r-")
        #cls.plotfunc(pdf, e1, e2, px, 1.0, "b--")
        plt.legend(loc=0)
        plt.title("{}-{}".format(e1, e2))

    @classmethod
    def run(cls):
        import core.calculation as calculation
        calc = calculation.Calculation("../results")
        filename = "../xyz/structure_c.xyz"
        resolution = 64
        frame = 9
        #filename = "../xyz/GST_111_128_bulk.xyz"
        #resolution = 256
        #frame = 0
        settings = calculation.CalculationSettings(
                {filename : [frame]},
                resolution, True, False, False)
        print "calculating..."
        res = calc.calculate(settings)[0][0]
        print "generating statistics..."
        pdf = PDF(res)
        #centers = cls.continuous_coordinates(res.domains.centers,
        #                                     res.atoms.volume,
        #                                     res.resolution)
        #pdf = PDF(res.atoms.positions, res.atoms.elements,
        #                  centers, res.atoms.volume)

        print "plotting..."
        #cls.plotpdf(pdf, "Ge", "Ge")
        #cls.plotpdf(pdf, "Ge", "Te")
        cls.plotpdf(pdf, "cav", "cav")
        plt.show()

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    _TestPDF.run()
    #x = np.linspace(-2, 2, 200)
    #plt.plot(x, Kernels.gauss(x))
    #plt.plot(x, Kernels.compact(x))
    #plt.plot(x, Kernels.triang(x))
    #plt.plot(x, Kernels.quad(x))
    #plt.plot(x, Kernels.epanechnikov(x))
    #plt.show()

max_stars_repo_name: sciapp/pyMolDyn

max_issues_repo_name: sciapp/pyMolDyn

max_forks_repo_name: sciapp/pyMolDyn


text:
# -*- coding: utf-8 -*-
"""Trend estimation for timeseries.

Estimate and remove trend in a timeseries.

In this package, trend removal is in service of isolating and
estimating periodic (non-trend) variation. "trend" is in the sense of
Cleveland's STL decomposition -- a lowpass smoothing of
the data that, when removed from the original series, preserves
original seasonal variation.  Detrending is accomplishd by a coarse
fitted spline, mean or median filters, or a fitted line.

Input samples are assumed evenly-spaced from an anomaly-free
continuous-time signal.

"""
from __future__ import division
import numpy as np
from scipy import stats
from scipy.interpolate import LSQUnivariateSpline
from .periodogram import periodogram_peaks

def fit_trend(data, kind="spline", period=None, ptimes=2):
    """Fit a trend for a possibly noisy, periodic timeseries.

    Trend may be modeled by a line, cubic spline, or mean or median
    filtered series.

    Parameters
    ----------
    data : ndarray
        list of observed values
    kind : string ("mean", "median", "line", "spline", None)
        if mean, apply a period-based mean filter
        if median, apply a period-based median filter
        if line, fit a slope to median-filtered data.
        if spline, fit a piecewise cubic spline to the data
        if None, return zeros
    period : number
        seasonal periodicity, for filtering the trend.
        if None, will be estimated.
    ptimes : number
        multiple of period to use as smoothing window size

    Returns
    -------
    trend : ndarray

    """
    if kind is None:
        return np.zeros(len(data)) + np.mean(data)
    if period is None:
        period = guess_trended_period(data)
    window = (int(period * ptimes) // 2) * 2 - 1 # odd window
    if kind == "median":
        filtered = aglet(median_filter(data, window), window)
    elif kind == "mean":
        filtered = aglet(mean_filter(data, window), window)
    elif kind == "line":
        filtered = line_filter(data, window)
    elif kind == "spline":
        nsegs = len(data) // (window * 2) + 1
        filtered = aglet(spline_filter(data, nsegs), window)
    else:
        raise Exception("adjust_trend: unknown filter type {}".format(kind))
    return filtered

def guess_trended_period(data):
    """return a rough estimate of the major period of trendful data.

    Periodogram wants detrended data to score periods reliably. To do
    that, apply a broad median filter based on a reasonable maximum
    period.  Return a weighted average of the plausible periodicities.

    Parameters
    ----------
    data : ndarray
        list of observed values, evenly spaced in time.

    Returns
    -------
    period : int

    """
    max_period = min(len(data) // 3, 512)
    broad = fit_trend(data, kind="median", period=max_period)
    peaks = periodogram_peaks(data - broad)
    if peaks is None:
        return max_period
    periods, scores, _, _ = zip(*peaks)
    period = int(round(np.average(periods, weights=scores)))
    return period

def aglet(src, window, dst=None):
    """straigten the ends of a windowed sequence.

    Replace the window/2 samples at each end of the sequence with
    lines fit to the full window at each end.  This boundary treatment
    for windowed smoothers is better behaved for detrending than
    decreasing window sizes at the ends.

    Parameters
    ----------
    src : ndarray
        list of observed values
    window : int
        odd integer window size (as would be provided to a windowed smoother)
    dst : ndarray
        if provided, write aglets into the boundaries of this array.
        if dst=src, overwrite ends of src in place. If None, allocate result.

    Returns
    -------
    dst : ndarray
        array composed of src's infield values with aglet ends.

    """
    if dst is None:
        dst = np.array(src)
    half = window // 2
    leftslope = stats.theilslopes(src[: window])[0]
    rightslope = stats.theilslopes(src[-window :])[0]
    dst[0:half] = np.arange(-half, 0) * leftslope + src[half]
    dst[-half:] = np.arange(1, half + 1) * rightslope + src[-half - 1]
    return dst

def median_filter(data, window):
    """Apply a median filter to the data.

    This implementation leaves partial windows at the ends untouched

    """
    filtered = np.copy(data)
    for i in range(window // 2, len(data) - window // 2):
        filtered[i] = np.median(data[max(0, i - window // 2) : i + window // 2 + 1])
    return filtered

def mean_filter(data, window):
    """Apply a windowed mean filter to the data.

    This implementation leaves partial windows at the ends untouched

    """
    filtered = np.copy(data)
    cum = np.concatenate(([0], np.cumsum(data)))
    half = window // 2
    filtered[half : -half] = (cum[window:] - cum[:-window]) / window
    return filtered

def line_filter(data, window):
    """fit a line to the data, after filtering"""
    # knock down seasonal variation with a median filter first
    half = window // 2
    coarse = median_filter(data, window)[half : -half] # discard crazy ends
    slope, _, lower, upper = stats.theilslopes(coarse)
    if lower <= 0.0 and upper >= 0.0:
        filtered = np.zeros(len(data)) + np.median(data)
    else:
        intercept = np.median(data) - (len(data) - 1) / 2 * slope
        filtered = slope * np.arange(len(data)) + intercept
    return filtered

def spline_filter(data, nsegs):
    """Detrend a possibly periodic timeseries by fitting a coarse piecewise
       smooth cubic spline

    Parameters
    ----------
    data : ndarray
        list of observed values
    nsegs : number
        number of spline segments

    Returns
    -------
    filtered : ndarray

    """
    index = np.arange(len(data))
    nknots = max(2, nsegs + 1)
    knots = np.linspace(index[0], index[-1], nknots + 2)[1:-2]
    return LSQUnivariateSpline(index, data, knots)(index)