text:
"""Orbital maneuvers.

"""
import numpy as np
from astropy import units as u

from poliastro.util import norm

class Maneuver(object):
    r"""Class to represent a Maneuver.

    Each ``Maneuver`` consists on a list of impulses \\(\\Delta v_i\\)
    (changes in velocity) each one applied at a certain instant \\(t_i\\).
    You can access them directly indexing the ``Maneuver`` object itself.

    >>> man = Maneuver((0 * u.s, [1, 0, 0] * u.km / u.s),
    ... (10 * u.s, [1, 0, 0] * u.km / u.s))
    >>> man[0]
    (<Quantity 0. s>, <Quantity [1., 0., 0.] km / s>)
    >>> man.impulses[1]
    (<Quantity 10. s>, <Quantity [1., 0., 0.] km / s>)

    """

    def __init__(self, *impulses):
        r"""Constructor.

        Parameters
        ----------
        impulses : list
            List of pairs (delta_time, delta_velocity)

        Notes
        -----
        TODO: Fix docstring, \*args convention

        """
        self.impulses = impulses
        # HACK: Change API or validation code
        _dts, _dvs = zip(*impulses)
        self._dts, self._dvs = self._initialize(
            [(_dt * u.one).value for _dt in _dts] * (_dts[0] * u.one).unit,
            [(_dv * u.one).value for _dv in _dvs] * (_dvs[0] * u.one).unit,
        )
        try:
            if not all(len(dv) == 3 for dv in self._dvs):
                raise TypeError
        except TypeError:
            raise ValueError("Delta-V must be three dimensions vectors")

    @u.quantity_input(dts=u.s, dvs=u.m / u.s)
    def _initialize(self, dts, dvs):
        return dts, dvs

    def __getitem__(self, key):
        return self.impulses[key]

    @classmethod
    def impulse(cls, dv):
        """Single impulse at current time.

        """
        return cls((0 * u.s, dv))

    @classmethod
    def hohmann(cls, orbit_i, r_f):
        """Compute a Hohmann transfer between two circular orbits.

        """
        # TODO: Check if circular?
        ss_i = orbit_i.state
        r_i = ss_i.a
        v_i = ss_i.v
        k = ss_i.attractor.k
        R = r_f / r_i
        dv_a = ((np.sqrt(2 * R / (1 + R)) - 1) * v_i).decompose()
        dv_b = (-(1 - np.sqrt(2 / (1 + R))) / np.sqrt(R) * v_i).decompose()
        t_trans = (np.pi * np.sqrt((r_i * (1 + R) / 2) ** 3 / k)).decompose()
        return cls((0 * u.s, dv_a), (t_trans, dv_b))

    @classmethod
    def bielliptic(cls, orbit_i, r_b, r_f):
        """Compute a bielliptic transfer between two circular orbits.

        """
        ss_i = orbit_i.state
        r_i = ss_i.a
        v_i = ss_i.v
        k = ss_i.attractor.k
        R = r_f / r_i
        Rs = r_b / r_i
        dv_a = ((np.sqrt(2 * Rs / (1 + Rs)) - 1) * v_i).decompose()
        dv_b = (
            -np.sqrt(2 / Rs) * (np.sqrt(1 / (1 + Rs / R)) - np.sqrt(1 / (1 + Rs))) * v_i
        ).decompose()
        dv_c = (-(np.sqrt(2 * Rs / (R + Rs)) - 1) / np.sqrt(R) * v_i).decompose()
        t_trans1 = (np.pi * np.sqrt((r_i * (1 + Rs) / 2) ** 3 / k)).decompose()
        t_trans2 = (np.pi * np.sqrt((r_i * (R + Rs) / 2) ** 3 / k)).decompose()
        return cls((0 * u.s, dv_a), (t_trans1, dv_b), (t_trans2, dv_c))

    def get_total_time(self):
        """Returns total time of the maneuver.

        """
        total_time = sum(self._dts, 0 * u.s)
        return total_time

    def get_total_cost(self):
        """Returns total cost of the maneuver.

        """
        dvs = [norm(dv) for dv in self._dvs]
        return sum(dvs, 0 * u.km / u.s)

max_stars_repo_name: aOrionis/poliastro

max_issues_repo_name: aOrionis/poliastro

max_forks_repo_name: aOrionis/poliastro


text:
import pickle
import numpy as np
import torch

def save_object(fname, obj):
    #the function is used to save some data in class or object in .pkl file
    with open(fname, 'wb') as out_file:
        pickle.dump(obj, out_file)
    out_file.close()

def load_object(fname):
    #the function is used to read the data in .pkl file
    with open(fname, 'rb') as in_file:
        return pickle.load(in_file)

class RandomLabel():
    #to retain the same training label
    def __init__(self, ratio, mode = 'create'):
        self.ratio = ratio
        self.mode = mode
        self.label = []

    def grab(self, index):
        if index > len(self.label):
            raise IndexError('Please check the iteration of dataloader.')

        self.mode = 'stop'
        return self.label[index]

    def generate(self, data):
        if self.mode == 'create':
            new = self._random_sample(data)
            self.label.append(new)

            return new
        else:
            raise RuntimeError('Please check the RandomLabel api usage.')

    def _random_sample(self, data):
        size = int(data.size(0) * self.ratio)
        rand_arr = torch.randint(0, 10, data.size())
        _, index = torch.rand(data.size(0)).sort()
        for i in range(size):
            data[index[i]] = rand_arr[index[i]]

        return data

ext: py

reason: import numpy


text:
import numpy as np
import matplotlib.pyplot as plt
import os
from edibles.utils.edibles_spectrum import EdiblesSpectrum
from edibles.utils.edibles_oracle import EdiblesOracle
from edibles.continuum import Continuum
from edibles import DATADIR

pythia = EdiblesOracle()
dirs = [f for f in os.listdir(DATADIR) if not f.startswith(".")]
figure_output = DATADIR.replace("/data", "/Na3302Plots")
if not os.path.isdir(figure_output):
    os.mkdir(figure_output)

for i, sightline in enumerate(dirs):
    print("="*20)
    print("Now working on Sightline " + sightline)
    print("Sightline {i}/{j}".format(i=i+1, j=len(dirs)))
    print("="*20)

    matched_files = pythia.getFilteredObsList(object=[sightline.replace("HD", "HD ")], OrdersOnly=True, Wave=3302)
    figure_output_sightline = os.path.join(figure_output, sightline)
    if not os.path.isdir(figure_output_sightline):
        os.mkdir(figure_output_sightline)

    for file in matched_files:
        sp = EdiblesSpectrum(file)
        sp.getSpectrum(xmin=3301, xmax=3305)
        date = sp.datetime.date()
        order = file.split("_")[-1].replace(".fits", "")

        # build a 4 anchor points spline
        cont = Continuum(sp, method="spline", n_anchors=4, plot=False, verbose=0)
        params = cont.model.guess(sp.flux, x=sp.wave)
        out = cont.model.eval(params=params, x=sp.wave)
        cont.add_to_csv(
            user="Na Worker", comments="1st order continuum for Na 3302"
        )

        # save plots
        filename = "_".join([sightline, str(date), order]) + ".png"
        filename = "/".join([figure_output_sightline, filename])
        plt.plot(sp.wave, sp.flux)
        plt.plot(sp.wave, out)
        plt.ylabel(sightline)
        plt.xlabel("  ".join([str(date), order]))
        plt.savefig(filename)
        plt.close()

max_stars_repo_name: jancami/edibles

max_issues_repo_name: jancami/edibles

max_forks_repo_name: jancami/edibles


text:
import numpy as np
import mglearn

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# Ignore Warnings
import warnings
warnings.filterwarnings("ignore", category=Warning)

X, y = mglearn.datasets.make_wave(n_samples=60)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

lr = LinearRegression().fit(X_train, y_train)

print("lr.coef_:{}".format(lr.coef_))
print("lr.intercept_:{}".format(lr.intercept_))

print("Training set score:{:.2f}".format(lr.score(X_train, y_train)))
print("Test set score:{:.2f}".format(lr.score(X_test, y_test)))

X, y = mglearn.datasets.load_extended_boston()

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
lr = LinearRegression().fit(X_train, y_train)

print("Training set score:{:.2f}".format(lr.score(X_train, y_train)))
print("Test set score:{:.2f}".format(lr.score(X_test, y_test)))

max_stars_repo_name: FDUJiaG/PyML-Course

max_issues_repo_name: FDUJiaG/PyML-Course

max_forks_repo_name: FDUJiaG/PyML-Course


text:
import numpy as np
from tqdm import tqdm
import scipy.io as scio

import utils.io as io
from data.hico.hico_constants import HicoConstants

class ConvertMat2Json():
    def __init__(self,const):
        self.const = const
        self.anno = scio.loadmat(self.const.anno_mat)
        self.anno_bbox = scio.loadmat(self.const.anno_bbox_mat)

    def create_hoi_list(self):
        num_hoi = self.anno['list_action'].shape[0]
        hoi_list = [None]*num_hoi
        for i in range(num_hoi):
            hoi_list[i] = {
                'id': str(i+1).zfill(3),
                'object': self.anno['list_action'][i,0][0][0],
                'verb': self.anno['list_action'][i,0][1][0],
            }

        return hoi_list

    def get_image_size(self,i,subset):
        W = self.anno_bbox[f'bbox_{subset}'][0,i][1][0,0][0][0,0]
        H = self.anno_bbox[f'bbox_{subset}'][0,i][1][0,0][1][0,0]
        C = self.anno_bbox[f'bbox_{subset}'][0,i][1][0,0][2][0,0]
        image_size = [int(v) for v in [H,W,C]]
        return image_size

    def get_hoi_bboxes(self,i,subset):
        num_hois = self.anno_bbox[f'bbox_{subset}'][0,i][2].shape[1]
        hois = [None]*num_hois
        for j in range(num_hois):
            hoi_data = self.anno_bbox[f'bbox_{subset}'][0,i][2][0,j]

            hoi_id = str(hoi_data[0][0,0]).zfill(3)

            num_boxes = hoi_data[1].shape[1]
            human_bboxes = [None]*num_boxes
            for b in range(num_boxes):
                human_bboxes[b] = \
                    [int(hoi_data[1][0,b][k][0,0]-1) for k in [0,2,1,3]]

            num_boxes = hoi_data[2].shape[1]
            object_bboxes = [None]*num_boxes
            for b in range(num_boxes):
                object_bboxes[b] = \
                    [int(hoi_data[2][0,b][k][0,0]-1) for k in [0,2,1,3]]

            connections = (hoi_data[3]-1).tolist()

            invis = int(hoi_data[4][0,0])

            hois[j] = {
                'id': hoi_id,
                'human_bboxes': human_bboxes,
                'object_bboxes': object_bboxes,
                'connections': connections,
                'invis': invis,
            }

        return hois

    def create_anno_list(self):
        anno_list = []
        for subset in ['train','test']:
            print(f'Adding {subset} data to anno list ...')
            num_samples = self.anno[f'anno_{subset}'].shape[1]
            for i in tqdm(range(num_samples)):
                image_jpg = self.anno[f'list_{subset}'][i][0][0]

                if image_jpg.endswith('.jpg'):
                    global_id = image_jpg[:-4]
                else:
                    assert(False), 'Image extension is not .jpg'

                anno = {
                    'global_id': global_id,
                    'image_path_postfix': f'{subset}2015/{image_jpg}',
                    'image_size': self.get_image_size(i,subset),
                    'hois': self.get_hoi_bboxes(i,subset)
                }

                anno['pos_hoi_ids'] = [str(k[0]+1).zfill(3) for k in \
                    np.argwhere(self.anno[f'anno_{subset}'][:,i]==1).tolist()]
                anno['neg_hoi_ids'] = [str(k[0]+1).zfill(3) for k in \
                    np.argwhere(self.anno[f'anno_{subset}'][:,i]==-1).tolist()]

                anno_list.append(anno)

        return anno_list

    def convert(self):
        print('Creating anno list ...')
        anno_list = self.create_anno_list()
        io.dump_json_object(anno_list,self.const.anno_list_json)

        print('Creating hoi list ...')
        hoi_list = self.create_hoi_list()
        io.dump_json_object(hoi_list,self.const.hoi_list_json)

        print('Creating object list ...')
        object_list = sorted(list(set([hoi['object'] for hoi in hoi_list])))
        for i,obj in enumerate(object_list):
            object_list[i] = {
                'id': str(i+1).zfill(3),
                'name': obj
            }

        io.dump_json_object(object_list,self.const.object_list_json)

        print('Creating verb list ...')
        verb_list = sorted(list(set([hoi['verb'] for hoi in hoi_list])))
        for i,verb in enumerate(verb_list):
            verb_list[i] = {
                'id': str(i+1).zfill(3),
                'name': verb
            }

        io.dump_json_object(verb_list,self.const.verb_list_json)

if __name__=='__main__':
    hico_const = HicoConstants()
    converter = ConvertMat2Json(hico_const)
    converter.convert()

reason: import numpy,import scipy


text:
from pywizard.userSettings import settings
from pywizard.Buffer import Buffer
import scipy as sp

class Segmenter(object):
    def __init__(self, buf, windowWidth):
        milliseconds = settings.frameRate
        self.size = int(sp.ceil(buf.sampleRate / 1e3 * milliseconds))
        self.buf = buf
        self.windowWidth = windowWidth

    def eachSegment(self):
        length = self.numberOfSegments()
        for i in range (0, length):
            samples = self.samplesForSegment(i)
            buf = Buffer(samples = samples, size=self.sizeForWindow, sampleRate=self.buf.sampleRate)
            yield (buf, i)

    def samplesForSegment(self, index):
        length = self.sizeForWindow()

        samples = self.buf.samples[ index * self.size : index * self.size + length ]
        samples = sp.append(samples, [0]*(length-len(samples)))

        return samples

    def sizeForWindow(self):
        return int(self.size * self.windowWidth)

    def numberOfSegments(self):
        return int(sp.ceil(float(self.buf.size) / float(self.size)))

reason: import scipy


text:
"""
    Dummy conftest.py for omlt.

    If you don't know what this is for, just leave it empty.
    Read more about conftest.py under:
    - https://docs.pytest.org/en/stable/fixture.html
    - https://docs.pytest.org/en/stable/writing_plugins.html
"""

# import pytest

import numpy as np

def get_neural_network_data(desc):
    """
    Return input and test data for a neural network.

    Parameters
    ----------
    desc : string
        model name. One of 131 or 2353.
    """
    if desc == "131":
        # build data with 1 input and 1 output and 500 data points
        x = np.random.uniform(-1, 1, 500)
        y = np.sin(x)
        x_test = np.random.uniform(-1, 1, 5)
        return x, y, x_test

    elif desc == "2353":
        # build data with 2 inputs, 3 outputs, and 500 data points
        np.random.seed(42)
        x = np.random.uniform([-1, 2], [1, 3], (500, 2))
        y1 = np.sin(x[:, 0] * x[:, 1])
        y2 = x[:, 0] + x[:, 1]
        y3 = np.cos(x[:, 0] / x[:, 1])
        y = np.column_stack((y1, y2, y3))
        x_test = np.random.uniform([-1, 2], [1, 3], (5, 2))
        return x, y, x_test

    return None

max_stars_repo_name: michaelbynum/OMLT

max_issues_repo_name: michaelbynum/OMLT

max_forks_repo_name: michaelbynum/OMLT


text:
#======================================================================
#
#     This routine solves an infinite horizon growth model
#     with dynamic programming and sparse grids
#
#     The model is described in Scheidegger & Bilionis (2017)
#     https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2927400
#
#     external libraries needed:
#     - IPOPT (https://projects.coin-or.org/Ipopt)
#     - PYIPOPT (https://github.com/xuy/pyipopt)
#     - scikit-learn GPR (https://scikit-learn.org)
#
#     Simon Scheidegger, 01/19
#
#     Cameron Gordon, 11/21 - conversion to Python3
#======================================================================

import nonlinear_solver_pre as solver     #solves opt. problems for terminal VF
import nonlinear_solver_post as solvpost   #solves opt. problems during VFI
from parameters import *                      #parameters of model
import iteration_pre as iter_pre              #interface to sparse grid library/terminal VF
import iteration_post as iter_post    #interface to sparse grid library/iteration
import postprocessing as post                 #computes the L2 and Linfinity error of the model
import numpy as np

#======================================================================
# Start with Value Function Iteration

for i in range(numstart, numits):
# terminal value function
    if (i==1):
        print("start with Value Function Iteration")
        iter_pre.GPR_init(i)

    else:
        print("Now, we are in Value Function Iteration step", i)
        iter_post.GPR_iter(i)

#======================================================================
print("===============================================================")
print(" ")
print(" Computation of a growth model of dimension ", n_agents ," finished after ", numits, " steps")
print(" ")
print("===============================================================")
#======================================================================

# compute errors
avg_err=post.ls_error(n_agents, numstart, numits, No_samples_postprocess)

#======================================================================
print("===============================================================")
print(" ")
#print " Errors are computed -- see error.txt"
print(" ")
print("===============================================================")
#======================================================================

max_stars_repo_name: uq-aibe/spir-oz

max_issues_repo_name: uq-aibe/spir-oz

max_forks_repo_name: uq-aibe/spir-oz


text:
def NAFAgent(self, nb_states, nb_actions):
    # Get the environment and extract the number of actions.

    # Build all necessary models: V, mu, and L networks.
    V_model = Sequential()
    V_model.add(Flatten(input_shape=(1, nb_states)))
    V_model.add(Dense(16))
    V_model.add(Activation('relu'))
    V_model.add(Dense(16))
    V_model.add(Activation('relu'))
    V_model.add(Dense(16))
    V_model.add(Activation('relu'))
    V_model.add(Dense(1))
    V_model.add(Activation('linear'))
    print(V_model.summary())

    mu_model = Sequential()
    mu_model.add(Flatten(input_shape=(1, nb_states)))
    mu_model.add(Dense(16))
    mu_model.add(Activation('relu'))
    mu_model.add(Dense(16))
    mu_model.add(Activation('relu'))
    mu_model.add(Dense(16))
    mu_model.add(Activation('relu'))
    mu_model.add(Dense(nb_actions))
    mu_model.add(Activation('linear'))
    print(mu_model.summary())

    action_input = Input(shape=(1, nb_actions,), name='action_input')
    observation_input = Input(shape=(1, nb_states), name='observation_input')
    x = Concatenate()([Flatten()(action_input), Flatten()(observation_input)])
    x = Dense(32)(x)
    x = Activation('relu')(x)
    x = Dense(32)(x)
    x = Activation('relu')(x)
    x = Dense(32)(x)
    x = Activation('relu')(x)
    # x = Dense((1))(x)
    x = Dense(((nb_actions * nb_actions + nb_actions) // 2))(x)
    x = Activation('linear')(x)
    L_model = Model(inputs=[action_input, observation_input], outputs=x)
    print(L_model.summary())

    # Finally, we configure and compile our agent. You can use every built-in Keras optimizer and
    # even the metrics!
    self.model = V_model
    processor = PendulumProcessor()
    memory = SequentialMemory(limit=100000, window_length=1)
    random_process = OrnsteinUhlenbeckProcess(theta=.15, mu=0., sigma=.3, size=nb_actions)
    agent = NAFAgent(nb_actions=nb_actions, V_model=self.model, L_model=L_model, mu_model=mu_model,
                     memory=memory, nb_steps_warmup=100, random_process=random_process,
                     gamma=.99, target_model_update=1e-3, processor=processor)
    return agent
    # agent.compile(Adam(lr=.001, clipnorm=1.), metrics=['mae'])

def agentDDP(self, nb_states, nb_actions):
    # Next, we build a very simple model.
    actor = Sequential()
    actor.add(Dense(16, activation='relu', input_shape=(1, nb_states)))
    actor.add(Dense(16, activation='relu'))
    actor.add(Dense(16, activation='relu'))
    actor.add(Dense(16, activation='relu'))
    actor.add(Dense(nb_actions, activation='linear'))
    print(actor.summary())

    action_input = Input(shape=(nb_actions,), name='action_input')
    observation_input = Input(shape=(nb_states,), name='observation_input')

    x1 = keras.layers.Dense(8)(action_input)
    x2 = keras.layers.Dense(8)(observation_input)
    # flattened_observation = Flatten()(observation_input)
    x = Concatenate()([x1, x2])
    x = Dense(32)(x)
    x = Activation('relu')(x)
    x = Dense(32)(x)
    x = Activation('relu')(x)
    x = Dense(32)(x)
    x = Activation('relu')(x)
    x = Dense(1)(x)
    x = Activation('linear')(x)
    critic = Model(inputs=[action_input, observation_input], outputs=x)
    print(critic.summary())

    # Finally, we configure and compile our agent. You can use every built-in Keras optimizer and
    # even the metrics!
    memory = SequentialMemory(limit=100000, window_length=1)
    random_process = OrnsteinUhlenbeckProcess(size=nb_actions, theta=.15, mu=0., sigma=.3)
    self.model = actor
    agent = DDPGAgent(nb_actions=nb_actions, actor=self.model, critic=critic, critic_action_input=action_input,
                      memory=memory, nb_steps_warmup_critic=100, nb_steps_warmup_actor=100,
                      random_process=random_process, gamma=.99, target_model_update=1e-3)
    agent.compile(Adam(lr=.001, clipnorm=1.), metrics=['mae'])
    return agent
    # Okay, now it's time to learn something! We visualize the training here for show, but this
    # slows down training quite a lot. You can always safely abort the training prematurely using
    # Ctrl + C.
    # agent.fit(env, nb_steps=50000, visualize=True, verbose=1, nb_max_episode_steps=200)
    #
    # # After training is done, we save the final weights.
    # agent.save_weights('ddpg_{}_weights.h5f'.format(ENV_NAME), overwrite=True)
    #
    # # Finally, evaluate our algorithm for 5 episodes.
    # agent.test(env, nb_episodes=5, visualize=True, nb_max_episode_steps=200)

import math
import gym
from gym import spaces
import numpy as np
import coordinateOperation
import fileOperation

class dobotGym(gym.Env):
    def __init__(self, plot = False, save= True, emulateOculus = True, episodeLength = 15):
        self.coordinateOperationInstance = coordinateOperation.coordinateOperation(plot = plot, save = save,emulateOculus = emulateOculus)
        self.coordinateOperationInstance.recording = False

        # Angle at which to fail the episode

        # Angle limit set to 2 * theta_threshold_radians so failing observation
        # is still within bounds.
        maxStep = 150
        #self.action_space = spaces.Box(low=np.array([-maxStep, -maxStep, -maxStep]), high=np.array([maxStep, maxStep, maxStep]),dtype=np.float32)

        self.action_space = spaces.Box(np.array([self.coordinateOperationInstance.minX,self.coordinateOperationInstance.minY,self.coordinateOperationInstance.minZ]),
                                       np.array([self.coordinateOperationInstance.maxX,self.coordinateOperationInstance.maxY,self.coordinateOperationInstance.maxZ])
                                            ,dtype=np.float32)

        self.observation_space = spaces.Box(np.array([self.coordinateOperationInstance.minX,self.coordinateOperationInstance.minY,self.coordinateOperationInstance.minZ]),
                                       np.array([self.coordinateOperationInstance.maxX,self.coordinateOperationInstance.maxY,self.coordinateOperationInstance.maxZ])
                                            ,dtype=np.float32)

        self.viewer = None
        self.state = None

        self.steps_beyond_done = None

        self.episodeLength = episodeLength -2
        self.episodeStep = 0

    def step(self, action):
        print("Action", action, self.action_space.low, self.action_space.high)
        action = np.clip(action, self.action_space.low, self.action_space.high)
        print(action)
        self.coordinateOperationInstance.coordinateFromOculusToDobotTranslation() #translating coordinates from oculus to dobot system
        self.coordinateOperationInstance.setDiffPositionToMove(action[0],action[1],action[2])
        self.state = self.coordinateOperationInstance.moveDobotToPreparedPosition()

        reward = 1/(self.coordinateOperationInstance.actualDiffXYZ+1)
        if self.coordinateOperationInstance.grip is True:
            done = False
        else:
            done = True

        if (self.episodeStep > self.episodeLength):
            done = True

        info = {
            'Diff': self.coordinateOperationInstance.actualDiffXYZ,
        }
        self.episodeStep += 1
        return np.array(self.state), reward, done, info

    def reset(self):
        self.coordinateOperationInstance.endOfMoving()
        self.episodeStep = 0
        self.coordinateOperationInstance.preparationForMoving()
        self.state = np.array([self.coordinateOperationInstance.rightXLastDobot,self.coordinateOperationInstance.rightYLastDobot,self.coordinateOperationInstance.rightZLastDobot])
        while(self.coordinateOperationInstance.grip is False):
            pass
        self.coordinateOperationInstance.startRecording()
        return np.array(self.state)

    def close(self):
        self.coordinateOperationInstance.endOfMoving()


text:
from nose.tools import *
import numpy as np

from nanosims_analysis.importer import Importer
from nanosims_analysis.data_structures import IsotopeData
from nanosims_analysis.data_structures import RatioData

class TestClass:

    @classmethod
    def setup_class(cls):
        cls.numerator_isotope = IsotopeData("numerator",
                                             np.array([[[782, 645, 275],
                                                        [733, 875, 590],
                                                        [902, 708, 739]],

                                                       [[789, 433, 716],
                                                        [ 72, 367, 390],
                                                        [465,  58, 995]]]))

        cls.denominator_isotope = IsotopeData("denominator",
                                              np.array([[[500, 475, 370],
                                                         [779,  37, 769],
                                                         [ 33, 796, 741]],

                                                        [[502, 618, 573],
                                                         [963, 473, 858],
                                                         [ 31, 734,   0]]]))

    def test_ratio_init(self):
        testRatio = RatioData("test_ratio",
                              self.numerator_isotope,
                              self.denominator_isotope)

        ans = np.array([[[ 1.564     ,  1.35789474,  0.74324324],
                         [ 0.94094994, 23.64864865,  0.76723017],
                         [27.33333333,  0.88944724,  0.99730094]],

                        [[ 1.57171315,  0.70064725,  1.2495637 ],
                         [ 0.07476636,  0.77589852,  0.45454545],
                         [15.        ,  0.07901907,         0]]])
        assert_true(np.allclose(ans, testRatio.get_data()))

    @raises(RuntimeError)
    def test_ratio_deadtime(self):
        testRatio = RatioData("test_ratio",
                              self.numerator_isotope,
                              self.denominator_isotope)
        testRatio.perform_deadtime_correction(2.0, 3.0)


text:
#This version of the code uses the simple blocking worker-controller
#But falls back to serial on 1 processor or if MPI isn't installed

try:
  from mpi4py import MPI as MPI
  has_mpi = True
except:
  has_mpi = False

import time
from math import sqrt, trunc
import numpy

primes = []

def check_prime(num):

  has_divisor = False

  #Check 2 specially
  if num%2 == 0:
    return False
  #Ignore all other as 2 was false
  for i in range(3, trunc(sqrt(num))+1, 2):
    if num%i == 0:
      has_divisor = True
      break
  return not has_divisor

def controller(lower, upper):

  #Set up the basic MPI stuff
  comm = MPI.COMM_WORLD
  nproc = comm.Get_size()
  rank = comm.Get_rank()

  #Setup values for array of flags
  length = upper - lower
  flags = numpy.zeros(length)
  #Offset of last dispatched value
  current_val = 0

  #Number of in-flight work packets
  inflight = 0

  #Arrays holding data per worker:
  #Value last sent to worker
  vals_in_use = numpy.zeros(nproc-1)

  #Workers stats - how many processed in how long
  processed = numpy.zeros(nproc-1)
  start_time = numpy.zeros(nproc-1)
  cum_time = numpy.zeros(nproc-1)
  end_time = numpy.zeros(nproc-1)

  #Some things need to have the correct type BEFORE the MPI calls
  info = MPI.Status()

  while True:
    #Wait for message from workers
    # On first pass we expect to get info.tag of zero signalling "ready"
    # After that we get a non-zero tag signalling "done"
    result = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=info)

    if info.tag > 0:
      # Capture stats
      end_time[info.source-1] = time.time()
      cum_time[info.source-1] = cum_time[info.source-1] + (end_time[info.source-1] - start_time[info.source-1])
      processed[info.source-1] = processed[info.source-1] + 1
      #Store result
      #Result is True (1) if prime, False (0) else
      offset = vals_in_use[info.source-1] - lower
      flags[int(offset)] = result
      inflight = inflight - 1

    if current_val < length:
      #If there is still work to do, reply with next package
      vals_in_use[info.source-1] = lower + current_val
      #print("Dispatching ", lower+current_val)

      start_time[info.source-1] = time.time()
      current_val = current_val + 1
      comm.send(vals_in_use[info.source-1], dest=info.source, tag=1)
      inflight = inflight + 1
    else:
      #No more work, shut down the worker
      comm.send(1, dest=info.source, tag=0)

    if inflight == 0:
      #Nothing is in flight, all done
      break

  #Summarize findings
  for i in range(0, nproc-1):
    print("Worker ", i, " processed ", int(processed[i-1]), " packets in ", cum_time[i-1], "s")

  print("Found ", int(numpy.sum(flags)), " primes")

def worker():

  comm = MPI.COMM_WORLD

  #Send initial message with dummy data and tag =0 for "ready"
  data = False
  comm.send(data, dest=0, tag=0)

  while True:
    #Wait for a message
    tag = 0
    info = MPI.Status()
    candidate = comm.recv(source=0, tag=MPI.ANY_TAG, status=info)
    tag = info.tag
    if(tag > 0):
      #Got number to check, check and return
      result = check_prime(candidate)
      comm.send(result, dest=0, tag=tag)
    else:
      #Shutdown - nothing more to do
      return

def main_parallel(lower, upper):
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()

  if rank == 0:
    controller(lower, upper)
  else:
    worker()

def main_serial(lower, upper):

  length = upper - lower
  flags = numpy.zeros(length)
  #Offset of last dispatched value
  start_time = time.time()
  for i in range(1, length):
    flags[i] = check_prime(lower+i)
  end_time = time.time()
  print("Found ", int(numpy.sum(flags)), " primes in", end_time-start_time, ' s')

if __name__ == "__main__":

  if has_mpi:
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    nproc = comm.Get_size()

     #Take two values describing range to check, [lower, upper]
    if(rank ==0):
      #Notify that multicore is available
      if(nproc == 1): print("Found mpi4py, can run on multiple cores")
      try:
        lower = int(input('Enter lower bnd: '))
      except:
        print("I didn't understand. I'll try 10000")
        lower = 10000
      try:
        upper = int(input('Enter upper bnd: '))
      except:
        print("I didn't understand. I'll try 20000")
        upper = 20000
    else:
      lower = 0
      upper = 0

    if nproc > 1:
      main_parallel(lower, upper)
    else:
      main_serial(lower, upper)

  else:
    try:
      lower = int(input('Enter lower bnd: '))
    except:
      print("I didn't understand. I'll try 10000")
      lower = 10000
    try:
      upper = int(input('Enter upper bnd: '))
    except:
      print("I didn't understand. I'll try 20000")
      upper = 20000

    main_serial(lower, upper)


text:
# transfer.py

import cv2
import numpy as np
from matplotlib import pyplot as plt

def readimg(name):
    img = cv2.imread(name)
    return (img)

def grayimg(img, name):
    img_gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    name = "./test/" + name
    cv2.imwrite(name, img_gray)
    return (img_gray)

def showdistribution(img):
    plt.hist(img.ravel(),256,[0,256])
    plt.savefig("hist.png")
    #plt.show()
    #hist = cv2.calcHist([img],[0],None,[256],[0,256])

def blurredimg(img, kernel_size=10):
    kernel = np.ones((kernel_size,kernel_size),np.float32)/(kernel_size*kernel_size)
    dst = cv2.filter2D(img,-1,kernel)
    cv2.imwrite('./test/blurred.png', dst)
    #cv2.imshow('dst', dst)
    return (dst)

def threshold(img, x1, x2):
    ret2, thresh = cv2.threshold(img,x1,x2,0) #127
    cv2.imwrite("./test/thrshold.png", thresh)
    return (ret2, thresh)

def kmeans(img, K = 3):
    Z = img.reshape((-1,3))
    Z = np.float32(Z)
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    ret,label,center=cv2.kmeans(Z,K,None,criteria,10,cv2.KMEANS_RANDOM_CENTERS)
    center = np.uint8(center)
    res = center[label.flatten()]
    res2 = res.reshape((img.shape))
    return (res2)

def drawContours(draw_img, thresh_img):
    contours, hierarchy = cv2.findContours(thresh_img,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
    img_contour = cv2.drawContours(draw_img, contours, -1, (0,255,0), 3)
    cv2.imwrite('./test/contour.png', img_contour)
    #cv2.imshow('contour',img_contour)
    return (img_contour)

def main():
    img = readimg('pikachu.jpg')
    hsv_img = cv2.cvtColor(img.copy(), cv2.COLOR_BGR2HSV)
    blurred_hsv_img = blurredimg(hsv_img, 15)
    gray_hsv_img = cv2.cvtColor(blurred_hsv_img, cv2.COLOR_BGR2GRAY)
    cv2.imwrite("./test/gray_hsv_img.png", gray_hsv_img)
    ret2, thresh = threshold(gray_hsv_img,140,255)
    hsv_img_contour = drawContours(img.copy(), thresh)
    cv2.imwrite('./test/contour_hsv.png', hsv_img_contour)
    """
    img_gray = grayimg(img, 'img_gray.png')
    showdistribution(img_gray)
    blurred_img = blurredimg(img, 20)
    k3_img = kmeans(blurred_img, K = 3)
    k3_gray_img = grayimg(k3_img, 'k3_img.png')
    blurred_gray_img = cv2.cvtColor(blurred_img, cv2.COLOR_BGR2GRAY)
    #ret2, thresh = threshold(blurred_gray_img,150,255) #127
    ret2, thresh = threshold(k3_gray_img,140,255)
    img3_contour = drawContours(k3_img, thresh)
    cv2.imwrite('./test/contour3.png', img3_contour)
    """

#cv2.waitKey(0)
#cv2.destroyAllWindows()

if __name__ == "__main__":
    main()


text:
import numpy as np

from arkouda.client import generic_msg, verbose
from arkouda.pdarrayclass import pdarray, create_pdarray

global verbose

__all__ = ["register_pda","attach_pda","unregister_pda"]

def register_pda(pda, user_defined_name):
    if not isinstance(user_defined_name, str):
        raise ValueError("user_defined_name must be a string")

    ret = None
    if isinstance(pda, pdarray):
        repMsg = generic_msg("register {} {}".format(pda.name, user_defined_name))
        ret = create_pdarray(repMsg)
    elif isinstance(pda, str):
        repMsg = generic_msg("register {} {}".format(pda, user_defined_name))
        ret = create_pdarray(repMsg)
    else:
        raise ValueError("pda must be pdarray or string")

    return ret

def attach_pda(name):
    if not isinstance(name, str):
        raise ValueError("user_defined_name must be a string")

    repMsg = generic_msg("attach {}".format(name))
    return create_pdarray(repMsg)

def unregister_pda(pda):
    if isinstance(pda, pdarray):
        repMsg = generic_msg("unregister {}".format(pda.name))
    elif isinstance(pda, str):
        repMsg = generic_msg("unregister {}".format(pda))
    else:
        raise ValueError("pda must be pdarray or string")

max_stars_repo_name: e-kayrakli/arkouda

max_issues_repo_name: e-kayrakli/arkouda

max_forks_repo_name: e-kayrakli/arkouda


text:
def asymptotic_variance( A, B, sigma, T=100 ):
    '''Computes asymptotic variance of the AR(1) process defined by:
    $X_t = A X_{t-1} + B \epsilon_t$
    where the $\epsilon_t$ follows a random law with covariance sigma.
    '''
    import numpy
    p = A.shape[0] # number of variables
    Q = numpy.zeros( (p,p) )
    for i in range(T):
        Q = numpy.dot( A, numpy.dot( Q, A.T)) + numpy.dot( B, numpy.dot( sigma, B.T))
    return Q

max_stars_repo_name: christophe-gouel/dolo

max_issues_repo_name: christophe-gouel/dolo

max_forks_repo_name: christophe-gouel/dolo


text:
from time import sleep
import numpy as np
import yaml
import time

class PTZ:
    """PTZ object and its basic functionality"""

    def __init__(self, onvif_device):
        print('PTZ initialized ...')
        # import params
        self.config = yaml.safe_load(open("config.yml"))
        self.pan_speed = self.config['pan_speed']
        self.tilt_speed = self.config['tilt_speed']
        self.zoom_speed = self.config['zoom_speed']
        self.threshold_speed = self.config['threshold_speed']

        # limit the speed values
        self.pan_speed = np.clip(self.pan_speed, -self.threshold_speed, self.threshold_speed)
        self.tilt_speed = np.clip(self.tilt_speed, -self.threshold_speed, self.threshold_speed)
        self.zoom_speed = np.clip(self.zoom_speed, -self.threshold_speed, self.threshold_speed)

        # Create ptz service object
        self.ptz = onvif_device.create_ptz_service()

        # Create media service object
        self.media = onvif_device.create_media_service()
        # Get target profile
        self.media_profile = self.media.GetProfiles()[0]
        # print('0 ', self.media_profile)

        # Create imaging service object
        self.imaging = onvif_device.create_imaging_service()

        # Get Status
        self.request_status = self.ptz.create_type('GetStatus')
        self.request_status.ProfileToken = self.media_profile.token

        # # Go to Home Position
        self.request_goto_home = self.ptz.create_type('GotoHomePosition')
        self.request_goto_home.ProfileToken = self.media_profile.token
        # self.ptz.GotoHomePosition(self.request_goto_home)

        # ABSOLUTE MOTION
        # create move request for ABSOLUTE motion
        self.move_request_abs = self.ptz.create_type('AbsoluteMove')
        self.move_request_abs.ProfileToken = self.media_profile.token
        if self.move_request_abs.Position is None:
            self.move_request_abs.Position = self.ptz.GetStatus({'ProfileToken': self.media_profile.token}).Position
        # if self.move_request_abs.Speed is None:
        #     self.move_request_abs.Speed = self.ptz.GetStatus({'ProfileToken': self.media_profile.token}).Position

        # # RELATIVE MOTION
        # This is not currently supported in Ulisse.
        # # # create move request for RELATIVE motion
        # self.move_request_rel = self.ptz.create_type('RelativeMove')
        # self.move_request_rel.ProfileToken = self.media_profile.token
        # if self.move_request_rel.Translation is None:
        #     print('self.move_request_rel.Position is None:')
        #     self.move_request_rel.Translation = self.ptz.GetStatus({'ProfileToken': self.media_profile.token}).Position

        # CONTINUOUS MOTION
        # create move request for continuous motion
        self.move_request_cont = self.ptz.create_type('ContinuousMove')
        self.move_request_cont.ProfileToken = self.media_profile.token
        if self.move_request_cont.Velocity is None:
            self.move_request_cont.Velocity = self.ptz.GetStatus({'ProfileToken': self.media_profile.token}).Position

        # Get PTZ configuration options for getting continuous move range
        self.request_config = self.ptz.create_type('GetConfigurationOptions')
        self.request_config.ConfigurationToken = self.media_profile.PTZConfiguration.token
        self.ptz_configuration_options = self.ptz.GetConfigurationOptions(self.request_config)

        # Get range of pan and tilt and zoom
        self.XMAX = self.ptz_configuration_options.Spaces.ContinuousPanTiltVelocitySpace[0].XRange.Max
        self.XMIN = self.ptz_configuration_options.Spaces.ContinuousPanTiltVelocitySpace[0].XRange.Min
        self.YMAX = self.ptz_configuration_options.Spaces.ContinuousPanTiltVelocitySpace[0].YRange.Max
        self.YMIN = self.ptz_configuration_options.Spaces.ContinuousPanTiltVelocitySpace[0].YRange.Min
        self.ZMIN = self.ptz_configuration_options.Spaces.ContinuousZoomVelocitySpace[0].XRange.Min
        self.ZMAX = self.ptz_configuration_options.Spaces.ContinuousZoomVelocitySpace[0].XRange.Max

        # print('max: ', self.ZMIN, self.ZMAX)

    # CONTINUOUS
    # basic motions
    def move_left_cont(self, timeout, speed):
        print('move left.. ')
        self.pan_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = -self.XMAX * self.pan_speed
        self.move_request_cont.Velocity.PanTilt.y = 0
        self.move_request_cont.Velocity.Zoom.x = 0
        self.perform_move_cont(self.move_request_cont, timeout)

    def move_right_cont(self, timeout, speed):
        print('move right.. ')
        self.pan_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = self.XMAX * self.pan_speed
        self.move_request_cont.Velocity.PanTilt.y = 0
        self.move_request_cont.Velocity.Zoom.x = 0
        # self.move_request_cont.Timeout = timeout

        self.perform_move_cont(self.move_request_cont, timeout)

    def move_up_cont(self, timeout, speed):
        print('move up.. ')
        self.tilt_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = 0
        self.move_request_cont.Velocity.PanTilt.y = self.YMAX * self.tilt_speed
        self.move_request_cont.Velocity.Zoom.x = 0
        self.perform_move_cont(self.move_request_cont, timeout)

    def move_down_cont(self, timeout, speed):
        print('move down.. ')
        self.tilt_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = 0
        self.move_request_cont.Velocity.PanTilt.y = -self.YMAX * self.tilt_speed
        self.move_request_cont.Velocity.Zoom.x = 0
        self.perform_move_cont(self.move_request_cont, timeout)

    # zoom
    def zoom_in_cont(self, timeout, speed):
        print('zoom in.. ')
        self.zoom_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = 0
        self.move_request_cont.Velocity.PanTilt.y = 0
        self.move_request_cont.Velocity.Zoom.x = self.ZMAX * self.zoom_speed
        # self.move_request_cont.Timeout = timeout

        self.perform_move_cont(self.move_request_cont, timeout)

    def zoom_out_cont(self, timeout, speed):
        print('zoom out.. ')
        self.zoom_speed = speed
        self.move_request_cont.Velocity.PanTilt.x = 0
        self.move_request_cont.Velocity.PanTilt.y = 0
        self.move_request_cont.Velocity.Zoom.x = -self.ZMAX * self.zoom_speed
        # self.move_request_cont.Timeout = timeout

        self.perform_move_cont(self.move_request_cont, timeout)

    def perform_move_cont(self, request, timeout):
        # Start continuous move
        self.ptz.ContinuousMove(request)
        # Wait a certain time
        sleep(timeout)
        # Stop continuous move
        self.ptz.Stop({'ProfileToken': request.ProfileToken})

    # ABSOLUTE MOTION
    def move_abs(self, angle_x, angle_y, angle_z):
        self.move_request_abs.Position.PanTilt.x = angle_x
        self.move_request_abs.Position.PanTilt.y = angle_y
        self.move_request_abs.Position.Zoom.x = angle_z             # zoom bug here

        # self.move_request_abs.Speed.PanTilt = [self.XMAX * self.pan_speed, self.YMAX * self.tilt_speed]
        # self.move_request_abs.Speed.Zoom = self.ZMAX * self.zoom_speed
        self.perform_move_abs(self.move_request_abs)

    def perform_move_abs(self, request):
        # Start absolute  move
        self.ptz.AbsoluteMove(request)
        # sleep(7)
        # Stop absolute  move
        # self.ptz.Stop({'ProfileToken': request.ProfileToken})

    # # # RELATIVE MOTION
    # def move_rel(self):
    #     print('move_abs .. ')
    #     self.move_request_rel.Translation.PanTilt.x = 0.1
    #     self.move_request_rel.Translation.PanTilt.y = 0.1
    #     #self.move_request_rel.Translation.Zoom.x = 0.1
    #
    #     # self.move_request_abs.Speed.PanTilt = [self.XMAX * self.pan_speed, self.YMAX * self.tilt_speed]
    #     # self.move_request_abs.Speed.Zoom = self.ZMAX * self.zoom_speed
    #     self.perform_move_rel(self.move_request_rel)
    #
    # def perform_move_rel(self, request):
    #     # Start absolute  move
    #     self.ptz.RelativeMove(request)
    #     # Stop absolute  move
    #     # self.ptz.Stop({'ProfileToken': request.ProfileToken})

    # GET POSITION
    def get_pose(self):
        start = time.time()
        pose = self.ptz.GetStatus(self.request_status).Position
        end = time.time()
        print('elapsed: ', end - start)
        pos_x = pose.PanTilt.x
        pos_y = pose.PanTilt.y
        pos_z = pose.Zoom.x
        print('x,y,z :', pos_x, pos_y, pos_z)
        return pos_x, pos_y, pos_z

    def get_move_status(self):
        start = time.time()
        pose = self.ptz.GetStatus(self.request_status)
        end = time.time()
        print('elapsed: ', end - start)
        return pose

    def get_time(self):
        ptz_utc_time = self.ptz.GetStatus(self.request_status).UtcTime
        ptz_time = {'year': ptz_utc_time.strftime("%Y"), 'month': ptz_utc_time.strftime("%m"),
                    'day': ptz_utc_time.strftime("%d"),
                    'hour': ptz_utc_time.strftime("%H"), 'minute': ptz_utc_time.strftime("%M"),
                    'second': ptz_utc_time.strftime("%S")}
        print('get time :', ptz_time)
        return ptz_time

    def move_home(self):
        print('Going to Home position ... ')
        self.ptz.GotoHomePosition(self.request_goto_home)
        # print('Home position: {}, {}, {}'.format(self.get_pose())

    # def get_dt(self):
    #     dt = self.device.devicemgmt.GetSystemDateAndTime()
    #     print('Dt: ', dt)

max_stars_repo_name: baimukashev/PTZ-Control

max_issues_repo_name: baimukashev/PTZ-Control

max_forks_repo_name: baimukashev/PTZ-Control


text:
import typing
from typing import Union
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf

def plot_residuals(residuals: Union[np.array, pd.Series]=None,
                   y_true: Union[np.array, pd.Series]=None,
                   y_pred: Union[np.array, pd.Series]=None,
                   fig: matplotlib.figure.Figure=None,
                   **kwargs) -> None:
    '''

    Parameters
    ----------
    residuals: np.array, default `None`.
        Values of residuals. If `None`, residuals are calculated internally using
        `y_true` and `y_true`.

    y_true: np.array, default `None`.
        Ground truth (correct) values. Ignored if residuals is not `None`.

    y_pred: np.array, default `None`.
        Values of predictions. Ignored if residuals is not `None`.

    fig: matplotlib.figure.Figure, default `None`.
        Pre-existing fig for the plot. Otherwise, call matplotlib.pyplot.figure()
        internally.

    kwargs
        Other keyword arguments are passed to matplotlib.pyplot.figure()

    Returns
    -------
        None

    '''

    if residuals is None and (y_true is None or y_pred is None):
        raise Exception(
            "If `residuals` aregument is None then, y_true and y_pred must be provided."
        )

    if residuals is None:
        residuals = y_pred - y_true

    if fig is None:
        fig = plt.figure(constrained_layout=True, **kwargs)

    gs  = matplotlib.gridspec.GridSpec(2, 2, figure=fig)
    ax1 = plt.subplot(gs[0, :])
    ax2 = plt.subplot(gs[1, 0])
    ax3 = plt.subplot(gs[1, 1])

    ax1.plot(residuals)
    sns.histplot(residuals, kde=True, bins=30, ax=ax2)
    plot_acf(residuals, ax=ax3, lags=60)

    ax1.set_title("Residuals")
    ax2.set_title("Distribution")
    ax3.set_title("Autocorrelation")

reason: import numpy,from statsmodels


text:
from spectral_cube import SpectralCube, VaryingResolutionSpectralCube
from radio_beam import Beam
import numpy as np
import astropy.units as u
import astropy.utils.console as console
import copy
import warnings

def ftconvolve(ImageIn, major = 1.0, minor = 1.0,
               angle = 0.0):
    NanMaskFlag = False
    nanmask = np.isnan(ImageIn)
    image = np.copy(ImageIn)
    if np.any(nanmask):
        wtimg = np.ones_like(ImageIn)
        NanMaskFlag = True
        image[nanmask] = 0.0
        wtimg[nanmask] = 0.0
        ftwtimg = np.fft.fftn(wtimg)

    ftimg = np.fft.fftn(image)

    if major == 0.0:
        sigmau = np.inf
    else:
        sigmau = 0.5/(np.pi * major)

    if minor == 0.0:
        sigmav = np.inf
    else:
        sigmav = 0.5/(np.pi * minor)

    FTPA = angle + np.pi

    a = 0.5 * (np.cos(FTPA)**2 / sigmau**2 +
               np.sin(FTPA)**2 / sigmav**2)
    c = 0.5 * (np.sin(FTPA)**2 / sigmau**2 +
               np.cos(FTPA)**2 / sigmav**2)
    b = 0.25 * np.sin(2 * FTPA) * (1.0 / sigmav**2 - 1.0 / sigmau**2)

    vv, uu = np.meshgrid(np.fft.fftfreq(ftimg.shape[0]),
                         np.fft.fftfreq(ftimg.shape[1]),
                         indexing='ij')

    FTkernel = np.exp(-a*(uu)**2 -c*(vv)**2 +2*b*(uu*vv))
    ConvolvedImage = (np.fft.ifftn(ftimg * FTkernel)).real

    if NanMaskFlag:
        ConvolvedMask = (np.fft.ifftn(ftwtimg * FTkernel)).real
        ConvolvedImage /= ConvolvedMask
        ConvolvedImage[nanmask] = np.nan
    return(ConvolvedImage)

def MakeRoundBeam(incube,
                  outfile=None,
                  overwrite=True):

    '''
    This takes a FITS file or a SpectralCube and outputs

    Parameters
    ----------
    filename : `string` or `SpectralCube`
       Input spectral cube

    Returns
    -------
    cube : `SpectralCube`

    '''
    if isinstance(incube,str):
        cube = SpectralCube.read(incube)

    if isinstance(incube, VaryingResolutionSpectralCube):
        cube = incube

    if not isinstance(cube, VaryingResolutionSpectralCube):
        warnings.warn("No information about multiple beams")
        return(None)

    beams = cube.beams
    major_axes = np.array([bm.major.to(u.deg).value for bm in beams])
    target_beamsize = np.array(major_axes.max())
    target_beam = Beam(major=target_beamsize*u.deg,
                       minor=target_beamsize*u.deg,
                       pa=0.0*u.deg)
    print "Target beam is :", target_beam

    # Let's assume square pixels
    pixsize = cube.wcs.pixel_scale_matrix[1,1]
    fwhm2sigma = np.sqrt(8*np.log(2))

    output = np.zeros(cube.shape)

    with console.ProgressBar(cube.shape[0]) as bar:

        for ii, plane in enumerate(cube.filled_data[:]):
            this_beam = beams[ii]
            conv_beam = target_beam - this_beam

            majpix = conv_beam.major.value / pixsize / fwhm2sigma
            minpix = conv_beam.minor.value / pixsize / fwhm2sigma

            output[ii,:,:] = ftconvolve(plane,
                                        major = majpix,
                                        minor = minpix,
                                        angle = conv_beam.pa.value)

            bar.update()

    hdr = copy.copy(cube.header)
    hdr['CASAMBM'] = False
    hdr['BMAJ'] = float(target_beam.major.value)
    hdr['BMIN'] = float(target_beam.major.value)
    hdr['BPA'] = 0.0
    outcube = SpectralCube(output, cube.wcs, header=hdr)
    if outfile:
        outcube.write(outfile,overwrite=overwrite)
        return None
    return(outcube)

reason: import numpy,import astropy


text:
# @Time   : 2021/12
# @Author : Yong Zheng
# @Notes  : added LabledDataSortEvalDataLoader for context-aware ranking evaluations

"""
deepcarskit.data.dataloader.general_dataloader
################################################
"""

import numpy as np
import torch

from recbole.data.dataloader.general_dataloader import FullSortEvalDataLoader
from recbole.data.interaction import Interaction, cat_interactions
from recbole.utils import InputType, ModelType
from collections import defaultdict
from logging import getLogger

class FullSortEvalDataLoader(FullSortEvalDataLoader):
    def __init__(self, config, dataset, sampler, shuffle=False, used_ids=None):
        super().__init__(config, dataset, sampler, shuffle=shuffle)

class LabledDataSortEvalDataLoader(FullSortEvalDataLoader):
    """:class:`FullSortEvalDataLoader` is a dataloader for full-sort evaluation. In order to speed up calculation,
    this dataloader would only return then user part of interactions, positive items and used items.
    It would not return negative items.

    Args:
        config (Config): The config of dataloader.
        dataset (Dataset): The dataset of dataloader.
        sampler (Sampler): The sampler of dataloader.
        shuffle (bool, optional): Whether the dataloader will be shuffle after a round. Defaults to ``False``.

        used_item = all items that users have interacted in the training and evaluation set.
        positve_item = all items that users have interacted in the evaluation set.
        history_item = all items that users have interacted in the training set.
    """

    def __init__(self, config, dataset, sampler, shuffle=False, used_ids=None):
        self.uid_field = dataset.uid_field
        self.iid_field = dataset.iid_field
        self.is_sequential = config['MODEL_TYPE'] == ModelType.SEQUENTIAL

        self.user_id = config['USER_ID_FIELD']
        self.item_id = config['ITEM_ID_FIELD']
        self.uc_id = config['USER_CONTEXT_FIELD']
        self.LABEL = config['LABEL_FIELD']

        if not self.is_sequential:
            multidict_uc_items = self._get_multidict(dataset) # uc and rated items

            '''
                uc_positve_item = all items that uc have rated and must be positive in the evaluation set.
                uc_history_item = all items that uc have rated in the training set.
            '''
            self.ucid_list = multidict_uc_items.keys()
            self.uc_num = max(self.ucid_list)+1
            self.ucid2items_num = np.zeros(self.uc_num, dtype=np.int64)
            self.ucid2positive_item = np.array([None] * self.uc_num)
            self.ucid2history_item = np.array([None] * self.uc_num)
            self.ucid_condidates={}

            # rated items (positive AND negative) for each uc in the training set
            ucid2used_item = used_ids

            for ucid in self.ucid_list:

                uc_positve_itemlist = set(multidict_uc_items[ucid])
                self.ucid2positive_item[ucid] = torch.tensor(list(uc_positve_itemlist), dtype=torch.int64)

                self.ucid2items_num[ucid] = len(uc_positve_itemlist)

                uc_history_itemlist = ucid2used_item[ucid]

                self.ucid2history_item[ucid] = torch.tensor(list(uc_history_itemlist), dtype=torch.int64)

            # get uid and context information from uc innerid
            context_fields = dataset._get_context_fields()
            uid_list = []
            dict_context = {}
            for context in context_fields:
                dict_context[context]=[]

            for ucid in self.ucid_list:
                uid = dataset._get_uid_from_usercontexts(ucid)
                uid_list.append(uid)
                tuple_context = dataset._get_context_tuple_from_usercontexts(ucid)
                for i in range(0,len(context_fields)):
                    context = context_fields[i]
                    dict_context[context].append(tuple_context[i])

            self.ucid_list = torch.tensor(list(self.ucid_list), dtype=torch.int64)
            uid_list = torch.tensor(list(uid_list), dtype=torch.int64)
            # add uc data into data for predictions
            self.uc_df = dataset.join(Interaction({self.uid_field: uid_list, self.uc_id: self.ucid_list}))
            for context in dict_context.keys():
                new_inter = dataset.join(Interaction({context: torch.tensor(list(dict_context[context]), dtype=torch.int64)}))
                self.uc_df.update(new_inter)

        self.config = config
        self.logger = getLogger()
        self.dataset = dataset
        self.sampler = sampler
        self.batch_size = self.step = None
        self.shuffle = shuffle
        self.pr = 0
        self._init_batch_size_and_step()

    def _get_multidict(self, dataset):
        matrix_uc_item = dataset._create_sparse_matrix(dataset.inter_feat, self.uc_id, self.item_id, 'coo',
                                                       self.LABEL)
        # multidict_u_uc = defaultdict(list)
        # key = userid, value = Dict (key = uc, value = decending ranked items with ratings)
        # will get a list of Dict, given a userid
        multidict_uc_items = defaultdict(list)
        multidict_uc_items_positives = defaultdict(list)

        rows, cols = matrix_uc_item.shape
        for uc_id in range(1, rows):
            # Index = 0 => [PAD]
            # uc_id == inner id for user_context
            uc_items = matrix_uc_item.getrow(uc_id)  # csr_matrix
            items = uc_items.indices  # a list of items
            rates = uc_items.data  # a list of ratings
            num_rates = len(rates)

            if num_rates == 0:
                continue

            dict_item_rating = {}

            for i in range(0, num_rates):
                key = items[i]
                value = rates[i]
                dict_item_rating[key] = value
            # sort items based on ratings
            dict_item_rating_decending = sorted(dict_item_rating.items(), key=lambda x: x[1], reverse=True)
            # add these items into dict which uses uc as key
            for items in dict_item_rating_decending:
                multidict_uc_items[uc_id].append(items[0])
        return multidict_uc_items

    @property
    def pr_end(self):
        if not self.is_sequential:
            return len(self.ucid_list)
        else:
            return len(self.dataset)

    def _next_batch_data(self):
        if not self.is_sequential:
            uc_df = self.uc_df[self.pr:self.pr + self.step]
            ucid_list = list(uc_df[self.uc_id])

            history_item = self.ucid2history_item[ucid_list]
            positive_item = self.ucid2positive_item[ucid_list]

            history_u = torch.cat([torch.full_like(hist_iid, i) for i, hist_iid in enumerate(history_item)])
            history_i = torch.cat(list(history_item))

            positive_u = torch.cat([torch.full_like(pos_iid, i) for i, pos_iid in enumerate(positive_item)])
            positive_i = torch.cat(list(positive_item))

            self.pr += self.step
            return uc_df, (history_u, history_i), positive_u, positive_i
        else:
            interaction = self.dataset[self.pr:self.pr + self.step]
            inter_num = len(interaction)
            positive_u = torch.arange(inter_num)
            positive_i = interaction[self.iid_field]

            self.pr += self.step
            return interaction, None, positive_u, positive_i

max_stars_repo_name: irecsys/DeepCARSKit

max_issues_repo_name: irecsys/DeepCARSKit

max_forks_repo_name: irecsys/DeepCARSKit


text:
import matplotlib
matplotlib.use('Agg')

import os
import sys
sys.path.insert(0, '/home/dmishkin/dev/irish-coffe/python')
if len (sys.argv) < 4:
  raise RuntimeError('Usage: python ' + sys.argv[0] + ' path_to_solver path_to_save_model mode')
import caffe

from pylab import *
import random
import numpy as np

gpu_runline = list()

solver_path = str(sys.argv[1])
init_path = str(sys.argv[2])
init_mode =  str(sys.argv[3])
margin = 0.1;
max_iter = 100;

#if len (sys.argv) >= 5:
#  neg_slope =  float(str(sys.argv[4]))
#else:
#  neg_slope=0;

#if len (sys.argv) >= 6:
#  margin =  float(str(sys.argv[5]))

mode_check=False;
if init_mode == 'Orthonormal':
  print init_mode
  mode_check=True
elif init_mode == 'IterNorm':
  print init_mode
  mode_check=True
elif init_mode == 'OrthonormalIterNorm':
  print init_mode
  mode_check=True
elif init_mode == 'OrthonormalIterNormMean':
  print init_mode
  mode_check=True
else:
  raise RuntimeError('Unknown mode. Try Orthonormal or IterNorm or  OrthonormalIterNorm')

caffe.set_mode_gpu()
solver = caffe.SGDSolver(solver_path)
if os.path.isfile(init_path):
  print "Loading"
  solver.net.copy_from(init_path)

def svd_orthonormal(shape):
        if len(shape) < 2:
            raise RuntimeError("Only shapes of length 2 or more are "
                               "supported.")

        flat_shape = (shape[0], np.prod(shape[1:]))
        a = standard_normal(flat_shape)
        a = standard_normal(flat_shape)
        a = standard_normal(flat_shape)
        a = standard_normal(flat_shape)
        u, _, v = np.linalg.svd(a, full_matrices=False)
        q = u if u.shape == flat_shape else v
        q = q.reshape(shape)
        return q

for k,v in solver.net.params.iteritems():
    print k, v[0].data.shape,v[1].data.shape
    if '_bn' in k:
      print 'BatchNorm layer'
      continue;
    if 'Orthonormal' in init_mode:
      weights=svd_orthonormal(v[0].data[:].shape)
      solver.net.params[k][0].data[:]=weights#* sqrt(2.0/(1.0+neg_slope*neg_slope));
      if 'Mean' in init_mode:
	biases=solver.net.params[k][1].data[:]
	biases[:]=0;
	solver.net.params[k][1].data[:]=biases;
    else:
      weights=solver.net.params[k][0].data[:]
      if 'Mean' in init_mode:
	biases=solver.net.params[k][1].data[:]
    if 'IterNorm' in init_mode:
      solver.net.forward()
      v = solver.net.blobs[k];
      var1=np.var(v.data[:]);
      mean1 = np.mean(v.data[:]);
      print k,'var = ', var1,'mean = ', mean1
      sys.stdout.flush()
      iter_num = 0;
      while ((abs(1.0 - var1) > margin) or (var1 < 1.0)):
	if 'Mean' in init_mode:
	  biases=solver.net.params[k][1].data[:]
	  solver.net.params[k][1].data[:] -= mean1 / len(biases[:]);
	weights=solver.net.params[k][0].data[:]
	solver.net.params[k][0].data[:] = weights / sqrt(var1);
        if (iter_num > 50) and (var1 < 1.0):
          solver.net.params[k][0].data[:]*= 2
        solver.net.forward()
        v = solver.net.blobs[k];
        var1=np.var(v.data[:]);
        mean1 = np.mean(v.data[:]);
        print k,'var = ', var1,'mean = ', mean1
        sys.stdout.flush()
        iter_num+=1;
        if iter_num > max_iter:
	  print 'Could not converge in ', iter_num, ' iterations, go to next layer'
	  break;

print "Initialization finished!"
solver.net.forward()
for k,vv in solver.net.blobs.iteritems():
    try:
      len1=len(vv.data[:])
      if len1==0:
	continue
    except:
      continue
    print k,vv.data[:].shape, ' var = ', np.var(vv.data[:]), ' mean = ', np.mean(vv.data[:])
print "Saving model..."
solver.net.save(init_path)
print "Model saved to", init_path

max_stars_repo_name: ducha-aiki/caffe

max_issues_repo_name: ducha-aiki/caffe

max_forks_repo_name: ducha-aiki/caffe


text:
# run this script if the underlying covid19 disease parameters have been modified
from ai4good.params.disease_params import covid_specific_parameters
import ai4good.utils.path_utils as pu
import pandas as pd
import numpy as np

if __name__=="__main__":
    generated_params = {}
    np.random.seed(42)
    generated_params['R0'] = np.random.normal(covid_specific_parameters["R0_medium"], 1, 1000)
    generated_params['LatentPeriod'] = np.random.normal(covid_specific_parameters["Latent_period"], 1, 1000)
    generated_params['RemovalPeriod'] = np.random.normal(covid_specific_parameters["Infectious_period"], 1, 1000)
    generated_params['HospPeriod'] = np.random.normal(covid_specific_parameters["Hosp_period"], 1, 1000)
    generated_params['DeathICUPeriod'] = np.random.normal(covid_specific_parameters["Death_period_withICU"], 1, 1000)
    generated_params['DeathNoICUPeriod'] = np.random.normal(covid_specific_parameters["Death_period"], 1, 1000)
    generated_params_df = pd.DataFrame(generated_params)
    generated_params_df[generated_params_df<0] = np.nan
    generated_params_df = generated_params_df.dropna()
    generated_params_df.to_csv(pu.params_path('generated_params.csv'))

max_stars_repo_name: AIforGoodSimulator/model-server

max_issues_repo_name: AIforGoodSimulator/model-server

max_forks_repo_name: AIforGoodSimulator/model-server


text:
import numpy as np
import matplotlib.pyplot as plt
from tools import get_params, get_filenames, sort_key

Tc = 2/np.log(1+np.sqrt(2))
lw = 3
c = ['indianred', 'orange', 'forestgreen', 'royalblue']

fig, ax = plt.subplots(figsize=(10,8))

for filename in get_filenames('test_scalar'):

    b, L, K = get_params(filename)
    results = np.loadtxt(filename, unpack=True)

    temp = results[1]
    mse = results[6]

    i = int(np.log(L)/np.log(b)-3)
    ax.semilogy(temp/Tc, mse, label="L="+str(L), color=c[i], linewidth=lw)

# sort legend
handles, labels = ax.get_legend_handles_labels()
labels, handles = zip(*sorted(zip(labels, handles), key=sort_key))
ax.legend(handles, labels)

# fig params
plt.ylabel('Mean Squared Error')
plt.xlim(temp[0]/Tc, temp[-1]/Tc)
plt.xlabel(r'$T/T_c$')
plt.title('Testing MSE of Outputs From Large and Small Lattices')
plt.grid(alpha=0.2)
plt.savefig('test_mse.pdf', format='pdf')

max_stars_repo_name: kim-jane/MCRG

max_issues_repo_name: kim-jane/MCRG

max_forks_repo_name: kim-jane/MCRG


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from copy import copy
import itertools
import json
import math
import os
import re
import sys

import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
import pandas as pd
import seaborn as sns

import click
from pylru import lrudecorator
import rasterio

import pdb

import projections.utils as utils

def replacer(source, target, replacement, replacements=None):
  return replacement.join(source.rsplit(target, replacements))

@lrudecorator(5)
def carea(bounds=None, height=None):
  with rasterio.open(utils.luh2_static('carea')) as ds:
    if bounds is None:
      return ds.read(1, masked=True)
    return ds.read(1, masked=True, window=ds.window(*bounds))

@lrudecorator(5)
def tarea(bounds=None, height=None):
  area = carea(bounds, height)
  ice_ds = rasterio.open(utils.luh2_static('icwtr'))
  if bounds is None:
    ice = ice_ds.read(1, masked=True)
  else:
    win = ice_ds.window(*bounds)
    if win[1][1] - win[0][1] > height:
      win = ((win[0][0], win[0][0] + height), win[1])
    ice = ice_ds.read(1, masked=True, window=win)
  return area * (1 - ice)

@lrudecorator(10)
def cnames_df():
  cnames = pd.read_csv(os.path.join(utils.data_root(), 'ssp-data',
                                    'country-names.csv'))
  return cnames

@lrudecorator(10)
def gdp_df():
  return pd.read_csv(utils.gdp_csv(), index_col=0).T

def gdp(year, fips=None):
  if fips:
    return gdp_df().loc[fips, year]
  else:
    return gdp_df().loc[:, year]

@lrudecorator(10)
def wjp_df():
  df = pd.read_excel(utils.wjp_xls(), sheet_name=5, index_col=0).T
  df = pd.merge(df, cnames_df().loc[:, ['fips', 'iso3c']],
                left_on='Country Code', right_on='iso3c')
  df.set_index('fips', inplace=True)
  return df

def wjp(attrs):
  return wjp_df().loc[:, attrs]

@lrudecorator(10)
def energy_c_df():
  df = pd.read_csv(utils.energy_c_csv())
  del df['Unnamed: 0']
  df[df == '--'] = np.nan
  return df

def cid_to_x(cid, x):
  if cid == 736:
    cid = 729
  df = cnames_df()
  row = df[df.un == cid]
  if not row.empty:
    return row[x].values[0]
  return str(int(cid))

def cid_to_fips(cid):
  return cid_to_x(cid, 'fips')

def cid_to_name(cid):
  return cid_to_x(cid, 'country.name.en')

def cid_to_ar5(cid):
  return cid_to_x(cid, 'ar5')

def sum_by(ccode, data):
  df = pd.DataFrame({'idx': ccode.reshape(-1),
                     'data': data.reshape(-1)}).dropna()
  agg = df.groupby(['idx'], sort=False).sum()
  return np.column_stack((agg.index.values.astype(int), agg.values))

def sum_by2(ccode, data, weights):
  dd = data * weights
  dd.mask = np.logical_or(data.mask, ccode.mask)
  return sum_by(ccode, dd)

def mean_by(idx, data):
  assert idx.shape == data.shape
  return np.bincount(idx, weights=data)

def weighted_mean_by_country(ccode, data, weights):
  dd = data * weights
  dd.mask = np.logical_or(data.mask, ccode.mask)
  save_mask = ccode.mask
  ccode.mask = ma.getmask(dd)
  ccode_idx = ccode.compressed().astype(int)
  ccode.mask = save_mask
  sums = mean_by(ccode_idx, dd.compressed())
  ncells = np.bincount(ccode_idx)
  idx = np.where(ncells > 0)
  carea = ncells[idx]
  return np.column_stack((idx[0].astype(int), carea, sums[idx] / carea))

def remap(what, table, nomatch=None):
  f = np.vectorize(lambda x: table.get(x, nomatch), otypes=[np.float32])
  shape = what.shape
  tmp = f(what.reshape(-1))
  return tmp.reshape(*shape)

def gen_mp4(oname, stack, ccode, title='', captions=None):
  from projections.mp4_utils import to_mp4

  all_max = stack[:, 2, :].max()
  all_min = stack[:, 2, :].min()
  cnorm = colors.Normalize(vmin=all_min, vmax=all_max)
  cmap = dict((k,v) for k,v in itertools.izip(stack[:, 0, 0],
                                              stack[:, 2, 0]))
  data = remap(ccode, cmap, ccode.fill_value)
  if captions is None:
    captions = tuple(itertools.repeat('', stack.shape[2]))
  elif len(captions) != stack.shape[2]:
    print("Error: not enough captions for all frames")
    return
  for idx, img, text in to_mp4(title, oname, stack.shape[2],
                               data, text='LogAbund', fps=10, cnorm=cnorm):
    cmap = dict((k,v) for k,v in itertools.izip(stack[:, 0, idx],
                                                stack[:, 2, idx]))
    img.set_array(remap(ccode, cmap, ccode.fill_value))
    text.set_text(captions[idx])

def parse_fname(fname):
  m = re.search(r'([a-zA-Z0-9.]+)-([a-zA-Z_]+)-([0-9]{4})\.tif$',
                os.path.basename(fname))
  if m:
    return int(m.group(3))
  return os.path.splitext(os.path.basename(fname))[0]

def parse_fname2(fname):
  return os.path.splitext(os.path.basename(fname))[0].rsplit('-', 2)

def printit(stacked):
  print("%4s %8s %8s %6s" % ('fips', 'start', 'end', '%'))
  for idx in range(stacked.shape[0]):
    print("%4s %8.2f %8.2f %6.2f%%" % (cid_to_fips(stacked[idx, 0, 0]),
                                       stacked[idx, 2, 0],
                                       stacked[idx, 2, -1],
                                       (100.0 *
                                        stacked[idx, 2, -1] /
                                        stacked[idx, 2, 0])))

def to_df(stacked, names):
  hs = {'fips': tuple(map(cid_to_fips, stacked[:, 0, 0])),
        'name': tuple(map(cid_to_name, stacked[:, 0, 0])),
        'ar5': tuple(map(cid_to_ar5, stacked[:, 0, 0])),
        'ratio': stacked[:, 2, -1] / stacked[:, 2, 0],
        'percent': (stacked[:, 2, -1] - stacked[:, 2, 0]) / stacked[:, 2, 0],
        'cells': stacked[:, 1, 0].astype(int)}
  assert len(names) == stacked.shape[2]
  for idx in range(stacked.shape[2]):
    hs[names[idx]] = stacked[:, 2, idx]
  df = pd.DataFrame(hs, index=stacked[:, 0, 0].astype(int))
  return df

def to_df2(stacked, names):
  hs = {'fips': tuple(map(cid_to_fips, stacked[:, 0, 0])),
        'name': tuple(map(cid_to_name, stacked[:, 0, 0])),
        'ar5': tuple(map(cid_to_ar5, stacked[:, 0, 0]))}
  assert len(names) == stacked.shape[2]
  for idx in range(stacked.shape[2]):
    hs[names[idx]] = stacked[:, 1, idx]
  df = pd.DataFrame(hs, index=stacked[:, 0, 0].astype(int))
  return df

@click.group(invoke_without_command=True)
@click.pass_context
def cli(ctx):
  if ctx.invoked_subcommand is None:
    click.echo('I was invoked without subcommand')
  else:
    click.echo('I am about to invoke %s' % ctx.invoked_subcommand)

@cli.command()
@click.argument('country-file', type=click.Path(dir_okay=False))
@click.argument('infiles', nargs=-1, type=click.Path(dir_okay=False))
@click.option('--npp', type=click.Path(dir_okay=False),
              help='Weight the abundance data with NPP per cell')
@click.option('--vsr', type=click.Path(dir_okay=False),
              help='Weight the richness data with vertebrate richness per cell')
@click.option('-b', '--band', type=click.INT, default=1,
              help='Index of band to process (default: 1)')
@click.option('--mp4', type=click.Path(dir_okay=False),
              help='Generate a video (in mp4 format) of the by country' +
              'weighted mean (default: False)')
@click.option('-l', '--log', is_flag=True, default=False,
              help='When set the data is in log scale and must be ' +
              'converted to linear scale (default: False)')
def countrify(infiles, band, country_file, npp, vsr, mp4, log):
  stack = []
  maps = []
  extent = None
  area = None
  if npp:
    npp_ds = rasterio.open(npp)
  with rasterio.open(country_file) as cc_ds:
    for arg in infiles:
      with rasterio.open(arg) as src:
        win = cc_ds.window(*src.bounds)
        if win[1][1] - win[0][1] > src.height:
          win = ((win[0][0], win[0][0] + src.height), win[1])
        ccode = cc_ds.read(1, masked=True, window=win)
        ccode = ma.masked_equal(ccode, -99)
        if extent is None:
          extent = (src.bounds.left, src.bounds.right,
                    src.bounds.bottom, src.bounds.top)
        data = src.read(band, masked=True)
        if log:
          data = ma.exp(data)
        if npp:
          npp_data = npp_ds.read(1, masked=True,
                                 window=npp_ds.window(*src.bounds))
          data *= npp_data
        if vsr:
          vsr_data = vsr_ds.read(1, masked=True,
                                 window=vsr_ds.window(*src.bounds))
          data *= vsr_data
        res = weighted_mean_by_country(ccode, data, carea(src.bounds,
                                                          src.height))
        if area is None:
          ice_ds = rasterio.open(utils.luh2_static('icwtr'))
          ice = ice_ds.read(1, window=ice_ds.window(*src.bounds))
          area = ma.MaskedArray(carea(src.bounds, src.height))
          area.mask = np.where(ice == 1, True, False)
          intercept = np.exp(3.477987) * area
          if npp:
            npp_data = npp_ds.read(1, masked=True,
                                   window=npp_ds.window(*src.bounds))
            intercept *= npp_data
          if vsr:
            vsr_data = vsr_ds.read(1, masked=True,
                                   window=vsr_ds.window(*src.bounds))
            intercept *= vsr_data

        #res = weighted_mean_by_country(ccode, data, 1)
        stack.append(res)
        maps.append(data)
        print('%40s: %8.2f / %8.2f' % (os.path.basename(arg),
                                       res[2, :].max(), res[2, :].min()))
    stacked = np.dstack(stack)
    names = tuple(map(parse_fname, infiles))
    df = to_df(stacked, names)

    ratio = maps[-1] / maps[0]
    a = ma.where(ratio > 1.05)
    b = ma.where(ratio < 0.85)
    w = ma.where((ratio >= 0.85) & (ratio <= 1.05))
    area.mask = ratio.mask

    above = ma.sum(area[a])
    below = ma.sum(area[b])
    within = ma.sum(area[w])
    total = ma.sum(area)
    unaccounted = ma.sum(area[ratio.mask != area.mask])
    print("Area: %6.4f / %6.4f / %6.4f" % (above / total, below / total, within / total))

    total1950 = ma.sum(ma.masked_invalid(maps[0] * area))
    # newbold-a intercept is 4.63955498
    pristine = ma.sum(intercept)

    npp_df = pd.DataFrame(weighted_mean_by_country(ccode, npp_data, area),
                          columns=['ID', 'Cells', 'npp_mean'])
    npp_df.index = npp_df.ID

    print("loss w.r.t. primary: %6.4f" % (total / pristine))
    print("loss w.r.t. 1950   : %6.4f" % (total / total1950))
    gdp_1950 = gdp([1950, 1970, 2010])
    gdp_1950.columns = ('gdp_1950', 'gdp_1970', 'gdp_2014')
    merged = pd.merge(df, gdp_1950, left_on='fips', right_index=True,
                      sort=False)#.sort_values(by=['gdp'])
    merged = merged.merge(npp_df, how='inner', left_index=True,
                          right_index=True)
    merged['ab_delta'] = np.log(merged[2014] / merged[1970])

    title = u'Abundance gain (loss) 1950 — 2014'
    idx = 0
    syms = ['x', 'o', '+', 'v', '^', '*']
    cols = ['r', 'g', 'blue', 'b', 'purple', 'y']

    plt.style.use('ggplot')
    fig1 = plt.figure(figsize=(6, 4))
    ax1 = plt.gca()
    for name, group in merged.groupby('ar5'):
      ax1.plot(group.pindex, group.ratio, label=name, marker=syms[idx],
               linestyle='', #ms=11,
               c=cols[idx])
      idx += 1
    ax1.plot([0, len(df)], [1.0, 1.0], color='k', linestyle='-', linewidth=2)

    ax1.set_title(title)
    ax1.set_ylabel('Mean area weighted abundance gain (%)')
    ax1.set_xlabel('Country sorted by GDP (1950)')
    ax1.xaxis.set_major_locator(plt.NullLocator())
    ax1.legend()
    fig1.savefig('ab-by-gdp.png')
    fig1.savefig('ab-by-gdp.pdf')
    plt.show()

    palette = copy(plt.cm.viridis)
    palette.set_over('w', 1.0)
    palette.set_under('r', 1.0)
    palette.set_bad('k', 1.0)

    fig2 = plt.figure(figsize=(6, 4))
    ax2 = plt.gca()
    title = u'Abundance ratio 2014 / 1950'
    ax2.set_title(title)
    ax2.axis('off')
    img = plt.imshow(maps[-1] / maps[0], cmap=palette, vmin=0.75, vmax=1.05,
                     extent=extent)
    plt.colorbar(orientation='horizontal')
    fig2.savefig('ab-1950-2010.png')
    fig2.savefig('ab-1950-2010.pdf')
    plt.show()

    if mp4:
      gen_mp4(mp4, stacked, ccode)

@cli.command()
@click.argument('country-file', type=click.Path(dir_okay=False))
@click.argument('infiles', nargs=-1, type=click.Path(dir_okay=False))
@click.option('--npp', type=click.Path(dir_okay=False),
              help='Weight the abundance data with NPP per cell')
@click.option('--vsr', type=click.Path(dir_okay=False),
              help='Weight the richness data with vertebrate richness per cell')
@click.option('-b', '--band', type=click.INT, default=1,
              help='Index of band to process (default: 1)')
@click.option('-l', '--log', is_flag=True, default=False,
              help='When set the data is in log scale and must be ' +
              'back-transformed to linear scale (default: False)')
@click.option('-o', '--out', type=click.File('w'),
              help='A file to write the merged data to')
def export(infiles, band, country_file, npp, vsr, log, out):
  stack = []
  maps = []
  area = None
  if npp:
    npp_ds = rasterio.open(npp)
  if vsr:
    vsr_ds = rasterio.open(vsr)
  types = list(set((x[0], x[1]) for x in
                   map(parse_fname2, infiles)))
  assert len(types) == 1
  scenario = types[0][0]
  metric = types[0][1]
  print('%s -- %s' % (scenario, metric))

  with rasterio.open(country_file) as cc_ds:
    for arg in infiles:
      with rasterio.open(arg) as src:
        win = cc_ds.window(*src.bounds)
        if win[1][1] - win[0][1] > src.height:
          win = ((win[0][0], win[0][0] + src.height), win[1])
        ccode = cc_ds.read(1, masked=True, window=win)
        ccode = ma.masked_equal(ccode, -99)
        data = src.read(band, masked=True)
        if log:
          data = ma.exp(data)
        if npp:
          npp_data = npp_ds.read(1, masked=True,
                                 window=npp_ds.window(*src.bounds))
          data *= npp_data
        if vsr:
          vsr_data = vsr_ds.read(1, masked=True,
                                 window=vsr_ds.window(*src.bounds))
          data *= vsr_data
        res = weighted_mean_by_country(ccode, data, carea(src.bounds))
        if area is None:
          ice_ds = rasterio.open(utils.luh2_static('icwtr'))
          ice = ice_ds.read(1, window=ice_ds.window(*src.bounds))
          area = ma.MaskedArray(carea(src.bounds) * (1 - ice))
          area.mask = np.where(ice == 1, True, False)
          if npp:
            npp_data = npp_ds.read(1, masked=True,
                                   window=npp_ds.window(*src.bounds))
            npp_res = weighted_mean_by_country(ccode, npp_data,
                                               carea(src.bounds))
          if vsr:
            vsr_data = vsr_ds.read(1, masked=True,
                                   window=vsr_ds.window(*src.bounds))
            vsr_res = weighted_mean_by_country(ccode, vsr_data,
                                               carea(src.bounds))
        stack.append(res)
        maps.append(data)
        print('%40s: %8.2f / %8.2f' % (os.path.basename(arg),
                                       res[:, 2].max(), res[:, 2].min()))
    stacked = np.dstack(stack)
    names = tuple(map(parse_fname, infiles))
    df = to_df(stacked, names)
    del df['percent']
    del df['ratio']

    df.rename(columns=dict((x, metric + '_' + str(x)) for x in
                           filter(lambda x: isinstance(x, int), df.columns)),
              inplace=True)
    if npp:
      df['npp_mean'] = npp_res[:, 2]
    if vsr:
      df['vsr_mean'] = vsr_res[:, 2]

    ##
    ## Fix a few mising FIPS codes.
    ##
    df.loc[df.name == 'South Sudan', 'fips'] = 'OD'
    df.loc[df.name == 'Curaçao', 'fips'] = 'UC'
    df.loc[df.name == 'Åland Islands', 'fips'] = 'AX'  # ISO 3166 code

    #
    # Compute the fraction of land that was primary in 1950
    #
    prim_fn = os.path.join(utils.outdir(),
                           'luh2', 'historical-primary-1950.tif')
    with rasterio.open(prim_fn) as src:
      ccode = cc_ds.read(1, masked=True, window=cc_ds.window(*src.bounds))
      ccode = ma.masked_equal(ccode, -99)
      with rasterio.open(utils.luh2_static('icwtr')) as ice_ds:
        ice = ice_ds.read(1, window=ice_ds.window(*src.bounds))
      prim = src.read(1, masked=True)
      area = ma.MaskedArray(carea(src.bounds) * (1 - ice),
                            mask=prim.mask)
    prim_by_cc = sum_by2(ccode, prim, area)
    area_by_cc = sum_by(ccode, area)
    prim_df = pd.DataFrame({'Primary': prim_by_cc[:, 1],
                            'Area': area_by_cc[:, 1]},
                             index=prim_by_cc[:, 0].astype(int)).sort_index()
    prim_df['prim_ratio'] = prim_df.Primary / prim_df.Area
    prim_df = prim_df.apply(pd.to_numeric)
    merged = df.merge(prim_df, how='left', left_index=True,
                      right_index=True)

    ##
    ## Read GDP data.
    ## Source: https://www.rug.nl/ggdc/historicaldevelopment/maddison/
    ## For this source I had to slightly massage the data (see the
    ## cleanup-maddison.py script) which extracts data since 1950 and
    ## adds a FIPS attribute to each country.
    ##
    merged2 = merged.merge(gdp_df().add_prefix('GDP_'), how='left',
                           left_on='fips', right_index=True)

    ##
    ## Read Rule of Law data (World Justice Project).
    ## Source: http://data.worldjusticeproject.org/#table
    ##
    wjp_attrs = ['WJP Rule of Law Index: Overall Score',
                 'Factor 1: Constraints on Government Powers',
                 'Factor 2: Absence of Corruption',
                 'Factor 3: Open Government ',
                 'Factor 4: Fundamental Rights',
                 'Factor 5: Order and Security',
                 'Factor 6: Regulatory Enforcement',
                 'Factor 7: Civil Justice',
                 'Factor 8: Criminal Justice']
    wjp_data = wjp(wjp_attrs)
    wjp_data[wjp_attrs] = wjp_data[wjp_attrs].apply(pd.to_numeric)
    merged3 = merged2.merge(wjp_data, how='left',
                            left_on='fips', right_index=True)

    ##
    ## Read human population density data and compute
    ## human population (per year).
    ## Source: projections
    ##
    hp_stk = []
    for fname in infiles:
      fname = replacer(fname, metric, 'hpd', 1)
      year = parse_fname(fname)
      print('hpd: ', year)
      with rasterio.open(fname) as src:
        hpd = src.read(1, masked=True)
        hp = hpd * carea(src.bounds)
        ccode = cc_ds.read(1, masked=True, window=cc_ds.window(*src.bounds))
        ccode = ma.masked_equal(ccode, -99)
        hp_stk.append(sum_by(ccode, hp))
    hp_df = to_df2(np.dstack(hp_stk), ['HP_' + str(x) for x in names])
    del hp_df['fips'], hp_df['name'], hp_df['ar5']
    merged4 = merged3.merge(hp_df,  how='left', left_index=True,
                            right_index=True)
    merged4['CID'] = merged4.index

    ##
    ## Read Energy consumption data.
    ## Source: https://www.eia.gov/beta/international/
    ##
    energy = energy_c_df()
    energy.rename(columns=dict((x, 'BTU_' + x) for x in
                               energy.columns[2:]),
                  inplace=True)
    del energy['Units']
    merged5 = merged4.merge(energy, how='left', left_on='name',
                            right_on='Country')

    ##
    ## Combine energy consumption and human population into
    ## energy consumption per capita.
    ##
    hp_cols = [col for col in merged5.columns if 'HP_' in col]
    hp_years = [int(x[-4:]) for x in hp_cols]
    btu_cols = [col for col in merged5.columns if 'BTU_' in col]
    btu_years = [int(x[-4:]) for x in btu_cols]
    years = sorted(set(hp_years).intersection(set(btu_years)))
    for year in years:
      yy = str(year)
      print('BTU / person: ' + yy)
      merged5['BTU_PC_' + yy] = (merged5['BTU_' + yy].astype(float) * 1e6 /
                                 merged5['HP_' + yy])

    ##
    ## Read economic complexity index (ECI) data and add it to DF.
    ## Source: https://atlas.media.mit.edu/en/
    ##
    cnames = cnames_df()
    eci = pd.read_csv(utils.eci_csv())
    eciw = eci.pivot(index='Country', values='ECI', columns='Year')
    eciw = eciw.add_prefix('ECI_')
    eci_plus = eciw.merge(cnames.loc[:, ['country.name.en', 'un']],
                          how='inner', left_index=True,
                          right_on='country.name.en')
    eci_plus.index = eci_plus.un.astype(int)
    eci_plus.index.name = None
    del eci_plus['un']
    del eci_plus['country.name.en']
    merged6 = merged5.merge(eci_plus, how='left',
                            left_on='CID', right_index=True)

    if out:
      merged6.to_csv(out.name, index=False, encoding='utf-8')

@cli.command()
@click.argument('infiles', nargs=-1, type=click.Path(dir_okay=False))
@click.option('--npp', type=click.Path(dir_okay=False),
              help='Weight the abundance data with NPP per cell')
@click.option('-b', '--band', type=click.INT, default=1,
              help='Index of band to process (default: 1)')
def timeline(infiles, npp, band):
  area = None
  if npp:
    npp_ds = rasterio.open(npp)
  parsed = map(lambda fname: parse_fname2(fname), infiles)
  scenarios, whats, years = zip(*parsed)
  yy = sorted(set(map(int, years)))
  assert len(set(whats)) == 1
  keys = tuple(set(scenarios))
  out = dict((key, [0.0] * len(yy)) for key in keys)
  out = [{'name': xx, 'data': [0.0] * len(yy)} for xx in keys]
  for scenario, year, arg in zip(scenarios, years, infiles):
    print(scenario, year)
    with rasterio.open(arg) as src:
      data = src.read(band, masked=True)
      if npp:
        npp_data = npp_ds.read(1, masked=True,
                               window=npp_ds.window(*src.bounds))
        data *= npp_data
      area = tarea(src.bounds, src.height)
      data *= area
      out[keys.index(scenario)]['data'][years.index(year)] = float(ma.sum(data))

  if 'historical' in keys:
    #ref = out[keys.index('historical')]['data'][0]
    ref = ma.sum(area)
    for jj, k in enumerate(out):
      for ii, v in enumerate(k['data']):
        out[jj]['data'][ii] /= ref
  print(json.dumps({'years': yy, 'data': out}))
  print('')

@cli.command()
@click.argument('country-file', type=click.Path(dir_okay=False))
@click.argument('infiles', nargs=-1, type=click.Path(dir_okay=False))
@click.option('--npp', type=click.Path(dir_okay=False),
              help='Weight the abundance data with NPP per cell')
@click.option('-b', '--band', type=click.INT, default=1,
              help='Index of band to process (default: 1)')
@click.option('-o', '--out', type=click.File('wb'))
def country_timeline(infiles, band, country_file, npp, out):
  stack = []
  maps = []
  extent = None
  extent_inset = None
  area = None
  if npp:
    npp_ds = rasterio.open(npp)
  with rasterio.open(country_file) as cc_ds:
    for arg in infiles:
      with rasterio.open(arg) as src:
        win = cc_ds.window(*src.bounds)
        if win[1][1] - win[0][1] > src.height:
          win = ((win[0][0], win[0][0] + src.height), win[1])
        ccode = cc_ds.read(1, masked=True, window=win)
        ccode = ma.masked_equal(ccode, -99)
        if extent is None:
          extent = (src.bounds.left, src.bounds.right,
                    src.bounds.bottom, src.bounds.top)
        data = src.read(band, masked=True)
        if npp:
          npp_data = npp_ds.read(1, masked=True,
                                 window=npp_ds.window(*src.bounds))
          if npp_data.shape != data.shape:
            import pdb; pdb.set_trace()
          data *= npp_data
        res = weighted_mean_by_country(ccode, data, carea(src.bounds,
                                                          src.height))
        if area is None:
          ice_ds = rasterio.open(utils.luh2_static('icwtr'))
          ice = ice_ds.read(1, window=ice_ds.window(*src.bounds))
          area = ma.MaskedArray(carea(src.bounds, src.height))
          area.mask = np.where(ice == 1, True, False)
          intercept = np.exp(4.63955498) * area
          if npp:
            npp_data = npp_ds.read(1, masked=True,
                                   window=npp_ds.window(*src.bounds))
            intercept *= npp_data

        #res = weighted_mean_by_country(ccode, data, 1)
        stack.append(res)
        maps.append(data)
        print('%40s: %8.2f / %8.2f' % (os.path.basename(arg),
                                       res[2, :].max(), res[2, :].min()))
    stacked = np.dstack(stack)
    names = tuple(map(parse_fname, infiles))
    df = to_df(stacked, names)
    if out:
      out.write(df.to_csv(index=False, encoding='utf-8').encode())
    print(df)

if __name__ == '__main__':
  cli()

max_stars_repo_name: ricardog/raster-project

max_issues_repo_name: NaturalHistoryMuseum/raster-project

max_forks_repo_name: NaturalHistoryMuseum/raster-project


text:
import os.path as pth
import json
import numpy as np

data_dir = 'data/'
inputJson = pth.join(data_dir, 'v2_vqa_info.json')

# load vqa info and find indices of train/val images
print(f'Loading vqa v2 info from: {inputJson}')
with open(inputJson, 'r') as f:
    v2_info = json.load(f)
questions = v2_info['imgs'] # yes, this is a list of questions
splits_by_ques = [ques['split'] for ques in questions]
val_ixs = [ix for ix, split in enumerate(splits_by_ques) if split == 'val']
train_ixs = [ix for ix, split in enumerate(splits_by_ques) if split == 'train']
assert len(val_ixs + train_ixs) == len(questions)

# Split the val set into val1 (validation) and val2 (test) so paired images
# are kept together.
np.random.seed(8)
# validation, val1 - 30%
# test, val2 - 70%
ques_id2ix = {questions[ix]['question_id']: ix for ix in val_ixs}
val_img_ids = set(questions[ix]['image_id'] for ix in val_ixs)
val1_percent = 0.3
val1_len = int(len(val_ixs) * val1_percent)
val1_ixs = set()
val2_ixs = set()
val1_img_ids = set()
val2_img_ids = set()
def is_assigned(ix):
    # Figure out if a question is assigned to a split.
    # This also tracks question ixs which have already been assigned by
    # another question with the same image id.
    img_id = questions[ix]['image_id']
    assigned = False
    if img_id in val1_img_ids:
        assigned = True
        val1_ixs.add(ix)
    elif img_id in val2_img_ids:
        assigned = True
        val2_ixs.add(ix)
    return assigned

#def add_neighbors(ques, ixs, img_ids):
#    for ques_id in ques['ques_with_diff_ans']:
#        alt_ix = ques_id2ix[ques_id]
#        if is_assigned(alt_ix):
#            continue
#        alt_img_id = questions[alt_ix]['image_id']
#        ixs.add(alt_ix)
#        img_ids.add(alt_img_id)

# If val1 is filled up first then it might get a better selection of
# which questions are available, potentially biasing val2's choice.
for i, ix in enumerate(val_ixs):
    if is_assigned(ix):
        continue
    if np.random.rand() < val1_percent:
        ixs = val1_ixs
        img_ids = val1_img_ids
    else:
        ixs = val2_ixs
        img_ids = val2_img_ids
    ques = questions[ix]
    img_id = ques['image_id']
    ixs.add(ix)
    img_ids.add(img_id)
    #add_neighbors(ques, ixs, img_ids)

# how did the sampling procedure turn out?
ratio = len(val1_ixs) / (len(val1_ixs) + len(val2_ixs))
print(f'val1 percent goal: {val1_percent}')
print(f'questions: val1 / (val1 + val2) = {ratio}')
ratio = len(val1_img_ids) / (len(val1_img_ids) + len(val2_img_ids))
print(f'val1 percent goal: {val1_percent}')
print(f'images: val1 / (val1 + val2) = {ratio}')
assert len(val1_ixs) + len(val2_ixs) == len(val_ixs)
assert len(val1_img_ids) + len(val2_img_ids) == len(val_img_ids)
# how many examples have a
def count_questions_with_diff_ans(ixs):
    import ipdb;ipdb.set_trace()
    N = len(ixs)
    count = 0
    for ix in ixs:
        ques = questions[ix]
        if len(ques['ques_with_diff_ans']) > 0:
            count += 1
    return f'{count} / {N} = {count / N}'
assert len(val1_img_ids.intersection(val2_img_ids)) == 0
#print(f'train diff ans questions: {count_questions_with_diff_ans(train_ixs)}')
#print(f'val1 diff ans questions: {count_questions_with_diff_ans(val1_ixs)}')
#print(f'val2 diff ans questions: {count_questions_with_diff_ans(val2_ixs)}')

# convert ixs to image_ids
train_img_ids = sorted(list(set([questions[ix]['image_id'] for ix in train_ixs])))
print(f'train: {len(train_img_ids)}')
val1_img_ids = sorted(list(val1_img_ids))
print(f'val1: {len(val1_img_ids)}')
val2_img_ids = sorted(list(val2_img_ids))
print(f'val2: {len(val2_img_ids)}')

# save splits to file
out_name = pth.join(data_dir, 'split_info_v3.json')
with open(out_name, 'w') as f:
    json.dump({
        'splits': ['train', 'val1', 'val2'],
        'split_img_ids': {
            'train': train_img_ids,
            'val1': val1_img_ids,
            'val2': val2_img_ids,
        }
    }, f)