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

"""NumPy implementation of the activation transformations used in DKS/TAT."""

import itertools
import os

from absl import logging

from autograd import elementwise_grad as egrad
from autograd import numpy as np
from dks.base.activation_getter import get_activation_function as _get_numpy_activation_function
import scipy.integrate as sp_int
import scipy.optimize as sp_opt
from scipy.special import roots_legendre

# pylint: disable=g-import-not-at-top
# This is a trick to achieve compatibility with multiple versions of SciPy:
try:
  from scipy.integrate._quadrature import _cached_roots_legendre
except ImportError:
  from scipy.integrate.quadrature import _cached_roots_legendre
# pylint: enable=g-import-not-at-top

# Random seed used to initialize activation function parameter searches after
# standard starting points fail (which almost never happens).
_RANDOM_SEED = 123

_CURRENT_DIR = os.path.dirname(os.path.abspath(__file__))

# Feel free to change the path of this file if needed. This file should ideally
# persist between experiments since it can take around 20 minutes to regenerate.
_ROOTS_CACHE_FILE = os.path.join(_CURRENT_DIR, "roots_{}.npy")

_QUADRATURE_ORDER = 100000

# These represent sensible values and ranges to initialize the activation
# function parameter search with. They are sensible because commonly used
# activation functions tend to have reasonable behavior for such values. If
# the searches fail for some activation function they may need to be
# tweaked/expanded.
_ALWAYS_TRY_VALUES = {"input_shift": (0.0, 1.0, -1.0),
                      "input_scale": (1.0, 0.1),
                      "output_shift": (0.0, 1.0, -1.0),
                      "negative_slope": (0.0, 0.5)}

_SAMPLE_RANGE_LOW = {"input_shift": -3.0,
                     "input_scale": 0.0,
                     "output_shift": -3.0,
                     "negative_slope": 0.0}

_SAMPLE_RANGE_HIGH = {"input_shift": 3.0,
                      "input_scale": 2.0,
                      "output_shift": 3.0,
                      "negative_slope": 1.0}

def _precompute_or_load_roots(order):
  """Compute or load the roots used by fixed_quad to save time."""

  if order not in _cached_roots_legendre.cache:

    roots_cache_file = _ROOTS_CACHE_FILE.format(order)

    if os.path.exists(roots_cache_file):

      with open(roots_cache_file, "rb") as fhandle:
        _cached_roots_legendre.cache[order] = np.load(fhandle,
                                                      allow_pickle=False)

    else:
      roots = roots_legendre(order)
      _cached_roots_legendre.cache[order] = roots

      with open(roots_cache_file, "wb") as fhandle:
        np.save(fhandle, roots, allow_pickle=False)

def _estimate_gaussian_mean(fn, order=_QUADRATURE_ORDER):
  """Estimate the mean of a function fn(x) where x ~ N(0,1)."""

  _precompute_or_load_roots(order)

  fn_weighted = lambda x: np.exp(-x**2 / 2) * fn(x)

  integral, _ = sp_int.fixed_quad(fn_weighted, -10., 10., n=order)

  return integral / np.sqrt(2*np.pi)

def _calc_c_map(activation, derivative_order=0, c=1.0, q_output=None):
  """Evaluate local C map value assuming an input q value of 1.

  Args:
    activation: A callable representing the activation function (applied
      elementwise) from which to define the local C map.
    derivative_order: An integer giving the order of the derivative of the C map
      to take before evaluating it. (Default: 0)
    c: A float giving the input point at which to evaluate the C map. Must
      be 0.0 or 1.0. (Default: 1.0)
    q_output: Float or None giving the output q value associated with
      ``activation``, if this is known. If None this will be computed from
      scratch. (Default: None)

  Returns:
     A float giving the value of the (derivative of) the local C map for the
     given activation function.
  """

  derivative = activation
  for _ in range(derivative_order):
    derivative = egrad(derivative)

  if c == 0.0:
    integral = _estimate_gaussian_mean(derivative)**2
  elif c == 1.0:
    integral = _estimate_gaussian_mean(lambda x: derivative(x)**2)
  else:
    raise NotImplementedError("Input c value must be 0.0 or 1.0.")

  if q_output is None:
    q_output = _estimate_gaussian_mean(lambda x: activation(x)**2)

  return integral / q_output

def _calc_c_slope(activation, c=1.0, q_output=None):
  """Evaluate local C map derivative assuming an input q value of 1."""
  return _calc_c_map(activation, derivative_order=1, c=c, q_output=q_output)

def _calc_c_curv(activation, c=1.0, q_output=None):
  """Evaluate local C map second derivative assuming an input q value of 1."""
  return _calc_c_map(activation, derivative_order=2, c=c, q_output=q_output)

def _calc_q_slope_1(activation):
  """Computes the derivative of a local Q map at q=1."""
  derivative = egrad(activation)
  return _estimate_gaussian_mean(lambda x: activation(x)*derivative(x)*x)

def _leaky_relu_cmap(c, negative_slope):
  """Evaluates the local C map for Leaky ReLU with the given negative slope."""
  return ((1 - negative_slope)**2 * (np.sqrt(1 - c**2)
                                     + (np.pi - np.arccos(c)) * c) / np.pi
          + 2 * negative_slope * c) / (1 + negative_slope**2)

def _compute_output_params(act, local_c_val_0_target):
  """Compute output params to achieve Q(1)=1 and C(0)=local_c_val_0_target."""

  if local_c_val_0_target is not None:
    output_shift = np.sqrt(local_c_val_0_target) - _estimate_gaussian_mean(act)
    act_shifted = lambda x: act(x) + output_shift
  else:
    output_shift = None
    act_shifted = act

  output_scale = 1. / np.sqrt(
      _estimate_gaussian_mean(lambda x: act_shifted(x)**2))

  return output_shift, output_scale

def _transform_activation(phi, params):
  """Transform an activation function phi using the given parameters."""

  params = params.copy()

  input_scale = params.pop("input_scale", None)
  input_shift = params.pop("input_shift", None)
  output_shift = params.pop("output_shift", None)
  output_scale = params.pop("output_scale", None)

  def activation(x):

    # Note: DO NOT use += and *= below! Bad things will happen.
    if input_scale is not None:
      x = x * input_scale
    if input_shift is not None:
      x = x + input_shift

    x = phi(x, **params)

    if output_shift is not None:
      x = x + output_shift
    if output_scale is not None:
      x = x * output_scale

    return x

  return activation

def _solve_for_activation_params(
    name, local_q_slope_target, local_c_val_0_target, local_c_slope_1_target,
    local_c_slope_0_target, local_c_curv_target, reject_condition=None,
    num_tries=50):
  """Computes activation function parameters to achieve the given targets."""

  # Making sure random starting points used in solvers will be the same for
  # each run, so that they find the same exact solutions (which is important
  # when using JAX with multiple processes).
  np.random.seed(_RANDOM_SEED)

  # RELU has 1 less degree of freedom so we use this special-case logic:
  if name == "relu":
    # The constant below is the maximum value of local_c_slope_1 that can be
    # achieved in ReLUs given a positive input shift value.
    if local_c_slope_1_target < 1.466942206924260361:
      input_shift = 1.0
    else:
      input_shift = -1.0

    constant_params = {"input_shift": input_shift}
    opt_params = ("input_scale",)

    local_q_slope_target = None  # we turn this off for ReLUs
    assert local_c_slope_0_target is None
    assert local_c_curv_target is None

  elif name == "leaky_relu":
    constant_params = {}
    opt_params = ("negative_slope", "input_shift")

  elif local_c_val_0_target is not None:
    constant_params = {}
    opt_params = ("input_scale", "input_shift")

  else:
    constant_params = {}
    opt_params = ("input_scale", "input_shift", "output_shift")

  def make_params_from_pvector(p):
    """Make a dictionary of activation parameters from a given vector."""

    params = constant_params.copy()

    for i, pname in enumerate(opt_params):
      assert pname not in params
      params[pname] = p[i]

    # We directly compute output_scale, and sometimes also output_shift, instead
    # of taking them from the parameter vector.
    output_shift, output_scale = _compute_output_params(
        _transform_activation(_get_numpy_activation_function(name), params),
        local_c_val_0_target)

    if output_shift is not None:
      assert "output_shift" not in params
      params["output_shift"] = output_shift

    if output_scale is not None:
      assert "output_scale" not in params
      params["output_scale"] = output_scale

    return params

  def compute_error_vector(p):
    """Computes vector of errors (value - target) over relevant quantities."""

    params = make_params_from_pvector(p)
    phi_hat = _transform_activation(_get_numpy_activation_function(name),
                                    params)

    residual_array = []

    if local_q_slope_target is not None:
      local_q_slope = _calc_q_slope_1(phi_hat)
      residual_array += [local_q_slope - local_q_slope_target]

    if local_c_slope_1_target is not None:
      local_c_1_slope = _calc_c_slope(phi_hat, c=1.0, q_output=1.0)
      residual_array += [local_c_1_slope - local_c_slope_1_target]

    if local_c_slope_0_target is not None:
      local_c_0_slope = _calc_c_slope(phi_hat, c=0.0, q_output=1.0)
      residual_array += [local_c_0_slope - local_c_slope_0_target]

    if local_c_curv_target is not None:
      local_c_curv = _calc_c_curv(phi_hat, c=1.0, q_output=1.0)
      residual_array += [local_c_curv - local_c_curv_target]

    return np.asarray(residual_array)

  # Make the starting points for the search:
  always_try = tuple(_ALWAYS_TRY_VALUES[pname] for pname in opt_params)
  starting_points = list(itertools.product(*always_try))
  for _ in range(num_tries - len(starting_points)):
    starting_points += [
        tuple(np.random.uniform(low=_SAMPLE_RANGE_LOW[pname],
                                high=_SAMPLE_RANGE_HIGH[pname])
              for pname in opt_params)]

  # For each starting point we run sp_opt.root to try and find a solution:
  for starting_point in starting_points:

    sol = sp_opt.root(compute_error_vector, np.asarray(starting_point),
                      method="hybr", jac=False, options=None)
    if sol.success:
      params = make_params_from_pvector(sol.x)
      if reject_condition is None or not reject_condition(name, params):
        break

    logging.debug("Failed to find parameters from starting point %s.",
                  starting_point)

  if not sol.success:
    raise ValueError(f"Failed to find parameters for '{name}'!")

  logging.info("Found parameters for '%s': %s", name, params)

  return params

def _solve_increasing(fn, target, input_, min_, max_, tol=1e-8, max_eval=100):
  """Solves for x in fn(x)=target, where fn is an increasing function.

  Args:
    fn: A callable which takes a scalar input x and produces a scalar output.
      Must compute an increasing function.
    target: The target output value of ``fn``.
    input_: The initial guess for x.
    min_: A lower bound on the possible value of x.
    max_: An upper bound on the possible value of x.
    tol: A float which giving the acceptable tolerance for ``|fn(x) - target|``.
      (Default: 1e-8)
    max_eval: An integer giving the maximum number of times to evaluate ``fn``
      before giving up. (Default: 100)

  Returns:
    A float giving a value of x such that ``|fn(x) - target| < tol``.
  """

  # The method used to find the solution is a simple binary search in the
  # inteval [min_, max_], where max_ or min_ will change each iteration. If
  # max_ is infinity we double our current guess instead of averaging it with
  # max_.
  for _ in range(max_eval):

    value = fn(input_)

    logging.debug("binary search vals: min = %f, input = %f, max = %f, "
                  "target = %f, value = %f", min_, input_, max_, target, value)

    if np.abs(value - target) < tol:
      return input_

    if value > target:
      max_ = input_
      input_ = 0.5 * (input_ + min_)

    elif value < target:
      min_ = input_

      if np.isinf(max_):
        input_ = input_ * 2
      else:
        input_ = 0.5 * (input_ + max_)

  raise ValueError(f"Maximum evaluations ({max_eval}) exceeded while searching "
                   "for solution. This is probably due the specified targets "
                   "being unachievable for the given architecture.")

def _compute_local_c_slope_1_target(max_slope_func, target_value):
  return _solve_increasing(max_slope_func, target_value, 1.1, 1., np.inf)

def _compute_local_c_slope_0_target(max_slope_func, target_value):
  return _solve_increasing(max_slope_func, target_value, 0.99, 0., 1.0)

def _compute_local_c_curv_target(max_curv_func, target_value):
  return _solve_increasing(max_curv_func, target_value, 0.1, 0., np.inf)

def _compute_negative_slope_param(max_lrelu_c0, target_value):
  return _solve_increasing(
      lambda a: 1.0 - max_lrelu_c0(a), 1.0 - target_value, 0.5, 0., 1.0)

def _verify_params_dict_and_set_defaults(params_dict, defaults):
  """Verify keys in parameter dict and set any missing ones to the defaults."""

  bad_keys = set(params_dict.keys()).difference(set(defaults.keys()))

  if bad_keys:
    raise ValueError(
        f"Parameter dictionary had unrecognized keys: '{bad_keys}'")

  for key in defaults.keys():
    if key not in params_dict:
      params_dict[key] = defaults[key]

def _get_activations_params(
    activation_names, method="DKS", dks_params=None, tat_params=None,
    max_slope_func=None, max_curv_func=None, subnet_max_func=None):
  """Get dict of optimized parameters for given named activation functions."""

  if not isinstance(activation_names, (list, tuple)):
    raise ValueError("activation_names argument must be a list or tuple of "
                     "strings.")

  # Note that using dictionaries as defaults in the function def is bad, hence
  # we do this instead:
  if dks_params is None:
    dks_params = {}
  if tat_params is None:
    tat_params = {}

  _verify_params_dict_and_set_defaults(
      dks_params, {"c_slope_1_target": 1.5, "local_q_slope_target": 1.0})
  _verify_params_dict_and_set_defaults(
      tat_params, {"c_val_0_target": 0.9, "c_curve_target": 0.3})

  local_q_slope_target = None
  local_c_val_0_target = None
  local_c_slope_1_target = None
  local_c_slope_0_target = None
  local_c_curv_target = None
  reject_condition = None

  if method == "DKS":

    if "relu" in activation_names:
      logging.warning("The use of ReLUs with DKS is *highly* discouraged. You "
                      "are advised to use Leaky ReLUs instead.")

    c_slope_1_target = dks_params["c_slope_1_target"]

    if c_slope_1_target <= 1.0:
      raise ValueError("Invalid value for DKS 'c_slope_1_target' parameter. "
                       "Must be a float greater than 1.0.")

    if max_slope_func is None:
      if subnet_max_func is None:
        raise ValueError("Must supply 'subnet_max_func' if using DKS and not "
                         "passing in 'max_slope_func'.")
      # We can compute the maximal slope function by replacing composition
      # with multiplication in the maximal c value function.
      max_slope_func = lambda x: subnet_max_func(1.0, lambda y: x * y)

    # Three of the four conditions used by DKS. The remaining condition Q(1) = 1
    # is implied.
    local_q_slope_target = dks_params["local_q_slope_target"]
    local_c_val_0_target = 0.0
    # We set the local slope to achieve C'_f(1) <= target over all subnetworks
    # f:
    local_c_slope_1_target = _compute_local_c_slope_1_target(
        max_slope_func, c_slope_1_target)

    logging.info("Found 'local_c_slope_1_target': %s", local_c_slope_1_target)

  elif method == "TAT" and "leaky_relu" not in activation_names:

    if "relu" in activation_names:
      raise ValueError("Standard ReLU is not supported with TAT. Use leaky "
                       "ReLU instead.")

    c_curve_target = tat_params["c_curve_target"]

    if c_curve_target <= 0.0:
      raise ValueError("Invalid value for TAT 'c_curve_target' parameter. Must "
                       "be a float greater than 0.0.")

    if max_curv_func is None:
      if subnet_max_func is None:
        raise ValueError("Must supply 'subnet_max_func' if using TAT with "
                         "smooth activations and not passing in "
                         "'max_curv_func'.")
      # We can compute the maximal curvature function by replacing composition
      # with addition in the maximal c value function.
      max_curv_func = lambda x: subnet_max_func(0.0, lambda y: x + y)

    # Three of the four conditions used by TAT in the smooth case. The remaining
    # condition Q(1) = 1 is implied.
    local_q_slope_target = 1.0
    local_c_slope_1_target = 1.0
    # We set the local second derivative to achieve C''_f(1) <= target over all
    # subnetworks f:
    local_c_curv_target = _compute_local_c_curv_target(
        max_curv_func, c_curve_target)

    logging.info("Found 'local_c_curv_target': %s", local_c_curv_target)

    # This is a hacky fix used to avoid certain 'bad' solutions we observed that
    # seem to have unstable Q maps and higher kernel approximation errors. It
    # should probably be replaced with something more principled.
    reject_condition = lambda name, params: (  # pylint: disable=g-long-lambda
        params["input_scale"] * params["output_scale"] >= 2.0)

  elif method == "TAT" and "leaky_relu" in activation_names:

    if len(activation_names) > 1:
      raise ValueError("When using Leaky ReLU with TAT it must be the only "
                       "activation function.")

    c_val_0_target = tat_params["c_val_0_target"]
    if c_val_0_target > 1.0 or c_val_0_target < 0.0:
      raise ValueError("Invalid value for TAT 'c_val_0_target' parameter. Must "
                       "be a float between 0.0 and 1.0.")

    if subnet_max_func is None:
      raise ValueError("Must supply 'subnet_max_func' if using TAT with Leaky "
                       "ReLU activation functions.")

    max_lrelu_c0 = lambda neg_slope: subnet_max_func(  # pylint: disable=g-long-lambda
        0.0, lambda c: _leaky_relu_cmap(c, neg_slope))

    # We set the negative slope parameter to achieve C_f(0) <= target over all
    # subnetworks f:
    negative_slope = _compute_negative_slope_param(
        max_lrelu_c0, c_val_0_target)

    # This is the value required to achieve Q(1) = 1 for Leaky ReLUs. See the
    # TAT paper for details.
    output_scale = np.sqrt(2.0 / (1.0 + negative_slope**2))

    logging.info("Found parameters for 'leaky_relu': negative_slope = %s, "
                 "output_scale = %s.", negative_slope, output_scale)

    return {"leaky_relu": {"output_scale": output_scale,
                           "negative_slope": negative_slope}}

  else:
    raise ValueError(f"Unrecognized value for argument 'method': {method}")

  params = {}
  for name in activation_names:

    params[name] = _solve_for_activation_params(
        name,
        local_q_slope_target=local_q_slope_target,
        local_c_val_0_target=local_c_val_0_target,
        local_c_slope_1_target=local_c_slope_1_target,
        local_c_slope_0_target=local_c_slope_0_target,
        local_c_curv_target=local_c_curv_target,
        reject_condition=reject_condition)

  return params

def get_transformed_activations(
    activation_names, method="DKS", dks_params=None, tat_params=None,
    max_slope_func=None, max_curv_func=None, subnet_max_func=None,
    activation_getter=_get_numpy_activation_function):
  """Gets transformed activation functions using the DKS or TAT method.

  See the DKS paper (https://arxiv.org/abs/2110.01765) and the TAT paper
  (https://openreview.net/forum?id=U0k7XNTiFEq) for details about what these
  are, how they are computed, and what their parameters mean. A procedure to
  compute the "maximal slope function" is given in the section titled "Summary
  of our method" of the DKS paper. Procedures to compute the "maximal curvature
  function", and the "subnet maximizing function", are given in the appendix
  section titled "Additional details and pseudocode for activation function
  transformations" of the TAT paper.

  Note that if you are using the JAX, PyTorch, or TensorFlow frameworks, you
  probably want to be using the version of get_transformed_activations() in the
  corresponding subpackage. (These are basically thin wrappers around this
  function that pass a framework-specific value to the ``activation_getter``
  argument.)

  Args:
    activation_names: An iterable of string names for the activation functions.
      Supported names are the intersection of those supported by
      dks.base.activation_getter.get_activation_function, and those supported
      by the getter passed to the ``activation_getter`` argument (which defaults
      to dks.base.activation_getter.get_activation_function). The built-in
      getters in this package (for each framework) currently support the
      following names: "tanh", "softplus", "leaky_relu", "relu" (not
      recommended; use "leaky_relu" instead), "selu", "elu", "swish", "sigmoid",
      "erf", "bentid", "atan", "asinh", "square", and "softsign".
    method: A string representing the method used to transform the activation
      functions. Can be "DKS", "TAT", or "untransformed". The latter choice
      will return activation functions without any transformations.
      (Default: "DKS")
    dks_params: A dictionary containing the parameters to use for DKS. Keys
      should be a subset of {"c_slope_1_target", "local_q_slope_target"}.
      "c_slope_1_target" gives the target maximal slope value for the network
      (corresponding to "zeta" from the paper), and defaults to 1.5.
      "local_q_slope_target" gives the target value for the local Q map slope
      of each nonlinear layer (which is kept at 1.0 in the paper -- except in
      ablation tests), and defaults to 1.0. If ``dks_params`` is passed as None,
      it defaults to the empty dictionary (so that the parameters will use their
      default values). (Default: None)
    tat_params: A dictionary containing the parameters to use for TAT. Keys
      should be a subset of {"c_val_0_target", "c_curve_target"}.
      "c_val_0_target" gives the maximum value of ``C_f(0)`` over subnetworks f,
      which is used when transforming Leaky ReLUs (and corresponds to "eta" from
      the paper), and defaults to 0.9. "c_curve_target" gives the maximum value
      of ``C''_f(1)``, which is used for all other activation functions (and
      corresponds to "tau" from the paper), and defaults to 0.3. If
      ``tat_params`` is passed as None, it defaults to the empty dictionary (so
      that the parameters will use their default values). (Default: None)
    max_slope_func: A callable which computes the "maximal slope function" of
      the network, as defined in the DKS paper. It should take a single argument
      representing the slope of each local C map at ``c=1``. If this is required
      (i.e. when using DKS) but passed as None, it will be generated using
      ``subnet_max_func`` when possible. (Default: None)
    max_curv_func: A callable which computes the "maximal curvature function" of
      the network, as defined in the TAT paper. It should take a single
      parameter representing the second derivative of each local C map at c=1.
      If this is required (i.e. when using TAT with smooth activation functions)
      but is passed as None, it will be generated using ``subnet_max_func`` when
      possible. (Default: None)
    subnet_max_func: A callable which computes the "subnetwork maximizing
      function" of the network, as defined in the TAT paper (and denoted
      ``M_{f,r}(x)``). It should take two arguments: the input value ``x``, and
      a callable ``r_fn`` which maps a float to a float. This is required when
      using TAT with Leaky ReLUs. (Default: None)
    activation_getter: A callable which takes a string name for an activation
      function and returns the (untransformed) activation function corresponding
      to this name. Defaults to one returning activation functions in NumPy
      (with autograd). Returned transformed activation functions will be based
      on the output of this callable. Other tensor frameworks can be supported
      by changing this argument. See the versions of
      get_transformed_activations() in the ``dks.jax``, ``dks.pytorch``, and
      ``dks.tensorflow`` subpackages.

  Returns:
    A dictionary mapping the activation function names to their corresponding
    transformed activation functions.
  """

  if method == "untransformed":
    return {name: activation_getter(name) for name in activation_names}

  params = _get_activations_params(
      activation_names, method=method, dks_params=dks_params,
      tat_params=tat_params, max_slope_func=max_slope_func,
      max_curv_func=max_curv_func, subnet_max_func=subnet_max_func)

  transformed_acts = {}

  for name in activation_names:
    transformed_acts[name] = _transform_activation(activation_getter(name),
                                                   params[name])
  return transformed_acts

max_stars_repo_name: deepmind/dks

max_issues_repo_name: deepmind/dks

max_forks_repo_name: deepmind/dks

reason: import numpy,import scipy,from scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 17 09:11:12 2020

@author: emc1977
"""

import numpy as np
from scipy import optimize

def f(x):   # The rosenbrock function
    return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2

def jacobian(x):
    return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))

solution = optimize.minimize(f, [2,-1], method="Newton-CG", jac=jacobian)
print(solution)

ext: py

max_stars_repo_path: Newton Scipy of Rosenbrock (1).py

max_issues_repo_path: Newton Scipy of Rosenbrock (1).py

max_forks_repo_path: Newton Scipy of Rosenbrock (1).py

reason: import numpy,from scipy


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

"""
    Graph from monitoring data
    Command line parameters
        1st csv input file
        2nd column number for point IDs (zero based, optional)
        3rd column number for x (datetime) values (optional)
        4th column number for y values (optional)
        5th point number to use in chart (optional) several point IDs can be given
            if no point IDs are given all points are plotted
"""

import sys
import numpy as np
import csv
import os
from datetime import datetime
import matplotlib.pyplot as plt

if len(sys.argv) < 2:
    print("""
usage: input_csv_file p_num_col x_colum y_column [p_num1 [pnum2 ...]]
    Command line parameters
        1st csv input file
        2nd column number for point IDs (zero based)
        3rd column number for x (datetime) values
        4th column number for y values
        5th point number to use in chart (optional) several point IDs can be given
            if no point IDs are given all points are plotted
        """)
    sys.exit()

if not os.path.exists(sys.argv[1]):
    print("File does not exist: " + sys.argv[1])
    sys.exit()
p_col = 0       # column number for point numbers values
if len(sys.argv) > 2:
    p_col = int(sys.argv[2])
x_col = 4       # column for x values
if len(sys.argv) > 3:
    x_col = int(sys.argv[3])
y_col = 3       # column for y values
if len(sys.argv) > 4:
    y_col = int(sys.argv[4])
p_nums = None
if len(sys.argv) > 5:
    p_nums = sys.argv[4:]
# load data into dict for point numbers
data = {}
with open(sys.argv[1], newline='') as f:
    for row in csv.reader(f, delimiter=";"):
        p_num = row[p_col]
        if p_nums is not None and p_num not in p_nums:
            continue    # skip point
        dt = datetime.strptime(row[x_col], "%Y-%m-%d %H:%M:%S")
        y = float(row[y_col])
        if not p_num in data:
            data[p_num] = [[], []]
        data[p_num][0].append(dt)
        data[p_num][1].append(y)
fig = plt.figure(1)
for key, val in data.items():
    plt.plot_date(val[0], val[1], '-', label=key)
#for key in data.keys():
#    plt.plot_date(data[key][0], data[key][1], '-', label=key)
plt.xticks(rotation=45)
plt.xlabel("date time")
plt.grid()
plt.legend()
plt.show()

max_stars_repo_name: hrutkabence/tutorials

max_issues_repo_name: hrutkabence/tutorials

max_forks_repo_name: hrutkabence/tutorials

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 15 11:54:23 2019

@author: smile
"""
import numpy as np
from sklearn import tree
from sklearn.datasets import load_iris
iris = load_iris()

print(iris.feature_names)  # like: ['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
print(iris.target_names)  # output labels like: ['setosa' 'versicolor' 'virginica']

print(iris.data[0])  # actual set of features like:  [5.1 3.5 1.4 0.2]
print(iris.target[0])  # current label

# iterate through all of the features:
for i in range(len(iris.target)):
    print("Example: %d, features: %s, label: %s" % (i, iris.data[i], iris.target[i]))

# save test index for later tests, we will not use it for training
test_idx = [0, 50, 100]

# training data
train_target = np.delete(iris.target, test_idx)  # removing 3 samples from the training labels set
train_data = np.delete(iris.data, test_idx, axis=0)  # removing 3 samples from the training data set

# testing data
test_target = iris.target[test_idx]
test_data = iris.data[test_idx]

# classifier training starts
clf = tree.DecisionTreeClassifier()
clf = clf.fit(train_data, train_target)

# out result
print(test_target)  # print this as example and compare to the next line
print(clf.predict(test_data))  # print out result and compare with previous line

# viz code
import matplotlib

vis_out = tree.plot_tree(clf.fit(iris.data, iris.target))
matplotlib.pyplot.show(vis_out)


text:
import covid19sim
import numpy as np
import os
import datetime
import matplotlib
import matplotlib.pyplot as plt
from covid19sim.utils.constants import SECONDS_PER_DAY
from pathlib import Path

def run(data, path, comparison_key):

    label2pkls = list()
    for method in data:
        for key in data[method]:
            label = f"{method}_{key}"
            pkls = [r["pkl"] for r in data[method][key].values()]
            label2pkls.append((label, pkls))

    for label, pkls in label2pkls:
        times = list()
        for data in pkls:
            data = data["infection_monitor"]
            for x in data:
                if x["from"]:
                    times.append((x["infection_timestamp"] - x["from_infection_timestamp"]).total_seconds() / SECONDS_PER_DAY)

        fig, ax = plt.subplots(figsize=(7.5, 7.5))
        ax.hist(times, color="royalblue", bins=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21])
        ax.set(title=f"Statistics of Generation Time {label}")
        ax.set_xlabel("Generation Time")
        ax.set_ylabel("Frequency")
        output_dir = os.path.join(path, "generation_times")
        Path(output_dir).mkdir(parents=True, exist_ok=True)
        output_file = os.path.join(output_dir, f"GenerationTime_{label}.png")
        print(f"printing: {output_file}")
        plt.savefig(output_file)

max_stars_repo_name: mila-iqia/COVI-AgentSim

max_issues_repo_name: mila-iqia/COVI-AgentSim

max_forks_repo_name: mila-iqia/COVI-AgentSim


text:
#!/usr/bin/env python3
# written by Toggi3 - <jeffj8086@gmail.com>

import discord
import asyncio
import subprocess
import json
import os.path
import numpy
from time import sleep

cfgData=json.loads(open("../learnpython/ToggBot_Working/config.json").read()) # Load config

client=discord.Client()

def runScript(args):
    botCmd=args[1:]
    cmd=botCmd.split(' ',1)[0]
    if (os.path.isfile("./scripts/%s/runScript" % cmd)==True):
        if (len(botCmd.split())==1):
            params=" " # So that params is not null
        else: # The below quotes the arguments passed for shell sanitation
            params=' '.join('"{}"'.format(word) for word in botCmd.split(' ',1)[1].split(' '))
        return(subprocess.check_output("./scripts/%s/runScript %s" % (cmd,params),shell=True).decode('utf-8').splitlines())
    else:
        return("$writeChannel$ No such command: \"%s\"" % cmd).splitlines()

def shiftList(data):
    list=numpy.roll(data,-1)
    return(list)

@client.event
async def on_ready():
    print('Logged in as')
    print(client.user.name)
    print(client.user.id)
    print('------')

@client.event
async def on_message(message):
    if message.content.startswith('!'):
        scriptReturn=runScript(message.content)
        for i in scriptReturn:
            action=scriptReturn[0].split(None, 1)[0]
            response=" ".join(scriptReturn[0].split()[1:])
            if (action=="$writeChannel$"):
                await client.send_message(message.channel,response)
            elif (action=="$tmpWriteChannel$"): # for writing messages to be edited later
                tmp=await client.send_message(message.channel,response)
            elif (action=="$editChannel$"):
                await client.edit_message(tmp,response)
                sleep(0.5) # for timing of silly things like animations
            else:
                print("Unrecognized command: %s %s" % (action, response))
            if len(scriptReturn)>1: # shift list of commands if multi-line.
                scriptReturn=shiftList(scriptReturn)
    else:
        return(0)
client.run(cfgData["discord_email"],cfgData["discord_password"])

max_stars_repo_name: Toggie/ToggBot

max_issues_repo_name: Toggie/ToggBot

max_forks_repo_name: Toggie/ToggBot


text:
#! /usr/bin/env python
"""Takes trajectory data from xyz and uses inverse kinematics algorithms to calculate the best state. """

import rospy
import numpy
from scipy.linalg import pinv
import time
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from math import pi
from math import cos
from math import sin
from sympy import *
from robot_functions import *	# Built functions for jaco model

if __name__ == '__main__':
  try:
    ti = time.time()
    rospy.init_node('move_robot_using_trajectory_msg')

    # Allow gazebo to launch
    rospy.sleep(1)

    ######################### EDIT THIS SECTION ###############################
    # Parameters
    numPts = 10		# Discretized path (Choose number of points)
    step   = 0.5	# Step size
    error  = 0.001	# Error between state and desired state (theshold value, in meters)

    # Initial state
    q10 = 0.0
    q20 = 2.9
    q30 = 1.3
    q40 = 4.2
    q50 = 1.4
    q60 = 0.0

    # Final state
    xf = 0.2
    yf = 0.6
    zf = 0.4
    ###########################################################################

    # Boundary conditions (DON'T WORRY ABOUT THESE, and DON'T UNCOMMENT)
    init_q  = Matrix([[-q10], [q20-pi/2], [q30+pi/2], [q40], [q50-pi], [q60+pi/2]])

    # Initial position
    _, init_X = systemKinematics(init_q)
    final_X = Matrix([[xf],[yf],[zf]])

    ######################### EDIT THIS SECTION ###############################
    # Discretize path (change this to change path shape, default is line)
    xPath = numpy.linspace(init_X[0], final_X[0], num=numPts)
    yPath = numpy.linspace(init_X[1], final_X[1], num=numPts)
    zPath = numpy.linspace(init_X[2], final_X[2], num=numPts)
    ###########################################################################

    # Move robot to initial state
    moveJoint([q10,q20,q30,q40,q50,q60])
    moveFingers([1.3,1.3,1.3])

    # Initialize variables
    delta_X = Matrix([[1000],[1000],[1000]])

    X = []
    Y = []
    Z = []

    # Matrix calculatations
    dynMatrix,_ = systemKinematics(init_q)
    row0 = dynMatrix.row(0)
    row1 = dynMatrix.row(1)
    row2 = dynMatrix.row(2)
    DynMatrix = Matrix([[row0],[row1],[row2]])
    J = robotJacobian(DynMatrix)

    # Connect all points in mesh path (START LOOP AT 1, NOT 0)
    for i in range(1,numPts):
      path_X = Matrix([[xPath[i]],[yPath[i]],[zPath[i]]])
      print(i)
      _, init_X = systemKinematics(init_q)
      delta_X = path_X - Matrix([[init_X[0]],[init_X[1]],[init_X[2]]])

      while (Abs(delta_X.norm()) > error):
        _, init_X = systemKinematics(init_q)
        delta_X = path_X - Matrix([[init_X[0]],[init_X[1]],[init_X[2]]])
        print(Abs(delta_X.norm()))

        # Linearize dynamics matrix
        useful_J = evalJacobian(J,init_q)

        # Calculate joint state change and new state
        delta_q = inverseKinematics(useful_J,step,delta_X)
        init_q = init_q + step*delta_q

        # Convert back from DH angles to JACO angles
        moveJoint([-init_q[0], init_q[1]+pi/2, init_q[2]-pi/2, init_q[3], init_q[4]+pi, init_q[5]-pi/2])

        # Trajectory plotting
        X.append(init_X[0])
        Y.append(init_X[1])
        Z.append(init_X[2])

    elapsed = time.time() - ti
    print(elapsed)

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(X, Y, Z, label='traveled')
    ax.plot(xPath, yPath, zPath, label='path')
    ax.legend()
    plt.show()

  except rospy.ROSInterruptException:
    print "program interrupted before completion"

reason: import numpy,from scipy,from sympy


text:
import os
import cv2 # pip install opencv-python
import time
import numpy
import mouse
import pyautogui

def get_pics(directory: str) -> list:
    """Grabs the images from /media/<directory> and returns their file names.

    Args:
        dir ([type]): The subfolder (/media/<directory>) for the images to load

    Returns:
        list[any]: A OpenCV image data list
    """

    loading_screens = []
    for entry in os.scandir(os.getcwd() + f'/media/{directory}'):
        if entry.is_file():
            file_path = os.getcwd().replace('\\', '/') + f'/media/{directory}/' + entry.name
            loading_screens.append({'data': cv2.imread(file_path, 0), 'file': file_path})
    return loading_screens

def check(sh: numpy.ndarray, template, threshold: float) -> bool:
    """Scans the screen and returns if it found the template

    Args:
        sh (numpy.ndarray): Image of the current screen
        template ([type]): Image of what to find
        threshold (float): Detection threshold

    Returns:
        bool: If it found the template given on screen
    """

    res = cv2.matchTemplate(sh, template, cv2.TM_CCOEFF_NORMED)
    loc = numpy.where(res >= threshold)
    if loc[0].size != 0:
        return True
    return False

def screen_loop(threshold: float, dir_content: dict) -> None:
    """Returns when found the screen

    Args:
        dir_name (str): Directory of training pictures
        threshold (float): [description]

    Returns:
        str: Path of the image that was found
    """

    while True:
        time.sleep(0.2)
        screenshot = numpy.asarray(pyautogui.screenshot().convert(mode='L'))

        for folder in [d for d in os.listdir() if d.isdir()]:
            for screen in get_pics(folder):
                if check(screenshot, screen['data'], threshold):
                    return screen['file']

def click(image: str, timeout, left=True):
    """Click the location of the image (if found)

    Args:
        image (str): Image path to click
        timeout (int or float): Time in seconds to hold the mouse down

    Returns:
        bool: If the location could be found
    """

    try:
        click_location = pyautogui.locateOnScreen(image)
        x, y = pyautogui.center(click_location)

    except:
        return False

    else:
        mouse.move(x, y)
        mouse.press(button=mouse.LEFT if left else mouse.RIGHT)
        time.sleep(timeout)
        return True

max_stars_repo_name: nsde/intellicraft

max_issues_repo_name: nsde/intellicraft

max_forks_repo_name: nsde/intellicraft


text:
"""
Copyright 2020 The OneFlow Authors. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

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

The raw ImageNet data set is expected to reside in JPEG files located in the
following directory structure.

  data_dir/n01440764/ILSVRC2012_val_00000293.JPEG
  data_dir/n01440764/ILSVRC2012_val_00000543.JPEG
  ...

where 'n01440764' is the unique synset label associated with
these images.

The training data set consists of 1000 sub-directories (i.e. labels)
each containing 1200 JPEG images for a total of 1.2M JPEG images.

The evaluation data set consists of 1000 sub-directories (i.e. labels)
each containing 50 JPEG images for a total of 50K JPEG images.

  train_directory/part-00000
  train_directory/part-00001
  ...
  train_directory/part-01023

and

  validation_directory/part-00000
  validation_directory/part-00001
  ...
  validation_directory/part-01023

Each record within the OFRecord file is a
serialized Example proto. The Example proto contains the following fields:

  image/encoded: string containing JPEG encoded image in RGB colorspace
  image/height: integer, image height in pixels
  image/width: integer, image width in pixels
  image/colorspace: string, specifying the colorspace, always 'RGB'
  image/channels: integer, specifying the number of channels, always 3
  image/format: string, specifying the format, always 'JPEG'

  image/filename: string containing the basename of the image file
            e.g. 'n01440764_10026.JPEG' or 'ILSVRC2012_val_00000293.JPEG'
  image/class/label: integer specifying the index in a classification layer.
    The label ranges from [1, 1000] where 0 is not used.
  image/class/synset: string specifying the unique ID of the label,
    e.g. 'n01440764'
  image/class/text: string specifying the human-readable version of the label
    e.g. 'red fox, Vulpes vulpes'

  image/object/bbox/xmin: list of integers specifying the 0+ human annotated
    bounding boxes
  image/object/bbox/xmax: list of integers specifying the 0+ human annotated
    bounding boxes
  image/object/bbox/ymin: list of integers specifying the 0+ human annotated
    bounding boxes
  image/object/bbox/ymax: list of integers specifying the 0+ human annotated
    bounding boxes
  image/object/bbox/label: integer specifying the index in a classification
    layer. The label ranges from [1, 1000] where 0 is not used. Note this is
    always identical to the image label.

Note that the length of xmin is identical to the length of xmax, ymin and ymax
for each example.
"""

from datetime import datetime
import os
import random
import sys
import threading
import re
import argparse
import glob

import numpy as np
import six
import cv2
import oneflow.core.record.record_pb2 as of_record
import struct

"""
# train dataset to ofrecord
python3 imagenet_ofrecord.py \
--train_directory data/imagenet/train  \
--output_directory data/imagenet/ofrecord/train   \
--label_file imagenet_lsvrc_2015_synsets.txt   \
--shards 1024  --num_threads 8 --name train  \
--bounding_box_file imagenet_2012_bounding_boxes.csv   \
--height 224 --width 224

# val dataset to ofrecord
python3 imagenet_ofrecord.py \
--validation_directory data/imagenet/validation  \
--output_directory data/imagenet/ofrecord/validation  \
--label_file imagenet_lsvrc_2015_synsets.txt --name validation  \
--shards 32 --num_threads 4 --name validation \
--bounding_box_file imagenet_2012_bounding_boxes.csv  \
--height 224 --width 224
"""

arg_parser = argparse.ArgumentParser(description =
            'The python script to resize pics ')

arg_parser.add_argument('--resize', dest = 'resize', default = False, help = 'resize image')

arg_parser.add_argument('--name', dest='name', default='train', \
help = 'data_file_type')

arg_parser.add_argument('--width', dest='width', default=0, \
type=int, help='fixed image width')
arg_parser.add_argument('--height', dest='height', default=0, \
type=int, help='fixed image height')

arg_parser.add_argument('--train_directory', dest = 'train_directory',\
default='/tmp/', help='Training data directory')
arg_parser.add_argument('--validation_directory', dest = 'validation_directory', \
default='/tmp/', help='Validation data directory')
arg_parser.add_argument('--output_directory', dest = 'output_directory', \
default='/tmp/', help = 'Output data directory')

arg_parser.add_argument('--shards', dest='shards',\
default=1024, type=int, help='Number of shards in making OFRecord files.')

arg_parser.add_argument('--num_threads', dest='num_threads', default = 8, \
type=int, help='Number of threads to preprocess the images.')

# The labels file contains a list of valid labels are held in this file.
# Assumes that the file contains entries as such:
#   n01440764
#   n01443537
#   n01484850
# where each line corresponds to a label expressed as a synset. We map
# each synset contained in the file to an integer (based on the alphabetical
# ordering). See below for details.
arg_parser.add_argument('--label_file', dest = 'labels_file', \
default = 'imagenet_lsvrc_2015_synsets.txt', help = 'Labels file')

# This file containing mapping from synset to human-readable label.
# Assumes each line of the file looks like:
#
#   n02119247    black fox
#   n02119359    silver fox
#   n02119477    red fox, Vulpes fulva
#
# where each line corresponds to a unique mapping. Note that each line is
# formatted as <synset>\t<human readable label>.
arg_parser.add_argument('--imagenet_metadata_file', dest = 'imagenet_metadata_file', \
default = 'imagenet_metadata.txt', help = 'ImageNet metadata file')

# This file is the output of process_bounding_box.py
# Assumes each line of the file looks like:
#
#   n00007846_64193.JPEG,0.0060,0.2620,0.7545,0.9940
#
# where each line corresponds to one bounding box annotation associated
# with an image. Each line can be parsed as:
#
#   <JPEG file name>, <xmin>, <ymin>, <xmax>, <ymax>
#
# Note that there might exist mulitple bounding box annotations associated
# with an image file.
arg_parser.add_argument('--bounding_box_file', dest = 'bounding_box_file', \
default = './imagenet_2012_bounding_boxes.csv', help = 'Bounding box file')

ARGS = arg_parser.parse_args()

def _int32_feature(value):
  """Wrapper for inserting int32 features into Example proto."""
  if not isinstance(value, list):
    value = [value]
  return of_record.Feature(int32_list=of_record.Int32List(value=value))

def _float_feature(value):
  """Wrapper for inserting float features into Example proto."""
  if not isinstance(value, list):
    value = [value]
  return of_record.Feature(float_list=of_record.FloatList(value=value))

def _double_feature(value):
  """Wrapper for inserting float features into Example proto."""
  if not isinstance(value, list):
    value = [value]
  return of_record.Feature(double_list=of_record.DoubleList(value=value))

def _bytes_feature(value):
  """Wrapper for inserting bytes features into Example proto."""
  # if isinstance(value, six.string_types):
  #  value = six.binary_type(value, encoding='utf-8')
  return of_record.Feature(bytes_list=of_record.BytesList(value=[value]))

def _convert_to_example(filename, image_buffer, label, index, synset, human, bbox,
                        height, width):
  """Build an Example proto for an example.

  Args:
    filename: string, path to an image file, e.g., '/path/to/example.JPG'
    image_buffer: string, JPEG encoding of RGB image
    label: integer, identifier for the ground truth for the network
    synset: string, unique WordNet ID specifying the label, e.g., 'n02323233'
    human: string, human-readable label, e.g., 'red fox, Vulpes vulpes'
    bbox: list of bounding boxes; each box is a list of integers
      specifying [xmin, ymin, xmax, ymax]. All boxes are assumed to belong to
      the same label as the image label.
    height: integer, image height in pixels
    width: integer, image width in pixels
  Returns:
    Example proto
  """
  xmin = []
  ymin = []
  xmax = []
  ymax = []
  for b in bbox:
    assert len(b) == 4
    # pylint: disable=expression-not-assigned
    [l.append(point) for l, point in zip([xmin, ymin, xmax, ymax], b)]
    # pylint: enable=expression-not-assigned

  colorspace = 'RGB'
  channels = 3
  image_format = 'JPEG'

  example = of_record.OFRecord(feature={
      'data_id': _bytes_feature(str(index).encode('utf-8')),
      'height': _int32_feature(height),
      'width': _int32_feature(width),
      'colorspace': _bytes_feature(colorspace.encode('utf-8')),
      'channels': _int32_feature(channels),
      'class/label': _int32_feature(label),
      'class/synset': _bytes_feature(synset.encode('utf-8')),
      'class/text': _bytes_feature(human.encode('utf-8')),
      'object/bbox/xmin': _float_feature(xmin),
      'object/bbox/xmax': _float_feature(xmax),
      'object/bbox/ymin': _float_feature(ymin),
      'object/bbox/ymax': _float_feature(ymax),
      'object/bbox/label': _int32_feature([label] * len(xmin)),
      'format': _bytes_feature(image_format.encode('utf-8')),
      'filename': _bytes_feature(os.path.basename(filename).encode('utf-8')),
      'encoded': _bytes_feature(image_buffer)})
  return example

class ImageCoder(object):
  """Helper class that provides image coding utilities."""

  def __init__(self,size = None):
    self.size = size

  def _resize(self, image_data):
    if self.size != None and image_data.shape[:2] != self.size:
      return cv2.resize(image_data, self.size)
    return image_data

  def image_to_jpeg(self, image_data, resize=ARGS.resize):
    # image_data = cv2.imdecode(np.fromstring(image_data, np.uint8), 1) # deprecated,
    image_data = cv2.imdecode(np.frombuffer(image_data, np.uint8), 1)
    if resize:
      image_data = self._resize(image_data)
    return cv2.imencode(".jpg", image_data)[1].tobytes(), image_data.shape[0], image_data.shape[1]

def _process_image(filename, coder):
  """Process a single image file.

  Args:
    filename: string, path to an image file e.g., '/path/to/example.JPG'.
    coder: instance of ImageCoder to provide image coding utils.
  Returns:
    image_buffer: string, JPEG encoding of RGB image.
    height: integer, image height in pixels.
    width: integer, image width in pixels.
  """
  # Read the image file.
  with open(filename, 'rb') as f:
    image_data = f.read()
    image_data,height, width =  coder.image_to_jpeg(image_data)
    #print(height, width)
    return image_data,height, width
    # Decode the RGB JPEG.
    # image_data = coder._resize(image_data)

def _process_image_files_batch(coder, thread_index, ranges, name, filenames,
                               synsets, labels, indexs, humans, bboxes, num_shards):
  """Processes and saves list of images as OFRecord in 1 thread.

  Args:
    coder: instance of ImageCoder to provide image coding utils.
    thread_index: integer, unique batch to run index is within [0, len(ranges)).
    ranges: list of pairs of integers specifying ranges of each batches to
      analyze in parallel.
    name: string, unique identifier specifying the data set
    filenames: list of strings; each string is a path to an image file
    synsets: list of strings; each string is a unique WordNet ID
    labels: list of integer; each integer identifies the ground truth
    humans: list of strings; each string is a human-readable label
    bboxes: list of bounding boxes for each image. Note that each entry in this
      list might contain from 0+ entries corresponding to the number of bounding
      box annotations for the image.
    num_shards: integer number of shards for this data set.
  """
  # Each thread produces N shards where N = int(num_shards / num_threads).
  # For instance, if num_shards = 128, and the num_threads = 2, then the first
  # thread would produce shards [0, 64).
  num_threads = len(ranges)
  assert not num_shards % num_threads
  num_shards_per_batch = int(num_shards / num_threads)

  shard_ranges = np.linspace(ranges[thread_index][0],
                             ranges[thread_index][1],
                             num_shards_per_batch + 1).astype(int)
  num_files_in_thread = ranges[thread_index][1] - ranges[thread_index][0]

  counter = 0

  for s in range(num_shards_per_batch):
    # Generate a sharded version of the file name, e.g. 'train-00002-of-00010'
    shard = thread_index * num_shards_per_batch + s
    # output_filename = '%s-%.5d-of-%.5d' % (name, shard, num_shards)
    output_filename = 'part-%.5d' % (shard)
    output_file = os.path.join(ARGS.output_directory, output_filename)

    f = open(output_file, 'wb')
    shard_counter = 0
    files_in_shard = np.arange(shard_ranges[s], shard_ranges[s + 1], dtype=int)
    for i in files_in_shard:
      filename = filenames[i]
      label = labels[i]
      synset = synsets[i]
      human = humans[i]
      bbox = bboxes[i]
      index = indexs[i]
      try:
        image_buffer, height, width = _process_image(filename, coder)
      except:
        print(filename)
        continue

      # print('filename, label, index,synset, human, bbox,height, width >>>>>>>>>>>>>>>>>>>>>>>>>>>>\n',
      #       filename , label, index,synset, human, bbox,height, width)
      example = _convert_to_example(filename, image_buffer, label, index,
                                    synset, human, bbox,
                                    height, width)
      l = example.ByteSize()
      f.write(struct.pack("q", l))
      f.write(example.SerializeToString())
      shard_counter += 1
      counter += 1

      if not counter % 1000:
        print('%s [thread %d]: Processed %d of %d images in thread batch.' %
              (datetime.now(), thread_index, counter, num_files_in_thread))
        sys.stdout.flush()

    f.close()
    print('%s [thread %d]: Wrote %d images to %s' %
          (datetime.now(), thread_index, shard_counter, output_file))
    sys.stdout.flush()
    shard_counter = 0
  print('%s [thread %d]: Wrote %d images to %d shards.' %
        (datetime.now(), thread_index, counter, num_files_in_thread))
  sys.stdout.flush()

def _process_image_files(name, filenames, synsets, labels, indexs, humans,
                         bboxes, num_shards):
  """Process and save list of images as OFRecord of Example protos.

  Args:
    name: string, unique identifier specifying the data set
    filenames: list of strings; each string is a path to an image file
    synsets: list of strings; each string is a unique WordNet ID
    labels: list of integer; each integer identifies the ground truth
    humans: list of strings; each string is a human-readable label
    bboxes: list of bounding boxes for each image. Note that each entry in this
      list might contain from 0+ entries corresponding to the number of bounding
      box annotations for the image.
    num_shards: integer number of shards for this data set.
  """
  assert len(filenames) == len(synsets)
  assert len(filenames) == len(labels)
  assert len(filenames) == len(humans)
  assert len(filenames) == len(bboxes)

  # Break all images into batches with a [ranges[i][0], ranges[i][1]].
  spacing = np.linspace(0, len(filenames), ARGS.num_threads + 1).astype(np.int)
  ranges = []
  threads = []
  for i in range(len(spacing) - 1):
    ranges.append([spacing[i], spacing[i + 1]])
  # Launch a thread for each batch.
  print('Launching %d threads for spacings: %s' % (ARGS.num_threads, ranges))
  sys.stdout.flush()

  # Create a generic utility for converting all image codings.
  if ARGS.width <= 0 or ARGS.height <= 0:
    coder = ImageCoder()
  else:
    coder = ImageCoder((ARGS.width, ARGS.height))

  threads = []
  for thread_index in range(len(ranges)):
    args = (coder, thread_index, ranges, name, filenames,
            synsets, labels, indexs, humans, bboxes, num_shards)
    t = threading.Thread(target=_process_image_files_batch, args=args)
    t.start()
    threads.append(t)

  # Wait for all the threads to terminate.
  for t in threads:
    t.join()
  print('%s: Finished writing all %d images in data set.' %
        (datetime.now(), len(filenames)))
  sys.stdout.flush()

def _find_image_files(data_dir, labels_file):
  """Build a list of all images files and labels in the data set.

  Args:
    data_dir: string, path to the root directory of images.

      Assumes that the ImageNet data set resides in JPEG files located in
      the following directory structure.

        data_dir/n01440764/ILSVRC2012_val_00000293.JPEG
        data_dir/n01440764/ILSVRC2012_val_00000543.JPEG

      where 'n01440764' is the unique synset label associated with these images.

    labels_file: string, path to the labels file.

      The list of valid labels are held in this file. Assumes that the file
      contains entries as such:
        n01440764
        n01443537
        n01484850
      where each line corresponds to a label expressed as a synset. We map
      each synset contained in the file to an integer (based on the alphabetical
      ordering) starting with the integer 1 corresponding to the synset
      contained in the first line.

      The reason we start the integer labels at 1 is to reserve label 0 as an
      unused background class.

  Returns:
    filenames: list of strings; each string is a path to an image file.
    synsets: list of strings; each string is a unique WordNet ID.
    labels: list of integer; each integer identifies the ground truth.
  """
  print('Determining list of input files and labels from %s.' % data_dir)
  challenge_synsets = [l.strip() for l in
                       open(labels_file, 'r').readlines()]

  labels = []
  filenames = []
  synsets = []

  # Leave label index 0 empty as a background class.
  label_index = 0# use to be 1

  # Construct the list of JPEG files and labels.
  for synset in challenge_synsets:
    if not os.path.exists(os.path.join(data_dir, synset)):
      continue

    jpeg_file_path = '%s/%s/*.JPEG' % (data_dir, synset)
    matching_files = glob.glob(jpeg_file_path)

    labels.extend([label_index] * len(matching_files))
    synsets.extend([synset] * len(matching_files))
    filenames.extend(matching_files)

    if not label_index % 100:
      print('Finished finding files in %d of %d classes.' % (
          label_index, len(challenge_synsets)))
    label_index += 1
  # Shuffle the ordering of all image files in order to guarantee
  # random ordering of the images with respect to label in the
  # saved OFRecord files. Make the randomization repeatable.
  shuffled_index = list(range(len(filenames)))
  random.seed(12345)
  random.shuffle(shuffled_index)

  filenames = [filenames[i] for i in shuffled_index]
  synsets = [synsets[i] for i in shuffled_index]
  labels = [labels[i] for i in shuffled_index]

  print('Found %d JPEG files across %d labels inside %s.' %
        (len(filenames), len(challenge_synsets), data_dir))
  return filenames, synsets, labels, shuffled_index

def _find_human_readable_labels(synsets, synset_to_human):
  """Build a list of human-readable labels.

  Args:
    synsets: list of strings; each string is a unique WordNet ID.
    synset_to_human: dict of synset to human labels, e.g.,
      'n02119022' --> 'red fox, Vulpes vulpes'

  Returns:
    List of human-readable strings corresponding to each synset.
  """
  humans = []
  for s in synsets:
    assert s in synset_to_human, ('Failed to find: %s' % s)
    humans.append(synset_to_human[s])
  return humans

def _find_image_bounding_boxes(filenames, image_to_bboxes):
  """Find the bounding boxes for a given image file.

  Args:
    filenames: list of strings; each string is a path to an image file.
    image_to_bboxes: dictionary mapping image file names to a list of
      bounding boxes. This list contains 0+ bounding boxes.
  Returns:
    List of bounding boxes for each image. Note that each entry in this
    list might contain from 0+ entries corresponding to the number of bounding
    box annotations for the image.
  """
  num_image_bbox = 0
  bboxes = []
  for f in filenames:
    basename = os.path.basename(f)
    if basename in image_to_bboxes:
      bboxes.append(image_to_bboxes[basename])
      num_image_bbox += 1
    else:
      bboxes.append([])
  print('Found %d images with bboxes out of %d images' % (
      num_image_bbox, len(filenames)))
  return bboxes

def _process_dataset(name, directory, num_shards, synset_to_human,
                     image_to_bboxes):
  """Process a complete data set and save it as a OFRecord.

  Args:
    name: string, unique identifier specifying the data set.
    directory: string, root path to the data set.
    num_shards: integer number of shards for this data set.
    synset_to_human: dict of synset to human labels, e.g.,
      'n02119022' --> 'red fox, Vulpes vulpes'
    image_to_bboxes: dictionary mapping image file names to a list of
      bounding boxes. This list contains 0+ bounding boxes.
  """
  filenames, synsets, labels, indexs = _find_image_files(directory, ARGS.labels_file)
  # ./train/n03085013/n03085013_23287.JPEG n03085013 508 652481

  humans = _find_human_readable_labels(synsets, synset_to_human)
  bboxes = _find_image_bounding_boxes(filenames, image_to_bboxes)

  _process_image_files(name, filenames, synsets, labels, indexs,
                       humans, bboxes, num_shards)

def _build_synset_lookup(imagenet_metadata_file):
  """Build lookup for synset to human-readable label.

  Args:
    imagenet_metadata_file: string, path to file containing mapping from
      synset to human-readable label.

      Assumes each line of the file looks like:

        n02119247    black fox
        n02119359    silver fox
        n02119477    red fox, Vulpes fulva

      where each line corresponds to a unique mapping. Note that each line is
      formatted as <synset>\t<human readable label>.

  Returns:
    Dictionary of synset to human labels, such as:
      'n02119022' --> 'red fox, Vulpes vulpes'
  """
  lines = open(imagenet_metadata_file, 'r').readlines()
  synset_to_human = {}
  for l in lines:
    if l:
      parts = l.strip().split('\t')
      assert len(parts) == 2
      synset = parts[0]
      human = parts[1]
      synset_to_human[synset] = human
  return synset_to_human

def _build_bounding_box_lookup(bounding_box_file):
  """Build a lookup from image file to bounding boxes.

  Args:
    bounding_box_file: string, path to file with bounding boxes annotations.

      Assumes each line of the file looks like:

        n00007846_64193.JPEG,0.0060,0.2620,0.7545,0.9940

      where each line corresponds to one bounding box annotation associated
      with an image. Each line can be parsed as:

        <JPEG file name>, <xmin>, <ymin>, <xmax>, <ymax>

      Note that there might exist mulitple bounding box annotations associated
      with an image file. This file is the output of process_bounding_boxes.py.

  Returns:
    Dictionary mapping image file names to a list of bounding boxes. This list
    contains 0+ bounding boxes.
  """
  lines = open(bounding_box_file, 'r').readlines()
  images_to_bboxes = {}
  num_bbox = 0
  num_image = 0
  for l in lines:
    if l:
      parts = l.split(',')
      assert len(parts) == 5, ('Failed to parse: %s' % l)
      filename = parts[0]
      xmin = float(parts[1])
      ymin = float(parts[2])
      xmax = float(parts[3])
      ymax = float(parts[4])
      box = [xmin, ymin, xmax, ymax]

      if filename not in images_to_bboxes:
        images_to_bboxes[filename] = []
        num_image += 1
      images_to_bboxes[filename].append(box)
      num_bbox += 1

  print('Successfully read %d bounding boxes '
        'across %d images.' % (num_bbox, num_image))
  return images_to_bboxes

def main():
  assert not ARGS.shards % ARGS.num_threads, (
      'Please make the ARGS.num_threads commensurate with ARGS.shards')

  print('Saving results to %s' % ARGS.output_directory)
  if not os.path.exists(ARGS.output_directory):
    os.makedirs(ARGS.output_directory)

  # Build a map from synset to human-readable label.
  synset_to_human = _build_synset_lookup(ARGS.imagenet_metadata_file)
  image_to_bboxes = _build_bounding_box_lookup(ARGS.bounding_box_file)

  # Run it!
  if ARGS.name == 'validation':
    _process_dataset('validation', ARGS.validation_directory,
                   ARGS.shards, synset_to_human, image_to_bboxes)
  else:
    _process_dataset('train', ARGS.train_directory, ARGS.shards,
                   synset_to_human, image_to_bboxes)

if __name__ == '__main__':
  main()


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 10 12:32:31 2019

@author: iagorosa
"""

#%%
import numpy as np
from PIL import Image as pil
#import matplotlib.pyplot as plt
from os import listdir
from os.path import isfile, join
from sklearn_extensions.extreme_learning_machines.elm import ELMRegressor, ELMClassifier

from sklearn.preprocessing import MinMaxScaler, MaxAbsScaler

import ClassELM as elm
import aux_functions as af

#%%

X = np.loadtxt("./Yale_32x32.csv", delimiter=',', skiprows=1)
# col 0: pessoa
# col 1: label

#%%

# Para reconstruir a imagem, precisa passar X para uint8

#for i in range(11):
#    m0 = X[10+i*11, 2:].reshape(32,32).T
#    img = pil.fromarray(m0, 'L')
#    img.show()

'''
for x in X[X[:, 1] == 8]:
#    m0 = X[10+i*11, 2:].reshape(32,32).T
    m0 = x[2:].reshape(32,32).T
    img = pil.fromarray(m0, 'L')
    img.show()
'''

'''
labels = {
         'centerlight':  0,
         'glasses':      1,
         'happy':        2,
         'leftlight':    3,
         'noglasses':    4,
         'normal':       5,
         'rightlight':   6,
         'sad':          7,
         'sleepy':       8,
         'surprised':    9,
         'wink':        10
         }
'''

#%%

#scaler = MinMaxScaler()
scaler = MaxAbsScaler()

X_ = scaler.fit_transform(X[:, 2:].T)

X[:, 2:] = X_.T

#%%

###### DEFINIÇÕES DE TREINO E TESTE ##############

v = 1

X_train = X[X[:, 0] != v][:, 2:]
y_train = X[X[:, 0] != v][:,  1]

X_test  = X[X[:, 0] == v][:, 2:]
y_test  = X[X[:, 0] == v][:,  1]

y_train_real = y_train.copy()

#%%
#### ELM DEFINIDO PELO SKLEARN

clsf = ELMClassifier(random_state=100)
#clsf.fit(X[:v, 2:], y[:v])
clsf.fit(X_train, y_train)

Htreino = clsf._get_weights()

print()

print(clsf.predict(X_test))

#print()

#print(y[v:])

print('Score:', clsf.score(X_test, y_test))

Hteste = clsf._get_weights()

#%%

################# ELM IMPLENTÇÃO DO KAGGLE  ################
# https://www.kaggle.com/robertbm/extreme-learning-machine-example #

# funções

def relu(a):
    return np.maximum(a, 0, a)

def input_to_hidden(x, W, B, bias = True): # produto pesos de entrada com neuronios de pesos aleatorios + função de ativação
    if bias:
        B_mat = np.tile(B, (x.shape[0], 1)) # repete a linha varias vezes
        a = np.dot(x, W) + B_mat
    else:
        a = np.dot(x, W)
    a = relu(a) # ReLU
    return a

def predict(x, W, beta, B, bias):
    x = input_to_hidden(x, W, B, bias)
    y = np.dot(x, beta)
    return y

#%%

try:
    y_train.shape[1]
except:
    labels = y_train.astype('int')
    CLASSES = len(np.unique(labels))
    y_train = np.zeros([labels.shape[0], CLASSES])
    for i in range(labels.shape[0]):
            y_train[i][labels[i]] = 1
#y_train.view(type=np.matrix)

#%%

INPUT_LENGHT = X_train.shape[1] # 784
HIDDEN_UNITS = 20

#np.random.seed(20)

W = np.random.normal(size=[INPUT_LENGHT, HIDDEN_UNITS],)
print('Input Weight shape: {shape}'.format(shape=W.shape))

B = np.random.normal(size=[HIDDEN_UNITS]) #Bias

#%%

### MODIFICACAO

def weight_func(X, W, d):
    v = np.average(X)
    const = ( 1.3 / (np.sqrt(1 + d*v**2)) )
    return const * W

#%%

### MODIFICACAO

Wuni = np.random.uniform(-1, 1, size=(INPUT_LENGHT, HIDDEN_UNITS))

W = weight_func(X_train, Wuni, HIDDEN_UNITS)

#%%

bias = False

H = input_to_hidden(X_train, W, B, bias)
Ht = np.transpose(H)
#beta = np.dot(np.linalg.inv(np.dot(Ht, H)), np.dot(Ht, y_train))
beta = np.dot(np.linalg.pinv(np.dot(Ht, H)), np.dot(Ht, y_train))

print('Output weights shape: {shape}'.format(shape=beta.shape))

#%%

y = predict(X_test, W, beta, B, bias)
correct = 0
total = y.shape[0]
r_test = []
for i in range(total):
    predicted = np.argmax(y[i])
#    test = np.argmax(y_test[i])
    test = y_test[i]
    r_test.append(predicted)
    correct = correct + (1 if predicted == test else 0)
print('Accuracy: {:f}'.format(correct/total))
print(r_test)

#%%

def inverse_ReLU(x):
    xx = np.where(x > 0, x, 0)
    pos = np.where(x <= 0, 1, 0)
    rand = np.random.uniform(-1, 0, size=x.shape)
    aux = pos * rand

    return pos, rand, aux + xx

##### MULTICAMADAS

# necessario calcular o beta na parte anterior

# H1 = T beta*
H1 = np.dot(y_train, np.linalg.pinv(beta))

# He = [H 1s], onde 1s eh vetor linha de 1
He = np.concatenate([H, np.ones(H.shape[0]).reshape(-1, 1)], axis=1)

#Whe = np.dot(np.linalg.pinv(He), inverse_ReLU(H1))

pos, rand, H1i = inverse_ReLU(H1)
Whe = np.dot(np.linalg.pinv(He), H1i)

Wh = Whe[:-1, :]
B1  = Whe[-1, :]

H2 = input_to_hidden(H, Wh, B1, bias)
beta2 = np.dot(np.linalg.pinv(H2), y_train)

#%%

#saida

p = np.dot(input_to_hidden(X_test, W, B), Wh) + B1
f = np.dot(relu(p), beta2)

#%%

#y = predict(X_test, W, f, B)
y = f
correct = 0
total = y.shape[0]
r_test = []
for i in range(total):
    predicted = np.argmax(y[i])
#    test = np.argmax(y_test[i])
    test = y_test[i]
    r_test.append(predicted)
    correct = correct + (1 if predicted == test else 0)
print('Accuracy: {:f}'.format(correct/total))
print(r_test)

#%%

celm = elm.ELMClassifier(activation_func='relu', func_hidden_layer = af.normal_random_layer, bias=True, random_state=10)

celm.fit(X_train, y_train)

celm.predict(X_test, y_test)

#%%

cmelm = elm.ELMMLPClassifier(activation_func='relu', func_hidden_layer = af.normal_random_layer, random_state=10)

cmelm.fit(X_train, y_train)

cmelm.predict(X_test, y_test)

#%%

max_stars_repo_path: em desuso/main.py

max_issues_repo_path: em desuso/main.py

max_forks_repo_path: em desuso/main.py


text:
# -*- coding: utf-8 -*-

import sys
import numpy as np
from multiprocessing import Pool
import os
import math
import argparse
from dateutil.relativedelta import relativedelta
from datetime import datetime
from configobj import ConfigObj

from lib.gif_lib_image import prj_gll, get_file_time
from lib.gif_lib_image import find_file, time_block, run_command
from lib.gif_lib_image import read_json, image_plot_title, image_plot_shp
# import matplotlib
# matplotlib.use('Agg')
import matplotlib.pyplot as plt

# 获取程序所在目录位置
g_path, _ = os.path.split(os.path.realpath(__file__))

# 模块的配置文件
g_file_cfg = os.path.join(g_path, 'cfg/gif.cfg')

# 读取模块配置文件内容
g_var_cfg = ConfigObj(g_file_cfg)
g_convert = '/bin/convert'
g_mogrify = '/bin/mogrify'

threads = int(g_var_cfg['CROND']['threads'])

def find_time_range_file(ipath, t1, t2, reg):

    date_start = datetime.strptime(t1, '%Y%m%d%H%M%S')
    date_end = datetime.strptime(t2, '%Y%m%d%H%M%S')

    date_s = date_start
    date_e = date_end
    file_list = []
    while date_s <= date_e:
        ymd = date_s.strftime('%Y%m%d')
        full_in_path = os.path.join(ipath, ymd)
        file_list_ymd = find_file(full_in_path, '.*')

        for in_file in file_list_ymd:
            file_name = os.path.basename(in_file)
            file_time = get_file_time(file_name, reg)

            if date_start <= file_time <= date_end:
                file_list.append(in_file)

        date_s = date_s + relativedelta(days=1)

    return file_list

def image_crop(file_list, lon_w, lon_e, lat_s, lat_n, opath):
    """
    抠图
    """
    latlon = g_var_cfg['PAIRS']['CLR_1000M']['latlon']
    rowcol = g_var_cfg['PAIRS']['CLR_1000M']['rowcol']
    row, col = [int(each) for each in rowcol]
    latlon = [float(each) for each in latlon]
    g_lon_w, g_lon_e, g_lat_s, g_lat_n = latlon

    proj = prj_gll(g_lat_n, g_lat_s, g_lon_w, g_lon_e, None, None, row, col)
    start_i, start_j = proj.lonslats2ij(lon_w, lat_n)
    end_i, end_j = proj.lonslats2ij(lon_e, lat_s)

    width = end_j - start_j
    height = end_i - start_i

    print 'width*height', width, height, start_j, start_i
    size = '%dx%d+%d+%d' % (width, height, start_j, start_i)

    cmd_list = []

    for in_file in file_list:
        file_name = os.path.basename(in_file)
        out_file = os.path.join(opath, file_name)
        cmd = '%s -crop %s %s %s' % (g_convert, size, in_file, out_file)
        print cmd
#         if os.path.isfile(out_file):
#             continue
        cmd_list.append(cmd)

    run_command(cmd_list, threads)

def imgage_resize(file_list, percent):
    """
    重置图像大小
    """

    cmd_list = []
    for in_file in file_list:
        cmd = '%s -resize %d%%x%d%% %s %s' % (g_convert,
                                              percent, percent, in_file, in_file)
        cmd_list.append(cmd)

    if len(cmd_list) > 0:
        run_command(cmd_list, threads)

def image_mogrify(file_list, x, y, out_shp_path):
    """
    叠加图片
    """

    # convert  -fill white -pointsize 24 -draw "text 10,15 'lifesinger 2006' "
    # 1.png  2.png
    # convert   -fill white -pointsize 24 -font
    # /home/gsics/Project/OM/hz/font/msyh.ttf  -draw "text 1950,50 '哈哈哈'"
    # H8XX_AHIXX_L2_GLL_20190421_0000_1000M_CLR_GeoColor.jpg 2.jpg

    # 获取字体文件

    # 获取数据所在位置 cd 过去
    os.chdir(out_shp_path)

    # 标题
    cmd_list = []
    for in_file in file_list:
        file_name = os.path.basename(in_file)
        file_name_shpae = file_name.split('.')[0] + '.png'
        # cmd = '''%s -fill white -pointsize 20  -stroke black -font %s
        # -draw "text %d,%d '%s'" -draw 'image SrcOver 0,0 %d,%d %s' %s''' % (
        # g_mogrify, font_file, title_y, title_x, new_title, img_y, img_x,
        # shp_file, file_name)

        cmd = '''%s -draw 'image SrcOver 0,0 %d,%d %s' %s''' % (
            g_mogrify,  y, x, file_name_shpae, in_file)
        cmd_list.append(cmd)

    run_command(cmd_list, threads)

    # shp边界
    cmd_list = []
    for in_file in file_list:
        cmd = '''%s -draw 'image SrcOver 0,0 %d,%d %s' %s''' % (
            g_mogrify,  y, x, 'overlap.png', in_file)
        cmd_list.append(cmd)

    run_command(cmd_list, threads)

def image_gif(file_list, out_gif_file):
    """
    制作动画
    """

    cmd_list = []
    out_file_list = []
    out_path = os.path.dirname(out_gif_file)
    if not os.path.isdir(out_path):
        os.makedirs(out_path)

    for in_file in file_list:
        file_name = os.path.basename(in_file)

        out_file = os.path.join(out_path, file_name.split('.')[0] + '.gif')
        out_file_list.append(out_file)
#         if os.path.isfile(out_file):
#             continue
        cmd = '%s -delay 25 %s  %s' % (g_convert, in_file, out_file)
        cmd_list.append(cmd)

    if len(cmd_list) > 0:
        run_command(cmd_list, threads)

    os.system('%s -loop 0 %s %s ' %
              (g_convert, ' '.join(out_file_list), out_gif_file))

if __name__ == '__main__':

    ipath = g_var_cfg['PATH']['IN']['h8_jpg']
    mpath = g_var_cfg['PATH']['MID']['h8_crop']
    opath = g_var_cfg['PATH']['OUT']['h8_gif']

    reg = g_var_cfg['PAIRS']['CLR_1000M']['reg']
    json_file = g_var_cfg['PAIRS']['CLR_1000M']['json']
    title = g_var_cfg['PAIRS']['CLR_1000M']['title']
    typhoon = g_var_cfg['PAIRS']['CLR_1000M']['typhoon']

    # 输入
    parser = argparse.ArgumentParser()
    parser.add_argument("-t1", "--time1", help="YYYYMMDDHHMMSS")
    parser.add_argument("-t2", "--time2", help="YYYYMMDDHHMMSS")
    parser.add_argument(
        "-A", "--area", help="lat,lon (Left upper) lat,lon (lower right)")
    parser.add_argument(
        "-N", "--name", help="area name")
    args = parser.parse_args()
    t1 = args.time1
    t2 = args.time2
    latlon = [float(each) for each in args.area.split(',')]
    lon_w, lon_e, lat_s, lat_n = latlon
    name = unicode(args.name)

    print 'args', t1, t2, latlon, name

    # 0.  找到符合条件的文件
    file_list = find_time_range_file(ipath, t1, t2, reg)

    # 1. 开始从原图把需要的区域扣图

    crop_opath = os.path.join(mpath, name, '%s_%s' % (t1, t2))
    if not os.path.isdir(crop_opath):
        os.makedirs(crop_opath)

    with time_block(u'抠图耗时', True):
        image_crop(file_list, lon_w, lon_e, lat_s, lat_n, crop_opath)
        pass

    # 2. 重新定义大小 判断是否缩小 100万像素是分界点

    file_list = find_file(crop_opath, '.*.jpg')
    in_file = file_list[0]
    img = plt.imread(in_file)
    img_x, img_y = img.shape[:2]

    dividing_point = 1000000.
    all_point = img_x * img_y

    if all_point > dividing_point:
        scale = math.ceil((dividing_point / all_point) * 100)
        if scale <= 30.:
            scale = 30
        with time_block(u'等比例缩小 %d' % scale, True):
            imgage_resize(file_list, scale)
            pass

    # 缩放完成重新取大小
    in_file = file_list[0]
    img = plt.imread(in_file)
    img_x, img_y = img.shape[:2]
    plt.close()

    print img_x, img_y

    # 2. 开始绘制边界
    shp_path = os.path.join(crop_opath, 'shp')
    if not os.path.isdir(shp_path):
        os.makedirs(shp_path)
    jcfg = read_json(json_file)
    pool = Pool(processes=int(threads))
    with time_block(u"制作标题png", switch=True):
        for in_file in file_list:
            file_name = os.path.basename(in_file)
            out_file = os.path.join(
                crop_opath, 'shp', file_name.split('.')[0] + '.png')
            tt = get_file_time(file_name, reg)
            tt = tt + relativedelta(hours=8)
            ymd_hm = tt.strftime('%Y-%m-%d %H:%M')
            format_title = title % ymd_hm
#             print img_x, img_y, latlon, format_title, out_file
            pool.apply_async(
                image_plot_title, (img_x, img_y, latlon, format_title, jcfg, out_file))

        pool.close()
        pool.join()

    # 2. 开始绘制边界

    with time_block(u"制作矢量边界", switch=True):
        out_shp_file = os.path.join(crop_opath, 'shp', 'overlap.png')
        image_plot_shp(img_x, img_y, latlon, jcfg, typhoon, out_shp_file)

    # 3. 贴图
    out_shp_path = os.path.join(crop_opath, 'shp')
    with time_block(u'贴图', True):
        image_mogrify(file_list, img_x, img_y, out_shp_path)
        pass

    # 4. 动画
    os.chdir(g_path)
    out_file = os.path.join(
        opath, name, '%s_%s' % (t1, t2), '%s_%s.gif' % (t1, t2))
    with time_block(u'动画', True):
        image_gif(file_list, out_file)
        pass


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

def filter_TDP(data_frame, thresh = 0.3):
    """[optional function that removes molecules that do not transition below threshold at some point]

    Args:
        data_frame ([dataframe]): [dataframe that has been cleaned to remove outliers by 'remove_outliers' function]
        thresh (float, optional): [FRET value to threshold - will filter molecules and keep those that go below the thresh]. Defaults to 0.3.

    Returns:
        [dataframe]: [contains only molecules of interest]
    """
    filtered_mol = []
    for treatment, df in data_frame.groupby("treatment_name"):
        mol_list = df[(df["FRET_before"] <= thresh)|(df["FRET_after"] <= thresh)].Molecule.unique().tolist()
        filtered = df[df["Molecule"].isin(mol_list)]
        filtered_mol.append(filtered)
    filtered_mol = pd.concat(filtered_mol)
    return filtered_mol

def remove_outliers(compiled, plot_type, data_type = "raw"):
    """[removes outliers from dataframe]

    Args:
        compiled ([dataframe]): [raw dataframe containing outliers to be removed]
        plot_type ([str]): [string can either be 'hist' for histogram data or 'TDP' for TDP data]
        data_type (str, optional): [removes either raw FRET values or 'idealized' FRET values]. Defaults to "raw".

    Returns:
        [dataframe]: [returns cleaned data without outliers]
    """
    if plot_type == 'hist':
        if data_type == "raw":
            rawFRET = compiled[(compiled[3] > -0.5) & (compiled[3] < 1.5)].copy()
            return rawFRET
        if data_type == "idealized":
            idealizedFRET = compiled[(compiled[4] > -0.5) & (compiled[4] < 1.5)].copy()
            return idealizedFRET
    elif plot_type == 'TDP':
        outliers = compiled[(compiled["FRET before transition"] < -0.5)|(compiled["FRET before transition"] > 1.5)|(compiled["FRET after transition"] < -0.5) | (compiled["FRET after transition"] > 1.5)].index
        compiled.drop(outliers, inplace = True)
        return compiled
    else:
        print('invalid plot type, please set plot_type as "hist" or "TDP" - you idiot')

def cleanup_dwell(data, fps, thresh, first_dwell = "delete"):
    """[Will convert the data frome frame number to unit of time (seconds) and then delete all dwell times
    that are smaller than the set threshold (defined previously) in seconds. Will also delete the first dwell
    state from each molecule]

    Args:
        data ([dataframe]): [raw data]
        first_dwell (str, optional): [Set to 'keep' to keep the first dwell state from each molecule otherwise
        Will delete the first dwell state from each molecule by default]. Defaults to "delete".

    Returns:
        [dataframe]: [Data is now cleaned and ready for subsequent processing]
    """
    if first_dwell == "delete":
        filtered = []
        for molecule, df in data.groupby("Molecule"):
            filtered.append(df.iloc[1:])
        filtered = pd.concat(filtered)  #####filtered = pd.concat([df.iloc[1:] for molecule, df in A.groupby("Molecule")]) ##code here is the same as the for loop but in a list comprehension format
        filtered["Time (s)"] = filtered["Time"]/fps
        filtered = filtered[filtered["Time (s)"] >= thresh]
        return filtered
    if first_dwell == "keep":
        data["Time (s)"] = data["Time"]/fps
        data = data[data["Time (s)"] >= thresh]
        return data

def filter_dwell(df, FRET_thresh, headers):
    """[Will take the cleaned TDP data and will filter it using a threshold (defined by FRET_thresh)
    into seperate types of transitions (e.g., < 0.5 to > 0.5 FRET if FRET_thresh is = 0.5 is one example
    of a type of transition).

    Args:
        df ([dataframe]): [contains cleaned data that has been processed using the 'cleanup_dwell' function]

    Returns:
        [dataframe]: [contains dwell time that has been categorized into each transition class]
    """
    filtered_lowtohigh = df[(df["FRET_before"] < FRET_thresh) & (df["FRET_after"] > FRET_thresh)].copy()
    filtered_lowtolow = df[(df["FRET_before"] < FRET_thresh) & (df["FRET_after"] < FRET_thresh)].copy()
    filtered_hightolow = df[(df["FRET_before"] > FRET_thresh) & (df["FRET_after"] < FRET_thresh)].copy()
    filtered_hightohigh = df[(df["FRET_before"] > FRET_thresh) & (df["FRET_after"] > FRET_thresh)].copy()
    dataf = [filtered_lowtolow["Time (s)"], filtered_lowtohigh["Time (s)"], filtered_hightohigh["Time (s)"], filtered_hightolow["Time (s)"]]
    df_col = pd.concat(dataf, axis = 1, keys = headers)
    df_col = df_col.apply(lambda x:pd.Series(x.dropna().values))  ## removes NaN values from each column in df_col
    return df_col

def transition_frequency(filt):
    """calculates the transition frequency (i.e., the number of transitions per transition class divided
    by the total number of transitions observed). For example if there are 40 transitions total, and a
    < 0.5 to > 0.5 transition occurs 10 times, then the transition probability for that transition type is
    0.25 or 25%.

    Args:
        filt (dataframe): contains the dataframe with filtered data (cleaned data has been filtered by
        'filter_dwell' function)

    Returns:
        dataframe: returns a dataframe containing the percentage for each transition type
    """
    count_df_col = pd.DataFrame(filt.count(axis = 0)).transpose()
    count_df_col["sum"] = count_df_col.sum(axis = 1)
    dwell_frequency = pd.DataFrame([(count_df_col[column]/count_df_col["sum"])*100 for column in count_df_col]).transpose()
    print(dwell_frequency)
    return dwell_frequency

def calculate_mean(filtered_data, treatment_name):
    """calculates the mean dwell time of each type of transition class

    Args:
        filtered_data (dataframe): dataframe generated after the 'cleanup_dwell' and 'filter_dwell' functions
        have been run
        treatment_name (str): not required, only present to receive input from for loop. set to treatment_name
    Returns:
        [dataframe]: returns dataframe containing the mean of each transition class
    """
    mean_dwell = pd.DataFrame([filtered_data.iloc[0:].mean()])
    mean_dwell["sample"] = treatment_name
    return mean_dwell

def float_generator(data_frame, treatment, FRET_thresh):
    """Will generate the float values used to scale the size of the arrows when plotting the summary heatmap

    Args:
        data_frame (dataframe): Takes dataframe containing the transition frequencies generated using the
    'transition_frequency' function
        treatment (str): will take 'treatment' from for loop. Leave as treatment.

    Returns:
        [dataframe]: returns values used to scale arrow -
    """
    transition_frequency_arrow = data_frame[data_frame["sample"]== treatment]
    normalised_number_lowtohigh = float(np.array(transition_frequency_arrow[f"< {FRET_thresh} to > {FRET_thresh}"])/1000)
    normalised_number_hightolow = float(np.array(transition_frequency_arrow[f"> {FRET_thresh} to < {FRET_thresh}"])/1000)
    normalised_number_hightohigh = float(np.array(transition_frequency_arrow[f"> {FRET_thresh} to > {FRET_thresh}"])/1000)
    normalised_number_lowtolow = float(np.array(transition_frequency_arrow[f"< {FRET_thresh} to < {FRET_thresh}"])/1000)
    arrow_list = [normalised_number_lowtohigh,normalised_number_hightolow,normalised_number_hightohigh,normalised_number_lowtolow]
    return arrow_list

def heatmap_prep(histogram_data, treatment, FRET_thresh):
    """Takes the data and calculates the number of data points total (total), how much data points are below
    a threshold (time below) and above a threshold (time above) and will use these values to calculate what proportion
    of time the FRET is below or above thresh. Will then feed into the heatmap when plotting

    Args:
        histogram_data (dataframe): data used to plot the histogram - will include all the cleaned FRET and
    idealized FRET values (time not a factor here, just number of frames)
        treatment (str): only used to be fed from for loop. do not change

    Returns:
        df: contains the data required to plot the heatmap. will be as a proportion of time spent below or above
        threshold
    """
    subset_data = histogram_data[histogram_data["treatment_name"]==treatment]
    total = len(subset_data[(subset_data["FRET"] < FRET_thresh) | (subset_data["FRET"] > FRET_thresh)])
    subset_data_largerthanthresh = len(subset_data[subset_data["FRET"] > FRET_thresh])
    subset_data_lessthanthresh = len(subset_data[subset_data["FRET"] < FRET_thresh])
    time_below = (subset_data_lessthanthresh/total)
    time_above = (subset_data_largerthanthresh/total)
    thresh_dicts = {f"< {FRET_thresh}":[time_below], f"> {FRET_thresh}":[time_above] }
    thresh_dfs = pd.DataFrame(thresh_dicts)
    #thresh_dfs["treatment"] = treatment
    return thresh_dfs

def mean_dwell_prep(mean_dwell_data, treatment, FRET_thresh):
    """calculates the mean values for each transition class and converts to float values for plotting in
    summary heatmap

    Args:
        mean_dwell_data (dataframe): dataframe containing mean values
        treatment (str): used for a for loop. do not change.

    Returns:
        [list]: contains list of float values
    """
    subset_mean_dwell = mean_dwell_data[mean_dwell_data["sample"]== treatment]
    meandwell_lowtohigh = float(np.array(subset_mean_dwell[f"< {FRET_thresh} to > {FRET_thresh}"]))
    meandwell_hightolow = float(np.array(subset_mean_dwell[f"> {FRET_thresh} to < {FRET_thresh}"]))
    meandwell_hightohigh = float(np.array(subset_mean_dwell[f"> {FRET_thresh} to > {FRET_thresh}"]))
    meandwell_lowtolow = float(np.array(subset_mean_dwell[f"< {FRET_thresh} to < {FRET_thresh}"]))
    mean_list = [meandwell_lowtohigh,meandwell_hightolow,meandwell_hightohigh,meandwell_lowtolow]
    return mean_list

def file_reader(input_folder, data_type, frame_rate = False, column_names = False):
    """will import data

    Args:
        input_folder (directory): where data is stored
        data_type (str): what data will be used to plot, needs to be either 'hist', 'TDP', 'transition_frequency'
        or 'other'.

    Returns:
        dataframe: dataframe with data to be used in subseqeunt codes
    """
    if data_type == 'hist':
        filenames = glob.glob(input_folder + "/*.dat")
        dfs = []
        for filename in filenames:
            molecule_number = filename.split('\\')[1].split('_')[0]
            hist_data = pd.read_table(filename, sep="\s+", header=None)
            hist_data['molecule number'] = molecule_number
            dfs.append(hist_data) ### will error if forward slash (e.g. "/s+")
        test = pd.concat(dfs)
        test_dfs = pd.DataFrame(test)
        return test_dfs
    elif data_type == 'TDP':
        filename = input_folder
        A = pd.read_table(filename, header = None, sep="\s+")
        A.columns = ['Molecule', 'Idealized_FRET']
        return A
    elif data_type == 'transition_frequency':
        if not column_names:
            print('no column_names found, specify list to use')
            return
        filenames = glob.glob(input_folder + "/*.csv")
        dfs = []
        for filename in filenames:
            dfs.append(pd.read_csv(filename, header=None))
        test = pd.concat(dfs, ignore_index=True)
        test_dfs = pd.DataFrame(test)
        test_dfs.columns = column_names
        return test_dfs
    elif data_type == 'heatmap':
        filenames = glob.glob(input_folder + "/*.dat")
        dfs = []
        for filename in filenames:
            name = filename.split('\\')[1].split('_')[0]
            mol_data = pd.read_table(filename, sep="\s+", header=None)
            mol_data.columns = ['frame', 'donor', 'acceptor', 'FRET', 'idealized_FRET']
            mol_data['molecule'] = name
            mol_data['time'] = mol_data['frame']/frame_rate
            dfs.append(mol_data) ### will error if forward slash (e.g. "/s+")
        test = pd.concat(dfs)
        test_dfs = pd.DataFrame(test)
        return test_dfs
    elif data_type == 'other':
        dfs = pd.read_csv(input_folder)
        return dfs
    else:
        print('invalid data_type, please set data_type as "hist", "TDP","transition_frequency" or "other" if using for violin or heatmap plots')

def count_filt_mol(df, thresh, dataname, order):
    """Will count the number of molecules in which the idealized FRET will go below a defined threshold at some point before photobleaching

    Args:
        df (dataframe): Contains all the data required to plot TDP and identify transitions below threshold (i.e., FRET, idealized FRET, molecule)
        thresh (float): Threshold to set. If set to 0.5, function will count the number of molecules that go below 0.5 at some point
        dataname (dict): Dictionary containing keys for each treatment - used to find mol count for each treatment

    Returns:
        dataframe: Will return dataframe with raw mol count and also corrected mol count. Corrected mol count is calculated as the Raw mol count subtracted
        by the molcount of another treatment. The treatment to subtract is defined by 'order', which is the index of the treatment you want to subtract
    """
    filtered_data = filter_TDP(df, thresh)
    data_paths = dataname
    percent_mol_concat = {}
    for data_name, data_path in data_paths.items():
        total_mol = len(df[df['treatment_name']==data_name].Molecule.unique())
        filt_mol = len(filtered_data[filtered_data['treatment_name']==data_name].Molecule.unique())
        percent_mol = (filt_mol/total_mol)*100
        percent_mol_concat[data_name] = percent_mol
    percent_mol_concat = pd.DataFrame(percent_mol_concat, index = ['percent_mol']).T.reset_index().rename(columns={'index':'treatment'})
    percent_mol_concat['norm_percent_mol'] = percent_mol_concat['percent_mol'] - percent_mol_concat.iloc[order,1]
    return percent_mol_concat

def fret_state_trans(dfs, thresh, fps_clean, thresh_clean, state = 'after'):
    """Prepares a dataset in which
    Will plot a violin plot of all the FRET states immediately prior to a transition to another FRET state that is below a defined threshold

    Args:
        dfs (dataframe): Contains all the data required to plot TDP and identify transitions below threshold (i.e., FRET, idealized FRET, molecule)
        thresh (float): Threshold that defines the FRET state that you want to look at. For example, if you want to look at the FRET state immediately priort
        to a transition below 0.3 FRET then you will set 'thresh' as 0.3
        fps_clean ([type]): Required for cleanup_dwell function. Needs this to convert frames to seconds
        thresh_clean ([type]): Required for cleanup_dwell function. Specifies the minimum residence time. All residence times less than thresh_clean will be deleted

    Returns:
        dataframe: Dataframe containing all the transitions in which the 'FRET after transition' is below 'thresh'
    """
    cleaned_df = []
    for treatment_name, df in dfs.groupby("treatment_name"):
        initial_data = df[df["treatment_name"] == treatment_name]
        cleaned = cleanup_dwell(initial_data, fps_clean, thresh_clean)
        cleaned_df.append(cleaned)
    cleaned_concat = pd.concat(cleaned_df)
    filt = []
    for treatment_name, df in cleaned_concat.groupby("treatment_name"):
        if state == 'after':
            filt.append(df[df['FRET_before'] <= thresh])
        elif state == 'before':
            filt.append(df[df['FRET_after'] <= thresh])
    filtered_f_state = pd.concat(filt)
    return filtered_f_state