text:
import math
import numpy as np
import time

from carla.sensor import Camera

from erdos.data_stream import DataStream
from erdos.message import Message
from erdos.op import Op
from erdos.timestamp import Timestamp
from erdos.utils import frequency, setup_logging

import agent_utils
from planner.map import CarlaMap
from planner.waypointer import Waypointer
from pid_controller.pid import PID
from utils import get_3d_world_position

class ERDOSAgentOperator(Op):
    def __init__(self,
                 name,
                 city_name,
                 goal_location,
                 goal_orientation,
                 depth_camera_name,
                 flags,
                 log_file_name=None):
        super(ERDOSAgentOperator, self).__init__(name)
        self._flags = flags
        self._logger = setup_logging(self.name, log_file_name)
        self._map = CarlaMap(city_name)
        self._goal_location = goal_location
        self._goal_orientation = goal_orientation
        self._wp_num_steer = 0.9  # Select WP - Reverse Order: 1 - closest, 0 - furthest
        self._wp_num_speed = 0.4  # Select WP - Reverse Order: 1 - closest, 0 - furthest
        self._waypointer = Waypointer(city_name)
        self._pid = PID(p=self._flags.pid_p,
                        i=self._flags.pid_i,
                        d=self._flags.pid_d)
        self._world_transform = []
        self._depth_imgs = []
        self._traffic_lights = []
        self._obstacles = []
        self._vehicle_pos = None
        self._vehicle_acc = None
        self._vehicle_speed = None
        (self._depth_intrinsic, self._depth_transform, self._depth_img_size) = self.__setup_camera_tranforms(
            name=depth_camera_name, postprocessing='SemanticSegmentation')

    def __setup_camera_tranforms(self,
                                 name,
                                 postprocessing,
                                 field_of_view=90.0,
                                 image_size=(800, 600),
                                 position=(0.3, 0, 1.3),
                                 rotation_pitch=0,
                                 rotation_roll=0,
                                 rotation_yaw=0):
        camera = Camera(
            name,
            PostProcessing=postprocessing,
            FOV=field_of_view,
            ImageSizeX=image_size[0],
            ImageSizeY=image_size[1],
            PositionX=position[0],
            PositionY=position[1],
            PositionZ=position[2],
            RotationPitch=rotation_pitch,
            RotationRoll=rotation_roll,
            RotationYaw=rotation_yaw)

        image_width = image_size[0]
        image_height = image_size[1]
        # (Intrinsic) K Matrix
        intrinsic_mat = np.identity(3)
        intrinsic_mat[0][2] = image_width / 2
        intrinsic_mat[1][2] = image_height / 2
        intrinsic_mat[0][0] = intrinsic_mat[1][1] = image_width / (2.0 * math.tan(90.0 * math.pi / 360.0))
        return (intrinsic_mat, camera.get_unreal_transform(), (image_width, image_height))

    @staticmethod
    def setup_streams(input_streams, depth_camera_name):
        def is_traffic_lights_stream(stream):
            return stream.labels.get('traffic_lights', '') == 'true'

        def is_segmented_frame_stream(stream):
            return stream.labels.get('segmented', '') == 'true'

        def is_obstacles_stream(stream):
            return stream.labels.get('obstacles', '') == 'true'

        input_streams.filter_name(depth_camera_name).add_callback(
            ERDOSAgentOperator.on_depth_camera_update)

        # XXX(ionel): We get the exact position from the simulator.
        input_streams.filter_name('world_transform').add_callback(
            ERDOSAgentOperator.on_world_transform_update)
        input_streams.filter_name('vehicle_pos').add_callback(
            ERDOSAgentOperator.on_vehicle_pos_update)
        input_streams.filter_name('acceleration').add_callback(
            ERDOSAgentOperator.on_vehicle_acceleration_update)
        input_streams.filter_name('forward_speed').add_callback(
            ERDOSAgentOperator.on_forward_speed_update)

        input_streams.filter(is_traffic_lights_stream).add_callback(
            ERDOSAgentOperator.on_traffic_lights_update)
        input_streams.filter(is_segmented_frame_stream).add_callback(
            ERDOSAgentOperator.on_segmented_frame)
        input_streams.filter(is_obstacles_stream).add_callback(
            ERDOSAgentOperator.on_obstacles_update)

        return [DataStream(name='action_stream')]

    # TODO(ionel): Set the frequency programmatically.
    @frequency(10)
    def run_step(self):
        self._logger.info("Running step")
        # Return if we haven't yet received all vehicle info data.
        if (self._vehicle_pos is None or
            self._vehicle_acc is None or
            self._vehicle_speed is None):
            return

        if (len(self._obstacles) == 0 or
            len(self._traffic_lights) == 0 or
            len(self._depth_imgs) == 0 or
            len(self._world_transform) == 0):
            return

        self._logger.info("Timestamps {} {} {} {}".format(
            self._obstacles[0].timestamp,
            self._traffic_lights[0].timestamp,
            self._depth_imgs[0].timestamp,
            self._world_transform[0].timestamp))

        world_transform = self._world_transform[0].data
        self._world_transform = self._world_transform[1:]

        depth_img = self._depth_imgs[0].data
        self._depth_imgs = self._depth_imgs[1:]

        traffic_lights = []
        for (corners, score, label) in self._traffic_lights[0].data:
            x = (corners[0] + corners[1]) / 2
            y = (corners[2] + corners[3]) / 2
            (x3d, y3d, z3d) = get_3d_world_position(
                x, y, self._depth_img_size, depth_img, self._depth_transform, world_transform)
            state = 0
            if label is not 'Green':
                state = 1
            traffic_lights.append((x3d, y3d, state))
        self._traffic_lights = self._traffic_lights[1:]

        vehicles = []
        pedestrians = []
        for (corners, score, label) in self._obstacles[0].data:
            x = (corners[0] + corners[1]) / 2
            y = (corners[2] + corners[3]) / 2
            if label == 'person':
                (x3d, y3d, z3d) = get_3d_world_position(x, y,
                                                        self._depth_img_size,
                                                        depth_img,
                                                        self._depth_transform,
                                                        world_transform)
                pedestrians.append((x3d, y3d))
            elif label == 'car' or label == 'bicycle' or label == 'motorcycle' or label == 'bus' or label == 'truck':
                (x3d, y3d, z3d) = get_3d_world_position(x, y,
                                                        self._depth_img_size,
                                                        depth_img,
                                                        self._depth_transform,
                                                        world_transform)
                vehicles.append((x3d, y3d))
        self._obstacles = self._obstacles[1:]

        start_time = time.time()

        wp_angle, wp_vector, wp_angle_speed, wp_vector_speed = agent_utils.get_waypoints(
            self._goal_location, self._goal_orientation, self._vehicle_pos,
            self._waypointer, self._wp_num_steer, self._wp_num_speed)

        runtime = time.time() - start_time
        self._logger.info('Waypointer {} runtime {}'.format(
            self.name, runtime))

        speed_factor, state = self.__stop_for_agents(
            wp_angle, wp_vector, vehicles, pedestrians, traffic_lights)

        control = self.get_control(wp_angle, wp_angle_speed, speed_factor,
                                   self._vehicle_speed * 3.6)
        output_msg = Message(control, Timestamp(coordinates=[0]))
        self.get_output_stream('action_stream').send(output_msg)

    def on_world_transform_update(self, msg):
        self._world_transform.append(msg)

    def on_depth_camera_update(self, msg):
        self._depth_imgs.append(msg)

    def on_segmented_frame(self, msg):
        self._logger.info("Received segmented frame update at {}".format(msg.timestamp))
        # TODO(ionel): Implement!

    def on_traffic_lights_update(self, msg):
        self._logger.info("Received traffic light update at {}".format(msg.timestamp))
        self._traffic_lights.append(msg)

    def on_traffic_signs_update(self, msg):
        # TODO(ionel): Implement!
        pass

    def on_obstacles_update(self, msg):
        self._logger.info("Received obstacle update at {}".format(msg.timestamp))
        self._obstacles.append(msg)

    def on_vehicle_pos_update(self, msg):
        self._logger.info("Received vehicle pos %s", msg)
        self._vehicle_pos = msg.data

    def on_vehicle_acceleration_update(self, msg):
        self._vehicle_acc = msg.data

    def on_forward_speed_update(self, msg):
        self._vehicle_speed = msg.data

    def execute(self):
        self.run_step()
        self.spin()

    def __stop_for_agents(
            self, wp_angle, wp_vector, vehicles, pedestrians, traffic_lights):
        speed_factor = 1
        speed_factor_tl = 1
        speed_factor_p = 1
        speed_factor_v = 1

        for vehicle in vehicles:
            if agent_utils.is_vehicle_on_same_lane(self._vehicle_pos, vehicle, self._map):
                new_speed_factor_v = agent_utils.stop_vehicle(
                    self._vehicle_pos, vehicle, wp_vector, speed_factor_v, self._flags)
                speed_factor_v = min(speed_factor_v, new_speed_factor_v)

        for pedestrian in pedestrians:
            if agent_utils.is_pedestrian_hitable(pedestrian, self._map):
                new_speed_factor_p = agent_utils.stop_pedestrian(
                    self._vehicle_pos,
                    pedestrian,
                    wp_vector,
                    speed_factor_p,
                    self._flags)
                speed_factor_p = min(speed_factor_p, new_speed_factor_p)

        for tl in traffic_lights:
            if (agent_utils.is_traffic_light_active(self._vehicle_pos, tl, self._map) and
                agent_utils.is_traffic_light_visible(self._vehicle_pos, tl, self._flags)):
                tl_state = tl[2]
                new_speed_factor_tl = agent_utils.stop_traffic_light(
                    self._vehicle_pos,
                    tl,
                    tl_state,
                    wp_vector,
                    wp_angle,
                    speed_factor_tl,
                    self._flags)
                speed_factor_tl = min(speed_factor_tl, new_speed_factor_tl)

        speed_factor = min(speed_factor_tl, speed_factor_p, speed_factor_v)
        state = {
            'stop_pedestrian': speed_factor_p,
            'stop_vehicle': speed_factor_v,
            'stop_traffic_lights': speed_factor_tl
        }

        return speed_factor, state

    def get_control(self, wp_angle, wp_angle_speed, speed_factor,
                    current_speed):
        current_speed = max(current_speed, 0)
        steer = self._flags.steer_gain * wp_angle
        if steer > 0:
            steer = min(steer, 1)
        else:
            steer = max(steer, -1)

        # Don't go to fast around corners
        if math.fabs(wp_angle_speed) < 0.1:
            target_speed_adjusted = self._flags.target_speed * speed_factor
        elif math.fabs(wp_angle_speed) < 0.5:
            target_speed_adjusted = 20 * speed_factor
        else:
            target_speed_adjusted = 15 * speed_factor

        self._pid.target = target_speed_adjusted
        pid_gain = self._pid(feedback=current_speed)
        throttle = min(
            max(self._flags.default_throttle - 1.3 * pid_gain, 0),
            self._flags.throttle_max)

        if pid_gain > 0.5:
            brake = min(0.35 * pid_gain * self._flags.brake_strength, 1)
        else:
            brake = 0

        return {
            'steer': steer,
            'throttle': throttle,
            'brake': brake,
            'hand_brake': False,
            'reverse': False
        }

max_stars_repo_name: yujialuo/erdos

max_issues_repo_name: yujialuo/erdos

max_forks_repo_name: yujialuo/erdos


text:
import numpy as np
from numpy.testing import assert_almost_equal
from pyscf import gto, dft
from mldftdat.dft.xc_models import NormGPFunctional, GPFunctional
import os
import joblib

def check_dm_rks(mf,i,j):
    dm0 = mf.make_rdm1()
    dm1 = dm0.copy()
    delta = np.maximum(dm0[i,j], 1) * 1e-8
    dm1[i,j] += delta / 2
    dm1[j,i] += delta / 2
    v0 = mf.get_veff(dm=dm0)
    e0 = v0.exc
    v0 -= v0.vj
    v1 = mf.get_veff(dm=dm1)
    e1 = v1.exc
    v1 -= v1.vj
    err = (e1-e0)/(delta) - v0[i,j]
    print('ERR', i,j, err)
    return err

BAS='def2-tzvp'
mol = gto.M(atom='Ne', basis=BAS, spin=0, verbose=4)

ks_ref = dft.RKS(mol)
ks_ref.xc = 'PBE0'
ks_ref.kernel()

from mldftdat.dft import numint

# Construct settings. 'xc' should compute all XC contributions EXCEPT CIDER.
# xmix should be the fraction of CIDER.
# The below settings represent PBE0 with 25% HF replaced with 25% CIDER.
settings = {
    'xc': 'GGA_C_PBE + 0.75*GGA_X_PBE',
    'xmix': 0.25
}

##### EXAMPLE 1: Simple GP Functional #####
# WARNING: Un-mapped GP functional provided as reference,
# but functional derivatives are buggy.
# Can load functionals using joblib
mlfunc = joblib.load('test_files/gpr_example.joblib')
mlfunc = GPFunctional(mlfunc)

ks = numint.setup_rks_calc(mol, mlfunc, **settings)
ks.kernel()

# This routine tests the functional derivatives.
# There are significant errors for the pure GP functional,
# but not for the spline functional.
for i in range(0,10):
    check_dm_rks(ks,i,i)

##### EXAMPLE 2: Simple Spline Functional #####
# Spline-based functionals are faster.
mlfunc = joblib.load('test_files/agpr_spline_example.joblib')
ks = numint.setup_rks_calc(mol, mlfunc, **settings)
ks.kernel()

for i in range(0,10):
    check_dm_rks(ks,i,i)

##### EXAMPLE 3: Accurate CIDER Functional #####
mlfunc = NormGPFunctional.load('functionals/CIDER_X_AHW.yaml')
ks = numint.setup_rks_calc(mol, mlfunc, **settings)
ks.kernel()

for i in range(0,10):
    check_dm_rks(ks,i,i)

ext: py

max_stars_repo_name: mir-group/CiderPress

max_issues_repo_name: mir-group/CiderPress

max_forks_repo_name: mir-group/CiderPress

reason: import numpy,from numpy


text:
import numpy as np
import json
import matplotlib.pyplot as plt
import multiprocessing
import threading
import queue
from inspect import signature
import copy
import types

from pycromanager import Acquisition, multi_d_acquisition_events, Bridge

event_list = multi_d_acquisition_events(num_time_points=10, z_start=0, z_end=10, z_step=0.5)
# first = event_list.pop(0)
# rest = event_list

#Version 1:
def hook_fn(event):
    # if np.random.randint(4) < 2:
    #     return event
    return event

# #Version 2:
# def hook_fn(event, bridge, event_queue):
#     return event

#Version 1:
def img_process_fn(image, metadata):
    image[250:350, 100:300] = np.random.randint(0, 4999)
    return image, metadata

# #Version 2:
# def img_process_fn(image, metadata, bridge, event_queue):
#     if not hasattr(img_process_fn, 'events'):
#         setattr(img_process_fn, 'events', multi_d_acquisition_events(
#             num_time_points=10, z_start=0, z_end=10, z_step=0.5))
#     events = getattr(img_process_fn, 'events')
#     if len(events) != 0:
#         event_queue.put(events.pop(0))
#     return image, metadata

# with Acquisition('/Users/henrypinkard/megllandump', 'pythonacqtest',
#                   image_process_fn=img_process_fn,
#                  post_hardware_hook_fn=hook_fn) as acq:
#
#     acq.acquire(event_list[0])

# acq.await_completion()

# #magellan example
with Acquisition(magellan_acq_index=0, post_hardware_hook_fn=hook_fn,
                  image_process_fn=img_process_fn, debug=True) as acq:
    pass
acq.await_completion()

from pycromanager import Acquisition, multi_d_acquisition_events

# with Acquisition(directory='/Users/henrypinkard/megllandump', name='tcz_acq') as acq:
#     # Generate the events for a single z-stack
#     events = multi_d_acquisition_events(
#         num_time_points=3, time_interval_s=0,
#         channel_group='channel', channels=['DAPI', 'FITC'],
#         z_start=0, z_end=6, z_step=0.4,
#         order='tcz')
#     acq.acquire(events)

with Acquisition('/Users/henrypinkard/megllandump', 'l_axis') as acq:
    #create one event for the image at each z-slice
    events = []
    for time in range(5):
        for index, z_um in enumerate(np.arange(start=0, stop=10, step=0.5)):
            evt = {
                    #'axes' is required. It is used by the image viewer and data storage to
                    #identify the acquired image
                    'axes': {'l': index, 'time': time},

                    #the 'z' field provides the z position in µm
                    'z': z_um}
            events.append(evt)

    acq.acquire(events)

max_stars_repo_name: nicost/pycro-manager

max_issues_repo_name: nicost/pycro-manager

max_forks_repo_name: nicost/pycro-manager

reason: import numpy


text:
import math
import statistics
import warnings

import numpy as np
from hmmlearn.hmm import GaussianHMM
from sklearn.model_selection import KFold
from asl_utils import combine_sequences

class ModelSelector(object):
    '''
    base class for model selection (strategy design pattern)
    '''

    def __init__(self, all_word_sequences: dict, all_word_Xlengths: dict, this_word: str,
                 n_constant=3,
                 min_n_components=2, max_n_components=10,
                 random_state=14, verbose=False):
        self.words = all_word_sequences
        self.hwords = all_word_Xlengths
        self.sequences = all_word_sequences[this_word]
        self.X, self.lengths = all_word_Xlengths[this_word]
        self.this_word = this_word
        self.n_constant = n_constant
        self.min_n_components = min_n_components
        self.max_n_components = max_n_components
        self.random_state = random_state
        self.verbose = verbose

    def select(self):
        raise NotImplementedError

    def base_model(self, num_states):
        # with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=DeprecationWarning)
        # warnings.filterwarnings("ignore", category=RuntimeWarning)
        try:
            hmm_model = GaussianHMM(n_components=num_states, covariance_type="diag", n_iter=1000,
                                    random_state=self.random_state, verbose=False).fit(self.X, self.lengths)
            if self.verbose:
                print("model created for {} with {} states".format(self.this_word, num_states))
            return hmm_model
        except:
            if self.verbose:
                print("failure on {} with {} states".format(self.this_word, num_states))
            return None

class SelectorConstant(ModelSelector):
    """ select the model with value self.n_constant

    """

    def select(self):
        """ select based on n_constant value

        :return: GaussianHMM object
        """
        best_num_components = self.n_constant
        return self.base_model(best_num_components)

class SelectorBIC(ModelSelector):
    """ select the model with the lowest Bayesian Information Criterion(BIC) score

    http://www2.imm.dtu.dk/courses/02433/doc/ch6_slides.pdf
    Bayesian information criteria: BIC = -2 * logL + p * logN
    """

    def select(self):
        """ select the best model for self.this_word based on
        BIC score for n between self.min_n_components and self.max_n_components

        :return: GaussianHMM object
        """
        #from office-hours:
        # d = features
        # n = HMM states
        # free parameters = n**2 + 2*d*n-1

        warnings.filterwarnings("ignore", category=DeprecationWarning)

        bestModel = None
        bestBIC = float("inf")
        best_num_components = self.min_n_components

        for n in range(self.min_n_components, self.max_n_components + 1):
            try:
                hmm_model = self.base_model(n)
                logL = hmm_model.score(self.X, self.lengths)
                free_params = n**2 + 2 * self.X.shape[1] * n - 1
                BIC = -2 * logL + free_params * np.log(sum(self.lengths))
                if BIC < bestBIC:
                    bestModel = hmm_model
                    bestBIC = BIC
            except:
                pass

        return bestModel

class SelectorDIC(ModelSelector):
    ''' select best model based on Discriminative Information Criterion

    Biem, Alain. "A model selection criterion for classification: Application to hmm topology optimization."
    Document Analysis and Recognition, 2003. Proceedings. Seventh International Conference on. IEEE, 2003.
    http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.58.6208&rep=rep1&type=pdf
    https://pdfs.semanticscholar.org/ed3d/7c4a5f607201f3848d4c02dd9ba17c791fc2.pdf
    DIC = log(P(X(i)) - 1/(M-1)SUM(log(P(X(all but i))
    '''

    def select(self):
        warnings.filterwarnings("ignore", category=DeprecationWarning)

        bestModel = None
        bestDIC = float("-inf")

        for n in range(self.min_n_components, self.max_n_components + 1):
            try:
                hmm_model = self.base_model(n)
                this_word_logL = hmm_model.score(self.X, self.lengths)

                sum_all_other_logL = 0
                for all_other_word in self.words:
                    if all_other_word != self.this_word:
                        wordX, wordLengths = self.hwords[all_other_word]
                        sum_all_other_logL += hmm_model.score(wordX, wordLengths)

                DIC = this_word_logL - (sum_all_other_logL / (len(self.words) - 1))
                if DIC > bestDIC:
                    bestDIC = DIC
                    bestModel = hmm_model

            except:
                pass

        return bestModel

class SelectorCV(ModelSelector):
    ''' select best model based on average log Likelihood of cross-validation folds

    '''

    def select(self):
        warnings.filterwarnings("ignore", category=DeprecationWarning)
        max_logL = float("-inf")
        n_splits = min(len(self.lengths), 3)
        best_num_components = self.min_n_components

        #print("N OF FOLDS>", n_splits)

        for n in range(self.min_n_components, self.max_n_components + 1):

            try:
                total_logL = 0

                if n_splits > 1:
                    split_method = KFold(n_splits)
                    for cv_train_idx, cv_test_idx in split_method.split(self.sequences):

                        # Create train and test data
                        X_train , X_lenghts_train = asarray(combine_sequences(cv_train_idx, self.sequences))
                        X_test  , X_lenghts_test  = asarray(combine_sequences(cv_test_idx, self.sequences))

                        #Train HMM
                        hmm_model = GaussianHMM(n_components=n, covariance_type="diag", n_iter=1000,
                                        random_state=self.random_state, verbose=False).fit(X_train, X_lenghts_train)

                        #Get LogL on the test data
                        logL_test = hmm_model.score(X_test, X_lenghts_test)
                        total_logL += logL_test

                    total_logL = total_logL / n_splits

                #pass kfolds
                else:
                    hmm_model = self.base_model(n)
                    total_logL = hmm_model.score(self.X, self.lengths)

                if total_logL > max_logL:
                        max_logL = total_logL
                        best_num_components = n

            except:
                pass

        return GaussianHMM(n_components=best_num_components, covariance_type="diag", n_iter=1000,
                                    random_state=self.random_state, verbose=False).fit(self.X, self.lengths)


text:
import numpy as np

import albumentations.augmentations.functional as af
from albumentations.core.transforms_interface import DualTransform

from allencv.data.transforms import _ImageTransformWrapper, ImageTransform

class CourtKeypointFlip(DualTransform):
    """Flip the input horizontally around the y-axis.

    Args:
        p (float): probability of applying the transform. Default: 0.5.

    Targets:
        image, mask, bboxes, keypoints

    Image types:
        uint8, float32
    """

    def apply(self, img, **params):
        if img.ndim == 3 and img.shape[2] > 1 and img.dtype == np.uint8:
            # Opencv is faster than numpy only in case of
            # non-gray scale 8bits images
            return af.hflip_cv2(img)
        else:
            return af.hflip(img)

    def apply_to_bbox(self, bbox, **params):
        return af.bbox_hflip(bbox, **params)

    def apply_to_keypoint(self, keypoint, **params):
        return af.keypoint_hflip(keypoint, **params)

    def apply_to_keypoints(self, keypoints, **params):
        keypoints = [list(keypoint) for keypoint in keypoints]
        kps = [self.apply_to_keypoint(keypoint[:4], **params) + keypoint[4:] for keypoint in keypoints]
        # print(kps)
        # sorted_x = sorted(kps, key=lambda x: x[0])
        # bottom_left = max(sorted_x[:2], key=lambda x: x[1])
        # top_left = min(sorted_x[:2], key=lambda x: x[1])
        # bottom_right = max(sorted_x[2:], key=lambda x: x[1])
        # top_right = min(sorted_x[2:], key=lambda x: x[1])
        # tmp = [bottom_left, bottom_right, top_right, top_left]
        # print(tmp)
        return kps

ImageTransform.register("keypoint_hflip")(_ImageTransformWrapper(CourtKeypointFlip))

max_stars_repo_name: sethah/deeptennis

max_issues_repo_name: sethah/deeptennis

max_forks_repo_name: sethah/deeptennis


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

from typing import Callable, Union

import numpy.testing as npt
import pytest
import tensorflow as tf

from tests.util.misc import TF_DEBUGGING_ERROR_TYPES, random_seed
from tests.util.models.gpflow.models import QuadraticMeanAndRBFKernel
from trieste.acquisition import (
    ExpectedImprovement,
    MinValueEntropySearch,
    PenalizationFunction,
    UpdatablePenalizationFunction,
)
from trieste.acquisition.function.local_penalization import (
    LocalPenalizationAcquisitionFunction,
    hard_local_penalizer,
    soft_local_penalizer,
)
from trieste.data import Dataset
from trieste.space import Box
from trieste.types import TensorType

def test_locally_penalized_expected_improvement_builder_raises_for_empty_data() -> None:
    data = Dataset(tf.zeros([0, 1]), tf.ones([0, 1]))
    space = Box([0, 0], [1, 1])
    with pytest.raises(tf.errors.InvalidArgumentError):
        LocalPenalizationAcquisitionFunction(search_space=space).prepare_acquisition_function(
            QuadraticMeanAndRBFKernel(),
            dataset=data,
        )
    with pytest.raises(tf.errors.InvalidArgumentError):
        LocalPenalizationAcquisitionFunction(search_space=space).prepare_acquisition_function(
            QuadraticMeanAndRBFKernel(),
        )

def test_locally_penalized_expected_improvement_builder_raises_for_invalid_num_samples() -> None:
    search_space = Box([0, 0], [1, 1])
    with pytest.raises(tf.errors.InvalidArgumentError):
        LocalPenalizationAcquisitionFunction(search_space, num_samples=-5)

@pytest.mark.parametrize("pending_points", [tf.constant([0.0]), tf.constant([[[0.0], [1.0]]])])
def test_locally_penalized_expected_improvement_builder_raises_for_invalid_pending_points_shape(
    pending_points: TensorType,
) -> None:
    data = Dataset(tf.zeros([3, 2], dtype=tf.float64), tf.ones([3, 2], dtype=tf.float64))
    space = Box([0, 0], [1, 1])
    builder = LocalPenalizationAcquisitionFunction(search_space=space)
    with pytest.raises(TF_DEBUGGING_ERROR_TYPES):
        builder.prepare_acquisition_function(QuadraticMeanAndRBFKernel(), data, pending_points)

@random_seed
@pytest.mark.parametrize(
    "base_builder",
    [
        ExpectedImprovement(),
        MinValueEntropySearch(Box([0, 0], [1, 1]), grid_size=10000, num_samples=10),
    ],
)
def test_locally_penalized_acquisitions_match_base_acquisition(
    base_builder: ExpectedImprovement | MinValueEntropySearch,
) -> None:
    data = Dataset(tf.zeros([3, 2], dtype=tf.float64), tf.ones([3, 2], dtype=tf.float64))
    search_space = Box([0, 0], [1, 1])
    model = QuadraticMeanAndRBFKernel()

    lp_acq_builder = LocalPenalizationAcquisitionFunction(
        search_space, base_acquisition_function_builder=base_builder
    )
    lp_acq = lp_acq_builder.prepare_acquisition_function(model, data, None)

    base_acq = base_builder.prepare_acquisition_function(model, dataset=data)

    x_range = tf.linspace(0.0, 1.0, 11)
    x_range = tf.cast(x_range, dtype=tf.float64)
    xs = tf.reshape(tf.stack(tf.meshgrid(x_range, x_range, indexing="ij"), axis=-1), (-1, 2))
    lp_acq_values = lp_acq(xs[..., None, :])
    base_acq_values = base_acq(xs[..., None, :])

    if isinstance(base_builder, ExpectedImprovement):
        npt.assert_array_equal(lp_acq_values, base_acq_values)
    else:  # check sampling-based acquisition functions are close
        npt.assert_allclose(lp_acq_values, base_acq_values, atol=0.001)

@random_seed
@pytest.mark.parametrize("penalizer", [soft_local_penalizer, hard_local_penalizer])
@pytest.mark.parametrize(
    "base_builder",
    [ExpectedImprovement(), MinValueEntropySearch(Box([0, 0], [1, 1]), grid_size=5000)],
)
def test_locally_penalized_acquisitions_combine_base_and_penalization_correctly(
    penalizer: Callable[..., Union[PenalizationFunction, UpdatablePenalizationFunction]],
    base_builder: ExpectedImprovement | MinValueEntropySearch,
) -> None:
    data = Dataset(tf.zeros([3, 2], dtype=tf.float64), tf.ones([3, 2], dtype=tf.float64))
    search_space = Box([0, 0], [1, 1])
    model = QuadraticMeanAndRBFKernel()
    pending_points = tf.zeros([2, 2], dtype=tf.float64)

    acq_builder = LocalPenalizationAcquisitionFunction(
        search_space, penalizer=penalizer, base_acquisition_function_builder=base_builder
    )
    lp_acq = acq_builder.prepare_acquisition_function(model, data, None)  # initialize
    lp_acq = acq_builder.update_acquisition_function(lp_acq, model, data, pending_points[:1], False)
    up_lp_acq = acq_builder.update_acquisition_function(lp_acq, model, data, pending_points, False)
    assert up_lp_acq == lp_acq  # in-place updates

    base_acq = base_builder.prepare_acquisition_function(model, dataset=data)

    best = acq_builder._eta
    lipshitz_constant = acq_builder._lipschitz_constant
    penalizer = penalizer(model, pending_points, lipshitz_constant, best)

    x_range = tf.linspace(0.0, 1.0, 11)
    x_range = tf.cast(x_range, dtype=tf.float64)
    xs = tf.reshape(tf.stack(tf.meshgrid(x_range, x_range, indexing="ij"), axis=-1), (-1, 2))

    lp_acq_values = lp_acq(xs[..., None, :])
    base_acq_values = base_acq(xs[..., None, :])
    penal_values = penalizer(xs[..., None, :])
    penalized_base_acq = tf.math.exp(tf.math.log(base_acq_values) + tf.math.log(penal_values))

    if isinstance(base_builder, ExpectedImprovement):
        npt.assert_array_equal(lp_acq_values, penalized_base_acq)
    else:  # check sampling-based acquisition functions are close
        npt.assert_allclose(lp_acq_values, penalized_base_acq, atol=0.001)

@pytest.mark.parametrize("penalizer", [soft_local_penalizer, hard_local_penalizer])
@pytest.mark.parametrize("at", [tf.constant([[0.0], [1.0]]), tf.constant([[[0.0], [1.0]]])])
def test_lipschitz_penalizers_raises_for_invalid_batch_size(
    at: TensorType,
    penalizer: Callable[..., PenalizationFunction],
) -> None:
    pending_points = tf.zeros([1, 2], dtype=tf.float64)
    best = tf.constant([0], dtype=tf.float64)
    lipshitz_constant = tf.constant([1], dtype=tf.float64)
    lp = penalizer(QuadraticMeanAndRBFKernel(), pending_points, lipshitz_constant, best)

    with pytest.raises(TF_DEBUGGING_ERROR_TYPES):
        lp(at)

@pytest.mark.parametrize("penalizer", [soft_local_penalizer, hard_local_penalizer])
@pytest.mark.parametrize("pending_points", [tf.constant([0.0]), tf.constant([[[0.0], [1.0]]])])
def test_lipschitz_penalizers_raises_for_invalid_pending_points_shape(
    pending_points: TensorType,
    penalizer: Callable[..., PenalizationFunction],
) -> None:
    best = tf.constant([0], dtype=tf.float64)
    lipshitz_constant = tf.constant([1], dtype=tf.float64)
    with pytest.raises(TF_DEBUGGING_ERROR_TYPES):
        soft_local_penalizer(QuadraticMeanAndRBFKernel(), pending_points, lipshitz_constant, best)

max_stars_repo_name: ChrisMorter/trieste

max_issues_repo_name: ChrisMorter/trieste

max_forks_repo_name: ChrisMorter/trieste


text:
#!/usr/bin/python3

# Copyright © 2020, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration.
# All rights reserved.
#
# The DELTA (Deep Earth Learning, Tools, and Analysis) platform is
# licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This script creates a pdf comparing the histograms of all tiff files input.

import sys

import matplotlib.pyplot as plt
import numpy as np

from matplotlib.backends.backend_pdf import PdfPages

from delta.extensions.sources.tiff import TiffImage

def plot_band(names, band):
    imgs = [TiffImage(n) for n in names]
    max_value = 2.0
    for img in imgs:
        data = np.ndarray.flatten(img.read(bands=band))
        data = data[data > 0.0]
        data[data > max_value] = max_value
        plt.hist(data, bins=200, alpha=0.5)
        plt.title('Band ' + str(band))
    pdf.savefig()
    plt.close()

assert len(sys.argv) > 1, 'No  input tiffs specified.'

with PdfPages('output.pdf') as pdf:
    a = TiffImage(sys.argv[1])
    for i in range(a.num_bands()):
        plot_band(sys.argv[1:], i)

max_stars_repo_name: bcoltin/delta

max_issues_repo_name: parkeraddison/delta

max_forks_repo_name: parkeraddison/delta


text:
"""
originally from: https://github.com/ml3705454/mapr2/blob/master/maci/environments/differential_game.py
"""

import numpy as np
from gym.spaces import Box
from rlkit.core.serializable import Serializable

class MultiDifferentialGame(Serializable):
    def __init__(self, game_name, agent_num=2, action_low=-10, action_high=10):
        Serializable.quick_init(self, locals())
        self.game = game_name
        self.agent_num = agent_num
        # self.action_num = action_num
        self.action_range = [action_low, action_high]
        self.payoff = {}

        if self.game == 'zero_sum':
            self.payoff[0] = lambda a1, a2: a1 * a2
            self.payoff[1] = lambda a1, a2: -a1 * a2
        elif self.game == 'cooperative':
            self.payoff[0] = lambda a1, a2: a1 * a2
            self.payoff[1] = lambda a1, a2: a1 * a2
        elif self.game == 'trigonometric':
            self.payoff[0] = lambda a1, a2: np.cos(a2) * a1
            self.payoff[1] = lambda a1, a2: np.sin(a1) * a2
        elif self.game == 'mataching_pennies':
            self.payoff[0] = lambda a1, a2: (a1-0.5)*(a2-0.5)
            self.payoff[1] = lambda a1, a2: (a1-0.5)*(a2-0.5)
        elif self.game == 'rotational':
            self.payoff[0] = lambda a1, a2: 0.5 * a1 * a1 + 10 * a1 * a2
            self.payoff[1] = lambda a1, a2: 0.5 * a2 * a2 - 10 * a1 * a2
        elif self.game == 'wolf':
            def V(alpha, beta, payoff):
                u = payoff[(0, 0)] - payoff[(0, 1)] - payoff[(1, 0)] + payoff[(1, 1)]
                return alpha * beta * u + alpha * (payoff[(0, 1)] - payoff[(1, 1)]) + beta * (
                            payoff[(1, 0)] - payoff[(1, 1)]) + payoff[(1, 1)]

            payoff_0 = np.array([[0, 3], [1, 2]])
            payoff_1 = np.array([[3, 2], [0, 1]])

            self.payoff[0] = lambda a1, a2: V(a1, a2, payoff_0)
            self.payoff[1] = lambda a1, a2: V(a1, a2, payoff_1)

        elif self.game == 'max2':
            assert self.agent_num == 2
            h1 = 0.8
            h2 = 1.
            s1 = 3.
            s2 = 1.
            x1 = -5.
            x2 = 5.
            y1 = -5.
            y2 = 5.
            c = 10.
            def max_f(a1, a2):
                f1 = h1 * (-(np.square(a1 - x1) / s1) - (np.square(a2 - y1) / s1))
                f2 = h2 * (-(np.square(a1 - x2) / s2) - (np.square(a2 - y2) / s2)) + c
                return max(f1, f2)
            self.payoff[0] = lambda a1, a2: max_f(a1, a2)
            self.payoff[1] = lambda a1, a2: max_f(a1, a2)
        else:
            raise NotImplementedError

    @property
    def observation_space(self):
        return Box(np.zeros(1), np.ones(1), dtype=np.float32)

    @property
    def action_space(self):
        return Box(-np.ones(1), np.ones(1), dtype=np.float32)

    def step(self, actions):
        assert len(actions) == self.agent_num
        a1 = np.mean(actions[:int(self.agent_num/2)])
        a2 = np.mean(actions[int(self.agent_num/2):])
        actions = np.array([a1,a2]) * self.action_range[1]
        # print('scaled', actions)
        reward_1 = self.payoff[0](*tuple(actions))
        reward_2 = self.payoff[1](*tuple(actions))
        reward_n = np.zeros((self.agent_num,))
        reward_n[:int(self.agent_num/2)] = reward_1
        reward_n[int(self.agent_num/2):] = reward_2

        # print(reward_n)
        state_n = np.array([[i] for i in range(self.agent_num)])/self.agent_num
        info = {}
        done_n = np.array([True] * self.agent_num)
        return state_n, reward_n, done_n, info

    def reset(self):
        return np.array([[i] for i in range(self.agent_num)])/self.agent_num

    def render(self, mode='human', close=False):
        if mode == 'human':
            print(self.__str__())

    def __str__(self):
        content = 'Game Name {}, Number of Agent {}, Action Range {}\n'.format(self.game, self.agent_num, self.action_range)
        return content

max_stars_repo_name: maxiaoba/rlk

max_issues_repo_name: maxiaoba/rlkit

max_forks_repo_name: maxiaoba/rlkit


text:
# Authors:
#     Loic Gouarin <loic.gouarin@polytechnique.edu>
#     Benjamin Graille <benjamin.graille@math.u-psud.fr>
#
# License: BSD 3 clause

"""
HDF5 module
"""
import os
import logging
from six.moves import range
import numpy as np
import h5py
import mpi4py.MPI as mpi

log = logging.getLogger(__name__) #pylint: disable=invalid-name

class H5File:
    """
    class to manage hfd5 and xdmf file.
    """
    def __init__(self, mpi_topo, filename, path='', timestep=None, init_xdmf=False):
        if timestep is not None:
            prefix = '_{}'.format(timestep)
        else:
            prefix = ''
        self.path = path
        name, ext = os.path.splitext(filename)
        self.filename = name + prefix
        self.h5filename = name + prefix + '.h5'

        self.origin = None
        self.dx = None
        self.dim = None
        self.n = None
        self.region = None
        self.global_size = None
        self.xdmf_file = None

        if mpi.COMM_WORLD.Get_rank() == 0:
            if not os.path.exists(path):
                os.mkdir(path)

            self.h5file = h5py.File(path + '/' + self.h5filename, "w")

        # All the processes wait for the creation of the output directory
        mpi.COMM_WORLD.Barrier()

        self.mpi_topo = mpi_topo
        self.scalars = {}
        self.vectors = {}
        self._init_grid = True
        self._init_xdmf = init_xdmf

    def set_grid(self, x, y=None, z=None):
        """
        create the hdf5 coordinate.

        Parameters
        ----------

        x : ndarray
            x-coordinate

        y : ndarray
            y-coordinate
            default is None

        z : ndarray
            z-coordinate
            default is None

        """
        if not self._init_grid:
            log.warning("h5 grid redefined.")

        self.origin = [x[0]]
        self.dx = [x[1] - x[0]]
        coords = [x]
        if y is not None:
            self.origin.append(y[0])
            self.dx.append(y[1] - y[0])
            coords.append(y)
        if z is not None:
            self.origin.append(z[0])
            self.dx.append(z[1] - z[0])
            coords.append(z)

        self.dim = len(self.origin)
        self.n = [0]*self.dim
        self.region = []
        self.global_size = []

        # get the region own by each processes
        # the global size
        # and the global coords
        for i in range(self.dim):
            sub = [False]*self.dim
            sub[i] = True
            comm = self.mpi_topo.cartcomm.Sub(sub)
            self.n[i] = comm.allreduce(coords[i].size, op=mpi.SUM)
            self.region.append(comm.gather(coords[i].size, root=0))
            self.global_size.append(np.sum(np.asarray(self.region[i])))
            coords[i] = np.asarray(comm.gather(coords[i], root=0))
            comm.Free()

        if mpi.COMM_WORLD.Get_rank() == 0:
            for i in range(len(self.region)):
                self.region[i].insert(0, 0)
                for j in range(1, len(self.region[i])):
                    self.region[i][j] += self.region[i][j-1]
            #self.region = np.asarray(self.region, dtype=np.int)
            #self.region = np.concatenate((np.zeros((self.dim, 1), dtype=np.int), np.cumsum(self.region, axis=1)), axis=1)
            #print(self.region)
            for i in range(self.dim):
                dset = self.h5file.create_dataset("x_{}".format(i), [self.global_size[i]], dtype=np.double)
                dset[:] = np.concatenate(coords[i])

    def _get_slice(self, rank):
        """
        get the part own by the processor rank and
        the size of the data in each direction
        """
        ind = []
        buffer_size = []
        mpi_coords = self.mpi_topo.cartcomm.Get_coords(rank)
        # print('coords', rank, mpi_coords)
        for i in range(self.dim):
            ind.append(slice(self.region[i][mpi_coords[i]], self.region[i][mpi_coords[i]+1]))
            buffer_size.append(self.region[i][mpi_coords[i]+1] - self.region[i][mpi_coords[i]])
        return ind[::-1], buffer_size[::-1]

    def _set_dset(self, dset, comm, data, index=0, with_index=False):
        """
        Merge data from multiple sub domains into a dataset.

        Parameters
        ----------

        dset : dataset
            hdf5 dataset where to store the data

        comm : comm
            mpi communicator

        data : array
            data on the sub-domain

        index : int
            use to store an index of a vector

        with_index : bool
            if we store a vector component

        """
        ind, buffer_size = self._get_slice(0)
        ind = tuple(ind)
        if with_index:
            ind = ind + (index,)
        dset[ind] = data.T

        for i in range(1, comm.Get_size()):
            ind, buffer_size = self._get_slice(i)
            #print(i, ind, buffer_size)
            ind = tuple(ind)
            if with_index:
                ind = ind + (index,)
            rcv_buffer = np.empty(buffer_size)
            comm.Recv([rcv_buffer, mpi.DOUBLE], source=i, tag=index)
            dset[ind] = rcv_buffer

    def add_scalar(self, name, f, *fargs):
        """
        store a scalar field.

        Parameters
        ----------

        name: string
            the name of the dataset entry

        f : array or function
            data to store

        fargs: tuple
            arguments of f if f is a function

        """
        if isinstance(f, np.ndarray):
            data = f
        else:
            data = f(*fargs)

        comm = self.mpi_topo.cartcomm
        if comm.Get_rank() == 0:
            dset = self.h5file.create_dataset(name, self.global_size[::-1], dtype=np.double)
            self._set_dset(dset, comm, data)
            self.scalars[name] = self.h5filename + ":/" + name
        else:
            comm.Send([np.ascontiguousarray(data.T, dtype=np.double), mpi.DOUBLE], dest=0, tag=0)

    def add_vector(self, name, f, *fargs):
        """
        store a vector field.

        Parameters
        ----------

        name: string
            the name of the dataset entry

        f : array or function
            data to store

        fargs: tuple
            arguments of f if f is a function

        """
        if isinstance(f, list):
            datas = f
        else:
            datas = f(*fargs)

        comm = self.mpi_topo.cartcomm
        if comm.Get_rank() == 0:
            dset = self.h5file.create_dataset(name, self.global_size[::-1] + [3], dtype=np.double)
            for i, data in enumerate(datas):
                self._set_dset(dset, comm, data, i, with_index=True)
            self.vectors[name] = self.h5filename + ":/" + name
        else:
            for i, data in enumerate(datas):
                comm.Send([np.ascontiguousarray(data.T), mpi.DOUBLE], dest=0, tag=i)

    def save(self):
        """
        save the hdf5 and the xdmf files.
        """
        comm = self.mpi_topo.cartcomm
        if comm.Get_rank() == 0:
            self.h5file.close()
            self.xdmf_file = open(self.path + '/' + self.filename + '.xdmf', "w")
            self.xdmf_file.write("""<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf>
 <Domain>
            """)
            if self.dim == 2:
                self.xdmf_file.write("""
                <Grid Name="Structured Grid" GridType="Uniform">
                    <Topology TopologyType="2DRectMesh" NumberOfElements="{0}"/>
                    <Geometry GeometryType="VXVY">
                """.format(' '.join(map(str, self.global_size))))
            else:
                self.xdmf_file.write("""
                <Grid Name="Structured Grid" GridType="Uniform">
                    <Topology TopologyType="3DRectMesh" NumberOfElements="{0}"/>
                    <Geometry GeometryType="VXVYVZ">
                """.format(' '.join(map(str, self.global_size))))
            for i in range(self.dim):
                self.xdmf_file.write("""
                <DataItem Format="HDF" Dimensions="{0}">
                    {1}:/x_{2}
                </DataItem>
                """.format(self.global_size[i], self.filename + '.h5', i))

            self.xdmf_file.write("</Geometry>\n")

            for k, v in self.scalars.items():
                self.xdmf_file.write("""
                <Attribute Name="{0}" AttributeType="Scalar" Center="Node">
                <DataItem Format="HDF" Dimensions="{1}">
                {2}
                </DataItem>
                </Attribute>
                """.format(k, ' '.join(map(str, self.global_size[::-1])), v))

            for k, v in self.vectors.items():
                self.xdmf_file.write("""
                <Attribute Name="{0}" AttributeType="Vector" Center="Node">
                <DataItem Format="HDF" Dimensions="{1} {2}">
                {3}
                </DataItem>
                </Attribute>
                """.format(k, ' '.join(map(str, self.global_size[::-1])), self.dim, v))

            self.xdmf_file.write("</Grid>\n</Domain>\n</Xdmf>\n")
            self.xdmf_file.close()

max_stars_repo_name: bgraille/pylbm

max_issues_repo_name: gouarin/pylbm

max_forks_repo_name: gouarin/pylbm


text:
import numpy as np
from numpy import mean, sqrt
import torch as torch

def average_correlation_coefficient(y_pred, y_true):
    """Calculate Average Correlation Coefficient

    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Returns:
        float: Average Relative Mean Squared Error

    Raises:
        ValueError : If Parameters are not both of type np.ndarray or torch.Tensor
    """
    if isinstance(y_true, np.ndarray) and isinstance(y_pred, np.ndarray):
        top = np.sum((y_true - mean(y_true, axis=0)) *
                     (y_pred - mean(y_pred, axis=0)), axis=0)
        bottom = np.sqrt(np.sum((y_true - mean(y_true, axis=0))**2, axis=0)
                         * np.sum((y_pred - mean(y_pred, axis=0))**2, axis=0))
        return np.sum(top / bottom) / len(y_true[0])
    elif isinstance(y_true, torch.Tensor) and isinstance(y_pred, torch.Tensor):
        top = torch.sum((y_true - torch.mean(y_true, dim=0))
                        * (y_pred - torch.mean(y_pred, dim=0)), dim=0)
        bottom = torch.sqrt(torch.sum((y_true - torch.mean(y_true, dim=0))**2, dim=0) *
                            torch.sum((y_pred - torch.mean(y_pred, dim=0))**2, dim=0))
        return torch.sum(top / bottom) / len(y_true[0])
    else:
        raise ValueError(
            'y_true and y_pred must be both of type numpy.ndarray or torch.Tensor')

def average_relative_error(y_pred, y_true):
    """Calculate Average Relative Error

    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Returns:
        float: Average Relative Mean Squared Error

    Raises:
        ValueError : If Parameters are not both of type np.ndarray or torch.Tensor
    """
    if isinstance(y_true, np.ndarray) and isinstance(y_pred, np.ndarray):
        return sum(sum(abs(y_true - y_pred) / y_true) /
                   len(y_true)) / len(y_true[0, :])
    elif isinstance(y_true, torch.Tensor) and isinstance(y_pred, torch.Tensor):
        return torch.sum(torch.sum(torch.abs(y_true - y_pred) / y_true, dim=0) /
                         len(y_true)) / len(y_true[0, :])
    else:
        raise ValueError(
            'y_true and y_pred must be both of type numpy.ndarray or torch.Tensor')

def average_relative_root_mean_squared_error(y_pred, y_true):
    """Calculate Average Relative Root Mean Squared Error (aRRMSE)

    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Returns:
        float : Average Relative Root Mean Squared Error

    Raises:
        ValueError : If Parameters are not both of type np.ndarray or torch.Tensor
    """
    if isinstance(y_true, np.ndarray) and isinstance(y_pred, np.ndarray):
        return sum(sqrt(sum((y_true - y_pred)**2) /
                        sum((y_true - mean(y_true, axis=0))**2))) / len(y_pred[0, :])
    elif isinstance(y_true, torch.Tensor) and isinstance(y_pred, torch.Tensor):
        return torch.sum(torch.sqrt(torch.sum((y_true - y_pred)**2, dim=0) /
                                    torch.sum(((y_true - torch.mean(y_true, dim=0))**2), dim=0))) / len(y_pred[0, :])
    else:
        raise ValueError(
            'y_true and y_pred must be both of type numpy.ndarray or torch.Tensor')

def mean_squared_error(y_pred, y_true):
    """Calculate Mean Squared Error (MSE)

    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Returns:
        float : Mean Squared Error

    Raises:
        ValueError : If Parameters are not both of type np.ndarray or torch.Tensor
     """
    if isinstance(y_true, np.ndarray) and isinstance(y_pred, np.ndarray):
        return sum((sum((y_true - y_pred)**2) /
                    len(y_true)))
    elif isinstance(y_true, torch.Tensor) and isinstance(y_pred, torch.Tensor):
        return torch.sum(torch.sum((y_true - y_pred)**2) /
                         len(y_true))
    else:
        raise ValueError(
            'y_true and y_pred must be both of type numpy.ndarray or torch.Tensor')

def average_root_mean_squared_error(y_pred, y_true):
    """Calculate Average Root Mean Squared Error (aRMSE)

    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Returns:
        float : Average Root Mean Squared Error

    Raises:
        ValueError : If Parameters are not both of type np.ndarray or torch.Tensor
    """
    if isinstance(y_true, np.ndarray) and isinstance(y_pred, np.ndarray):
        return sum(sqrt((sum((y_true - y_pred)**2) /
                         len(y_true))))/len(y_true[0, :])
    elif isinstance(y_true, torch.Tensor) and isinstance(y_pred, torch.Tensor):
        return torch.sum(torch.sqrt(torch.sum((y_true - y_pred)**2, dim=0)
                                    / len(y_true) ))/len(y_true[0])
    else:
        raise ValueError(
            'y_true and y_pred must be both of type numpy.ndarray or torch.Tensor')

def __validate_dimensions(y_pred, y_true):
    """Validates dimensions of the two input parameters
    Args:
        y_true (array-like): np.ndarray or torch.Tensor of dimension N x d with actual values
        y_pred (array-like): np.ndarray or torch.Tensor of dimension N x d with predicted values

    Raises:
        ValueError: If dimensions are not identical
    """
    if len(y_true) is not len(y_pred) and len(y_true[0]) is not len(y_pred[0]):
        raise ValueError('Dimensions of y_true and y_pred do not match.')

max_stars_repo_name: DSAAR/amorf

max_issues_repo_name: DSAAR/amorf

max_forks_repo_name: DSAAR/amorf


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

from collections import OrderedDict
from typing import Union

import os
import csv
import numpy as np
import torch

class CSVSaver:
    """
    Save the data in a dictionary format cache, and write to a CSV file finally.
    Typically, the data can be classification predictions, call `save` for single data
    or call `save_batch` to save a batch of data together, and call `finalize` to write
    the cached data into CSV file. If no meta data provided, use index from 0 to save data.
    """

    def __init__(self, output_dir: str = "./", filename: str = "predictions.csv", overwrite: bool = True):
        """
        Args:
            output_dir: output CSV file directory.
            filename: name of the saved CSV file name.
            overwrite: whether to overwriting existing CSV file content. If we are not overwriting,
                then we check if the results have been previously saved, and load them to the prediction_dict.

        """
        self.output_dir = output_dir
        self._cache_dict: OrderedDict = OrderedDict()
        assert isinstance(filename, str) and filename[-4:] == ".csv", "filename must be a string with CSV format."
        self._filepath = os.path.join(output_dir, filename)
        self.overwrite = overwrite
        self._data_index = 0

    def finalize(self) -> None:
        """
        Writes the cached dict to a csv

        """
        if not self.overwrite and os.path.exists(self._filepath):
            with open(self._filepath, "r") as f:
                reader = csv.reader(f)
                for row in reader:
                    self._cache_dict[row[0]] = np.array(row[1:]).astype(np.float32)

        if not os.path.exists(self.output_dir):
            os.makedirs(self.output_dir)
        with open(self._filepath, "w") as f:
            for k, v in self._cache_dict.items():
                f.write(k)
                for result in v.flatten():
                    f.write("," + str(result))
                f.write("\n")

    def save(self, data: np.ndarray, meta_data=None):
        """Save data into the cache dictionary. The metadata should have the following key:
            - ``'filename_or_obj'`` -- save the data corresponding to file name or object.
        If meta_data is None, use the default index from 0 to save data instead.

        Args:
            data (Tensor or ndarray): target data content that save into cache.
            meta_data (dict): the meta data information corresponding to the data.

        """
        save_key = meta_data["filename_or_obj"] if meta_data else str(self._data_index)
        self._data_index += 1
        if torch.is_tensor(data):
            data = data.detach().cpu().numpy()
        self._cache_dict[save_key] = data.astype(np.float32)

    def save_batch(self, batch_data: Union[torch.Tensor, np.ndarray], meta_data=None):
        """Save a batch of data into the cache dictionary.

        Args:
            batch_data (Tensor or ndarray): target batch data content that save into cache.
            meta_data (dict): every key-value in the meta_data is corresponding to 1 batch of data.

        """
        for i, data in enumerate(batch_data):  # save a batch of files
            self.save(data, {k: meta_data[k][i] for k in meta_data} if meta_data else None)


text:
# -*- coding: utf-8 -*-
"""
Created on Mon May 26 15:51:17 2014

@author:
    Georg Ramer
    georg@ramer.at

"""

from __future__ import division
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from COS import cdos, testdata, COS_plot

# create testdata
data, wn =  testdata(200,400)

plt.figure("spectra")
plt.plot(wn, data)
plt.xlim((wn[0],wn[-1]))
plt.figure()
plt.contourf( range(400),wn,data)

#calculate 2D-CDS
result = cdos(data =data, epsilon = 10**-10)

sync = result.sync
async = result.async

#display 2D-CDS results
axsync = plt.figure("sync").add_subplot(111)
COS_plot(sync,result.means, ax=axsync, wn = wn, lines = 10)

axasync = plt.figure("async").add_subplot(111)
COS_plot(async,result.means, ax=axasync, wn = wn, lines=10)

plt.show()


text:
#
# Copyright 2018 Antoine Sanner
#
# ### MIT license
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#

import pytest

import numpy as np
from NuMPI.Tools import  Reduction
from runtests.mpi import MPITest
from tests.MPI_minimization_problems import MPI_Objective_Interface
import scipy.optimize
import  time
import tests.minimization_problems as mp

from NuMPI.Optimization.MPI_LBFGS_Matrix_H import steepest_descent_wolfe2, LBFGS

def timer(fun, *args, **kwargs):
    start = time.perf_counter()
    res = fun(*args, **kwargs)
    delay = time.perf_counter() - start
    return res, delay

def test_linesearch():
    pass

def test_MPI_Parallel_Interface(comm): # TODO: Duplicate with test_LBFGS_interface ?
    """ Tests the parallel wrapping of objectives
    Test if parallel Version gives the same as the serial Version
    :param comm:
    :return:
    """

    def printMPI(msg):
        for i in range(comm.Get_size()):
            comm.barrier()
            if comm.Get_rank() == i:
                print("Proc {}: {}".format(i, msg))
    n = 10

    par = MPI_Objective_Interface(mp.Extended_Rosenbrock,nb_domain_grid_pts=n,comm=comm)

    printMPI(par.counts)

    ref = mp.Extended_Rosenbrock

    np.testing.assert_array_equal(mp.Extended_Rosenbrock.startpoint(n)[par.subdomain_slices], par.startpoint())
    np.testing.assert_almost_equal(mp.Extended_Rosenbrock.f(mp.Extended_Rosenbrock.startpoint(n)),par.f(par.startpoint()),err_msg="Different Function Value at startpoint")
    np.testing.assert_allclose(mp.Extended_Rosenbrock.grad(mp.Extended_Rosenbrock.startpoint(n))[par.subdomain_slices],
                                   par.grad(par.startpoint()), err_msg="Different Gradient Value at startpoint")

#@pytest.mark.parametrize("PObjective",[mp.Extended_Rosenbrock]) # Objective should support parallelization !
#@pytest.mark.parametrize("n",[10,20])

def test_analytical_min(comm):
    """
    Compares the result with the analyticaly known posistion of the minimum
    :return:
    """

    def printMPI(msg):
        for i in range(comm.Get_size()):
            comm.barrier()
            if comm.Get_rank() == i:
                print("Proc {}: {}".format(i, msg))

    n = 50

    Objective = mp.Extended_Rosenbrock

    PObjective=MPI_Objective_Interface(Objective, nb_domain_grid_pts=n, comm=comm)

    x0 = PObjective.startpoint()

    res = LBFGS(PObjective.f_grad, x0, jac=True, maxcor=5, maxiter=100,
                gtol=1e-12, ftol=0, pnp= Reduction(comm))
    #                        ^ only terminates if gradient condition is satisfied
    assert res.success
    assert res.message=="CONVERGENCE: NORM_OF_GRADIENT_<=_GTOL"
    np.testing.assert_allclose(res.x,PObjective.xmin(),atol=1e-16,rtol = 1e-7)

    assert np.abs(res.fun-Objective.minVal(n))< 1e-7

def test_ftol(comm):
    pass

def test_gtol(comm):
    pass

def test_g2tol(comm):
    pass

def test_alltol(comm):
    pass

@pytest.mark.skip(reason="just plotting")
def test_time_complexity(comm):

    maxcor = 5
    Objective = mp.Extended_Rosenbrock
    n= np.array([10,100,1000,1e4,1e5,1e6],dtype = int)
    t = np.zeros(len(n),dtype=float)
    res = [None] * len(n)
    pnp = Reduction(comm)
    for i in range(len(n)):
        PObjective = MPI_Objective_Interface(Objective, nb_domain_grid_pts=n[i], comm=comm)
        x0 = PObjective.startpoint()

        res[i],t[i] = timer(LBFGS,PObjective.f_grad, x0, jac=True, maxcor=maxcor, maxiter=100000, gtol=(1e-5), pnp=pnp)

        assert res[i].success

    if True:
        import matplotlib.pyplot as plt

        fig,ax = plt.subplots()
        ax.plot(n,t/n,'+-', label = "time / DOF")
        ax.plot(n, [t[i] / n[i] / res[i].nit for i in range(len(n))], '+-',label = "time per DOF per iteration")
        ax2 = plt.twinx(ax)
        ax2.plot(n,[res[i].nit for i in range(len(n))],'o',label = "nit")
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_ylabel("t/n (s)")
        ax.set_xlabel("DOF")
        ax.legend()
        ax2.legend()
        #ax.plot(n,n,c='gray')
        plt.show(block = True)

reason: import numpy,import scipy


text:
import numpy as np

np.linspace(0, 1, 5) # array([0., 0.25, 0.5, 0.75, 1.0])
np.arange(0, 1, 0.3) # array([0.0, 0.3, 0.6, 0.9])

max_stars_repo_name: arifmudi/Applying-Math-with-Python

max_issues_repo_name: arifmudi/Applying-Math-with-Python

max_forks_repo_name: arifmudi/Applying-Math-with-Python


text:
""" Class for IGM Absorption sightline
 Uses AbsSightline from linetools
"""
from __future__ import print_function, absolute_import, division, unicode_literals

import pdb
import warnings

from astropy import units as u
from astropy.coordinates import SkyCoord

from linetools.isgm.abssightline import AbsSightline
from linetools.isgm.abssystem import add_comps_from_dict
from linetools.isgm.utils import build_systems_from_components
from linetools import utils as ltu

class IGMSightline(AbsSightline):
    """Class for IGM Absorption Sightline
    """

    @classmethod
    def from_json(cls, jsonfile, **kwargs):
        """ Instantiate from a JSON file

        Parameters
        ----------
        jsonfile : str
          Filename
          See from_dict for required keys
        kwargs : passed to from_dict

        Returns
        -------

        """
        jdict = ltu.loadjson(jsonfile)
        slf = cls.from_dict(jdict, **kwargs)
        # Return
        return slf

    @classmethod
    def from_igmguesses(cls, igmgfile, name=None, radec=None, zem=None,
                        buildsys = False, **kwargs):
        """ Instantiate from a JSON file from IGMGuesses

        Parameters
        ----------
        igmgfile : str
          Filename
        name : str, optional
          Name of the IGMSightline
        radec : RA/DEC input, optional
          See ltu.radec_to_coord for options on format
          If given, it will overwrite the RA/DEC in IGMGuesses file (if any)
        zem : float, optional
          Emission redshift of sightline
          If given, it will overwrite the zem in IGMGuesses file (if any)
        buildsys : bool, optional
          If True, instantiate AbsSystem objects

        Returns
        -------

        """
        from linetools.isgm import abscomponent
        # Read
        jdict = ltu.loadjson(igmgfile)  # cmps, specfile, meta
        # Add in additional keys
        # Coords
        if radec is not None:
            coord = ltu.radec_to_coord(radec)
        else:
            coord = SkyCoord(jdict['meta']['RA'], jdict['meta']['DEC'], unit='deg')
        jdict['RA'] = coord.icrs.ra.deg
        jdict['DEC'] = coord.icrs.dec.deg
        # zem
        if zem is not None:
            jdict['zem'] = zem
        else:
            jdict['zem'] = jdict['meta']['zem']
            if jdict['zem'] == 0.:  # conforming IGMGuesses current rule zem = 0. means zem not set
                jdict['zem'] = None
        # Name
        if name is None:
            if jdict['zem'] is None:
                zem_name = 'z-unknown'
            else:
                zem_name = 'z{:0.3f}'.format(jdict['zem'])
            name = 'J{:s}{:s}_{:s}'.format(
                coord.icrs.ra.to_string(unit=u.hour, sep='', pad=True)[0:4],
                coord.icrs.dec.to_string(sep='', pad=True, alwayssign=True)[0:5], zem_name)
        jdict['name'] = name
        # Components
        jdict['components'] = jdict.pop('cmps')

        kwargs['use_coord'] = True
        slf = cls.from_dict(jdict, **kwargs)

        # Separate IGMGuesses attributes from AbsComponent
        acomp_keys = list(abscomponent.init_attrib.keys())
        for comp in slf._components:
            comp.igmg_attrib = {}
            akeys = list(comp.attrib.keys())
            for key in akeys:
                if key in acomp_keys:
                    pass
                else:
                    comp.igmg_attrib[key] = comp.attrib.pop(key)
            # Add other IGMGuesses specific attributes
            comp.igmg_attrib['top_level'] = {}
            for key in ['Nfit', 'bfit', 'zfit', 'mask_abslines', 'wrest', 'vlim']:  # I added vlim because these aren't always consistent. Must be something bad in igmguesses
                if key in jdict['components'][comp.name].keys():
                    comp.igmg_attrib['top_level'][key] = jdict['components'][comp.name][key]
        # Slurp a few other IGMGuesses things
        slf.igmg_dict = {}
        for key in ['spec_file', 'meta', 'fwhm']:
            slf.igmg_dict[key] = jdict[key]
        # Return

        if buildsys:
            slf.make_igmsystems()
        return slf

    @classmethod
    def from_dict(cls, idict, **kwargs):
        """ Instantiate from a dict

        Parameters
        ----------
        idict : dict
          Required keys are:
           'RA' -- float (deg)
           'DEC' -- float(deg)
           'zem' -- float
           'name' -- str
           'components' -- list
         Other keys are added as attributes to the IgmSightline object

        Returns
        -------

        """
        if 'linelist' not in kwargs.keys():
            from linetools.lists.linelist import LineList
            ism = LineList('ISM')
            kwargs['linelist'] = ism
        from pyigm.abssys.utils import class_by_type
        ism = LineList('ISM')
        kwargs['linelist'] = ism
        # Load ISM to speed things up
        if 'meta' in idict.keys():
            meta = idict['meta']
        else:
            meta = idict
        # Components -- backwards compatability
        if 'cmps' in idict.keys():
            idict['components'] = idict['cmps'].copy()
        # Instantiate
        try:
            slf = cls(SkyCoord(ra=meta['RA'], dec=meta['DEC'], unit='deg'),
                  zem=meta['zem'], name=meta['JNAME'], **kwargs)
        except KeyError:
            slf = cls(SkyCoord(ra=meta['RA'], dec=meta['DEC'], unit='deg'),
                      zem=meta['zem'], name=meta['name'], **kwargs)
        # Other
        for key in idict.keys():
            if key in ['RA', 'DEC', 'zem', 'name', 'components', 'meta', 'cmps']:
                continue
            else:
                setattr(slf, key, idict[key])
        add_comps_from_dict(slf, idict, **kwargs)

        # Systems
        if 'systems' in idict.keys():
            for key in idict['systems'].keys():
                asys = class_by_type(idict['systems'][key]['abs_type']).from_dict(idict['systems'][key], **kwargs)
                slf._abssystems.append(asys)
        # Return
        return slf

    def __init__(self, radec, zem=None, **kwargs):
        AbsSightline.__init__(self, radec, sl_type='IGM', **kwargs)
        if zem is None:
            self.zem = None
            warnings.warn("You really should set zem to create IGMSightline")
        else:
            self.zem = zem

    def make_igmsystems(self, igmsystem=None, **kwargs):
        """ Use the component list to generate a list of IGMSystems

        Returns
        -------
        igm_systems : list
          list of IGMSystem objects

        """
        if igmsystem is None:
            from pyigm.abssys.igmsys import IGMSystem
            igmsystem = IGMSystem
        # Main call
        igm_sys = build_systems_from_components(self._components, systype=igmsystem, **kwargs)
        # Return
        self._abssystems = igm_sys
        return igm_sys

    def write_to_json(self, outfile):
        # Generate the dict
        igms_dict = self.to_dict()
        # Jsonify
        clean_dict = ltu.jsonify(igms_dict)
        # Write
        ltu.savejson(outfile, clean_dict, overwrite=True)

    def write_to_igmguesses(self, outfile, fwhm=3., specfilename=None, creator=None,
                            instrument='unknown',
                            altname='unknown', date=None, overwrite=False):
        import json
        """
        Writes an IGMGuesses formatted JSON file

        Parameters
        ----------
        outfile : str
            Name of the IGMGuesses JSON file to write to.
        fwhm : int
            FWHM for IGMguesses
        specfilename : str
            Name of the spectrum file these guesses are associated to.
        altname : str
            Alternative name for the sightline. e.g. 3C273
        creator : str
            Name of the person who is creating the file. Important for tracking.
        instrument : str
            String indicating the instrument and its configuration associated to the specfilename.
            e.g. HST/COS/G130M+G160M/LP2
        overwrite : bool
            Wthether to overwrite output

        Returns
        -------

        """
        import datetime
        # Slurp IGMGuesses component attributes
        from pyigm.guis.igmguesses import comp_init_attrib
        from linetools.isgm.abscomponent import AbsComponent
        tmp = AbsComponent((10.0*u.deg, 45*u.deg), (14,2), 1.0, [-300,300]*u.km/u.s)
        comp_init_attrib(tmp)
        igmg_akeys = list(tmp.attrib.keys())
        # components
        comp_list = self._components
        coord_ref = comp_list[0].coord
        if date is None:
            date = str(datetime.date.today().strftime('%Y-%b-%d'))

        # spec_file, meta
        if hasattr(self, 'igmg_dict'):
            spec_file = self.igmg_dict['spec_file']
            meta = self.igmg_dict['meta']
            # Updates
            meta['Date'] = date
            if creator is not None:
                meta['Creator'] = creator
            #
            fwhm = self.igmg_dict['fwhm']
        else:
            spec_file = specfilename
            # coordinates and meta
            RA = coord_ref.ra.to('deg').value
            DEC = coord_ref.dec.to('deg').value
            jname = ltu.name_from_coord(coord_ref, precision=(2, 1))
            if self.zem is None:
                zem = 0.  # IGMGuesses rules
            else:
                zem = self.zem
            if creator is None:
                creator = 'unknown'
            meta = {'RA': RA, 'DEC': DEC, 'ALTNAME': altname,
                    'zem': zem, 'Creator': creator,
                    'Instrument': instrument, 'Date': date, 'JNAME': jname}

        # Create dict of the components
        out_dict = dict(cmps={},
                        spec_file=spec_file,
                        fwhm=fwhm, bad_pixels=[],
                        meta=meta)

        for comp in comp_list:
            key = comp.name
            out_dict['cmps'][key] = comp.to_dict()
            # import pdb; pdb.set_trace()
            # check coordinate
            if comp.coord != coord_ref:
                raise ValueError("All AbsComponent objects must have the same coordinates!")
            out_dict['cmps'][key]['zcomp'] = comp.zcomp
            # IGMGuesses specific component attr
            for igm_key in ['zfit','Nfit','bfit', 'wrest', 'mask_abslines', 'vlim']:
                out_dict['cmps'][key][igm_key] = comp.igmg_attrib['top_level'][igm_key]
            # IGMGuesses attribute dict
            out_dict['cmps'][key]['attrib'] = {}
            for igm_key in igmg_akeys:
                try:
                    out_dict['cmps'][key]['attrib'][igm_key] = comp.igmg_attrib[igm_key]
                except KeyError:
                    out_dict['cmps'][key]['attrib'][igm_key] = comp.attrib[igm_key]
            #out_dict['cmps'][key]['vlim'] = list(comp.vlim.value)
            out_dict['cmps'][key]['reliability'] = str(comp.reliability)
            out_dict['cmps'][key]['comment'] = str(comp.comment)
            # Compatability on sig_logN
            out_dict['cmps'][key]['attrib']['sig_logN'] = comp.attrib['sig_logN'][0]

        # JSONify
        gd_dict = ltu.jsonify(out_dict)

        # Write file
        ltu.savejson(outfile, gd_dict, overwrite=overwrite, sort_keys=True, indent=4, separators=(',', ': '))
        print('Wrote: {:s}'.format(outfile))

    def __repr__(self):
        if self.zem is None:
            zem = 'zem = unknown'
        else:
            zem = 'zem={:f}'.format(self.zem)
        txt = '<{:s}: {:s} {:s}, {:s}'.format(
                self.__class__.__name__, self.coord.icrs.ra.to_string(unit=u.hour,sep=':', pad=True),
                self.coord.icrs.dec.to_string(sep=':',pad=True,alwayssign=True), zem)

        # Type?
        if self.em_type is not None:
            txt = txt + ', em_type={:s}'.format(self.em_type)

        # Finish
        txt += '>'
        return (txt)

max_stars_repo_name: pyigm/pyigm

max_issues_repo_name: pyigm/pyigm

max_forks_repo_name: pyigm/pyigm

reason: from astropy


text:
from __future__ import print_function
import pytest

import numpy as np

# FIDIMAG libraries
from fidimag.common import CuboidMesh as Mesh
from fidimag.atomistic import Sim, UniformExchange, Anisotropy
import fidimag.common.constant as const

mu0 = 4 * np.pi * 1e-7

# Initial State, a rough DW in a 1D chain
def init_m(pos):

    z = pos[2]

    if z < 450:
        return (0, 0, 1)
    elif 450 <= z < 550:
        return (1, 1, 0)
    else:
        return (0, 0, -1)

def relax_system(mesh):

    # Only relaxation
    sim = Sim(mesh, name='relax')

    # Simulation parameters
    sim.driver.set_tols(rtol=1e-8, atol=1e-10)
    sim.alpha = 0.5
    sim.driver.gamma = 2.211e5 / mu0
    sim.mu_s = 1e-27 / mu0
    sim.driver.do_precession = False

    # The initial state passed as a function
    sim.set_m(init_m)
    # sim.set_m(np.load('m0.npy'))

    # Energies
    exch = UniformExchange(J=2e-20)
    sim.add(exch)

    anis = Anisotropy(0.01*2e-20, axis=(0, 0, 1))
    sim.add(anis)

    # dmi = DMI(D=8e-4)
    # sim.add(dmi)

    # Start relaxation and save the state in m0.npy
    sim.relax(dt=1e-14, stopping_dmdt=1e4, max_steps=5000,
              save_m_steps=None, save_vtk_steps=None)

    np.save('m0.npy', sim.spin)
    #sim.save_vtk()

# This function excites the system with a
# current in the x-direction using Spin Transfer Torque
# formalism
def excite_system(mesh, time=0.1, snaps=11):

    # Specify the stt dynamics in the simulation
    sim = Sim(mesh, name='dyn', driver='llg_stt')

    # Set the simulation parameters
    sim.driver.set_tols(rtol=1e-12, atol=1e-12)
    sim.driver.gamma = 2.211e5 / mu0
    sim.mu_s = 1e-27 / mu0
    sim.alpha = 0.05

    sim.set_m(np.load('m0.npy'))

    # Energies
    exch = UniformExchange(J=2e-20)
    sim.add(exch)

    anis = Anisotropy(0.01*2e-20, axis=(0, 0, 1))
    sim.add(anis)
    # dmi = DMI(D=8e-4)
    # sim.add(dmi)

    # Set the current in the x direction, in A / m
    # beta is the parameter in the STT torque
    sim.driver.jz = -1e12
    sim.driver.beta = 0.1

    # The simulation will run for x ns and save
    # 'snaps' snapshots of the system in the process
    ts = np.linspace(0, time * 1e-9, snaps)

    for t in ts:
        print('time', t)
        sim.driver.run_until(t)
        #sim.save_vtk()
    np.save('m1.npy', sim.spin)

    print(np.load('m1.npy')[:100])

def test_stt_dw_atomistic():
    mesh = Mesh(nx=1, ny=1, nz=1000,
                dx=1.0, dy=1.0, dz=1.0,
                unit_length=1e-9)

    # Relax the initial state. It will save the last state
    # to the m0.npy file
    relax_system(mesh)

    m0 = np.load('m0.npy')
    m0.shape = (-1, 3)
    mx0 = np.average(m0[:, 2])

    excite_system(mesh)

    m1 = np.load('m1.npy')
    m1.shape = (-1, 3)
    mx1 = np.average(m1[:, 2])

    mu_s = (1e-27 / mu0)
    v = 1e-27
    p = 0.5
    jx = 1e12
    alpha = 0.05
    beta = 0.1
    time = 0.1

    u = const.g_e * const.mu_B / (2 * const.c_e) * v / mu_s * p * jx
    distance = (1 + alpha * beta) / (1 + alpha * alpha) * u * time
    dmx = distance / 1000 * 2

    print(mx0, mx1, dmx, abs(mx1 - mx0 - dmx))
    assert(mx1 - mx0 > 0)
    assert(abs(mx1 - mx0 - dmx) < 5e-5)

if __name__ == '__main__':
    test_stt_dw_atomistic()

max_stars_repo_name: computationalmodelling/fidimag

max_issues_repo_name: computationalmodelling/fidimag

max_forks_repo_name: computationalmodelling/fidimag


text:
import os
import errno
import argparse

import numpy as np
import tensorflow as tf

# Data loader/manager
from data_manager import DataManager

# Models
from logistic_regression import LogisticRegression
from feed_forward_neural_network import FeedForwardNN
from lstm_network import LSTMNN

def parse_arguments():
    """
        Parse the command line arguments of the program.
    """

    parser = argparse.ArgumentParser(description='Train a model on data.')
    parser.add_argument(
        "data_dir",
        type=str,
        nargs="?",
        help="The data directory",
    )
    parser.add_argument(
        "output_dir",
        type=str,
        nargs="?",
        help="The output directory",
        default="out/",
    )
    parser.add_argument(
        "-b",
        "--batch_size",
        type=int,
        nargs="?",
        help="The number of examples in a batch",
        default=128
    )
    parser.add_argument(
        "-c",
        "--epoch_count",
        type=int,
        nargs="?",
        help="Number of training steps.",
        default=100
    )
    parser.add_argument(
        "-l",
        "--learning_rate",
        type=float,
        nargs="?",
        help="The initial learning rate.",
        default=0.01
    )
    parser.add_argument(
        "-r",
        "--train_test_ratio",
        type=float,
        nargs="?",
        help="The ratio between train and test examples in the data.",
        default=0.85
    )
    parser.add_argument(
        "-m",
        "--model",
        type=int,
        nargs="?",
        help="The model to use:\n \tLogReg (1)\n \tFeedForwardNN (2)\n \tLSTMNN (3)"
    )

    return parser.parse_args()

def main():
    """
        Description: Train the model, outputs the performance, saves the .pb
        Usage: python run.py [data_dir] [ratio] [out_dir] [epoch_count]
    """

    # Parse the arguments
    args = parse_arguments()

    # Create the directory if it does not exist.
    try:
        os.makedirs(args.output_dir)
    except OSError as e:
        if e.errno != errno.EEXIST:
            raise e

    data_manager = DataManager(args.data_dir, args.train_test_ratio, args.batch_size)

    model = None

    if args.model == 1:
        model = LogisticRegression(
            args.learning_rate,
            args.epoch_count,
            data_manager,
            args.output_dir
        )
    elif args.model == 2:
        model = FeedForwardNN(
            args.learning_rate,
            args.epoch_count,
            data_manager,
            args.output_dir
        )
    elif args.model == 3:
        model = LSTMNN(
            args.learning_rate,
            args.epoch_count,
            data_manager,
            args.output_dir
        )
    else:
        raise Exception("Unknown model")

    try:
        raise Exception()
        model.load(
            os.path.join(
                args.output_dir,
                'graph.pb'
            )
        )
    except:
        model.train()

    model.test()

    model.save(
        os.path.join(
            args.output_dir,
            'graph.pb'
        )
    )

if __name__ == '__main__':
    main()

max_issues_repo_name: keaideii/ML-IDS

max_forks_repo_name: keaideii/ML-IDS


text:
import os
import warnings
import tempfile

import pandas as pd
import numpy as np
from scipy.stats import pearsonr

import tensorflow.keras as keras
from keras import backend as K
from keras.models import Model,model_from_json
from keras.layers import Dense,Dropout,Input
from keras.callbacks import EarlyStopping
import keras.losses

import tensorflow as tf

os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'

def get_distance_matrix(raw, n_pred=None):

    VMR = raw.std() / raw.mean()
    VMR[np.isinf(VMR)] = 0

    if n_pred is None:
        potential_pred = raw.columns[VMR > 0]
    else:
        print("Using {} predictors".format(n_pred))
        potential_pred = VMR.sort_values(ascending=False).index[:n_pred]

    covariance_matrix = pd.DataFrame(np.abs(np.corrcoef(raw.T.loc[potential_pred])),
                                     index=potential_pred,
                                     columns=potential_pred).fillna(0)
    return covariance_matrix

def wMSE(y_true, y_pred, binary=False):
    if binary:
        weights = tf.cast(y_true>0, tf.float32)
    else:
        weights = y_true
    return tf.reduce_mean(weights*tf.square(y_true-y_pred))

def inspect_data(data):
    # Check if there area any duplicated cell/gene labels

    if sum(data.index.duplicated()):
        print("ERROR: duplicated cell labels. Please provide unique cell labels.")
        exit(1)

    if sum(data.columns.duplicated()):
        print("ERROR: duplicated gene labels. Please provide unique gene labels.")
        exit(1)

    max_value = np.max(data.values)
    if max_value < 10:
        print("ERROR: max value = {}. Is your data log-transformed? Please provide raw counts"
              .format(max_value))
        exit(1)

    print("Input dataset is {} cells (rows) and {} genes (columns)"
          .format(*data.shape))
    print("First 3 rows and columns:")
    print(data.iloc[:3,:3])

class MultiNet:

    def __init__(self,
                 learning_rate=1e-4,
                 batch_size=64,
                 max_epochs=500,
                 patience=5,
                 ncores=-1,
                 loss="wMSE",
                 output_prefix=tempfile.mkdtemp(),
                 sub_outputdim=512,
                 verbose=1,
                 seed=1234,
                 architecture=None
    ):
        self.NN_parameters = {"learning_rate": learning_rate,
                              "batch_size": batch_size,
                              "loss": loss,
                              "architecture": architecture,
                              "max_epochs": max_epochs,
                              "patience": patience}
        self.sub_outputdim = sub_outputdim
        self.outputdir = output_prefix
        self.verbose = verbose
        self.seed = seed
        self.setCores(ncores)

    def setCores(self, ncores):
        if ncores > 0:
            self.ncores = ncores
        else:
            self.ncores = os.cpu_count()
            print("Using all the cores ({})".format(self.ncores))

    def loadDefaultArchitecture(self):
        self.NN_parameters['architecture'] = [
                {"type": "dense", "neurons": self.sub_outputdim//2, "activation": "relu"},
                {"type": "dropout", "rate": 0.2},
            ]

    def save(self, model):
        os.system("mkdir -p {}".format(self.outputdir))

        model_json = model.to_json()

        with open("{}/model.json".format(self.outputdir), "w") as json_file:
            json_file.write(model_json)

        # serialize weights to HDF5
        model.save_weights("{}/model.h5".format(self.outputdir))
        print("Saved model to disk in {}".format(self.outputdir))

    def load(self):
        json_file = open('{}/model.json'.format(self.outputdir), 'r')
        loaded_model_json = json_file.read()
        json_file.close()
        model = model_from_json(loaded_model_json)
        model.load_weights('{}/model.h5'.format(self.outputdir))

        return model

    def build(self, inputdims):
        if self.NN_parameters['architecture'] is None:
            self.loadDefaultArchitecture()

        print(self.NN_parameters['architecture'])

        inputs = [ Input(shape=(inputdim,)) for inputdim in inputdims ]
        outputs = inputs

        for layer in self.NN_parameters['architecture']:
            if layer['type'].lower() == 'dense':
                outputs = [ Dense(layer['neurons'], activation=layer['activation'])(output)
                            for output in outputs ]
            elif layer['type'].lower() == 'dropout':
                outputs = [ Dropout(layer['rate'], seed=self.seed)(output)
                            for output in outputs]
            else:
                print("Unknown layer type.")

        outputs = [Dense(self.sub_outputdim, activation="softplus")(output)
                   for output in outputs]

        model = Model(inputs=inputs, outputs=outputs)

        loss = self.NN_parameters['loss']

        if loss in [k for k, v in globals().items() if callable(v)]:
            # if loss is a defined function
            loss = eval(self.NN_parameters['loss'])

        if not callable(loss):
            # it is defined in Keras
            if hasattr(keras.losses, loss):
                loss = getattr(keras.losses, loss)
            else:
                print('Unknown loss: {}. Aborting.'.format(loss))
                exit(1)

        model.compile(optimizer=keras.optimizer_v2.adam.Adam(lr=self.NN_parameters['learning_rate']),
                      loss=loss)

        return model

    def fit(self,
            raw,
            cell_subset=1,
            NN_lim=None,
            genes_to_impute=None,
            n_pred=None,
            ntop=5,
            minVMR=0.5,
            mode='random',
    ):

        inspect_data(raw)

        if self.seed is not None:
            np.random.seed(self.seed)

        if cell_subset != 1:
            if cell_subset < 1:
                raw = raw.sample(frac=cell_subset)
            else:
                raw = raw.sample(cell_subset)

        gene_metric = (raw.var()/(1+raw.mean())).sort_values(ascending=False)
        gene_metric = gene_metric[gene_metric > 0]

        if genes_to_impute is None:
            genes_to_impute = self.filter_genes(gene_metric, minVMR, NN_lim=NN_lim)
        else:
            # Make the number of genes to impute a multiple of the network output dim
            n_genes = len(genes_to_impute)
            if n_genes % self.sub_outputdim != 0:
                print("The number of input genes is not a multiple of {}. Filling with other genes.".format(n_genes))
                fill_genes = gene_metric.index[:self.sub_outputdim-n_genes]

                if len(fill_genes) < self.sub_outputdim-n_genes:
                    # Not enough genes in gene_metric. Sample with replacement
                    rest = self.sub_outputdim - n_genes - len(fill_genes)
                    fill_genes = np.concatenate([fill_genes,
                                                 np.random.choice(gene_metric.index, rest, replace=True)])

                genes_to_impute = np.concatenate([genes_to_impute, fill_genes])

        covariance_matrix = get_distance_matrix(raw, n_pred=n_pred)

        self.setTargets(raw.reindex(columns=genes_to_impute), mode=mode)
        self.setPredictors(covariance_matrix, ntop=ntop)

        print("Normalization")
        norm_data = np.log1p(raw).astype(np.float32) # normalizer.transform(raw)

        np.random.seed(self.seed)
        tf.random.set_seed(self.seed)

        tf.config.threading.set_inter_op_parallelism_threads(self.ncores)
        tf.config.threading.set_intra_op_parallelism_threads(self.ncores)

        print("Building network")
        model = self.build([len(genes) for genes in self.predictors])

        test_cells = np.random.choice(norm_data.index, int(0.05 * norm_data.shape[0]), replace=False)
        train_cells = np.setdiff1d(norm_data.index, test_cells)

        X_train = [norm_data.loc[train_cells, inputgenes].values for inputgenes in self.predictors]
        Y_train = [norm_data.loc[train_cells, targetgenes].values for targetgenes in self.targets]

        X_test = [norm_data.loc[test_cells, inputgenes].values for inputgenes in self.predictors]
        Y_test = [norm_data.loc[test_cells, targetgenes].values for targetgenes in self.targets]

        print("Fitting with {} cells".format(norm_data.shape[0]))
        result = model.fit(X_train, Y_train,
                           validation_data=(X_test,Y_test),
                           epochs=self.NN_parameters["max_epochs"],
                           batch_size=self.NN_parameters["batch_size"],
                           callbacks=[EarlyStopping(monitor='val_loss',
                                                    patience=self.NN_parameters["patience"])],
                           verbose=self.verbose)

        self.trained_epochs = len(result.history['loss'])
        print("Stopped fitting after {} epochs".format(self.trained_epochs))

        self.save(model)

        # Save some metrics on test data
        Y_test_raw = np.hstack(Y_test).flatten()
        Y_test_imputed = np.hstack(model.predict(X_test)).flatten()

        # Keep only positive values (since negative values could be dropouts)
        Y_test_imputed = Y_test_imputed[Y_test_raw>0]
        Y_test_raw = Y_test_raw[Y_test_raw>0]

        self.test_metrics = {
            'correlation': pearsonr(Y_test_raw,Y_test_imputed)[0],
            'MSE': np.sum((Y_test_raw-Y_test_imputed)**2)/len(Y_test_raw)
        }

        return self

    def predict(self,
                raw,
                imputed_only=False,
                policy="restore"):

        norm_raw = np.log1p(raw)

        inputs = [ norm_raw.loc[:,predictors].values.astype(np.float32)
                   for predictors in self.predictors ]

        model = self.load()

        predicted = model.predict(inputs)
        if len(inputs)>1:
            predicted = np.hstack(predicted)

        predicted = pd.DataFrame(predicted, index=raw.index, columns=self.targets.flatten())

        predicted = predicted.groupby(by=predicted.columns, axis=1).mean()
        not_predicted = norm_raw.drop(self.targets.flatten(), axis=1)

        imputed = (pd.concat([predicted,not_predicted],axis=1)
                   .loc[raw.index, raw.columns]
                   .values)

        # To prevent overflow
        imputed[ (imputed > 2*norm_raw.values.max()) | (np.isnan(imputed)) ] = 0
        # Convert back to counts
        imputed = np.expm1(imputed)

        if policy == "restore":
            print("Filling zeros")
            mask = (raw.values > 0)
            imputed[mask] = raw.values[mask]
        elif policy == "max":
            print("Imputing data with 'max' policy")
            mask = (raw.values > imputed)
            imputed[mask] = raw.values[mask]

        imputed = pd.DataFrame(imputed, index=raw.index, columns=raw.columns)

        if imputed_only:
            return imputed.loc[:, predicted.columns]
        else:
            return imputed

    def filter_genes(self,
                    gene_metric, # assumes gene_metric is sorted
                    threshold,
                    NN_lim=None
    ):
        if not str(NN_lim).isdigit():
            NN_lim = (gene_metric > threshold).sum()

        n_subsets = int(np.ceil(NN_lim / self.sub_outputdim))
        genes_to_impute = gene_metric.index[:n_subsets*self.sub_outputdim]

        rest = self.sub_outputdim - (len(genes_to_impute) % self.sub_outputdim)

        if rest > 0:
            fill_genes = np.random.choice(gene_metric.index, rest)
            genes_to_impute = np.concatenate([genes_to_impute, fill_genes])

        print("{} genes selected for imputation".format(len(genes_to_impute)))

        return genes_to_impute

    def setTargets(self,data, mode='random'):

        n_subsets = int(data.shape[1]/self.sub_outputdim)

        if mode == 'progressive':
            self.targets = data.columns.values.reshape([n_subsets, self.sub_outputdim])
        else:
            self.targets = np.random.choice(data.columns,
                                            [n_subsets, self.sub_outputdim],
                                            replace=False)

    def setPredictors(self, covariance_matrix, ntop=5):
        self.predictors = []

        for i,targets in enumerate(self.targets):

            genes_not_in_target = np.setdiff1d(covariance_matrix.columns, targets)

            if genes_not_in_target.size == 0:
                warnings.warn('Warning: number of target genes lower than output dim. Consider lowering down the sub_outputdim parameter',
                              UserWarning)
                genes_not_in_target = covariance_matrix.columns

            subMatrix = ( covariance_matrix
                          .loc[targets, genes_not_in_target]
                          )
            sorted_idx = np.argsort(-subMatrix.values, axis=1)
            predictors = subMatrix.columns[sorted_idx[:,:ntop].flatten()]

            self.predictors.append(predictors.unique())

            print("Net {}: {} predictors, {} targets"
                  .format(i,len(np.unique(predictors)),len(targets)))

    def score(self, data, policy=None):
        warnings.warn(
            "This method is deprecated. Please use model.test_metrics to measure model accuracy instead",
            DeprecationWarning)
        Y_hat = self.predict(data, policy=policy)
        Y = data.loc[Y_hat.index, Y_hat.columns]

        return pearsonr(Y_hat.values.reshape(-1), Y.values.reshape(-1))

max_stars_repo_name: mestagtx/deepimpute

max_issues_repo_name: mestagtx/deepimpute

max_forks_repo_name: mestagtx/deepimpute

reason: import numpy,from scipy