text:
# Copyright (C) 2020 Arm Limited or its affiliates. All rights reserved.
#
# SPDX-License-Identifier: Apache-2.0
#
# 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
#
# 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.
#
# Description:
# Contains external APIs
from enum import auto
from enum import Enum
from typing import List
from typing import NamedTuple
from typing import Optional
from typing import Tuple

import numpy

API_version_major = 1
API_version_minor = 0
api_version = f"{API_version_major}.{API_version_minor}"

class NpuAccelerator(Enum):
    """
    Supported accelerators
    """

    Ethos_U55_32 = auto()
    Ethos_U55_64 = auto()
    Ethos_U55_128 = auto()
    Ethos_U55_256 = auto()
    Ethos_U65_256 = auto()
    Ethos_U65_512 = auto()

class NpuElementWiseOp(Enum):
    """
    Elementwise operation
    """

    ADD = auto()
    SUB = auto()
    MUL = auto()
    ABS = auto()
    MIN = auto()
    MAX = auto()
    LRELU = auto()  # Leaky relu
    CLZ = auto()  # Number leading zeros
    SHR = auto()  # Rounded right-shift
    SHL = auto()  # Bitwise shift-left

class NpuPoolingOp(Enum):
    """
    Pooling operation
    """

    MAX = auto()
    AVERAGE = auto()
    REDUCE_SUM = auto()

class NpuActivationOp(Enum):
    """
    Activation function
    """

    NONE_OR_RELU = auto()  # Clamps output using min/max
    TANH = auto()
    SIGMOID = auto()
    TABLE_LOOKUP = auto()  # Performs table look-up, using the provided table lookup index

class NpuRoundingMode(Enum):
    """
    Available rounding modes
    """

    TFL = auto()  # TensorFlow Lite rounding
    TRUNCATE = auto()  # Truncate towards zero
    NATURAL = auto()  # Round to nearest with x.5 rounded up, towards +infinity

class NpuLayout(Enum):
    """
    Tensor layout of feature maps
    """

    NHWC = auto()
    NHCWB16 = auto()

    def __str__(self):
        return self.name

class NpuResamplingMode(Enum):
    """
    Resampling mode
    """

    NONE = auto()  # No resampling is performed
    NEAREST = auto()  # 2x2 insert nearest
    TRANSPOSE = auto()  # 2x2 transpose

class NpuBlockTraversal(Enum):
    """
    Block-traversal of weights
    """

    DEPTH_FIRST = auto()
    PART_KERNEL_FIRST = auto()

class NpuDataType(Enum):
    """
    Supported data types in feature maps
    """

    UINT8 = 8, False, auto()
    INT8 = 8, True, auto()
    UINT16 = 16, False, auto()
    INT16 = 16, True, auto()
    INT32 = 32, True, auto()

    def is_signed(self) -> bool:
        """Checks if this data type is signed or unsigned"""
        return self.value[1]

    def size_in_bits(self) -> int:
        """ Size of the data type in bits"""
        return self.value[0]

    def size_in_bytes(self) -> int:
        """ Size of the data type in bytes"""
        return self.value[0] // 8

    def min_value(self) -> int:
        """Minimum value of this type"""
        if self.is_signed():
            return -(1 << (self.size_in_bits() - 1))
        else:
            return 0

    def max_value(self) -> int:
        """Maximum value of this type"""
        if self.is_signed():
            return (1 << (self.size_in_bits() - 1)) - 1
        else:
            return (1 << self.size_in_bits()) - 1

    def __str__(self):
        return self.name

    __repr__ = __str__

class NpuAddressRange(NamedTuple):
    """
    Address range
    """

    region: int  # Memory region, a value between 0 and 7
    address: int  # Address, offset from the region's base address
    length: int  # The length of the range, in bytes

    def __str__(self):
        return f"(region={self.region}, address={hex(self.address)}, length={self.length})"

class NpuTileBox(NamedTuple):
    """
    Specifies the addresses and dimensions of the tiles of a feature map.
    A feature map can use 1 to 4 tiles
    """

    height_0: int  # The height of tile 0
    height_1: int  # The height of tile 1, 0 if unused
    width_0: int  # the width of tile 0, and tile 2 (if used)
    addresses: List[int]  # A list of 4 addresses, set unused addresses to 0

class NpuShape3D(NamedTuple):
    """
    Shape of (part of) a feature map
    """

    height: int
    width: int
    depth: int

class NpuQuantization(NamedTuple):
    """
    Quantization parameters
    """

    scale_f32: Optional[float]
    zero_point: int

class NpuPadding(NamedTuple):
    """
    Padding to be applied to a convolution operation
    """

    top: int
    left: int
    bottom: int
    right: int

class NpuActivation:
    """
    Activation function, fused with NPU operations
    """

    def __init__(self, op_type: NpuActivationOp):
        self.op_type = op_type  # The activation operation to be performed
        # min/max are optional
        self.min: Optional[float] = None  # E.g. set to 0.0 for RELU
        self.max: Optional[float] = None  # E.g. set to 6.0 for RELU6
        # Table lookup index, only applicable for TABLE_LOOKUP activation, 0-7
        self.lookup_table_index: int = 0

class NpuFeatureMap:
    """
    Basic information about IFM, IFM2, OFM
    """

    def __init__(self):
        self.data_type: NpuDataType = NpuDataType.UINT8
        # The memory region, a value 0-7
        self.region: int = 0
        # Shape of the feature map
        self.shape: NpuShape3D = NpuShape3D(height=0, width=0, depth=0)
        # The tiles that comprise the feature map. In the normal case when only 1 tile is used,
        # height_0 == self.shape.height, height_1 is 0, width_0 == self.shape.width, addresses[1:] are set to 0
        self.tiles: NpuTileBox = NpuTileBox(height_0=0, height_1=0, width_0=0, addresses=[0, 0, 0, 0])
        self.quantization: Optional[NpuQuantization]
        self.layout: NpuLayout = NpuLayout.NHWC
        # x/y/c strides used by the NPU when traversing the feature map, if None, vela will use default strides
        self.strides: Optional[NpuShape3D] = None

class NpuKernel:
    """
    Kernel information for NPU operations
    """

    def __init__(self, w: int, h: int, stride_x: int = 1, stride_y: int = 1, dilation_x: int = 1, dilation_y: int = 1):
        assert stride_x > 0 and stride_y > 0
        assert dilation_x > 0 and dilation_y > 0
        self.width = w
        self.height = h
        self.stride_x = stride_x
        self.stride_y = stride_y
        self.dilation_x = dilation_x
        self.dilation_y = dilation_y

class NpuOperationType(Enum):
    """
    Type of NPU operation
    """

    Dma = auto()
    Conv2D = auto()
    ConvDepthWise = auto()
    Pooling = auto()
    ElementWise = auto()

class NpuOperation:
    """
    Base class for all NPU operations
    """

    def __init__(self, op_type: NpuOperationType):
        self.op_type = op_type

class NpuDmaOperation(NpuOperation):
    """
    DMA operation
    """

    def __init__(self, src: NpuAddressRange, dest: NpuAddressRange):
        super().__init__(NpuOperationType.Dma)
        self.src = src
        self.dest = dest
        # DMA channel, usually 0 (user channel)
        self.channel: int = 0
        # Channel mode, 0 = external, 1 = internal (should usually be 0)
        self.mode: int = 0

class NpuBlockOperation(NpuOperation):
    """
    Base class for operations which produce an OFM
    """

    def __init__(self, op_type: NpuOperationType):
        super().__init__(op_type)
        self.ifm: Optional[NpuFeatureMap] = None
        self.ifm2: Optional[NpuFeatureMap] = None
        # The non-quantized scalar value in a binary elementwise operation. Only set if IFM2 is scalar
        self.ifm2_scalar: Optional[float] = None
        self.ofm: Optional[NpuFeatureMap] = None
        self.kernel: Optional[NpuKernel] = None
        # Weights, one element for each NPU core, empty if no weights are used.
        # Must have been compressed using weight_compressor.encode_weights()
        self.weights: List[NpuAddressRange] = []
        # Biases, one element for each NPU core, empty if no bias is used.
        # Must have been encoded using weight_compressor.encode_bias()
        self.biases: List[NpuAddressRange] = []
        self.padding: Optional[NpuPadding] = None
        # Optional activation function to be applied
        self.activation: Optional[NpuActivation] = None
        # The block config is the unit of work in which the NPU generates the OFM.
        # If the operation has weights, the depth of the block config must be the same as
        # the ofm depth used in the call to weight_compressor.encode_weights()
        # If set to None, vela will determine a suitable block size (can only be used if there are no weights)
        # If block_config.width and height are set to -1, vela will determine suitable width/height
        self.block_config: Optional[NpuShape3D] = None  # OFM_BLK parameters
        self.rounding_mode: NpuRoundingMode = NpuRoundingMode.TFL
        # Set to True if the operations is fused with a Quantize operation (affects scaling)
        self.fused_quantize: bool = False
        # IFM upscaling to be applied
        self.ifm_upscale: NpuResamplingMode = NpuResamplingMode.NONE

class NpuConv2DOperation(NpuBlockOperation):
    """
    NPU_OP_CONV operation
    """

    def __init__(self):
        super().__init__(NpuOperationType.Conv2D)
        # Block traversal must be consistent with the block_traversal parameter specified in
        # weight_compressor.encode_weights()
        self.block_traversal: NpuBlockTraversal = NpuBlockTraversal.PART_KERNEL_FIRST

class NpuConvDepthWiseOperation(NpuBlockOperation):
    """
    NPU_OP_DEPTHWISE operation
    """

    def __init__(self):
        super().__init__(NpuOperationType.ConvDepthWise)

class NpuPoolingOperation(NpuBlockOperation):
    """
    NPU_OP_POOL operation
    """

    def __init__(self, pooling_op_type: NpuPoolingOp):
        super().__init__(NpuOperationType.Pooling)
        self.sub_op_type: NpuPoolingOp = pooling_op_type
        # Set to a float value for ResizeBilinear operations (affects scaling), else to None
        self.rescale: Optional[float] = None

class NpuElementWiseOperation(NpuBlockOperation):
    """
    NPU_OP_ELEMENTWISE operation
    """

    def __init__(self, elementwise_op_type: NpuElementWiseOp):
        super().__init__(NpuOperationType.ElementWise)
        self.sub_op_type: NpuElementWiseOp = elementwise_op_type
        # Set to True for binary operators where IFM2 should be used as first operand
        self.reversed_operands: bool = False
        # Set to a tuple (scale, shift) for explicit rescale, else to None
        self.rescale: Optional[Tuple] = None

def npu_get_API_version():
    """
    Public facing API to get the API version
    :return: int, the 16 most significant bits, corresponding to major version
            the 16 least significant bits, corresponding to minor version
    """
    version = (API_version_major << 16) | (API_version_minor & 0xFFFF)
    return version

def npu_encode_weights(
    accelerator: NpuAccelerator,
    weights_volume: numpy.ndarray,
    dilation_xy: Tuple[int, int],
    ifm_bitdepth: int,
    ofm_block_depth: int,
    is_depthwise: bool,
    block_traversal: NpuBlockTraversal,
):
    """
    Public facing API to use the Ethos-U weight encoding.

    :param accelerator: NpuAccelerator enum to pick the correct accelerator
    :param weights_volume: numpy.ndarray in OHWI layout with a shape of four
    :param dilation_xy: a two element tuple of dilation attributes in x,y dimension
    :param ifm_bitdepth: the bitdepth of input feature map
    :param ofm_block_depth: the depth of blocks for processing
    :param is_depthwise: a boolean indicating these weights are used for a depthwise traversal
    :param block_traversal: indicates how these weights are traversed on sub-kernel basis
    :return: a bytearray of compressed weights
    """
    from .architecture_features import Accelerator
    from . import weight_compressor

    acc = Accelerator.from_npu_accelerator(accelerator)
    return weight_compressor.encode_weights(
        acc, weights_volume, dilation_xy, ifm_bitdepth, ofm_block_depth, is_depthwise, block_traversal
    )

def npu_encode_bias(bias: numpy.int64, scale: int, shift: int):
    """
    Public facing API to pack bias and scale values as required by the hardware
    :param bias: 64-bit signed number that includes 40-bit signed bias
    :param scale: 32-bit scale value
    :param shift: 6-bit shift value
    :return: packed 80-bit [0(2-bits),shift(6-bits),scale(32-bits),bias(40-bits)]
    """
    from . import weight_compressor

    return weight_compressor.encode_bias(bias, scale, shift)

def npu_generate_register_command_stream(npu_op_list: List[NpuOperation], accelerator: NpuAccelerator) -> List[int]:
    """
    Public facing API for generating an Ethos-U register command stream.
    Calculates dependencies between commands and inserts wait operations if needed.

    :param npu_op_list: List[NpuOperation] list of high level NPU operations
    :param accelerator: NpuAccelerator enum to pick the correct accelerator
    :return register commands, as a list of 32-bit integers
    """
    from . import register_command_stream_generator

    return register_command_stream_generator.generate_register_command_stream(npu_op_list, accelerator)

def npu_create_driver_payload(register_command_stream: List[int], accelerator: NpuAccelerator) -> bytes:
    """
    Public facing API for generating driver payload, containing a driver header
    and the given Ethos-U register command stream.
    Returns the payload, in little endian format, which must be placed in memory on a 16-byte aligned
    address.

    :param register_command_stream: List[int] register commands, as a list of 32-bit integers
    :param accelerator: NpuAccelerator enum to pick the correct accelerator
    :return driver payload, as a byte array
    """
    from . import driver_actions

    return driver_actions.npu_create_driver_payload(register_command_stream, accelerator)

max_stars_repo_name: DemonGiggle/vela

max_issues_repo_name: DemonGiggle/vela

max_forks_repo_name: DemonGiggle/vela


text:
import math
from math import sqrt, sin, cos, tan, degrees, radians, acos
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd

# IMPORT OF DATASET
dt = 0.1  # FRAME RATE
pixel12nm = 0.117  # PIXEL CONVERSION FACTOR (um)
df = pd.read_csv('1.csv')
frames_val = df.iloc[:, 0].values
x_val = df.iloc[:, 2].values
y_val = df.iloc[:, 3].values
frames = [i for i in range(0, len(frames_val))]
frames = np.array(frames)
x = np.empty(len(frames_val), dtype=float)
y = np.empty(len(frames_val), dtype=float)
for i in range(0, len(frames)):
    x[i] = x_val[i] * pixel12nm
    y[i] = y_val[i] * pixel12nm

# DIRECT MOTION TRAJECTORY ANALYSIS
# Number of trajectory points grouped together
points_grouped = 4
# Threshold ratio to distinguish between direct(>difTH) or non-direct (<difTH)
threshold_ratio = 0.7
# Minimum number of points to consider a pause [3 points for 0.1s fr rate -> 0.2s pausetime]
min_points = 5

# Displacement between adjacent frames (um)
disp = []
for i in range(0, len(frames_val) - 1):
    d = sqrt((x[i + 1] - x[i])**2 + (y[i + 1] - y[i])**2)
    disp.append(d)
disp = np.array(disp)
# print((disp))

# Displacement every "points_grouped" frames (um)
grouped_disp = []
for i in range(0, (len(frames_val) - points_grouped + 1)):
    d = sqrt((x[i + (points_grouped - 1)] - x[i])**2 +
             (y[i + (points_grouped - 1)] - y[i])**2)
    grouped_disp.append(d)
grouped_disp = np.array(grouped_disp)
# print(grouped_disp)

# Sum of displacements for first "points_grouped" points (um)
sum_disp = []
for i in range(0, len(disp) - (points_grouped - 2)):
    sum_disp.append(sum(disp[i:(i + (points_grouped - 1))]))
sum_disp = np.array(sum_disp)
# print((sum_disp))
ratio = np.true_divide(grouped_disp, sum_disp)
# print(ratio.size)

# Assign an average ratio to each data point
p = []
for i in range(0, len(ratio) + (points_grouped - 1)):
    if i == 0:
        pc = np.mean(ratio[0])
    elif i < points_grouped and i > 0:
        pc = np.mean(ratio[:i + 1])
    elif i > len(ratio):
        pc = np.mean(ratio[i - (points_grouped - 1):])
    else:
        pc = np.mean(ratio[i - (points_grouped - 1):i + 1])
    p.append(pc)
p = np.array(p)
# print(p)

# First screening to identify truly direct regions based on displacement ratio condition
direct_1 = []
for i in range(0, len(p)):
    if p[i] >= threshold_ratio:
        direct_1.append(1)
    else:
        direct_1.append(0)
direct_1 = np.array(direct_1)
# print(direct_1.shape)

# 2. Direct motion identification based on angle between adjacent tracsvectors

# Threshold angle to distinguish between direct(<threshold_angle) or non-direct (>threshold_angle)
threshold_angle = 90

disp_vector = []
for i in range(0, len(frames) - 1):
    d = ((x[i + 1] - x[i]), (y[i + 1] - y[i]))
    disp_vector.append(d)
disp_vector = np.array(disp_vector)

nv = []
for i in range(0, len(disp_vector)):
    d = np.linalg.norm(disp_vector[i, :])
    nv.append(d)
nv = np.array(nv)
# print(nv)

cos_theta = []
sign = []
for i in range(0, len(disp_vector) - 1):
    dot_mul = np.sum(disp_vector[i, :].conj(
    ) * disp_vector[i + 1, :], axis=0) / (nv[i] * nv[i + 1])
    cos_theta.append(dot_mul)

    mul = (disp_vector[i, 0] * disp_vector[i, 1])
    sign.append(np.sign(mul))

cos_theta = np.array(cos_theta)
sign = np.array(sign)
# print(cos_theta, sign)

theta_in_degrees = []
theta_in_degrees_sign = []
for i in range(0, len(cos_theta)):
    theta = acos(cos_theta[i]) * 180 / math.pi
    theta_in_degrees.append(theta)

    theta_in_degrees_sign.append(
        sign[i] * (acos(cos_theta[i]) * 180 / math.pi))

# Angle in degrees (without sign - mod value)
theta_in_degrees = np.array(theta_in_degrees)
# Angle in degrees (with sigh - Actual Angle)
theta_in_degrees_sign = np.array(theta_in_degrees_sign)
# print(theta_in_degrees_sign)

# Second screening to identify truly direct regions based on angle condition
direct_2 = []
direct_2.append(direct_1[0])

for i in range(0, len(theta_in_degrees)):
    if theta_in_degrees[i] < threshold_angle:
        direct_2.append(1)
    else:
        direct_2.append(0)

direct_2.append(direct_1[len(theta_in_degrees) + 1])
direct_2 = np.array(direct_2)
# print(direct_2)

# 3. Combine both screening methods (Angle & Displacement)
Angle_condition = True
direct_trajectory = []
if Angle_condition:
    for i in range(0, len(direct_2)):
        if direct_1[i] == 1 and direct_2[i] == 1:
            direct_trajectory.append(1)
        else:
            direct_trajectory.append(0)
else:
    direct_trajectory = direct_1

direct_trajectory = np.array(direct_trajectory)
# print(direct_trajectory)

# 4. Filter out short nondirect segments <min_points number of points (sequences of less than 5 zeros)
i = 0
idx = []

while i < len(direct_trajectory):
    if direct_trajectory[i] == 0:
        idx.append(i)
        if len(direct_trajectory) == i and not (not idx):
            if len(idx) < min_points:
                direct_trajectory[idx] = 1

        i = i + 1
    else:
        if not (not idx):
            if len(idx) + 1 < min_points:
                direct_trajectory[idx] = 1

            idx = []
            i = i + 1
        else:
            i = i + 1
# print((direct_trajectory[-1]))

# 5. Filter out direct segments < 0.5pixel max confinement radius
i = 0
s = 0
idx = []
D = []
while i < len(direct_trajectory):
    if direct_trajectory[i] == 1:
        idx.append(i)
        i = i + 1
        if len(direct_trajectory) - 1 == i and not (not idx):
            if idx[0] != 1:
                idx.extend((idx[0] - 1, idx[-1] + 1))
            else:
                idx.append(idx[-1] + 1)
            for s in range(1, len(idx)):
                D.append(sqrt((x[idx[s]] - x[idx[0]])**2 +
                              (y[idx[s]] - y[idx[0]])**2))
            R_max = max(D)

            if (len(idx) < 12 and R_max <= 0.5 * pixel12nm):
                direct_trajectory[idx] = 0
    else:
        if not (not idx):
            if idx[0] != 1:
                idx.extend((idx[0] - 1, idx[-1] + 1))
            else:
                idx.append(idx[-1] + 1)
            for s in range(0, len(idx)):
                D.append(sqrt((x[idx[s]] - x[idx[0]])**2 +
                              (y[idx[s]] - y[idx[0]])**2))
            R_max = max(D)

            if (len(idx) < 12 and R_max <= 0.5 * pixel12nm):
                direct_trajectory[idx] = 0
            idx = []
            D.clear()
            R_max = None
            i = i + 1
        else:
            i = i + 1

# print(direct_trajectory)

# 6. Find non-direct indices (diffusive or confined motion)

non_direct_idx = np.where(direct_trajectory == 0)
# +1 For Python Syntax correction as Python numbering starts from 0
non_direct_idx = np.array(non_direct_idx[0] + 1)

non_direct_start = []
non_direct_end = []

if not(not non_direct_idx.all):
    non_direct_start.append(non_direct_idx[0])
    # j = 0
    # k = 1
# range starting from 1 as range does not include the last number in Python
# so taking values from range(1, len(non_direct_idx)) rather than range(0, len(non_direct_idx)-1)
    for i in range(1, len(non_direct_idx)):
        if i != len(non_direct_idx) - 1:
            if non_direct_idx[i] - non_direct_idx[i - 1] > 1:
                non_direct_end.append(non_direct_idx[i - 1])
                # j = j + 1
                non_direct_start.append(non_direct_idx[i])
                # k = k + 1

        else:
            if non_direct_idx[i] - non_direct_idx[i - 1] > 1:
                non_direct_end.append(non_direct_idx[i - 1])
            else:
                non_direct_end.append(non_direct_idx[i])

else:
    non_direct_start = []
    non_direct_end = []

non_direct_start = np.array(non_direct_start)
non_direct_end = np.array(non_direct_end)
# print(non_direct_idx)
# print(type(non_direct_start))
# print(non_direct_end)

# 7. Indices for direct parts
# direct_start = []
# direct_end=[]

if not non_direct_idx.all:
    direct_start = 0
    direct_end = len(x)
else:
    direct_start = []
    direct_end = []

    for i in range(0, len(non_direct_start)):
        if i == 0:
            if non_direct_start[i] != 1:
                direct_start.append(1)
                direct_end.append(non_direct_start[i])
        else:
            direct_start.append(non_direct_end[i - 1])
            direct_end.append(non_direct_start[i])

    if non_direct_end[len(non_direct_end) - 1] < len(x) - 3:
        direct_start.append(non_direct_end[len(non_direct_end) - 1])
        direct_end.append(len(x))

direct_start = np.array(direct_start)
direct_end = np.array(direct_end)
# print(direct_start)
# print(direct_end)

# CALCULATIONS

# Time difference between frames(Delta Time)
delta_time = []
for i in range(0, len(frames) - 1):
    delta_time.append(np.multiply((frames[i + 1] - frames[i]), dt))

# Instantaneous Speed (um/s)
velocity = np.true_divide(disp, delta_time)

Rc = []
Rc_max = []

track_info = [None] * 5
track_info[0] = ['f']
DATASUBHEADER_keys = ['Phase type', 'idx', 'xy(pixels)', 'Instantaneous Runlength(um)',
                      'Total Runlength(um)', 'Instantaneous speed(um/s)',
                      'Average speed(um/s)', 'Duration(s)', 'D (um2/s)', 'Alpha',
                      'Displacement start-end(nm)', 'Maximum confinement radius(nm)']
DATASUBHEADER = dict.fromkeys(DATASUBHEADER_keys, None)
track_info[1] = DATASUBHEADER
# print(track_info[1]['idx'])
s = 0
m = 1
idx_start = []
idx_end = []
# print(non_direct_start, direct_start)

if not(not non_direct_start.all()) and not (not direct_start.all()):
    idx_start = np.sort(np.concatenate((non_direct_start, direct_start)))
    idx_end = np.sort(np.concatenate((non_direct_end, direct_end)))
elif not non_direct_start:
    idx_start = direct_start
    idx_end = direct_end
else:
    idx_start = non_direct_start
    idx_end = non_direct_end

idx_start = np.append(idx_start, idx_start[0])
idx_end = np.append(idx_end, idx_end[-1])
length_idx = idx_end - idx_start + 1
# print(idx_start)
# arr_1 = []
# for a in range(idx_start[1], idx_end[1] + 1):
#
#     arr_1.append(a)
# print(arr_1)

track_info[1]['Alpha'] = []
for i in range(0, len(idx_start)):
    # PARAMETER TRANSPORT

    # Indices to distinguish direct vs diffusive or Frame ('idx')
    arr_1 = []
    # farr_1 = []
    for a in range(idx_start[i], idx_end[i] + 1):
        arr_1 = np.append(arr_1, a)
    #     farr_1.append(arr_1)
    track_info[1]['idx'] = arr_1
    # print(track_info[1]['idx'])

    # xy coordinates ('xy(pixels)')
    arr_2_x = []
    arr_2_y = []
    arr_2 = []
    for b in range(idx_start[i], idx_end[i] + 1):
        arr_2_x = np.append(arr_2_x, x[b - 1])
        arr_2_y = np.append(arr_2_y, y[b - 1])
        arr_2 = np.append(arr_2_x, arr_2_y)
    track_info[1]['xy(pixels)'] = arr_2
    # print(track_info[1]['xy(pixels)'])

    # Instantaneous runlength (um) ('Instantaneous Runlength(um)')
    arr_3 = []
    for c in range(idx_start[i], idx_end[i]):
        arr_3 = np.append(arr_3, disp[c - 1])
    track_info[1]['Instantaneous Runlength(um)'] = arr_3
    # print(track_info[1]['Instantaneous Runlength(um)'])

    # Total runlength per segment ('Total Runlength(um)')
    track_info[1]['Total Runlength(um)'] = sum(
        disp[idx_start[i] - 1:idx_end[i]])
    # print(track_info[1]['Total Runlength(um)'])

    # Instantaneous speed ('Instantaneous speed(um/s)')
    track_info[1]['Instantaneous speed(um/s)'] = velocity[idx_start[i] -
                                                          1:idx_end[i] - 1]
    # print(track_info[1]['Instantaneous speed(um/s)'])

    # Average speed per segment ('Average speed(um/s)')
    track_info[1]['Average speed(um/s)'] = np.mean(
        velocity[idx_start[i] - 1:idx_end[i] - 1])
    # print(track_info[1]['Average speed(um/s)'])

    # Time spent (PAUSING for diffusive & PROCESSIVITY for direct) ('Duration(s)')
    track_info[1]['Duration(s)'] = dt * ((length_idx[i]) - 1)
    # print(track_info[1]['Duration(s)'])

    x_c = x[idx_start[i] - 1:idx_end[i]]
    y_c = y[idx_start[i] - 1:idx_end[i]]
    D = np.empty(len(x_c))
    for s in range(0, len(x_c)):
        D[s] = (sqrt((x_c[s] - x_c[0])**2 +
                     (y_c[s] - y_c[0])**2))

    R_max = max(D) * 1000
    Rc_max.append(R_max)  # Distance between most separated points
    # print(Rc_max)

    D_2 = sqrt((x[idx_end[i] - 1] - x[idx_start[i] - 1])**2 +  # (-1) in both to compensate for the '0' starting point of Python
               (y[idx_end[i] - 1] - y[idx_start[i] - 1])**2)
    Rc.append(D_2 * 1000)  # Distance between start and end points
    # print(Rc)

    # Confinement diameter (nm) start-end ('Displacement start-end(nm)')
    track_info[1]['Displacement start-end(nm)'] = np.array(Rc)

    # Confinement diameter (nm) maximum ('Maximum confinement radius(nm)')
    track_info[1]['Maximum confinement radius(nm)'] = np.array(Rc_max)

    # MSD ANALYSIS
    msd = []
    tau = []
    coef = []
    coef_pl = []
    if length_idx[i] >= 12:  # Minimum number of points to calculate MSD
        for j in range(1, length_idx[i]):
            tau.append(j * dt)
            # x_msd_start = x[idx_start[i] + j:idx_end[i] + 1]
            # y_msd_start = y[idx_start[i] + j:idx_end[i] + 1]
            # x_msd_end = x[idx_start[i]:idx_end[i] - j + 1]4
            # y_msd_end = y[idx_start[i]:idx_end[i] - j + 1]
            # for h in range(0, len(x_msd_start)):  # np.diff can also be used
            #     xdata = (x_msd_start[h] - x_msd_end[h])
            #     ydata = (y_msd_start[h] - y_msd_end[h])
            # msd.append(np.mean(xdata**2 + ydata**2))
            x = x.reshape((len(x), 1))
            y = y.reshape((len(y), 1))
            t1 = np.power(x[idx_start[i] + j: idx_end[i], :] -
                          x[idx_start[i]: idx_end[i] - j, :], 2)
            t2 = np.power(y[idx_start[i] + j: idx_end[i], :] -
                          y[idx_start[i]: idx_end[i] - j, :], 2)
            msd.append(np.mean(t1 + t2, axis=0)[0])

        tau = np.transpose(tau)
        msd = np.transpose(msd)
        thr = 5
        # assp = np.empty((5, 1), float)
        tau_thr = tau[0:thr].reshape(thr, 1)
        msd_thr = msd[0:thr].reshape(thr, 1)
        mtx = np.ones(len(tau[0:thr])).reshape(thr, 1)
        # for g in range(0, len(tau_thr)):
        #     mtx.append([tau_thr[g], 1], axis=1)
        assp = np.append(tau_thr, mtx, axis=1)
        log_assp = np.append(np.log10(tau_thr), mtx, axis=1)
        coef.append(np.linalg.lstsq(assp, msd_thr, rcond=None)[0])
        coef_pl.append(np.linalg.lstsq(
            log_assp, np.log10(msd_thr), rcond=None)[0][0])
        # print(coef)
        if coef_pl[0][0] > 0:  # Alpha coefficient > 0
            track_info[1]['D (um2/s)'] = coef[0] / 4
            track_info[1]['Alpha'].append(coef_pl[0][0])
        else:
            track_info[1]['D (um2/s)'] = None
            track_info[1]['Alpha'].append(np.nan)
    else:
        track_info[1]['D (um2/s)'] = None
        track_info[1]['Alpha'].append(np.nan)

    # print(round(np.nan * 100) / 100)
# print(idx_start[i])

# Assign Phase Type
    if i == len(idx_start):
        track_info[1]['Phase type'] = 'Full Track'
    elif direct_start.any() == idx_start[i]:
        if not np.isnan(track_info[1]['Alpha'][-1]):
            if round(track_info[1]['Alpha'][-1] * 100) / 100 >= 1.45:
                if track_info[1]['Average speed(um/s)'] > 0.25:
                    track_info[1]['Phase type'] = 'Direct-Fast'
                else:
                    track_info[1]['Phase type'] = 'Direct-Slow'

            elif (track_info[1]['Alpha'][-1] * 100) / 100 >= 1.0 and (track_info[1]['Alpha'][-1] * 100) / 100 < 1.45:
                track_info[1]['Phase type'] = 'Diffusive'
                if track_info[1]['Average speed(um/s)'] > 0.25:
                    track_info[1]['Phase type'] = 'Diffusive CHECK'
                else:
                    track_info[1]['Phase type'] = 'Diffusive'

        else:
            if track_info[1]['Average speed(um/s)'] > 0.25:
                track_info[1]['Phase type'] = 'Direct-Short-Fast'
            else:
                track_info[1]['Phase type'] = 'Direct-Short-Slow'

    else:
        if not np.isnan(track_info[1]['Alpha'][-1]):
            if (track_info[1]['Alpha'][-1] * 100) / 100 < 1:
                track_info[1]['Phase type'] = 'Confined'
            elif (track_info[1]['Alpha'][-1] * 100) / 100 >= 1.0 and (track_info[1]['Alpha'][-1] * 100) / 100 < 1.45:
                track_info[1]['Phase type'] = 'Diffusive'
            else:
                track_info[1]['Phase type'] = 'Recategorize to Dir'

        else:
            track_info[1]['Phase type'] = 'Confined-Short'

    if track_info[1]['Phase type'] == 'Direct-Fast' or track_info[1]['Phase type'] == 'Direct-Short-Fast':
        plt.plot(x[idx_start[i]: idx_end[i]],
                 y[idx_start[i]: idx_end[i]], '-*', color='blue', label='Direct Fast')
    elif track_info[1]['Phase type'] == 'Direct-Slow' or track_info[1]['Phase type'] == 'Direct-Short-Slow':
        plt.plot(x[idx_start[i]: idx_end[i]],
                 y[idx_start[i]: idx_end[i]], '-*', color='cyan', label='Direct Slow')
    elif track_info[1]['Phase type'] == 'Confined' or track_info[1]['Phase type'] == 'Confined-Short':
        plt.plot(x[idx_start[i]: idx_end[i]],
                 y[idx_start[i]: idx_end[i]], '-*', color='red', label='Confined')
    elif track_info[1]['Phase type'] == 'Diffusive':
        plt.plot(x[idx_start[i]: idx_end[i]],
                 y[idx_start[i]: idx_end[i]], '-*', color='magenta', label='Diffusive')

plt.xlabel('X - coordinates(um)')
plt.ylabel('Y - coordinates(um)')
plt.title('Particle Trajectory')
plt.legend()
# plt.colorbar(mappable=np.array(frames))
plt.grid()
plt.show()

ext: py

reason: import numpy


text:
"""Calib file from KITTI object detection benchmark, with modifications for KITTI road benchmark"""

import os

import numpy as np

root_dir = os.path.dirname(os.path.abspath(__file__))
data_dir = os.path.join(root_dir, 'data')
image_shape = 375, 1242

def get_drive_dir(drive, date='2011_09_26'):
    return os.path.join(data_dir, date, date + '_drive_%04d_sync' % drive)

def get_inds(path, ext='.png'):
    inds = [int(os.path.splitext(name)[0]) for name in os.listdir(path)
            if os.path.splitext(name)[1] == ext]
    inds.sort()
    return inds

def read_calib_file(path):
    float_chars = set("0123456789.e+- ")

    data = {}

    # print(path)

    with open(path, 'r') as f:
        for line in f.readlines():
            key, value = line.split(':', 1)
            value = value.strip()
            data[key] = value
            if float_chars.issuperset(value):
                # try to cast to float array
                try:
                    data[key] = np.array(list(map(float, value.split(' '))))
                except ValueError:
                    pass  # casting error: data[key] already eq. value, so pass

    # print(data)

    return data

def homogeneous_transform(points, transform, keep_last=False):
    """
    Parameters
    ----------
    points : (n_points, M) array-like
        The points to transform. If `points` is shape (n_points, M-1), a unit
        homogeneous coordinate will be added to make it (n_points, M).
    transform : (M, N) array-like
        The right-multiplying transformation to apply.
    """

    transform = transform.T

    points = np.asarray(points)
    transform = np.asarray(transform)
    n_points, D = points.shape
    M, N = transform.shape

    # do transformation in homogeneous coordinates
    if D == M - 1:
        points = np.hstack([points, np.ones((n_points, 1), dtype=points.dtype)])
    elif D != M:
        raise ValueError("Number of dimensions of points (%d) does not match"
                         "input dimensions of transform (%d)." % (D, M))
    n_points, D = points.shape

    new_points = np.dot(points, transform)

    # normalize homogeneous coordinates
    if not keep_last:
        new_points = new_points[:, :-1] / new_points[:, [-1]]
    else:
        new_points = np.hstack(
            [new_points[:, :-1] / new_points[:, [-1]], new_points[:, [-1]]])

    return new_points

def filter_disps(xyd, shape, max_disp=255, return_mask=False):
    x, y, d = xyd.T
    mask = ((x >= 0) & (x <= shape[1] - 1) &
            (y >= 0) & (y <= shape[0] - 1) &
            (d >= 0) & (d <= max_disp))
    xyd = xyd[mask]
    return (xyd, mask) if return_mask else xyd

def filter_depths(xyd, shape):
    x, y, d = xyd.T
    mask = ((x >= 0) & (x <= shape[1] - 1) &
            (y >= 0) & (y <= shape[0] - 1))
    xyd = xyd[mask]
    return xyd

class Calib(object):
    """Convert between coordinate frames.
    This class loads the calibration data from file, and creates the
    corresponding transformations to convert between various coordinate
    frames.
    Each `get_*` function returns a 3D transformation in homogeneous
    coordinates between two frames. All transformations are right-multiplying,
    and can be applied with `homogeneous_transform`.
    """

    def __init__(self, calib_path, color=False):
        self.calib_path = calib_path
        self.data = read_calib_file(self.calib_path)

    def get_imu2velo(self):
        RT_imu2velo = np.eye(4)
        RT_imu2velo[:3, :3] = self.data['R0_rect'].reshape(3, 3)
        RT_imu2velo[:3, 3] = self.data['Tr_imu_to_velo']
        return RT_imu2velo.T

    def get_velo2rect(self):
        RT_velo2cam = np.eye(4)
        RT_velo2cam[:3, :] = self.data['Tr_velo_to_cam'].reshape(3, 4)
        R_rect00 = np.eye(4)
        R_rect00[:3, :3] = self.data['R0_rect'].reshape(3, 3)

        RT_velo2rect = np.dot(R_rect00, RT_velo2cam)
        return RT_velo2rect

    def get_rect2disp(self):
        cam0, cam1 = (0, 1) if not self.color else (2, 3)
        P_rect0 = self.data['P%d' % cam0].reshape(3, 4)
        P_rect1 = self.data['P%d' % cam1].reshape(3, 4)

        P0, P1, P2 = P_rect0
        Q0, Q1, Q2 = P_rect1

        # create disp transform
        T = np.array([P0, P1, P0 - Q0, P2])
        return T.T

    def get_imu2rect(self):
        return np.dot(self.get_imu2velo(), self.get_velo2rect())

    def get_imu2disp(self):
        return np.dot(self.get_imu2rect(), self.get_rect2disp())

    def get_velo2disp(self):
        return np.dot(self.get_velo2rect(), self.get_rect2disp())

    def get_disp2rect(self):
        return np.linalg.inv(self.get_rect2disp())

    def get_disp2imu(self):
        return np.linalg.inv(self.get_imu2disp())

    def rect2disp(self, points):
        return homogeneous_transform(points, self.get_rect2disp())

    def disp2rect(self, xyd):
        return homogeneous_transform(xyd, self.get_disp2rect())

    def velo2rect(self, points):
        return homogeneous_transform(points, self.get_velo2rect())

    def velo2disp(self, points):
        return homogeneous_transform(points, self.get_velo2disp())

    def imu2rect(self, points):
        return homogeneous_transform(points, self.get_imu2rect())

    def rect2imu(self, points):
        return homogeneous_transform(points, self.get_rect2imu())

    def filter_disps(self, xyd, max_disp=255, return_mask=False):
        return filter_disps(
            xyd, image_shape, max_disp=max_disp, return_mask=return_mask)

    def get_proj(self, cam_idx):
        P = self.data['P%d' % cam_idx].reshape(3, 4)
        return P

    def get_velo2depth(self, cam_idx):
        return np.dot(self.get_proj(cam_idx), self.get_velo2rect())

    def velo2depth(self, points, cam_idx):
        return homogeneous_transform(points, self.get_velo2depth(cam_idx),
                                     keep_last=True)

    def velo2img(self, points, cam_idx):
        return homogeneous_transform(points, self.get_velo2depth(cam_idx))

    def filter_depths(self, xyd, image_shape):
        return filter_depths(xyd, image_shape)


text:
###display temperature on sense###
from sense_hat import SenseHat
sense=SenseHat()
from pixelcolours import *
import snumbers as snum
import time
import numpy as np
import os

totaltime=float(input("Enter total run time in seconds: ")) #total time for the experiment
filename=input("Enter a file name: ")
printQ=" "
while printQ not in ["y","n"]:
    printQ=input("Would you like to enable hat display? y or n: ")
    printQ=printQ.lower()

filepath="/home/pi/"+filename+".txt"

tsign=[
    e,e,e,e,w,e,e,e,
    w,w,w,e,e,e,e,e,
    e,w,e,e,e,w,w,w,
    e,w,e,e,e,w,e,e,
    e,w,e,e,e,w,e,e,
    e,w,e,e,e,w,w,w,
    e,e,e,e,e,e,e,e,
    e,e,e,e,e,e,e,e
    ]
initialtime=time.time()
np_Pitemp=np.array([])
np_temp=np.array([])
np_ptemp=np.array([])
np_humidity=np.array([])
np_pressure=np.array([])
np_time=np.array([])

def getpitemp():
    """Gets the temperature of the Pi chip"""
    resline=os.popen("vcgencmd measure_temp").readline()
    res=resline.replace("temp=","").replace("'C\n","")
    return(float(res))

def getenvirodata():
    """Gets temp from humidity and pressure, and humidity and pressure"""
    ht=sense.get_temperature()
    pt=sense.get_temperature_from_pressure()
    hu=sense.get_humidity()
    pr=sense.get_pressure() #in millibar
    if pr<1: #sometimes first attempt =0
        pr=sense.get_pressure() #in millibar
        pt=sense.get_temperature_from_pressure()

    return ht,pt,hu,pr

def printenvirodata(Tem,Hum,Pres):
    if Tem!=None:
        if Tem>=25:
            sense.show_message("T= "+str(round(Tem,1)),scroll_speed=0.07,text_colour=r)
        if 19<Tem<25:
            sense.show_message("T= "+str(round(Tem,1)),scroll_speed=0.07,text_colour=g)
        if Tem<=19:
            sense.show_message("T= "+str(round(Tem,1)),scroll_speed=0.07,text_colour=b)
    if Hum!=None:
        sense.show_message("H= "+str(round(Hum,1)),scroll_speed=0.07)
    if Pres!=None:
        sense.show_message("P= "+str(round(Pres,1)),scroll_speed=0.07)
    return
i=3
while(time.time()<(initialtime+totaltime+5)):
    starttime=time.time()

    htemp,ptemp,humidity,pressure=getenvirodata()
    chiptemp=getpitemp()

    np_Pitemp=np.append(np_Pitemp,chiptemp)
    np_time=np.append(np_time,int(time.time()-initialtime))
    np_temp=np.append(np_temp,htemp)
    np_ptemp=np.append(np_ptemp,ptemp)
    np_humidity=np.append(np_humidity,humidity)
    np_pressure=np.append(np_pressure,pressure)
    avgPiT=(np_Pitemp.sum()/len(np_Pitemp))
    TcalcH=np_temp+(0.201*avgPiT)-13.965
    TcalcP=np_ptemp+(0.145*avgPiT)-9.812
    Tcalc=((TcalcH+TcalcP)/2)-1 #i don't know why it seems to be 1 out

    if i==3 and printQ=="y": #every 4th measurement it prints
        printenvirodata(Tcalc[-1],humidity,pressure)
        i=0
    else:
        i+=1

    print("t="+str(np_time[-1])+
          " Th="+str(round(TcalcH[-1],1))+
          " Tp="+str(round(TcalcP[-1],1))+
          " Tcalc="+str(round(Tcalc[-1],1))+
          " H="+str(round(humidity,1))+
          " P="+str(round(pressure,0)))

    medtime=time.time()
    waittime=30-(medtime-starttime)
    time.sleep(waittime)

finaldata=np.array([np_time,Tcalc,np_Pitemp,np_temp,np_ptemp,np_humidity,np_pressure])

np.savetxt(filepath,finaldata,delimiter=",")

print("measurement finished")

max_stars_repo_name: scotwheeler/piHaT

max_issues_repo_name: scotwheeler/piHaT

max_forks_repo_name: scotwheeler/piHaT


text:
# Identify location
import socket
location = socket.gethostname()
if location == 'saruman':
     dropbox = '/home/herdata/spx7cjc/Dropbox/'

# Import smorgasbord
import os
import sys
import numpy as np
import scipy.ndimage
import multiprocessing as mp
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
import pdb
import gc
import time
import datetime
import astropy.io.fits
import astropy.wcs
import astropy.io.votable
import astropy.modeling.models
import astropy.modeling.fitting
import aplpy
import ChrisFuncs
plt.ioff()

# Define function to provide appropriate Padding for FITS header entires
def Padding(entry):
    return ( ' ' * np.max([ 0, 19-len( str(entry) ) ]) ) + ' / '

# Define function to finalise Spitzer image of a given source in a given band
def Spitzer_Generator(name, ra, dec, in_dir, band_dict, thumbnails):
    band_short = band_dict['band_short']
    band_long = band_dict['band_long']
    wavelength = band_dict['wavelength']
    pix_size = band_dict['pix_size']
    colourmap = band_dict['colourmap']

    # Check if source has coverage in this band
    if not os.path.exists(in_dir+name+'_Spitzer_'+band_short+'.fits.gz'):
        print name+' has no coverage at '+wavelength
    else:

        # Read in image map
        image_in_fitsdata = astropy.io.fits.open(in_dir+name+'_Spitzer_'+band_short+'.fits.gz')
        image_in_cutout = image_in_fitsdata[0].data
        image_in_header = image_in_fitsdata[0].header
        image_in_wcs = astropy.wcs.WCS(image_in_header)
        image_out_cutout = image_in_cutout.copy()

        # Read in error map
        error_in_fitsdata = astropy.io.fits.open(in_dir+name+'_Spitzer_'+band_short+'_Error.fits.gz')
        error_in_cutout = error_in_fitsdata[0].data
        error_out_cutout = error_in_cutout.copy()

        # Set non-coverage pixels to be nan
        image_out_cutout[ np.where(image_out_cutout==0) ] = np.NaN
        error_out_cutout[ np.where(image_out_cutout==0) ] = np.NaN

        # Double-recheck that there is actual coverage at region of source
        i_target, j_target = np.round( image_in_wcs.wcs_world2pix( [[float(ra),float(dec)]], 0 )[0] ).astype(int).tolist()
        slice_target = image_out_cutout[ i_target-50:i_target+50, j_target-5:j_target+5 ]
        if np.where( np.isnan(slice_target)==False )[0].shape[0]==0:
            print name+' has no coverage at '+wavelength
        else:

            # Convert pixels from MJy/sr to Jy/sr
            image_out_cutout *= 1E6
            error_out_cutout *= 1E6

            # Calculate sr/pixel
            sr_per_sqarcsec = 2.3504E-11
            sr_per_pixels = sr_per_sqarcsec * pix_size**2

            # Convert pixels from Jy/sr to Jy/pix
            image_out_cutout *= sr_per_pixels
            error_out_cutout *= sr_per_pixels

            # Subtract minimum pixel value from entire map, to prevent issues with negative pixels later on
            image_out_cutout -= np.nanmin(image_out_cutout)

            # Only include coverage region that contains target source (ie, remove other, discontiguous regions)
            cont_binary = np.zeros([(image_out_cutout.shape)[0], (image_out_cutout.shape)[1]])
            cont_binary[ np.where( np.isnan(image_out_cutout)==False ) ] = 1
            cont_label = np.zeros([(image_out_cutout.shape)[0], (image_out_cutout.shape)[1]])
            cont_structure = np.array([[1,1,1], [1,1,1], [1,1,1]])
            scipy.ndimage.measurements.label(cont_binary, structure=cont_structure, output=cont_label)
            target_label = cont_label[i_target, j_target]
            image_out_cutout[ np.where( cont_label!=target_label ) ] = np.NaN
            error_out_cutout[ np.where( cont_label!=target_label ) ] = np.NaN

            # Create standard header
            out_header = astropy.io.fits.Header()
            date = datetime.datetime.now().isoformat()

            # Populate standard header entries
            out_header.set('TARGET', name , 'Target source of this map')
            out_header.set('COORDSYS', 'IRCS', 'Coordinate reference frame for the RA and DEC')
            out_header.set('SIGUNIT', 'Jy/pix', 'Unit of the map')
            out_header.set('TELESCOP', 'Spitzer', 'Telescope that made this observation')
            out_header.set('INSTMNT', band_dict['instrument'], 'Instrument used for this observation')
            out_header.set('FILTER', band_long, 'Filter used for this observation')
            out_header.set('WVLNGTH', wavelength, 'Wavelength of observation')
            out_header.set('DATABASE', 'Map produced using data acquired from the Spitzer Heritage Archive (SHA)', 'Database from which the reduced data was acquired')
            out_header.set('MAPDATE', date, 'Date this map was produced by CJR Clark (Cardiff) from the reduced data')

            # Construct WCS for output header
            out_wcs = astropy.wcs.WCS(naxis=2)
            out_wcs.wcs.crpix = [image_in_header['CRPIX1'], image_in_header['CRPIX2']]
            out_wcs.wcs.cdelt = [image_in_header['CDELT1'], image_in_header['CDELT2']]
            out_wcs.wcs.crval = [image_in_header['CRVAL1'], image_in_header['CRVAL2']]
            out_wcs.wcs.ctype = [image_in_header['CTYPE1'], image_in_header['CTYPE2']]
            out_wcs_header = out_wcs.to_header()
            for card in out_wcs_header.cards:
                out_header.append(card)

            # Create image-specific header and HDU
            image_out_header = out_header.copy()
            image_out_header.set('EXTNAME', 'image')
            image_out_hdu = astropy.io.fits.PrimaryHDU(data=image_out_cutout, header=image_out_header)

            # Create error-specific header, with appropriate WCS
            error_out_header = out_header.copy()
            error_out_header.set('EXTNAME', 'error')
            error_out_hdu = astropy.io.fits.ImageHDU(data=error_out_cutout, header=error_out_header)

            # Create hdulist and save to file
            image_cutout_hdulist = astropy.io.fits.HDUList([image_out_hdu, error_out_hdu])
            image_cutout_hdulist.writeto(our_dir+name+'_Spitzer_'+band_short+'.fits.gz', clobber=True)

            # Make thumbnail image of cutout
            if thumbnails:
                thumb_out = aplpy.FITSFigure(our_dir+name+'_Spitzer_'+band_short+'.fits.gz')
                thumb_out.show_colorscale(cmap=colourmap, stretch='linear')#, pmin=7.5, pmax=99.5)
                thumb_out.axis_labels.hide()
                thumb_out.tick_labels.hide()
                thumb_out.ticks.hide()
                thumb_out.set_nan_color('black')
                thumb_out.show_markers(np.array([float(ra)]), np.array([float(dec)]), marker='+', s=200, lw=2.0, edgecolor='#01DF3A', facecolor='#01DF3A' )
                thumb_out.save(our_dir+'Thumbnails/'+name+'_Spitzer_'+band_short+'.png')
                thumb_out.close()

            # Close file, and clean memory
            image_in_fitsdata.close()
            error_in_fitsdata.close()
            gc.collect()

# Commence main task
if __name__ == "__main__":

    # State if thumbnail images are wanted
    thumbnails = False

    # Decide what intrument to work for
    instrument = 'IRAC'

    # State band information
    if instrument=='IRAC':
        bands_dict = {'3.6um':{'instrument':'IRAC','band_short':'3.6','band_long':'IRAC 1','wavelength':'3.6um','pix_size':0.6,'colourmap':'pink'},
                      '4.5um':{'instrument':'IRAC','band_short':'4.5','band_long':'IRAC 2','wavelength':'4.5um','pix_size':0.6,'colourmap':'pink'},
                      '5.8um':{'instrument':'IRAC','band_short':'5.8','band_long':'IRAC 3','wavelength':'5.8um','pix_size':0.6,'colourmap':'pink'},
                      '8.0um':{'instrument':'IRAC','band_short':'8.0','band_long':'IRAC 4','wavelength':'8.0um','pix_size':0.6,'colourmap':'pink'}}
    elif instrument=='MIPS':
        bands_dict = {'24um':{'instrument':'MIPS','band_short':'24','band_long':'MIPS 1','wavelength':'24um','pix_size':2.45,'colourmap':'gist_heat'},
                      '70um':{'instrument':'MIPS','band_short':'70','band_long':'MIPS 2','wavelength':'70um','pix_size':4.0,'colourmap':'gist_heat'},
                      '160um':{'instrument':'MIPS','band_short':'160','band_long':'MIPS 3','wavelength':'160um','pix_size':8.0,'colourmap':'gist_heat'}}

    # Read in source catalogue
    ness_cat = np.genfromtxt(dropbox+'Work/Tables/NESS/NESS_Sample.csv', delimiter=',', names=True, dtype=None)
    name_list = ness_cat['name']
    ra_list = ness_cat['ra']
    dec_list = ness_cat['dec']

    # Give paths
    in_dir = '/home/sarumandata2/spx7cjc/NESS/Ancillary_Data/Spitzer/Mosaics_'+instrument+'/'
    our_dir = '/home/sarumandata2/spx7cjc/NESS/Ancillary_Data/Spitzer/Cutouts/'

    # Record time taken
    time_total = 0.0
    source_counter = 0.0
    source_total = name_list.shape[0]
    time_source_list = []

    # Loop over each source
    for i in range(0, name_list.shape[0]):#np.where(name_list=='NGC3031')[0]:
        source_counter += 1.0
        name = name_list[i].replace(' ','_')
        ra = ra_list[i]
        dec = dec_list[i]
        name.replace('[','')
        name.replace(']','')

        # In parallel, process each band
        print 'Processing '+str(name)
        pool = mp.Pool(processes=5)
        for key in bands_dict.keys():
            band_dict = bands_dict[key]
            #pool.apply_async( Spitzer_Generator, args=(name, in_dir, band_dict, thumbnails,) )#
            Spitzer_Generator(name, ra, dec, in_dir, band_dict, thumbnails)
        pool.close()
        pool.join()

    # Jubilate
    print 'All done!'

max_stars_repo_name: drtobybrown/AncillaryDataButton

max_issues_repo_name: drtobybrown/AncillaryDataButton

max_forks_repo_name: drtobybrown/AncillaryDataButton

reason: import numpy,import scipy,import astropy


text:
import numpy as np
from spn.structure.StatisticalTypes import MetaType
from spn.structure.leaves.histogram.Histograms import Histogram, create_histogram_leaf

# class SPMNLeaf(Histogram):
#
#     def __init__(self, scope):

def create_spmn_leaf(data, ds_context, scope):
    assert len(scope) == 1, "scope of univariate histogram for more than one variable?"
    assert data.shape[1] == 1, "data has more than one feature?"

    # data = data[~np.isnan(data)]

    idx = scope[0]
    meta_type = ds_context.meta_types[idx]

    if meta_type == MetaType.UTILITY:
        hist = create_histogram_leaf(data, ds_context, scope)
        return Utility(hist.breaks, hist.densities, hist.bin_repr_points, scope=idx)
    else:
        return create_histogram_leaf(data, ds_context, scope)

class Utility(Histogram):

    def __init__(self, breaks, densities, bin_repr_points, scope=None, type_=None, meta_type=MetaType.UTILITY):

        # has same member variables as histogram
        Histogram.__init__(self, breaks, densities, bin_repr_points, scope,
                           type_=None, meta_type=MetaType.UTILITY)


text:
# copyright 2003-2013 LOGILAB S.A. (Paris, FRANCE), all rights reserved.
# contact http://www.logilab.fr/ -- mailto:contact@logilab.fr
#
# This file is part of logilab-astng.
#
# logilab-astng is free software: you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 2.1 of the License, or (at your
# option) any later version.
#
# logilab-astng is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
# for more details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with logilab-astng. If not, see <http://www.gnu.org/licenses/>.

from logilab.common.testlib import unittest_main, TestCase

from logilab.astng import ResolveError, MANAGER, Instance, nodes, YES, InferenceError
from logilab.astng.builder import ASTNGBuilder
from logilab.astng.raw_building import build_module
from logilab.astng.manager import ASTNGManager

import sys
from os.path import join, abspath, dirname

class NonRegressionTC(TestCase):

    def setUp(self):
        sys.path.insert(0, join(dirname(abspath(__file__)), 'regrtest_data'))

    def tearDown(self):
        sys.path.pop(0)

    def brainless_manager(self):
        manager = ASTNGManager()
        # avoid caching into the ASTNGManager borg since we get problems
        # with other tests :
        manager.__dict__ = {}
        manager.astng_cache = {}
        manager._mod_file_cache = {}
        manager.transformers = {}
        return manager

    def test_module_path(self):
        man = self.brainless_manager()
        mod = man.astng_from_module_name('package.import_package_subpackage_module')
        package = next(mod.igetattr('package'))
        self.assertEqual(package.name, 'package')
        subpackage = next(package.igetattr('subpackage'))
        self.assertIsInstance(subpackage, nodes.Module)
        self.assertTrue(subpackage.package)
        self.assertEqual(subpackage.name, 'package.subpackage')
        module = next(subpackage.igetattr('module'))
        self.assertEqual(module.name, 'package.subpackage.module')

    def test_package_sidepackage(self):
        manager = self.brainless_manager()
        assert 'package.sidepackage' not in MANAGER.astng_cache
        package = manager.astng_from_module_name('absimp')
        self.assertIsInstance(package, nodes.Module)
        self.assertTrue(package.package)
        subpackage = next(package.getattr('sidepackage')[0].infer())
        self.assertIsInstance(subpackage, nodes.Module)
        self.assertTrue(subpackage.package)
        self.assertEqual(subpackage.name, 'absimp.sidepackage')

    def test_living_property(self):
        builder = ASTNGBuilder()
        builder._done = {}
        builder._module = sys.modules[__name__]
        builder.object_build(build_module('module_name', ''), Whatever)

    def test_new_style_class_detection(self):
        try:
            import pygtk
        except ImportError:
            self.skipTest('test skipped: pygtk is not available')
        # XXX may fail on some pygtk version, because objects in
        # gobject._gobject have __module__ set to gobject :(
        builder = ASTNGBuilder()
        data = """
import pygtk
pygtk.require("2.6")
import gobject

class A(gobject.GObject):
    pass
"""
        astng = builder.string_build(data, __name__, __file__)
        a = astng['A']
        self.assertTrue(a.newstyle)

    def test_pylint_config_attr(self):
        try:
            from pylint import lint
        except ImportError:
            self.skipTest('pylint not available')
        mod = MANAGER.astng_from_module_name('pylint.lint')
        pylinter = mod['PyLinter']
        expect = ['OptionsManagerMixIn', 'object', 'MessagesHandlerMixIn',
                  'ReportsHandlerMixIn', 'BaseRawChecker', 'BaseChecker',
                  'OptionsProviderMixIn', 'ASTWalker']
        self.assertListEqual([c.name for c in pylinter.ancestors()],
                             expect)
        self.assertTrue(list(Instance(pylinter).getattr('config')))
        infered = list(Instance(pylinter).igetattr('config'))
        self.assertEqual(len(infered), 1)
        self.assertEqual(infered[0].root().name, 'optparse')
        self.assertEqual(infered[0].name, 'Values')

    def test_numpy_crash(self):
        """test don't crash on numpy"""
        #a crash occured somewhere in the past, and an
        # InferenceError instead of a crash was better, but now we even infer!
        try:
            import numpy
        except ImportError:
            self.skipTest('test skipped: numpy is not available')
        builder = ASTNGBuilder()
        data = """
from numpy import multiply

multiply(1, 2, 3)
"""
        astng = builder.string_build(data, __name__, __file__)
        callfunc = astng.body[1].value.func
        infered = callfunc.infered()
        self.assertEqual(len(infered), 1)
        self.assertIsInstance(infered[0], Instance)

class Whatever(object):
    a = property(lambda x: x, lambda x: x)

if __name__ == '__main__':
    unittest_main()

max_stars_repo_name: tlksio/tlksio

max_issues_repo_name: tlksio/tlksio

max_forks_repo_name: tlksio/tlksio

reason: import numpy,from numpy


text:
import numpy as np
import joblib

# --------------------------------------------------------
# ---------------- Save and read collection dumps --------
# --------------------------------------------------------
def read_pkl_spectra(file_path, read_rows=None, read_cols=None):
    data_read = joblib.load(file_path)
    if read_cols is not None:
        data_read = data_read[:, read_cols]
    if read_rows is not None:
        data_read = data_read[read_rows, :]
    return data_read

def save_pkl_spectra(data, file_path):
    joblib.dump(data, file_path)

# --------------------------------------------------------
# ---------------- Spectra collection class --------------
# --------------------------------------------------------
class CollectionParameters:
    def __init__(self, filename):
        self.filename_full = filename
        self.filename_split = '.'.join(filename.split('.')[:-1]).split('_')

    def __get_str_pos__(self, search_str):
        try:
            idx = self.filename_split.index(search_str)
        except ValueError:
            idx = None
        return idx

    def __get_value__(self, search_str):
        idx_value = self.__get_str_pos__(search_str)
        if idx_value is not None:
            value_str = self.filename_split[idx_value + 1]
            try:
                return int(value_str)
            except ValueError:
                return float(value_str)
        else:
            return None

    def get_ccd(self):
        return self.filename_split[2][3]  # the fourth character is the number of the ccd

    def get_wvl_range(self):
        wvl_start = float(self.filename_split[3])
        wvl_end = float(self.filename_split[4])
        return wvl_start, wvl_end

    def get_wvl_values(self):
        min, max = self.get_wvl_range()
        return np.arange(min, max, self.get_wvl_step())

    def get_interp_method(self):
        return self.filename_split[5]

    def get_wvl_step(self):
        return self.__get_value__('wvlstep')

    def get_snr_limit(self):
        return self.__get_value__('snr')

    def get_teff_step(self):
        return self.__get_value__('teff')

    def get_logg_step(self):
        return self.__get_value__('logg')

    def get_feh_step(self):
        return self.__get_value__('feh')

max_stars_repo_name: kcotar/Priprava-GALAH-podatkov-studenti

max_issues_repo_name: kcotar/Priprava-GALAH-podatkov-studenti

max_forks_repo_name: kcotar/Priprava-GALAH-podatkov-studenti


text:
from typing import Optional, Union

import numpy as np
import torch
from torch.distributions.categorical import Categorical
from scipy.sparse import csr_matrix

from ..models import BaseEmbeddingModel

class BaseSampler:
    def __init__(
        self,
        train_set: torch.Tensor,
        n_user: Optional[int] = None,
        n_item: Optional[int] = None,
        pos_weight: Optional[np.ndarray] = None,
        neg_weight: Union[np.ndarray, BaseEmbeddingModel] = None,
        device: Optional[torch.device] = None,
        batch_size: int = 256,
        n_neg_samples: int = 10,
        strict_negative: bool = False,
        neutral: torch.Tensor = torch.Tensor([]),
    ):
        """Class of Base Sampler for get positive and negative batch.
        Args:
            train_set (torch.Tensor): training interaction data which columns are [user_id, item_id]
            n_user (Optional[int], optional): A number of user considered. Defaults to None.
            n_item (Optional[int], optional): A number of item considered. Defaults to None.
            pos_weight (Optional[np.ndarray], optional): Sampling weight for positive pair. Defaults to None.
            neg_weight (Optional[np.ndarray], optional): Sampling weight for negative item. Defaults to None.
            device (Optional[torch.device], optional): Device name. Defaults to None.
            batch_size (int, optional): Length of mini-batch. Defaults to 256.
            n_neg_samples (int, optional): A number of negative samples. Defaults to 10.
            strict_negative (bool, optional): If removing positive items from negative samples or not. Defaults to False.

        Raises:
            NotImplementedError: [description]
            NotImplementedError: [description]
        """
        # set positive pairs
        self.train_set = train_set

        # set some hyperparameters
        self.n_neg_samples = n_neg_samples
        self.batch_size = batch_size
        self.strict_negative = strict_negative
        self.two_stage = False

        if n_user is None:
            self.n_user = np.unique(train_set[:, 0].cpu()).shape[0]
        else:
            self.n_user = n_user

        if n_user is None:
            self.n_item = np.unique(train_set[:, 1].cpu()).shape[0]
        else:
            self.n_item = n_item

        # set flag matrix whose element indicates the pair is not negative.
        train_set_cpu = train_set.cpu()
        neutral_cpu = neutral.cpu()
        not_negative = torch.cat([train_set_cpu, neutral_cpu])
        self.not_negative_flag = csr_matrix(
            (np.ones(not_negative.shape[0]), (not_negative[:, 0], not_negative[:, 1])),
            [n_user, n_item],
        )
        self.not_negative_flag.sum_duplicates()
        self.not_negative_flag.data[:] = 1

        # device
        if device is None:
            self.device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
        else:
            self.device = device

        # set pos weight
        if pos_weight is not None:  # weighted
            if len(pos_weight) == len(train_set):
                pos_weight_pair = pos_weight

            elif len(pos_weight) == self.n_item:
                pos_weight_pair = pos_weight[train_set[:, 1].cpu()]

            elif len(pos_weight) == self.n_user:
                pos_weight_pair = pos_weight[train_set[:, 0].cpu()]

            else:
                raise NotImplementedError

        else:  # uniform
            pos_weight_pair = torch.ones(train_set.shape[0])

        self.pos_weight_pair = torch.Tensor(pos_weight_pair).to(self.device)
        self.pos_sampler = Categorical(probs=self.pos_weight_pair)

        # set neg weight
        if neg_weight is None:  # uniform
            self.negative_weighted_by_model = False
            self.neg_item_weight = torch.ones(self.n_item).to(self.device)

        elif isinstance(neg_weight, BaseEmbeddingModel):  # user-item weighted
            self.negative_weighted_by_model = True
            self.neg_weight_model = neg_weight

        elif len(neg_weight) == self.n_item:  # item weighted
            self.negative_weighted_by_model = False
            self.neg_item_weight = torch.Tensor(neg_weight).to(self.device)

        else:
            raise NotImplementedError

    def get_pos_batch(self) -> torch.Tensor:
        """Method for positive sampling.

        Returns:
            torch.Tensor: positive batch.
        """
        batch_indices = self.pos_sampler.sample([self.batch_size])
        batch = self.train_set[batch_indices]
        return batch

    def get_neg_batch(self, users: torch.Tensor) -> torch.Tensor:
        """Method of negative sampling

        Args:
            users (torch.Tensor): indices of users in pos pairs.

        Returns:
            torch.Tensor: negative samples.
        """

        if self.negative_weighted_by_model and self.strict_negative:
            flag = torch.Tensor(self.not_negative_flag[users.to("cpu")].A)
            mask = 1 - flag.to(self.device)
            weight = self.neg_weight_model.get_item_weight(users)
            weight *= mask
            neg_sampler = Categorical(probs=weight)
            neg_samples = neg_sampler.sample([self.n_neg_samples]).T

        elif self.negative_weighted_by_model and not self.strict_negative:
            weight = self.neg_weight_model.get_item_weight(users)
            neg_sampler = Categorical(probs=weight)
            neg_samples = neg_sampler.sample([self.n_neg_samples]).T

        elif not self.negative_weighted_by_model and self.strict_negative:
            flag = torch.Tensor(self.not_negative_flag[users.to("cpu")].A)
            mask = 1 - flag.to(self.device)
            weight = mask * self.neg_item_weight
            neg_sampler = Categorical(probs=weight)
            neg_samples = neg_sampler.sample([self.n_neg_samples]).T

        else:
            neg_sampler = Categorical(probs=self.neg_item_weight)
            neg_samples = neg_sampler.sample([self.batch_size, self.n_neg_samples])

        return neg_samples

reason: import numpy,from scipy


text:
import torch
import torch.nn as nn
import os
import numpy as np
from scipy.interpolate import interp1d
from glob import glob
from lib.utils import as_numpy
import csv
import pickle

# temp til better solution
fc_feature = None
pred_feature = None

class Feature:
    def __init__(self, fc_feature=fc_feature, pred_feature=pred_feature, device=0, num_classes=150, save_path=None, threshold=0.1):
        self.device = torch.device('cuda:{}'.format(device))
        self.num_classes = num_classes
        self.save_path = save_path
        self.threshhold = threshold

        self.fc_feat = fc_feature.to(self.device).type(
                torch.cuda.HalfTensor)
        self.pred_feat = pred_feature.to(self.device).type(
                torch.cuda.HalfTensor)

        self.fc_feat_size = fc_feature.shape
        self.fc_dim = fc_feature.shape[-1]

        self.pred_feat_size = pred_feature.shape

    def compute_image_rep(self):
        interp_feat = self.interpolate_feature(
            self.pred_feat, (self.fc_feat_size[1], self.fc_feat_size[2]))
        accepted_classes = self.get_accepted_classes(interp_feat)
        computed_features = self.pool_feature(accepted_classes, interp_feat)

    def pool_feature(self, accepted_classes, interp_feat):
        for idx, acc in enumerate(accepted_classes):
            locations = np.argwhere(
                interp_feat.numpy() == accepted_classes[accepted_classes.index(acc)])

            img_rep = torch.zeros((len(accepted_classes), self.fc_dim), device=self.device).type(
                torch.cuda.HalfTensor)
            for loc in locations:
                img_rep = img_rep[idx].add(self.fc_feat[0][loc[0]][loc[1]])
            img_rep[idx] = torch.div(img_rep[idx], len(locations))
        return img_rep

    def get_accepted_classes(self, feature):
        total_pixels = self.fc_feat_size[1] * self.fc_feat_size[2]

        accepted_classes = []
        for i in range(self.num_classes):
            if (len(np.argwhere(feature.numpy() == i+1)) / total_pixels) > self.threshhold:
                accepted_classes.append(i+1)
        return accepted_classes

    def interpolate_feature(self, feature, size):
        _, feat_150_pred = torch.max(feature, dim=1)
        feat_150_pred = as_numpy(feat_150_pred.squeeze(0))

        interp_feature = nn.functional.interpolate(torch.from_numpy(feat_150_pred).unsqueeze(
            0).unsqueeze(0).type(torch.FloatTensor), size=size, mode='bilinear', align_corners=False)
        interp_feature = interp_feature.squeeze().squeeze().type(torch.int64)

        return interp_feature

    # TODO: SAVE THE FEATURES
    def save_pooled_features(self, image_rep):
        pass


text:
"""
A simple version of Distributed Proximal Policy Optimization (DPPO) using distributed TensorFlow.

Based on:
1. Emergence of Locomotion Behaviours in Rich Environments (Google Deepmind): [https://arxiv.org/abs/1707.02286]
2. Proximal Policy Optimization Algorithms (OpenAI): [https://arxiv.org/abs/1707.06347]
3. Generalized Advantage Estimation [https://arxiv.org/abs/1506.02438]

Thanks to OpenAI and morvanzhou for their examples
"""
import sys
sys.path.append("/home/anjum/PycharmProjects/Aivot/")
import argparse
import tensorflow as tf
import numpy as np
import gym
import os
import scipy.signal
from time import time
from gym import wrappers
from tensorflow.python.training.summary_io import SummaryWriterCache
from utils import RunningStats, discount, add_histogram
OUTPUT_RESULTS_DIR = "./"

EP_MAX = 10000
GAMMA = 0.99
LAMBDA = 0.95
ENTROPY_BETA = 0.01  # 0.01 for discrete, 0.0 for continuous
LR = 0.0001
BATCH = 128  # 128 for discrete, 8192 for continuous
MINIBATCH = 32
EPOCHS = 10
EPSILON = 0.1
VF_COEFF = 1.0
L2_REG = 0.001
SIGMA_FLOOR = 0.0
PORT_OFFSET = 20  # number to offset the port from 2222. Useful for multiple runs

class PPO(object):
    def __init__(self, environment, wid, greyscale=True):
        if len(environment.action_space.shape) > 0:
            self.discrete = False
            self.s_dim, self.a_dim = environment.observation_space.shape, environment.action_space.shape[0]
            self.a_bound = (environment.action_space.high - environment.action_space.low) / 2
            self.actions = tf.placeholder(tf.float32, [None, self.a_dim], 'action')
        else:
            self.discrete = True
            self.s_dim, self.a_dim = environment.observation_space.shape, environment.action_space.n
            self.actions = tf.placeholder(tf.int32, [None, 1], 'action')
        self.cnn = len(self.s_dim) == 3
        self.greyscale = greyscale  # If not greyscale and using RGB, make sure to divide the images by 255
        self.wid = wid
        is_chief = wid == 0

        self.state = tf.placeholder(tf.float32, [None] + list(self.s_dim), 'state')
        self.advantage = tf.placeholder(tf.float32, [None, 1], 'advantage')
        self.rewards = tf.placeholder(tf.float32, [None, 1], 'discounted_r')

        self.dataset = tf.data.Dataset.from_tensor_slices({"state": self.state, "actions": self.actions,
                                                           "rewards": self.rewards, "advantage": self.advantage})
        self.dataset = self.dataset.shuffle(buffer_size=10000)
        self.dataset = self.dataset.batch(MINIBATCH)
        self.dataset = self.dataset.cache()
        self.dataset = self.dataset.repeat(EPOCHS)
        self.iterator = self.dataset.make_initializable_iterator()
        batch = self.iterator.get_next()

        pi_old, pi_old_params = self._build_anet(batch["state"], 'oldpi')
        pi, pi_params = self._build_anet(batch["state"], 'pi')
        pi_eval, _ = self._build_anet(self.state, 'pi', reuse=True)

        vf_old, vf_old_params = self._build_cnet(batch["state"], "oldvf")
        self.v, vf_params = self._build_cnet(batch["state"], "vf")
        self.vf_eval, _ = self._build_cnet(self.state, 'vf', reuse=True)

        self.sample_op = tf.squeeze(pi_eval.sample(1), axis=0, name="sample_action")
        self.eval_action = pi_eval.mode()  # Used mode for discrete case. Mode should equal mean in continuous
        self.global_step = tf.train.get_or_create_global_step()

        with tf.variable_scope('loss'):
            epsilon_decay = tf.train.polynomial_decay(EPSILON, self.global_step, 1e5, 0.01, power=0.0)

            with tf.variable_scope('policy'):
                # Use floor functions for the probabilities to prevent NaNs when prob = 0
                ratio = tf.maximum(pi.prob(batch["actions"]), 1e-6) / tf.maximum(pi_old.prob(batch["actions"]), 1e-6)
                ratio = tf.clip_by_value(ratio, 0, 10)
                surr1 = batch["advantage"] * ratio
                surr2 = batch["advantage"] * tf.clip_by_value(ratio, 1 - epsilon_decay, 1 + epsilon_decay)
                loss_pi = -tf.reduce_mean(tf.minimum(surr1, surr2))
                tf.summary.scalar("loss", loss_pi)

            with tf.variable_scope('value_function'):
                # Sometimes values clipping helps, sometimes just using raw residuals is better ¯\_(ツ)_/¯
                clipped_value_estimate = vf_old + tf.clip_by_value(self.v - vf_old, -epsilon_decay, epsilon_decay)
                loss_vf1 = tf.squared_difference(clipped_value_estimate, batch["rewards"])
                loss_vf2 = tf.squared_difference(self.v, batch["rewards"])
                loss_vf = tf.reduce_mean(tf.maximum(loss_vf1, loss_vf2)) * 0.5
                # loss_vf = tf.reduce_mean(tf.square(self.v - batch["rewards"])) * 0.5
                tf.summary.scalar("loss", loss_vf)

            with tf.variable_scope('entropy'):
                entropy = pi.entropy()
                pol_entpen = -ENTROPY_BETA * tf.reduce_mean(entropy)

            loss = loss_pi + loss_vf * VF_COEFF + pol_entpen
            tf.summary.scalar("total", loss)
            # tf.summary.scalar("epsilon", epsilon_decay)

        with tf.variable_scope('train'):
            self.global_step = tf.train.get_or_create_global_step()
            opt = tf.train.AdamOptimizer(LR)
            opt = tf.train.SyncReplicasOptimizer(opt, replicas_to_aggregate=N_AGG, total_num_replicas=N_WORKER)

            self.train_op = opt.minimize(loss, global_step=self.global_step, var_list=pi_params + vf_params)

            # grads, vs = zip(*opt.compute_gradients(loss, var_list=pi_params + vf_params))
            # Need to split the two networks so that clip_by_global_norm works properly
            # pi_grads, pi_vs = grads[:len(pi_params)], vs[:len(pi_params)]
            # vf_grads, vf_vs = grads[len(pi_params):], vs[len(pi_params):]
            # pi_grads, _ = tf.clip_by_global_norm(pi_grads, 10.0)
            # vf_grads, _ = tf.clip_by_global_norm(vf_grads, 10.0)
            # self.train_op = opt.apply_gradients(zip(pi_grads + vf_grads, pi_vs + vf_vs), global_step=self.global_step)

            self.sync_replicas_hook = opt.make_session_run_hook(is_chief)

        with tf.variable_scope('update_old'):
            self.update_pi_old_op = [oldp.assign(p) for p, oldp in zip(pi_params, pi_old_params)]
            self.update_vf_old_op = [oldp.assign(p) for p, oldp in zip(vf_params, vf_old_params)]

        tf.summary.scalar("value", tf.reduce_mean(self.v))
        tf.summary.scalar("policy_entropy", tf.reduce_mean(entropy))
        if not self.discrete:
            tf.summary.scalar("sigma", tf.reduce_mean(pi.stddev()))
        self.summarise = tf.summary.merge(tf.get_collection(tf.GraphKeys.SUMMARIES))

    def _build_anet(self, state_in, name, reuse=False):
        w_reg = tf.contrib.layers.l2_regularizer(L2_REG)

        with tf.variable_scope(name, reuse=reuse):
            if self.cnn:
                if self.greyscale:
                    state_in = tf.image.rgb_to_grayscale(state_in)
                conv1 = tf.layers.conv2d(inputs=state_in, filters=32, kernel_size=8, strides=4, activation=tf.nn.relu)
                conv2 = tf.layers.conv2d(inputs=conv1, filters=64, kernel_size=4, strides=2, activation=tf.nn.relu)
                conv3 = tf.layers.conv2d(inputs=conv2, filters=64, kernel_size=3, strides=1, activation=tf.nn.relu)
                state_in = tf.layers.flatten(conv3)

            layer_1 = tf.layers.dense(state_in, 400, tf.nn.relu, kernel_regularizer=w_reg, name="pi_l1")
            layer_2 = tf.layers.dense(layer_1, 400, tf.nn.relu, kernel_regularizer=w_reg, name="pi_l2")

            if self.discrete:
                a_logits = tf.layers.dense(layer_2, self.a_dim, kernel_regularizer=w_reg, name="pi_logits")
                dist = tf.distributions.Categorical(logits=a_logits)
            else:
                mu = tf.layers.dense(layer_2, self.a_dim, tf.nn.tanh, kernel_regularizer=w_reg, name="pi_mu")
                log_sigma = tf.get_variable(name="pi_sigma", shape=self.a_dim, initializer=tf.zeros_initializer())
                dist = tf.distributions.Normal(loc=mu * self.a_bound, scale=tf.maximum(tf.exp(log_sigma), SIGMA_FLOOR))
        params = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope=name)
        return dist, params

    def _build_cnet(self, state_in, name, reuse=False):
        w_reg = tf.contrib.layers.l2_regularizer(L2_REG)

        with tf.variable_scope(name, reuse=reuse):
            if self.cnn:
                if self.greyscale:
                    state_in = tf.image.rgb_to_grayscale(state_in)
                conv1 = tf.layers.conv2d(inputs=state_in, filters=32, kernel_size=8, strides=4, activation=tf.nn.relu)
                conv2 = tf.layers.conv2d(inputs=conv1, filters=64, kernel_size=4, strides=2, activation=tf.nn.relu)
                conv3 = tf.layers.conv2d(inputs=conv2, filters=64, kernel_size=3, strides=1, activation=tf.nn.relu)
                state_in = tf.layers.flatten(conv3)

            l1 = tf.layers.dense(state_in, 400, tf.nn.relu, kernel_regularizer=w_reg, name="vf_l1")
            l2 = tf.layers.dense(l1, 400, tf.nn.relu, kernel_regularizer=w_reg, name="vf_l2")
            vf = tf.layers.dense(l2, 1, kernel_regularizer=w_reg, name="vf_output")

        params = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope=name)
        return vf, params

    def update(self, s, a, r, adv, sess):
        start = time()
        e_time = []

        sess.run([self.update_pi_old_op, self.update_vf_old_op, self.iterator.initializer],
                 feed_dict={self.state: s, self.actions: a, self.rewards: r, self.advantage: adv})

        while True:
            try:
                e_start = time()
                summary, step, _ = sess.run([self.summarise, self.global_step, self.train_op])
                e_time.append(time() - e_start)
            except tf.errors.OutOfRangeError:
                break
        print("Worker_%i Trained in %.3fs at %.3fs/batch. Global step %i" % (self.wid, time() - start, np.mean(e_time), step))
        return summary

    def evaluate_state(self, state, sess, stochastic=True):
        if stochastic:
            action, value = sess.run([self.sample_op, self.vf_eval], {self.state: state[np.newaxis, :]})
        else:
            action, value = sess.run([self.eval_action, self.vf_eval], {self.state: state[np.newaxis, :]})
        return action[0], np.squeeze(value)

def start_parameter_server(pid, spec):
    cluster = tf.train.ClusterSpec(spec)
    server = tf.train.Server(cluster, job_name="ps", task_index=pid)
    print("Starting PS #{}".format(pid))
    server.join()

# class StopAtEpisodeHook(tf.train.SessionRunHook):
#     def __init__(self, max_episodes, worker):
#         self.ep_max = max_episodes
#         self.wid = worker
#
#     def after_run(self, run_context, run_values):
#         global episode
#         if episode >= self.ep_max and self.wid == 0:
#             print("Reached %i episodes" % episode)
#             run_context.request_stop()

class Worker(object):
    def __init__(self, wid, spec):
        self.wid = wid
        self.env = gym.make(ENVIRONMENT)

        print("Starting Worker #{}".format(wid))
        cluster = tf.train.ClusterSpec(spec)
        self.server = tf.train.Server(cluster, job_name="worker", task_index=wid)

        with tf.device(tf.train.replica_device_setter(worker_device="/job:worker/task:%d" % wid, cluster=cluster)):
            if self.wid == 0:
                self.env = wrappers.Monitor(self.env, os.path.join(SUMMARY_DIR, ENVIRONMENT), video_callable=None)
            self.ppo = PPO(self.env, self.wid)

    def work(self):
        hooks = [self.ppo.sync_replicas_hook]
        sess = tf.train.MonitoredTrainingSession(master=self.server.target, is_chief=(self.wid == 0),
                                                 checkpoint_dir=SUMMARY_DIR,
                                                 save_summaries_steps=None, save_summaries_secs=None, hooks=hooks)
        if self.wid == 0:
            writer = SummaryWriterCache.get(SUMMARY_DIR)

        t, episode, terminal = 0, 0, False
        buffer_s, buffer_a, buffer_r, buffer_v, buffer_terminal = [], [], [], [], []
        rolling_r = RunningStats()

        while not sess.should_stop() and not (episode > EP_MAX and self.wid == 0):

            s = self.env.reset()
            ep_r, ep_t, ep_a = 0, 0, []

            while True:
                a, v = self.ppo.evaluate_state(s, sess)

                # Update ppo
                if t == BATCH:  # or (terminal and t < BATCH):
                    # Normalise rewards
                    rewards = np.array(buffer_r)
                    rolling_r.update(rewards)
                    rewards = np.clip(rewards / rolling_r.std, -10, 10)

                    v_final = [v * (1 - terminal)]  # v = 0 if terminal, otherwise use the predicted v
                    values = np.array(buffer_v + v_final)
                    terminals = np.array(buffer_terminal + [terminal])

                    # Generalized Advantage Estimation - https://arxiv.org/abs/1506.02438
                    delta = rewards + GAMMA * values[1:] * (1 - terminals[1:]) - values[:-1]
                    advantage = discount(delta, GAMMA * LAMBDA, terminals)
                    returns = advantage + np.array(buffer_v)
                    advantage = (advantage - advantage.mean()) / np.maximum(advantage.std(), 1e-6)

                    bs, ba, br, badv = np.reshape(buffer_s, (t,) + self.ppo.s_dim), np.vstack(buffer_a), \
                                       np.vstack(returns), np.vstack(advantage)

                    graph_summary = self.ppo.update(bs, ba, br, badv, sess)
                    buffer_s, buffer_a, buffer_r, buffer_v, buffer_terminal = [], [], [], [], []
                    t = 0

                buffer_s.append(s)
                buffer_a.append(a)
                buffer_v.append(v)
                buffer_terminal.append(terminal)
                ep_a.append(a)

                if not self.ppo.discrete:
                    a = np.clip(a, self.env.action_space.low, self.env.action_space.high)
                s, r, terminal, _ = self.env.step(a)
                buffer_r.append(r)

                ep_r += r
                ep_t += 1
                t += 1

                if terminal:
                    # End of episode summary
                    print('Worker_%i' % self.wid,
                          '| Episode: %i' % episode, "| Reward: %.2f" % ep_r, '| Steps: %i' % ep_t)

                    if self.wid == 0:
                        worker_summary = tf.Summary()
                        worker_summary.value.add(tag="Reward", simple_value=ep_r)

                        # Create Action histograms for each dimension
                        actions = np.array(ep_a)
                        if self.ppo.discrete:
                            add_histogram(writer, "Action", actions, episode, bins=self.ppo.a_dim)
                        else:
                            for a in range(self.ppo.a_dim):
                                add_histogram(writer, "Action/Dim" + str(a), actions[:, a], episode)

                        try:
                            writer.add_summary(graph_summary, episode)
                        except NameError:
                            pass
                        writer.add_summary(worker_summary, episode)
                        writer.flush()

                    episode += 1
                    break

        self.env.close()
        print("Worker_%i finished" % self.wid)

def main(_):
    pass

if __name__ == '__main__':
    # Discrete environments
    # ENVIRONMENT = 'CartPole-v1'
    # ENVIRONMENT = 'MountainCar-v0'
    # ENVIRONMENT = 'LunarLander-v2'
    ENVIRONMENT = 'Pong-v0'

    # Continuous environments
    # ENVIRONMENT = 'Pendulum-v0'
    # ENVIRONMENT = 'MountainCarContinuous-v0'
    # ENVIRONMENT = 'LunarLanderContinuous-v2'
    # ENVIRONMENT = 'BipedalWalker-v2'
    # ENVIRONMENT = 'BipedalWalkerHardcore-v2'
    # ENVIRONMENT = 'CarRacing-v0'

    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
    tf.logging.set_verbosity(tf.logging.ERROR)
    parser = argparse.ArgumentParser()
    parser.add_argument('--job_name', action='store', dest='job_name', help='Either "ps" or "worker"')
    parser.add_argument('--task_index', action='store', dest='task_index', help='ID number of the job')
    parser.add_argument('--timestamp', action='store', dest='timestamp', help='Timestamp for output directory')
    parser.add_argument('--workers', action='store', dest='n_workers', help='Number of workers')
    parser.add_argument('--agg', action='store', dest='n_agg', help='Number of gradients to aggregate')
    parser.add_argument('--ps', action='store', dest='ps', help='Number of parameter servers')
    args = parser.parse_args()

    N_WORKER = int(args.n_workers)
    N_AGG = int(args.n_agg)
    PS = int(args.ps)
    TIMESTAMP = args.timestamp
    SUMMARY_DIR = os.path.join(OUTPUT_RESULTS_DIR, "DPPO", ENVIRONMENT, TIMESTAMP)

    # This could be defined in a parameter file
    if PS == 0:
        spec = {"worker": ["localhost:" + str(2222 + PS + i + PORT_OFFSET) for i in range(N_WORKER)]}
    else:
        spec = {"ps": ["localhost:" + str(2222 + i + PORT_OFFSET) for i in range(PS)],
                "worker": ["localhost:" + str(2222 + PS + i + PORT_OFFSET) for i in range(N_WORKER)]}

    if args.job_name == "ps":
        tf.app.run(start_parameter_server(int(args.task_index), spec))
    elif args.job_name == "worker":
        w = Worker(int(args.task_index), spec)
        tf.app.run(w.work())

reason: import numpy,import scipy


text:
import logging
logging.disable(logging.CRITICAL)

from tabulate import tabulate
from mjrl.utils.make_train_plots import make_train_plots
from mjrl.utils.gym_env import GymEnv
from mjrl.samplers.core import sample_paths
import numpy as np
import pickle
import time as timer
import os
import copy

def train_agent(job_name, agent,
                seed = 0,
                niter = 101,
                gamma = 0.995,
                gae_lambda = None,
                num_cpu = 1,
                sample_mode = 'trajectories',
                num_traj = 50,
                num_samples = 50000, # has precedence, used with sample_mode = 'samples'
                save_freq = 10,
                evaluation_rollouts = None,
                plot_keys = ['stoc_pol_mean'],
                env_kwargs=None
                ):

    np.random.seed(seed)
    if os.path.isdir(job_name) == False:
        os.mkdir(job_name)
    previous_dir = os.getcwd()
    os.chdir(job_name) # important! we are now in the directory to save data
    if os.path.isdir('iterations') == False: os.mkdir('iterations')
    if os.path.isdir('logs') == False and agent.save_logs == True: os.mkdir('logs')
    best_policy = copy.deepcopy(agent.policy)
    best_perf = -1e8
    train_curve = best_perf*np.ones(niter)
    mean_pol_perf = 0.0
    e = GymEnv(agent.env.env_id, env_kwargs=env_kwargs)

    for i in range(niter):
        print("......................................................................................")
        print("ITERATION : %i " % i)

        if train_curve[i-1] > best_perf:
            best_policy = copy.deepcopy(agent.policy)
            best_perf = train_curve[i-1]

        N = num_traj if sample_mode == 'trajectories' else num_samples
        args = dict(N=N, sample_mode=sample_mode, gamma=gamma, gae_lambda=gae_lambda, num_cpu=num_cpu, env=e.env_id, env_kwargs=env_kwargs)
        stats = agent.train_step(**args)
        train_curve[i] = stats[0]

        if evaluation_rollouts is not None and evaluation_rollouts > 0:
            print("Performing evaluation rollouts ........")
            eval_paths = sample_paths(num_traj=evaluation_rollouts, policy=agent.policy, num_cpu=num_cpu,
                                      env=e.env_id, eval_mode=True, base_seed=seed,env_kwargs=env_kwargs)
            mean_pol_perf = np.mean([np.sum(path['rewards']) for path in eval_paths])
            if agent.save_logs:
                agent.logger.log_kv('eval_score', mean_pol_perf)

        if i % save_freq == 0 and i > 0:
            if agent.save_logs:
                agent.logger.save_log('logs/')
                make_train_plots(log=agent.logger.log, keys=plot_keys, save_loc='logs/')
            policy_file = 'policy_%i.pickle' % i
            baseline_file = 'baseline_%i.pickle' % i
            pickle.dump(agent.policy, open('iterations/' + policy_file, 'wb'))
            pickle.dump(agent.baseline, open('iterations/' + baseline_file, 'wb'))
            pickle.dump(best_policy, open('iterations/best_policy.pickle', 'wb'))

        # print results to console
        if i == 0:
            result_file = open('results.txt', 'w')
            print("Iter | Stoc Pol | Mean Pol | Best (Stoc) \n")
            result_file.write("Iter | Sampling Pol | Evaluation Pol | Best (Sampled) \n")
            result_file.close()
        print("[ %s ] %4i %5.2f %5.2f %5.2f " % (timer.asctime(timer.localtime(timer.time())),
                                                 i, train_curve[i], mean_pol_perf, best_perf))
        result_file = open('results.txt', 'a')
        result_file.write("%4i %5.2f %5.2f %5.2f \n" % (i, train_curve[i], mean_pol_perf, best_perf))
        result_file.close()
        if agent.save_logs:
            print_data = sorted(filter(lambda v: np.asarray(v[1]).size == 1,
                                       agent.logger.get_current_log().items()))
            print(tabulate(print_data))

    # final save
    pickle.dump(best_policy, open('iterations/best_policy.pickle', 'wb'))
    if agent.save_logs:
        agent.logger.save_log('logs/')
        make_train_plots(log=agent.logger.log, keys=plot_keys, save_loc='logs/')
    os.chdir(previous_dir)


text:
#: -*- coding: utf-8 -*-
"""
Helper wrapper class for mathplotlib
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

mpl.rc('font', **{'sans-serif' : 'Arial', 'family' : 'sans-serif'})
majorFmt = mpl.dates.DateFormatter('%Y-%m-%d')

class Plotter:
	"""
	Helper class designed for drawing prediction model results.
	Can draw interval-valued data as bars and scalar data as point graph.
	Just create an instance: plotter = Plotter()
	And add dataframes: plotter.add(dataframe)
	Show graphs: plotter.show()
	"""

	def __init__(self, rows=1, cols=1, total=None):
		"""
		:param title: common title of graph
		:param xlabel: x-axis label, printed below graph
		:param ylabel: y-axis label, printed at left of graph
		"""
		self.axs = []
		self.fig = plt.figure(facecolor='white')

		width_ratios = [1 if x < 2 else 1 for x in range(cols)]
		gs = mpl.gridspec.GridSpec(rows, cols, width_ratios=width_ratios)

		for x in range(rows * cols):
			ax = plt.subplot(gs[x])
			ax.grid(alpha=0.4)
			ax.grid(which='minor', alpha=0.1)
			ax.grid(True)
			self.axs.append(ax)
		plt.tight_layout(w_pad=0.5, h_pad=3.0)

	def add(self, data, **kwargs):
		"""
		Analyzes data and draws suitable graph
		:param data: interval-valued or scalar dataframe
		:param label: name of dataframe
		"""
		pos = kwargs.get('pos', 0)
		self.axs[pos].xaxis.set_major_formatter(majorFmt)

		if isinstance(data['values'].iloc[0], (np.ndarray, list)):
			self._draw_intervals(data, **kwargs)
		else:
			self._draw_points(data, **kwargs)

	def _draw_points(self, data, **kwargs):
		"""
		Add scalar values graph to plot
		:param data: dataframe containing interval-valued data
		:param label: name of dataseries
		:param color: matplotlib color identifier , default is autoselected
		:param pos: position in graphs
		"""
		label = kwargs.get('label', '')
		color = kwargs.get('color', '')
		pos = kwargs.get('pos', 0)
		self.axs[pos].plot(data, label=label, marker='o', markersize=3,
		                   alpha=0.6, color=color)

	def _draw_intervals(self, data, **kwargs):
		"""Plots interval-valued data as bars
		:param data: dataframe containing interval-valued data
		:param label: name of interval-valued data
		:param color: matplotlib color identifier , default is autoselected
		:param pos: position in graphs
		"""
		label = kwargs.get('label', '')
		color = kwargs.get('color', '')
		pos = kwargs.get('pos', 0)
		heights = [x[1] for x in data['values']]
		bars = [x[0]-x[1] for x in data['values']]
		self.axs[pos].bar(data.index, bars, 0.015, heights, alpha=0.5,
		                  label=label, color=color, edgecolor=color)

	def add_bars(self, bars, labels, pos=0, **kwargs):
		"""
		Draws bar with given labels
		"""
		x = range(len(bars))
		color = kwargs.get('color', ['r', 'g', 'b', 'c'])

		#self.axs[pos].bar(x, bars, color=color, align='center')
		#self.axs[pos].set_xticks(x, minor=False)
		#self.axs[pos].set_xticklabels(labels)

		self.axs[pos].barh(x, bars, color=color, align='center')
		self.axs[pos].set_yticks(x, minor=False)
		self.axs[pos].set_yticklabels(labels)

	def prepare(self, pos=0, title='', xlabel='', ylabel=''):
		"""
		Prepares customized graph
		"""
		self.axs[pos].legend(loc=3, prop={'size':10})
		#self.fig.autofmt_xdate()
		self.axs[pos].set_xlabel(xlabel)
		self.axs[pos].set_ylabel(ylabel)
		self.axs[pos].set_title(title)

	def show(self):
		"""Shows prepared figure"""
		plt.show()

	def save(self, name):
		"""Saves prepared figure
		:param name: imange name and path where image is saved
		"""
		plt.ioff()
		self.fig.set_size_inches(20, 10)
		self.fig.savefig(name, dpi=200, bbox_inches='tight')

	def add_table(self, cols, labels, pos=0):
		"""Adds one row information table on top side of figure
		:param cols: table columns data
		:labels: table header labels
		"""
		#: no need in title when table used
		self.axs[pos].set_title('')
		self.axs[pos].table(cellText=[cols], colLabels=labels, loc='top')

	def __del__(self):
		"""Clean figure"""
		plt.close(self.fig)

max_stars_repo_name: zenio/ipredictor

max_issues_repo_name: zenio/ipredictor

max_forks_repo_name: zenio/ipredictor


text:
import tifffile
import os
import numpy as np
import sys
from scandir import scandir
from multiprocessing import Pool

# First nbg frames of input: Estimate median bg, generate sample
# f >= nbg < total-nbg:
def crop_tiff(tiffinput, tiffoutput, xymin, xymax):
    print (tiffoutput)
    with tifffile.TiffWriter(tiffoutput) as out_tif:
        with tifffile.TiffFile(tiffinput) as in_tif:
            total = len(in_tif.pages)

            for f in range(total):
                img = in_tif.pages[f].asarray()
                img = img[xymin[1]:xymax[1], xymin[0]:xymax[0]]
                out_tif.save(img)
                sys.stdout.write(f"\rframe {f}/{total} ({f/total*100:.2f}%)")
    print()

def _process(args):
    path, outdir, xymin, xymax = args
    print(f"pid={os.getpid()}: {path}")
    filename = os.path.split(path)[1]
    outfile = outdir + filename
    crop_tiff(path,outfile,xymin,xymax)

def crop_dir(inputdir, outputdir, xymin, xymax):

    def cb(fn):
        args=[ fn, outputdir, xymin, xymax]
        #params.append(args)
        _process(args)

    os.makedirs(outputdir,exist_ok=True)

    scandir(inputdir, "*.tif", cb)

#    p = Pool(8)
#    p.map(_process, params)

if __name__ == "__main__":
    crop_dir('O:/mod/', 'O:/mod-crop2/',  [50,50],[350,350])

max_stars_repo_name: qnano/simflux

max_issues_repo_name: qnano/simflux

max_forks_repo_name: qnano/simflux