text:
# Copyright (c) 2021 PaddlePaddle 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.

import paddle
import os
import paddle.nn as nn
import numpy as np
from numpy import *
import time
import logging
import sys
import threading
from threading import Thread
from importlib import import_module
__dir__ = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.abspath(os.path.join(__dir__, '../..')))
from utils.utils_single import load_yaml, load_dy_model_class, get_abs_model
from utils.save_load import save_model, load_model
from paddle.io import DistributedBatchSampler, DataLoader
import argparse
from paddle.inference import Config
from paddle.inference import create_predictor
#import pynvml
#import psutil
#import GPUtil
import multiprocessing
from queue import Queue
import paddle.distributed.fleet as fleet

logging.basicConfig(
    format='%(asctime)s - %(levelname)s - %(message)s', level=logging.INFO)
logger = logging.getLogger(__name__)

q = Queue()

class Reader(fleet.MultiSlotDataGenerator):
    def init(self, args):
        self.args = args
        self.padding = 0
        self.slots = self.slot_reader()
        self.fea_dict = {}
        self.visit = {}
        for slot in self.slots:
            self.visit[slot] = False

    def slot_reader(self):
        slot_num = 300
        slots = []
        if slot_num > 0:
            for i in range(2, slot_num + 2):
                slots.append(str(i))
        else:
            with open("./slot", "r") as rf:
                for line in rf.readlines():
                    slots.append(line.strip())
        return slots

    def batch_process(self, lines):
        output = {}
        lod = {}
        for slot in self.slots:
            output[slot] = []
            lod[slot] = [[0]]

        for line in lines:
            line = line.strip().split(" ")
            for i in range(len(line)):
                if i == 0:
                    continue
                slot_feasign = line[i].split(":")
                if len(slot_feasign) < 2:
                    print(i)
                slot = slot_feasign[1]
                if slot not in self.slots:
                    continue
                feasign = int(slot_feasign[0])
                if feasign not in self.fea_dict:
                    self.fea_dict[feasign] = len(self.fea_dict)
                output[slot].append(self.fea_dict[feasign])
                self.visit[slot] = True
            for slot in self.visit:
                if not self.visit[slot]:
                    output[slot].extend([self.padding])
                else:
                    self.visit[slot] = False
            for slot in self.slots:
                lod[slot][0].append(len(output[slot]))

        return output, lod

    def generate_batchSample(self):
        with open(self.args.reader_file, "r") as rf:
            cnt = 0
            lines = []
            for line in rf.readlines():
                cnt += 1
                lines.append(line)
                if (cnt == args.batchsize):
                    output, lod = self.batch_process(lines)
                    cnt = 0
                    i = 0
                    while True:
                        if i >= args.iteration_num:
                            break
                        q.put([lod, output])
                        i += 1
        return q.qsize()

def init_predictor(args):
    if args.model_dir:
        config = Config(args.model_dir)
    else:
        config = Config(args.model_file, args.params_file)

    if args.use_gpu:
        config.enable_use_gpu(1000, 0)
        if args.enable_tensorRT:
            config.enable_tensorrt_engine(
                max_batch_size=args.batchsize,
                min_subgraph_size=1,
                precision_mode=paddle.inference.PrecisionType.Float32)
    else:
        config.disable_gpu()
        config.switch_ir_optim()
        config.delete_pass("repeated_fc_relu_fuse_pass")
        config.set_cpu_math_library_num_threads(args.cpu_threads)
        if args.enable_mkldnn:
            config.enable_mkldnn()
    predictor = create_predictor(config)
    return predictor, config

class Times(object):
    def __init__(self):
        self.time = 0.
        self.st = 0.
        self.et = 0.

    def start(self):
        self.st = time.time()

    def end(self, accumulative=True):
        self.et = time.time()
        if accumulative:
            self.time += self.et - self.st
        else:
            self.time = self.et - self.st

    def reset(self):
        self.time = 0.
        self.st = 0.
        self.et = 0.

    def value(self):
        return round(self.time, 4)

def infer_main(predictor_cost_time, i, args):
    tid = threading.currentThread()
    print('Thread id : %d' % tid.ident)
    predictor, pred_config = init_predictor(args)
    place = paddle.set_device('gpu' if args.use_gpu else 'cpu')
    args.place = place
    input_names = predictor.get_input_names()
    results = {}
    output_names = predictor.get_output_names()
    for output_name in output_names:
        results[output_name] = []

    global q
    fin = open(args.log_file, 'w')
    batchCnt = 0
    t = Times()
    while True:
        if q.empty():
            break
        blob = q.get()
        batchCnt += 1
        for name in input_names:
            input_tensor = predictor.get_input_handle(name)
            input_tensor.set_lod(blob[0][name])
            nd = np.array(blob[1][name]).reshape((len(blob[1][name]), 1))
            input_tensor.copy_from_cpu(nd)
            '''
            fin.write("lod: ")
            for i in blob[0][name]:
                fin.write(str(i) + " ")
            fin.write("data: ")
            for i in blob[1][name]:
                fin.write(str(i) + " ")
            fin.write("\n")
            '''
        if args.test_predictor:
            t.start()
        predictor.run()
        if args.test_predictor:
            t.end(accumulative=True)
        '''
        for name in output_names:
            output_tensor = predictor.get_output_handle(name)
            output_data = output_tensor.copy_to_cpu()
            results[name].append(output_data)
        '''
        #fin.write("thread: " + str(tid.ident) + ", ")
        #fin.write("result: " + "\n")
        #fin.write(str(output_data))
    logger.info("processed batch num: {}".format(batchCnt))
    logger.info("predictor time cost: {}".format(t.value()))
    if args.test_predictor:
        predictor_cost_time[i] = t.value()
    fin.close()
    return results

class WrapperThread(Thread):
    def __init__(self, func, cost_time, thread_id, args):
        super(WrapperThread, self).__init__()
        self.func = func
        self.cost_time = cost_time
        self.thread_id = thread_id
        self.args = args

    def run(self):
        self.result = self.func(self.cost_time, self.thread_id, self.args)

    def get_result(self):
        return self.result

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--test_predictor",
        type=str,
        default="False",
        help="test paddle predictor")
    parser.add_argument("--thread_num", type=int, default=2, help="thread num")
    parser.add_argument("--batchsize", type=int, default=5, help="batch size")
    parser.add_argument(
        "--iteration_num", type=int, default=10, help="iteration num")
    parser.add_argument(
        "--reader_file", default="../data/out_test.1", type=str)
    parser.add_argument(
        "--params_file", default="../model/rec_inference.pdiparams", type=str)
    parser.add_argument(
        "--model_file", default="../model/rec_inference.pdmodel", type=str)
    parser.add_argument("--log_file", default="./slot_dnn_infer.log", type=str)
    parser.add_argument(
        "--performance_file", default="./performance.txt", type=str)
    parser.add_argument("--model_dir", type=str)
    parser.add_argument("--use_gpu", type=str, default="False")
    parser.add_argument("--data_dir", type=str)
    parser.add_argument("--model_name", type=str, default="rec_model")
    parser.add_argument("--cpu_threads", type=int, default=1)
    parser.add_argument("--enable_mkldnn", type=str, default="False")
    parser.add_argument("--enable_tensorRT", type=str, default="False")
    args = parser.parse_args()
    args.use_gpu = (True if args.use_gpu.lower() == "true" else False)
    args.enable_mkldnn = (True
                          if args.enable_mkldnn.lower() == "true" else False)
    args.enable_tensorRT = (True if args.enable_tensorRT.lower() == "true" else
                            False)
    args.test_predictor = (True
                           if args.test_predictor.lower() == "true" else False)
    return args

if __name__ == '__main__':
    args = parse_args()

    reader = Reader()
    reader.init(args)
    batchCnt = reader.generate_batchSample()
    samplesCnt = batchCnt * args.batchsize
    thread_num = args.thread_num
    predictor_cost_time = list(range(thread_num))
    threads = []
    ff = open(args.performance_file, 'a')
    inference_time = Times()
    inference_time.start()
    for i in range(thread_num):
        t = WrapperThread(infer_main, predictor_cost_time, i, args=args)
        threads.append(t)
        t.start()
    for i in threads:
        t.join()
    inference_time.end(accumulative=True)
    cost_time = inference_time.value()
    ff.write("thread num: " + str(args.thread_num) + "\n")
    ff.write("batch size: " + str(args.batchsize) + "\n")
    ff.write("total sample num: " + str(samplesCnt) + "\n")
    qps = samplesCnt / cost_time
    ff.write("qps: " + str(qps) + "\n")
    latency = 1 / qps
    ff.write("latency: " + str(latency) + "\n")
    if args.test_predictor:
        predictor_cost_time_mean = mean(predictor_cost_time)
        predictor_qps = samplesCnt / predictor_cost_time_mean
        ff.write("paddle predictor time cost: " + str(predictor_cost_time_mean)
                 + "\n")
        ff.write("paddle predictor qps: " + str(predictor_qps) + "\n")
    ff.write("++++++++++++++++++++++++++++++++\n")

max_issues_repo_name: jiangcongxu/PaddleRec

max_forks_repo_name: jiangcongxu/PaddleRec

reason: import numpy,from numpy


text:
# coding=utf-8
# Copyright 2021 The Google Research Authors.
#
# 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.

"""Script for training and evaluating RPC agents."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import functools
import json
import os
import time

from absl import app
from absl import flags
from absl import logging
import gin
import numpy as np
import rpc_agent
import rpc_utils
from six.moves import range
import tensorflow as tf
import tensorflow_probability as tfp
from tf_agents.agents import data_converter
from tf_agents.drivers import dynamic_step_driver
from tf_agents.environments import suite_gym
from tf_agents.environments import tf_py_environment
from tf_agents.eval import metric_utils
from tf_agents.metrics import tf_metrics
from tf_agents.policies import greedy_policy
from tf_agents.policies import random_tf_policy
from tf_agents.replay_buffers import tf_uniform_replay_buffer
from tf_agents.utils import common

flags.DEFINE_string('root_dir', os.getenv('TEST_UNDECLARED_OUTPUTS_DIR'),
                    'Root directory for writing logs/summaries/checkpoints.')
flags.DEFINE_multi_string('gin_file', None, 'Path to the trainer config files.')
flags.DEFINE_multi_string('gin_bindings', None, 'Gin binding to pass through.')

FLAGS = flags.FLAGS

DEFAULT_KL_CONSTRAINT = 1.0
DEFAULT_DUAL_LR = 1.0

@gin.configurable
def train_eval(
    root_dir,
    env_name='HalfCheetah-v2',
    num_iterations=3000000,
    actor_fc_layers=(),
    critic_obs_fc_layers=None,
    critic_action_fc_layers=None,
    critic_joint_fc_layers=(256, 256),
    initial_collect_steps=10000,
    collect_steps_per_iteration=1,
    replay_buffer_capacity=1000000,
    # Params for target update
    target_update_tau=0.005,
    target_update_period=1,
    # Params for train
    train_steps_per_iteration=1,
    batch_size=256,
    actor_learning_rate=3e-4,
    critic_learning_rate=3e-4,
    alpha_learning_rate=3e-4,
    dual_learning_rate=3e-4,
    td_errors_loss_fn=tf.math.squared_difference,
    gamma=0.99,
    reward_scale_factor=0.1,
    gradient_clipping=None,
    use_tf_functions=True,
    # Params for eval
    num_eval_episodes=30,
    eval_interval=10000,
    # Params for summaries and logging
    train_checkpoint_interval=50000,
    policy_checkpoint_interval=50000,
    rb_checkpoint_interval=50000,
    log_interval=1000,
    summary_interval=1000,
    summaries_flush_secs=10,
    debug_summaries=False,
    summarize_grads_and_vars=False,
    eval_metrics_callback=None,
    latent_dim=10,
    log_prob_reward_scale=0.0,
    predictor_updates_encoder=False,
    predict_prior=True,
    use_recurrent_actor=False,
    rnn_sequence_length=20,
    clip_max_stddev=10.0,
    clip_min_stddev=0.1,
    clip_mean=30.0,
    predictor_num_layers=2,
    use_identity_encoder=False,
    identity_encoder_single_stddev=False,
    kl_constraint=1.0,
    eval_dropout=(),
    use_residual_predictor=True,
    gym_kwargs=None,
    predict_prior_std=True,
    random_seed=0,):
  """A simple train and eval for SAC."""
  np.random.seed(random_seed)
  tf.random.set_seed(random_seed)
  if use_recurrent_actor:
    batch_size = batch_size // rnn_sequence_length
  root_dir = os.path.expanduser(root_dir)
  train_dir = os.path.join(root_dir, 'train')
  eval_dir = os.path.join(root_dir, 'eval')

  train_summary_writer = tf.compat.v2.summary.create_file_writer(
      train_dir, flush_millis=summaries_flush_secs * 1000)
  train_summary_writer.set_as_default()

  eval_summary_writer = tf.compat.v2.summary.create_file_writer(
      eval_dir, flush_millis=summaries_flush_secs * 1000)

  global_step = tf.compat.v1.train.get_or_create_global_step()
  with tf.compat.v2.summary.record_if(
      lambda: tf.math.equal(global_step % summary_interval, 0)):

    _build_env = functools.partial(suite_gym.load, environment_name=env_name,  # pylint: disable=invalid-name
                                   gym_env_wrappers=(), gym_kwargs=gym_kwargs)

    tf_env = tf_py_environment.TFPyEnvironment(_build_env())
    eval_vec = []  # (name, env, metrics)
    eval_metrics = [
        tf_metrics.AverageReturnMetric(buffer_size=num_eval_episodes),
        tf_metrics.AverageEpisodeLengthMetric(buffer_size=num_eval_episodes)
    ]
    eval_tf_env = tf_py_environment.TFPyEnvironment(_build_env())
    name = ''
    eval_vec.append((name, eval_tf_env, eval_metrics))

    time_step_spec = tf_env.time_step_spec()
    observation_spec = time_step_spec.observation
    action_spec = tf_env.action_spec()
    if latent_dim == 'obs':
      latent_dim = observation_spec.shape[0]

    def _activation(t):
      t1, t2 = tf.split(t, 2, axis=1)
      low = -np.inf if clip_mean is None else -clip_mean
      high = np.inf if clip_mean is None else clip_mean
      t1 = rpc_utils.squash_to_range(t1, low, high)

      if clip_min_stddev is None:
        low = -np.inf
      else:
        low = tf.math.log(tf.exp(clip_min_stddev) - 1.0)
      if clip_max_stddev is None:
        high = np.inf
      else:
        high = tf.math.log(tf.exp(clip_max_stddev) - 1.0)
      t2 = rpc_utils.squash_to_range(t2, low, high)
      return tf.concat([t1, t2], axis=1)

    if use_identity_encoder:
      assert latent_dim == observation_spec.shape[0]
      obs_input = tf.keras.layers.Input(observation_spec.shape)
      zeros = 0.0 * obs_input[:, :1]
      stddev_dim = 1 if identity_encoder_single_stddev else latent_dim
      pre_stddev = tf.keras.layers.Dense(stddev_dim, activation=None)(zeros)
      ones = zeros + tf.ones((1, latent_dim))
      pre_stddev = pre_stddev * ones  # Multiply to broadcast to latent_dim.
      pre_mean_stddev = tf.concat([obs_input, pre_stddev], axis=1)
      output = tfp.layers.IndependentNormal(latent_dim)(pre_mean_stddev)
      encoder_net = tf.keras.Model(inputs=obs_input,
                                   outputs=output)
    else:
      encoder_net = tf.keras.Sequential([
          tf.keras.layers.Dense(256, activation='relu'),
          tf.keras.layers.Dense(256, activation='relu'),
          tf.keras.layers.Dense(
              tfp.layers.IndependentNormal.params_size(latent_dim),
              activation=_activation,
              kernel_initializer='glorot_uniform'),
          tfp.layers.IndependentNormal(latent_dim),
      ])

    # Build the predictor net
    obs_input = tf.keras.layers.Input(observation_spec.shape)
    action_input = tf.keras.layers.Input(action_spec.shape)

    class ConstantIndependentNormal(tfp.layers.IndependentNormal):
      """A keras layer that always returns N(0, 1) distribution."""

      def call(self, inputs):
        loc_scale = tf.concat([
            tf.zeros((latent_dim,)),
            tf.fill((latent_dim,), tf.math.log(tf.exp(1.0) - 1))
        ],
                              axis=0)
        # Multiple by [B x 1] tensor to broadcast batch dimension.
        loc_scale = loc_scale * tf.ones_like(inputs[:, :1])
        return super(ConstantIndependentNormal, self).call(loc_scale)

    if predict_prior:
      z = encoder_net(obs_input)
      if not predictor_updates_encoder:
        z = tf.stop_gradient(z)
      za = tf.concat([z, action_input], axis=1)
      if use_residual_predictor:
        za_input = tf.keras.layers.Input(za.shape[1])
        loc_scale = tf.keras.Sequential(
            predictor_num_layers * [tf.keras.layers.Dense(256, activation='relu')] + [  # pylint: disable=line-too-long
                tf.keras.layers.Dense(
                    tfp.layers.IndependentNormal.params_size(latent_dim),
                    activation=_activation,
                    kernel_initializer='zeros'),
            ])(za_input)
        if predict_prior_std:
          combined_loc_scale = tf.concat([
              loc_scale[:, :latent_dim] + za_input[:, :latent_dim],
              loc_scale[:, latent_dim:]
          ],
                                         axis=1)
        else:
          # Note that softplus(log(e - 1)) = 1.
          combined_loc_scale = tf.concat([
              loc_scale[:, :latent_dim] + za_input[:, :latent_dim],
              tf.math.log(np.e - 1) * tf.ones_like(loc_scale[:, latent_dim:])
          ],
                                         axis=1)
        dist = tfp.layers.IndependentNormal(latent_dim)(combined_loc_scale)
        output = tf.keras.Model(inputs=za_input, outputs=dist)(za)
      else:
        assert predict_prior_std
        output = tf.keras.Sequential(
            predictor_num_layers * [tf.keras.layers.Dense(256, activation='relu')] +  # pylint: disable=line-too-long
            [tf.keras.layers.Dense(
                tfp.layers.IndependentNormal.params_size(latent_dim),
                activation=_activation,
                kernel_initializer='zeros'),
             tfp.layers.IndependentNormal(latent_dim),
             ])(za)
    else:
      # scale is chosen by inverting the softplus function to equal 1.
      if len(obs_input.shape) > 2:
        input_reshaped = tf.reshape(
            obs_input, [-1, tf.math.reduce_prod(obs_input.shape[1:])])
        #  Multiply by [B x 1] tensor to broadcast batch dimension.
        za = tf.zeros(latent_dim + action_spec.shape[0],) * tf.ones_like(input_reshaped[:, :1])  # pylint: disable=line-too-long
      else:
        #  Multiple by [B x 1] tensor to broadcast batch dimension.
        za = tf.zeros(latent_dim + action_spec.shape[0],) * tf.ones_like(obs_input[:, :1])  # pylint: disable=line-too-long
      output = tf.keras.Sequential([
          ConstantIndependentNormal(latent_dim),
      ])(
          za)
    predictor_net = tf.keras.Model(inputs=(obs_input, action_input),
                                   outputs=output)
    if use_recurrent_actor:
      ActorClass = rpc_utils.RecurrentActorNet  # pylint: disable=invalid-name
    else:
      ActorClass = rpc_utils.ActorNet  # pylint: disable=invalid-name
    actor_net = ActorClass(
        input_tensor_spec=observation_spec,
        output_tensor_spec=action_spec,
        encoder=encoder_net,
        predictor=predictor_net,
        fc_layers=actor_fc_layers)

    critic_net = rpc_utils.CriticNet(
        (observation_spec, action_spec),
        observation_fc_layer_params=critic_obs_fc_layers,
        action_fc_layer_params=critic_action_fc_layers,
        joint_fc_layer_params=critic_joint_fc_layers,
        kernel_initializer='glorot_uniform',
        last_kernel_initializer='glorot_uniform')
    critic_net_2 = None
    target_critic_net_1 = None
    target_critic_net_2 = None

    tf_agent = rpc_agent.RpAgent(
        time_step_spec,
        action_spec,
        actor_network=actor_net,
        critic_network=critic_net,
        critic_network_2=critic_net_2,
        target_critic_network=target_critic_net_1,
        target_critic_network_2=target_critic_net_2,
        actor_optimizer=tf.compat.v1.train.AdamOptimizer(
            learning_rate=actor_learning_rate),
        critic_optimizer=tf.compat.v1.train.AdamOptimizer(
            learning_rate=critic_learning_rate),
        alpha_optimizer=tf.compat.v1.train.AdamOptimizer(
            learning_rate=alpha_learning_rate),
        target_update_tau=target_update_tau,
        target_update_period=target_update_period,
        td_errors_loss_fn=td_errors_loss_fn,
        gamma=gamma,
        reward_scale_factor=reward_scale_factor,
        gradient_clipping=gradient_clipping,
        debug_summaries=debug_summaries,
        summarize_grads_and_vars=summarize_grads_and_vars,
        train_step_counter=global_step)
    dual_optimizer = tf.compat.v1.train.AdamOptimizer(
        learning_rate=dual_learning_rate)
    tf_agent.initialize()

    # Make the replay buffer.
    replay_buffer = tf_uniform_replay_buffer.TFUniformReplayBuffer(
        data_spec=tf_agent.collect_data_spec,
        batch_size=tf_env.batch_size,
        max_length=replay_buffer_capacity)
    replay_observer = [replay_buffer.add_batch]

    train_metrics = [
        tf_metrics.NumberOfEpisodes(),
        tf_metrics.EnvironmentSteps(),
        tf_metrics.AverageReturnMetric(
            buffer_size=num_eval_episodes, batch_size=tf_env.batch_size),
        tf_metrics.AverageEpisodeLengthMetric(
            buffer_size=num_eval_episodes, batch_size=tf_env.batch_size),
    ]
    kl_metric = rpc_utils.AverageKLMetric(
        encoder=encoder_net,
        predictor=predictor_net,
        batch_size=tf_env.batch_size)
    eval_policy = greedy_policy.GreedyPolicy(tf_agent.policy)
    initial_collect_policy = random_tf_policy.RandomTFPolicy(
        tf_env.time_step_spec(), tf_env.action_spec())
    collect_policy = tf_agent.collect_policy

    checkpoint_items = {
        'ckpt_dir': train_dir,
        'agent': tf_agent,
        'global_step': global_step,
        'metrics': metric_utils.MetricsGroup(train_metrics, 'train_metrics'),
        'dual_optimizer': dual_optimizer,
    }
    train_checkpointer = common.Checkpointer(**checkpoint_items)

    policy_checkpointer = common.Checkpointer(
        ckpt_dir=os.path.join(train_dir, 'policy'),
        policy=eval_policy,
        global_step=global_step)
    rb_checkpointer = common.Checkpointer(
        ckpt_dir=os.path.join(train_dir, 'replay_buffer'),
        max_to_keep=1,
        replay_buffer=replay_buffer)

    train_checkpointer.initialize_or_restore()
    rb_checkpointer.initialize_or_restore()

    initial_collect_driver = dynamic_step_driver.DynamicStepDriver(
        tf_env,
        initial_collect_policy,
        observers=replay_observer + train_metrics,
        num_steps=initial_collect_steps,
        transition_observers=[kl_metric])

    collect_driver = dynamic_step_driver.DynamicStepDriver(
        tf_env,
        collect_policy,
        observers=replay_observer + train_metrics,
        num_steps=collect_steps_per_iteration,
        transition_observers=[kl_metric])

    if use_tf_functions:
      initial_collect_driver.run = common.function(initial_collect_driver.run)
      collect_driver.run = common.function(collect_driver.run)
      tf_agent.train = common.function(tf_agent.train)

    if replay_buffer.num_frames() == 0:
      # Collect initial replay data.
      logging.info(
          'Initializing replay buffer by collecting experience for %d steps '
          'with a random policy.', initial_collect_steps)
      initial_collect_driver.run()

    for name, eval_tf_env, eval_metrics in eval_vec:
      results = metric_utils.eager_compute(
          eval_metrics,
          eval_tf_env,
          eval_policy,
          num_episodes=num_eval_episodes,
          train_step=global_step,
          summary_writer=eval_summary_writer,
          summary_prefix='Metrics-%s' % name,
      )
      if eval_metrics_callback is not None:
        eval_metrics_callback(results, global_step.numpy())
      metric_utils.log_metrics(eval_metrics, prefix=name)

    time_step = None
    policy_state = collect_policy.get_initial_state(tf_env.batch_size)

    timed_at_step = global_step.numpy()
    time_acc = 0
    train_time_acc = 0
    env_time_acc = 0

    if use_recurrent_actor:  # default from sac/train_eval_rnn.py
      num_steps = rnn_sequence_length + 1
      def _filter_invalid_transition(trajectories, unused_arg1):
        return tf.reduce_all(~trajectories.is_boundary()[:-1])

      tf_agent._as_transition = data_converter.AsTransition(  # pylint: disable=protected-access
          tf_agent.data_context, squeeze_time_dim=False)
    else:
      num_steps = 2

      def _filter_invalid_transition(trajectories, unused_arg1):
        return ~trajectories.is_boundary()[0]
    dataset = replay_buffer.as_dataset(
        sample_batch_size=batch_size,
        num_steps=num_steps).unbatch().filter(_filter_invalid_transition)

    dataset = dataset.batch(batch_size).prefetch(5)
    # Dataset generates trajectories with shape [Bx2x...]
    iterator = iter(dataset)

    @tf.function
    def train_step():
      experience, _ = next(iterator)

      prior = predictor_net((experience.observation[:, 0],
                             experience.action[:, 0]), training=False)
      z_next = encoder_net(experience.observation[:, 1], training=False)
      # predictor_kl is a vector of size batch_size.
      predictor_kl = tfp.distributions.kl_divergence(z_next, prior)

      with tf.GradientTape() as tape:
        tape.watch(actor_net._log_kl_coefficient)  # pylint: disable=protected-access
        dual_loss = -1.0 * actor_net._log_kl_coefficient * (  # pylint: disable=protected-access
            tf.stop_gradient(tf.reduce_mean(predictor_kl)) - kl_constraint)
      dual_grads = tape.gradient(dual_loss, [actor_net._log_kl_coefficient])  # pylint: disable=protected-access
      grads_and_vars = list(zip(dual_grads, [actor_net._log_kl_coefficient]))  # pylint: disable=protected-access
      dual_optimizer.apply_gradients(grads_and_vars)

      # Clip the dual variable so exp(log_kl_coef) <= 1e6.
      log_kl_coef = tf.clip_by_value(actor_net._log_kl_coefficient,  # pylint: disable=protected-access
                                     -1.0 * np.log(1e6), np.log(1e6))
      actor_net._log_kl_coefficient.assign(log_kl_coef)  # pylint: disable=protected-access

      with tf.name_scope('dual_loss'):
        tf.compat.v2.summary.scalar(name='dual_loss',
                                    data=tf.reduce_mean(dual_loss),
                                    step=global_step)
        tf.compat.v2.summary.scalar(name='log_kl_coefficient',
                                    data=actor_net._log_kl_coefficient,  # pylint: disable=protected-access
                                    step=global_step)

      z_entropy = z_next.entropy()
      log_prob = prior.log_prob(z_next.sample())
      with tf.name_scope('rp-metrics'):
        common.generate_tensor_summaries('predictor_kl', predictor_kl,
                                         global_step)
        common.generate_tensor_summaries('z_entropy', z_entropy, global_step)
        common.generate_tensor_summaries('log_prob', log_prob, global_step)
        common.generate_tensor_summaries('z_mean', z_next.mean(), global_step)
        common.generate_tensor_summaries('z_stddev', z_next.stddev(),
                                         global_step)
        common.generate_tensor_summaries('prior_mean', prior.mean(),
                                         global_step)
        common.generate_tensor_summaries('prior_stddev', prior.stddev(),
                                         global_step)

      if log_prob_reward_scale == 'auto':
        coef = tf.stop_gradient(tf.exp(actor_net._log_kl_coefficient))  # pylint: disable=protected-access
      else:
        coef = log_prob_reward_scale
      tf.debugging.check_numerics(
          tf.reduce_mean(predictor_kl), 'predictor_kl is inf or nan.')
      tf.debugging.check_numerics(coef, 'coef is inf or nan.')
      new_reward = experience.reward - coef * predictor_kl[:, None]

      experience = experience._replace(reward=new_reward)
      return tf_agent.train(experience)

    if use_tf_functions:
      train_step = common.function(train_step)

    # Save the hyperparameters
    operative_filename = os.path.join(root_dir, 'operative.gin')
    with tf.compat.v1.gfile.Open(operative_filename, 'w') as f:
      f.write(gin.operative_config_str())
      print(gin.operative_config_str())

    global_step_val = global_step.numpy()
    while global_step_val < num_iterations:
      start_time = time.time()
      time_step, policy_state = collect_driver.run(
          time_step=time_step,
          policy_state=policy_state,
      )
      env_time_acc += time.time() - start_time
      train_start_time = time.time()
      for _ in range(train_steps_per_iteration):
        train_loss = train_step()
      train_time_acc += time.time() - train_start_time
      time_acc += time.time() - start_time

      global_step_val = global_step.numpy()

      if global_step_val % log_interval == 0:
        logging.info('step = %d, loss = %f', global_step_val,
                     train_loss.loss)
        steps_per_sec = (global_step_val - timed_at_step) / time_acc
        logging.info('%.3f steps/sec', steps_per_sec)
        tf.compat.v2.summary.scalar(
            name='global_steps_per_sec', data=steps_per_sec, step=global_step)
        train_steps_per_sec = (global_step_val - timed_at_step) / train_time_acc
        logging.info('Train: %.3f steps/sec', train_steps_per_sec)
        tf.compat.v2.summary.scalar(
            name='train_steps_per_sec',
            data=train_steps_per_sec,
            step=global_step)
        env_steps_per_sec = (global_step_val - timed_at_step) / env_time_acc
        logging.info('Env: %.3f steps/sec', env_steps_per_sec)
        tf.compat.v2.summary.scalar(
            name='env_steps_per_sec', data=env_steps_per_sec, step=global_step)
        timed_at_step = global_step_val
        time_acc = 0
        train_time_acc = 0
        env_time_acc = 0

      for train_metric in train_metrics + [kl_metric]:
        train_metric.tf_summaries(
            train_step=global_step, step_metrics=train_metrics[:2])

      if global_step_val % eval_interval == 0:
        start_time = time.time()
        for name, eval_tf_env, eval_metrics in eval_vec:
          results = metric_utils.eager_compute(
              eval_metrics,
              eval_tf_env,
              eval_policy,
              num_episodes=num_eval_episodes,
              train_step=global_step,
              summary_writer=eval_summary_writer,
              summary_prefix='Metrics-%s' % name,
          )
          if eval_metrics_callback is not None:
            eval_metrics_callback(results, global_step_val)
          metric_utils.log_metrics(eval_metrics, prefix=name)
        logging.info('Evaluation: %d min', (time.time() - start_time) / 60)
        for prob_dropout in eval_dropout:
          rpc_utils.eval_dropout_fn(
              eval_tf_env, actor_net, global_step, prob_dropout=prob_dropout)

      if global_step_val % train_checkpoint_interval == 0:
        train_checkpointer.save(global_step=global_step_val)

      if global_step_val % policy_checkpoint_interval == 0:
        policy_checkpointer.save(global_step=global_step_val)

      if global_step_val % rb_checkpoint_interval == 0:
        rb_checkpointer.save(global_step=global_step_val)

def main(_):
  os.environ['TF_GPU_THREAD_MODE'] = 'gpu_private'

  tf.compat.v1.enable_v2_behavior()
  logging.set_verbosity(logging.INFO)

  gin.parse_config_files_and_bindings(FLAGS.gin_file, FLAGS.gin_bindings)

  if 'xm_parameters' in FLAGS and FLAGS.xm_parameters:
    hparams = json.loads(FLAGS.xm_parameters)
    with gin.unlock_config():
      for (key, value) in hparams.items():
        print('Setting: %s = %s' % (key, value))
        gin.bind_parameter(key, value)

  root_dir = FLAGS.root_dir
  train_eval(root_dir)

if __name__ == '__main__':
  flags.mark_flag_as_required('root_dir')
  app.run(main)

ext: py

max_stars_repo_name: xxdreck/google-research

max_issues_repo_name: xxdreck/google-research

max_forks_repo_name: admariner/google-research

reason: import numpy


text:
from abc import ABC

import numpy as np
import scipy.sparse as sps

from src.model.HybridRecommender.AbstractHybridRecommender import AbstractHybridRecommender

# --------- UTILITY ---------
def get_ui_matrix_from_ranking(ranking, scores, n_users, n_items):
    ranking_ui_matrix = sps.coo_matrix((n_users, n_items))
    item_per_user = ranking.shape[1]
    ranking_ui_matrix.row = np.repeat(np.arange(0, n_users), item_per_user)
    ranking_ui_matrix.col = np.reshape(ranking, newshape=ranking_ui_matrix.row.shape)
    ranking_ui_matrix.data = np.reshape(scores, newshape=ranking_ui_matrix.row.shape)
    return ranking_ui_matrix.tocsr()

class RerankingStrategyInterface(ABC):

    def get_sub_scores_related_to_main_ranking(self, main_ranking, sub_ranking, sub_scores, weight):
        """
        Get new rankings based on the recommendation strategy implemented

        :param main_ranking: main recommendations of the main recommender
        :param sub_ranking: sub recommendations of the sub recommender
        :param sub_scores: sub scores of the sub recommender
        :param weight: weight of the sub recommender
        :return: the new sub scores to be added on the main scores
        """
        raise NotImplementedError("This method has to be implemented by another class")

class VotingStrategy(RerankingStrategyInterface):

    def get_sub_scores_related_to_main_ranking(self, main_ranking: np.ndarray, sub_ranking, sub_scores, weight):
        # Set sub scores all to 1 since it is a voting and each recommender votes the items in the same way
        sub_scores[:, :] = 1

        main_max_item_id = np.max(np.max(main_ranking, axis=-1))
        sub_max_item_id = np.max(np.max(sub_ranking, axis=-1))
        max_item_id = np.maximum(main_max_item_id, sub_max_item_id)
        n_items = max_item_id + 1

        main_ranking_ui_matrix = get_ui_matrix_from_ranking(main_ranking,
                                                            np.ones(shape=main_ranking.shape),
                                                            main_ranking.shape[0], n_items)
        sub_ranking_ui_matrix = get_ui_matrix_from_ranking(sub_ranking, sub_scores * weight,
                                                           main_ranking.shape[0],
                                                           n_items)
        new_sub_scores_ui_matrix = main_ranking_ui_matrix.multiply(sub_ranking_ui_matrix)

        # Reshape the new sub scores matrix into a main_ranking shape
        new_sub_scores = new_sub_scores_ui_matrix[np.arange(main_ranking.shape[0]), main_ranking.T].T.todense()
        return np.array(new_sub_scores)

class WeightedAverageStrategy(RerankingStrategyInterface):

    def get_sub_scores_related_to_main_ranking(self, main_ranking, sub_ranking, sub_scores, weight):
        main_max_item_id = np.max(np.max(main_ranking, axis=-1))
        sub_max_item_id = np.max(np.max(sub_ranking, axis=-1))
        max_item_id = np.maximum(main_max_item_id, sub_max_item_id)
        n_items = max_item_id + 1

        main_ranking_ui_matrix = get_ui_matrix_from_ranking(main_ranking,
                                                            np.ones(shape=main_ranking.shape),
                                                            main_ranking.shape[0], n_items)
        sub_ranking_ui_matrix = get_ui_matrix_from_ranking(sub_ranking, sub_scores * weight,
                                                           main_ranking.shape[0],
                                                           n_items)
        new_sub_scores_ui_matrix = main_ranking_ui_matrix.multiply(sub_ranking_ui_matrix)

        # Reshape the new sub scores matrix into a main_ranking shape
        new_sub_scores = new_sub_scores_ui_matrix[np.arange(main_ranking.shape[0]), main_ranking.T].T.todense()
        return np.array(new_sub_scores)

class HybridRerankingRecommender(AbstractHybridRecommender):
    """
    Hybrid reranking recommender
    """

    RECOMMENDER_NAME = "HybridRerankingRecommender"

    def __init__(self, URM_train, main_recommender):
        super().__init__(URM_train)
        self.hybrid_strategy = None
        self.main_recommender = main_recommender
        self.weights = None
        self.STRATEGY_MAPPER = {"weighted_avg": WeightedAverageStrategy(),
                                "norm_weighted_avg": WeightedAverageStrategy(),
                                "voting": VotingStrategy()}

    def fit(self, main_cutoff=20, sub_cutoff=5, bias=True, strategy="voting", **weights):
        """
        Fit the hybrid model by setting the weight of each recommender

        :param main_cutoff: cutoff of the main recommendations that has to be reranked
        :param bias: True if the main_recommender contributes to the reranking, False otherwise
        :param strategy: an object of strategy pattern that handles the hybrid core functioning of the recommender
        system
        :param sub_cutoff: the multiplier used for multiplying the cutoff for handling more recommended items
        to average
        :return: None
        """
        if strategy not in self.STRATEGY_MAPPER:
            raise ValueError("The strategy name passed does not correspond to the implemented strategy: "
                             "{}".format(self.STRATEGY_MAPPER))
        self.weights = weights
        self.main_cutoff = main_cutoff
        self.sub_cutoff = sub_cutoff
        self.bias = bias
        self.hybrid_strategy = self.STRATEGY_MAPPER[strategy]
        self.normalize = False
        if strategy == "norm_weighted_avg":
            self.normalize = True

    def recommend(self, user_id_array, cutoff=None, remove_seen_flag=True, items_to_compute=None,
                  remove_top_pop_flag=False, remove_custom_items_flag=False, return_scores=False):
        """
        Recommend "number of cutoff" items to a given user_id_array

        REMARK: The RMSE measure got from this hybrid model is incorrect because the score returned is the one from
                last model in the for
        :param user_id_array:
        :param cutoff:
        :param remove_seen_flag:
        :param items_to_compute:
        :param remove_top_pop_flag:
        :param remove_custom_items_flag:
        :param return_scores:
        :return: for each user in user_id_array, a list of recommended items
        """
        # handle user_id_array in case it is scalar
        if np.isscalar(user_id_array):
            user_id_array = np.atleast_1d(user_id_array)
            single_user = True
        else:
            single_user = False

        main_ranking, main_scores_all = self.main_recommender.recommend(user_id_array, cutoff=self.main_cutoff,
                                                                        remove_seen_flag=remove_seen_flag,
                                                                        items_to_compute=items_to_compute,
                                                                        remove_top_pop_flag=remove_top_pop_flag,
                                                                        remove_custom_items_flag=remove_custom_items_flag,
                                                                        return_scores=True)
        main_ranking = np.array(main_ranking, dtype=np.int)
        main_scores = np.zeros(shape=main_ranking.shape)
        if self.bias:
            main_scores = main_scores_all[np.arange(len(main_ranking)), main_ranking.T].T

        for idx, recommender_key_value in enumerate(self.models.items()):
            recommender_name, recommender_model = recommender_key_value
            ranking, scores_all = recommender_model.recommend(user_id_array, cutoff=self.sub_cutoff,
                                                              remove_seen_flag=remove_seen_flag,
                                                              items_to_compute=items_to_compute,
                                                              remove_top_pop_flag=remove_top_pop_flag,
                                                              remove_custom_items_flag=remove_custom_items_flag,
                                                              return_scores=True)
            if self.normalize:
                maximum = np.max(scores_all, axis=1).reshape((len(user_id_array), 1))

                scores_batch_for_minimum = scores_all.copy()
                scores_batch_for_minimum[scores_batch_for_minimum == float('-inf')] = np.inf
                minimum = np.min(scores_batch_for_minimum, axis=1).reshape((len(user_id_array), 1))

                denominator = maximum - minimum
                denominator[denominator == 0] = 1.0

                scores_all = (scores_all - minimum) / denominator

            ranking = np.array(ranking, dtype=np.int)
            ranking_scores = scores_all[np.arange(len(ranking)), ranking.T].T
            sub_scores = self.hybrid_strategy.get_sub_scores_related_to_main_ranking(main_ranking=main_ranking,
                                                                                sub_ranking=ranking,
                                                                                sub_scores=ranking_scores,
                                                                                weight=self.weights[
                                                                                recommender_name])
            main_scores += sub_scores
        main_scores_argsort_indices = np.argsort(-main_scores, axis=1)
        rankings = main_ranking[np.arange(len(main_ranking)), main_scores_argsort_indices.T].T
        rankings = rankings[:, :cutoff]
        if single_user:
            rankings = rankings[0]

        if return_scores:
            return rankings, np.full(shape=(len(user_id_array), self.n_items), fill_value=1)
        else:
            return rankings

    def _compute_item_score(self, user_id_array, items_to_compute=None):
        # Useless for this recommender
        pass

    def copy(self):
        copy = HybridRerankingRecommender(URM_train=self.URM_train, main_recommender=self.main_recommender)
        copy.models = self.models
        return copy

    def save_model(self, folder_path, file_name=None):
        pass

    @classmethod
    def get_possible_strategies(cls):
        return ["voting", "weighted_avg", "norm_weighted_avg"]

reason: import numpy,import scipy


text:
#import pandas as pd
import numpy as np
#import matplotlib.pyplot as plt
import os
from random import randint
from skimage import io, transform
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader
from torch.utils.data.sampler import SubsetRandomSampler

def save_checkpoint(state, filename='/home/jamal/Downloads/GTSRB_Final_Training_Images/GTSRB/Final_Training/checkpoint.pth.tar'):
    torch.save(state, filename)

file_path ='/home/jamal/Downloads/GTSRB_Final_Training_Images/GTSRB/Final_Training/path_labels_all_train.txt'

load_net = 0
lr = 0.01

import torch
import torch.nn.functional as F
import torch.nn as nn

class Net(nn.Module):

    def __init__(self):
        super(Net,self).__init__()
        self.conv11 = nn.Conv2d(3,64,3)
        self.n11 = nn.BatchNorm2d(64)
        self.conv21 = nn.Conv2d(64,96,3)
        self.n21 = nn.BatchNorm2d(96)
        self.conv31 = nn.Conv2d(96,128,4)

        self.ip1 = nn.Linear(128,128)
        self.ip2 = nn.Linear(128,43)

    def forward(self,x):
        x = self.conv11(x)
        x = self.n11(x)
        x = F.relu(x)

        x = F.max_pool2d(x,3,stride = 3)

        x = self.conv21(x)

        x = self.n21(x)
        x = F.relu(x)
        x = F.max_pool2d(x,3,stride = 3)
        x = self.conv31(x)
        #
        # x = self.n31(x)
        # x = F.relu(x)
        # x = F.max_pool2d(x,3)
        x = x.view(-1,128)
        x = self.ip1(x)
        x = F.relu(x)
        x = self.ip2(x)
        return x

class trafic_sign_dataset(Dataset):

    def __init__(self,images_labels_txt,shape,transform=None):
        self.txt_file = open(images_labels_txt,'r')
        # self.root_dir = root_dir
        self.transform = transform
        self.shape = shape
        lst = self.txt_file.readlines()
        self.txt_file.close()
        paths = [item.rstrip().split('\t')[0] for item in lst]
        labels = [int(item.rstrip().split('\t')[1])  for item in lst]
        self.paths = paths
        self.labels = labels

    def __len__(self):

        return len(self.labels)

    def __getitem__(self,idx):
        img1_name = self.paths[idx]
        image1 = io.imread(img1_name)
        image1 = transform.resize(image1,shape)
        image1 = torch.Tensor.float(torch.from_numpy(image1))
        image1 = (torch.Tensor.permute(image1,(2,0,1)))
        image1 = image1/256.0

        annot = self.labels[idx]
        # annot = annot.astype('float').reshape(-1,2)
        sample = {'image1':image1,'lab' : annot}
        return sample

use_cuda = torch.cuda.is_available()
device = torch.device( "cpu")
shape = (48,48)

dataset = trafic_sign_dataset(file_path,shape)
if load_net:
    net = Net().to(device)
    checkpoint = torch.load('/home/jamal/Downloads/GTSRB_Final_Training_Images/GTSRB/Final_Training\checkpoint.pth.tar')
    net.load_state_dict ( checkpoint['state_dict'])
    optimizer = optim.Adam(net.parameters(),lr = lr, weight_decay = 0.0005)
    optimizer.load_state_dict = checkpoint['optimizer']
else:

    net = Net().to(device)
    optimizer = optim.Adam(net.parameters(),lr = lr, weight_decay = 0.0005)
    checkpoint= {'epoch':0}

optimizer = optim.Adam(net.parameters(),lr = 0.01, weight_decay = 0.0005)

train_dataloader = DataLoader(dataset,batch_size=32,shuffle = True)

err = []
acc = []

criterion = nn.CrossEntropyLoss()
for ep in range(checkpoint['epoch'],10):
    for i,data in enumerate(train_dataloader):
        optimizer.zero_grad()
        input1 , label = data.items()
        input1 = input1[1].to(device)
        label1 = label[1].to(device)
        label1 = torch.Tensor.long(label1)
        output = net(input1)
        loss = criterion (output,label1)
        loss.backward()
        optimizer.step()
        err.append(loss.item())
        if (i%5==0):
            print (loss.item(),i,ep)

    save_checkpoint({
            'epoch': ep + 1,
            'state_dict': net.state_dict(),

            'optimizer' : optimizer.state_dict(),
        })


text:
from functions import defun, defun_wrapped

@defun
def j0(ctx, x):
    """Computes the Bessel function `J_0(x)`. See :func:`~mpmath.besselj`."""
    return ctx.besselj(0, x)

@defun
def j1(ctx, x):
    """Computes the Bessel function `J_1(x)`.  See :func:`~mpmath.besselj`."""
    return ctx.besselj(1, x)

@defun
def besselj(ctx, n, z, derivative=0, **kwargs):
    if type(n) is int:
        n_isint = True
    else:
        n = ctx.convert(n)
        n_isint = ctx.isint(n)
        if n_isint:
            n = int(ctx._re(n))
    if n_isint and n < 0:
        return (-1)**n * ctx.besselj(-n, z, derivative, **kwargs)
    z = ctx.convert(z)
    M = ctx.mag(z)
    if derivative:
        d = ctx.convert(derivative)
        # TODO: the integer special-casing shouldn't be necessary.
        # However, the hypergeometric series gets inaccurate for large d
        # because of inaccurate pole cancellation at a pole far from
        # zero (needs to be fixed in hypercomb or hypsum)
        if ctx.isint(d) and d >= 0:
            d = int(d)
            orig = ctx.prec
            try:
                ctx.prec += 15
                v = ctx.fsum((-1)**k * ctx.binomial(d,k) * ctx.besselj(2*k+n-d,z)
                    for k in range(d+1))
            finally:
                ctx.prec = orig
            v *= ctx.mpf(2)**(-d)
        else:
            def h(n,d):
                r = ctx.fmul(ctx.fmul(z, z, prec=ctx.prec+M), -0.25, exact=True)
                B = [0.5*(n-d+1), 0.5*(n-d+2)]
                T = [([2,ctx.pi,z],[d-2*n,0.5,n-d],[],B,[(n+1)*0.5,(n+2)*0.5],B+[n+1],r)]
                return T
            v = ctx.hypercomb(h, [n,d], **kwargs)
    else:
    # Fast case: J_n(x), n int, appropriate magnitude for fixed-point calculation
        if (not derivative) and n_isint and abs(M) < 10 and abs(n) < 20:
            try:
                return ctx._besselj(n, z)
            except NotImplementedError:
                pass
        if not z:
            if not n:
                v = ctx.one + n+z
            elif ctx.re(n) > 0:
                v = n*z
            else:
                v = ctx.inf + z + n
        else:
            #v = 0
            orig = ctx.prec
            try:
                # XXX: workaround for accuracy in low level hypergeometric series
                # when alternating, large arguments
                ctx.prec += min(3*abs(M), ctx.prec)
                w = ctx.fmul(z, 0.5, exact=True)
                def h(n):
                    r = ctx.fneg(ctx.fmul(w, w, prec=max(0,ctx.prec+M)), exact=True)
                    return [([w], [n], [], [n+1], [], [n+1], r)]
                v = ctx.hypercomb(h, [n], **kwargs)
            finally:
                ctx.prec = orig
        v = +v
    return v

@defun
def besseli(ctx, n, z, derivative=0, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    if not z:
        if derivative:
            raise ValueError
        if not n:
            # I(0,0) = 1
            return 1+n+z
        if ctx.isint(n):
            return 0*(n+z)
        r = ctx.re(n)
        if r == 0:
            return ctx.nan*(n+z)
        elif r > 0:
            return 0*(n+z)
        else:
            return ctx.inf+(n+z)
    M = ctx.mag(z)
    if derivative:
        d = ctx.convert(derivative)
        def h(n,d):
            r = ctx.fmul(ctx.fmul(z, z, prec=ctx.prec+M), 0.25, exact=True)
            B = [0.5*(n-d+1), 0.5*(n-d+2), n+1]
            T = [([2,ctx.pi,z],[d-2*n,0.5,n-d],[n+1],B,[(n+1)*0.5,(n+2)*0.5],B,r)]
            return T
        v = ctx.hypercomb(h, [n,d], **kwargs)
    else:
        def h(n):
            w = ctx.fmul(z, 0.5, exact=True)
            r = ctx.fmul(w, w, prec=max(0,ctx.prec+M))
            return [([w], [n], [], [n+1], [], [n+1], r)]
        v = ctx.hypercomb(h, [n], **kwargs)
    return v

@defun_wrapped
def bessely(ctx, n, z, derivative=0, **kwargs):
    if not z:
        if derivative:
            # Not implemented
            raise ValueError
        if not n:
            # ~ log(z/2)
            return -ctx.inf + (n+z)
        if ctx.im(n):
            return nan * (n+z)
        r = ctx.re(n)
        q = n+0.5
        if ctx.isint(q):
            if n > 0:
                return -ctx.inf + (n+z)
            else:
                return 0 * (n+z)
        if r < 0 and int(ctx.floor(q)) % 2:
            return ctx.inf + (n+z)
        else:
            return ctx.ninf + (n+z)
    # XXX: use hypercomb
    ctx.prec += 10
    m, d = ctx.nint_distance(n)
    if d < -ctx.prec:
        h = +ctx.eps
        ctx.prec *= 2
        n += h
    elif d < 0:
        ctx.prec -= d
    # TODO: avoid cancellation for imaginary arguments
    cos, sin = ctx.cospi_sinpi(n)
    return (ctx.besselj(n,z,derivative,**kwargs)*cos - \
        ctx.besselj(-n,z,derivative,**kwargs))/sin

@defun_wrapped
def besselk(ctx, n, z, **kwargs):
    if not z:
        return ctx.inf
    M = ctx.mag(z)
    if M < 1:
        # Represent as limit definition
        def h(n):
            r = (z/2)**2
            T1 = [z, 2], [-n, n-1], [n], [], [], [1-n], r
            T2 = [z, 2], [n, -n-1], [-n], [], [], [1+n], r
            return T1, T2
    # We could use the limit definition always, but it leads
    # to very bad cancellation (of exponentially large terms)
    # for large real z
    # Instead represent in terms of 2F0
    else:
        ctx.prec += M
        def h(n):
            return [([ctx.pi/2, z, ctx.exp(-z)], [0.5,-0.5,1], [], [], \
                [n+0.5, 0.5-n], [], -1/(2*z))]
    return ctx.hypercomb(h, [n], **kwargs)

@defun_wrapped
def hankel1(ctx,n,x,**kwargs):
    return ctx.besselj(n,x,**kwargs) + ctx.j*ctx.bessely(n,x,**kwargs)

@defun_wrapped
def hankel2(ctx,n,x,**kwargs):
    return ctx.besselj(n,x,**kwargs) - ctx.j*ctx.bessely(n,x,**kwargs)

@defun_wrapped
def whitm(ctx,k,m,z,**kwargs):
    if z == 0:
        # M(k,m,z) = 0^(1/2+m)
        if ctx.re(m) > -0.5:
            return z
        elif ctx.re(m) < -0.5:
            return ctx.inf + z
        else:
            return ctx.nan * z
    x = ctx.fmul(-0.5, z, exact=True)
    y = 0.5+m
    return ctx.exp(x) * z**y * ctx.hyp1f1(y-k, 1+2*m, z, **kwargs)

@defun_wrapped
def whitw(ctx,k,m,z,**kwargs):
    if z == 0:
        g = abs(ctx.re(m))
        if g < 0.5:
            return z
        elif g > 0.5:
            return ctx.inf + z
        else:
            return ctx.nan * z
    x = ctx.fmul(-0.5, z, exact=True)
    y = 0.5+m
    return ctx.exp(x) * z**y * ctx.hyperu(y-k, 1+2*m, z, **kwargs)

@defun
def hyperu(ctx, a, b, z, **kwargs):
    a, atype = ctx._convert_param(a)
    b, btype = ctx._convert_param(b)
    z = ctx.convert(z)
    if not z:
        if ctx.re(b) <= 1:
            return ctx.gammaprod([1-b],[a-b+1])
        else:
            return ctx.inf + z
    bb = 1+a-b
    bb, bbtype = ctx._convert_param(bb)
    try:
        orig = ctx.prec
        try:
            ctx.prec += 10
            v = ctx.hypsum(2, 0, (atype, bbtype), [a, bb], -1/z, maxterms=ctx.prec)
            return v / z**a
        finally:
            ctx.prec = orig
    except ctx.NoConvergence:
        pass
    def h(a,b):
        w = ctx.sinpi(b)
        T1 = ([ctx.pi,w],[1,-1],[],[a-b+1,b],[a],[b],z)
        T2 = ([-ctx.pi,w,z],[1,-1,1-b],[],[a,2-b],[a-b+1],[2-b],z)
        return T1, T2
    return ctx.hypercomb(h, [a,b], **kwargs)

@defun
def struveh(ctx,n,z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/StruveH/26/01/02/
    def h(n):
        return [([z/2, 0.5*ctx.sqrt(ctx.pi)], [n+1, -1], [], [n+1.5], [1], [1.5, n+1.5], -(z/2)**2)]
    return ctx.hypercomb(h, [n], **kwargs)

@defun
def struvel(ctx,n,z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/StruveL/26/01/02/
    def h(n):
        return [([z/2, 0.5*ctx.sqrt(ctx.pi)], [n+1, -1], [], [n+1.5], [1], [1.5, n+1.5], (z/2)**2)]
    return ctx.hypercomb(h, [n], **kwargs)

def _anger(ctx,which,v,z,**kwargs):
    v = ctx._convert_param(v)[0]
    z = ctx.convert(z)
    def h(v):
        b = ctx.mpq_1_2
        u = v*b
        m = b*3
        a1,a2,b1,b2 = m-u, m+u, 1-u, 1+u
        c, s = ctx.cospi_sinpi(u)
        if which == 0:
            A, B = [b*z, s], [c]
        if which == 1:
            A, B = [b*z, -c], [s]
        w = ctx.square_exp_arg(z, mult=-0.25)
        T1 = A, [1, 1], [], [a1,a2], [1], [a1,a2], w
        T2 = B, [1], [], [b1,b2], [1], [b1,b2], w
        return T1, T2
    return ctx.hypercomb(h, [v], **kwargs)

@defun
def angerj(ctx, v, z, **kwargs):
    return _anger(ctx, 0, v, z, **kwargs)

@defun
def webere(ctx, v, z, **kwargs):
    return _anger(ctx, 1, v, z, **kwargs)

@defun
def lommels1(ctx, u, v, z, **kwargs):
    u = ctx._convert_param(u)[0]
    v = ctx._convert_param(v)[0]
    z = ctx.convert(z)
    def h(u,v):
        b = ctx.mpq_1_2
        w = ctx.square_exp_arg(z, mult=-0.25)
        return ([u-v+1, u+v+1, z], [-1, -1, u+1], [], [], [1], \
            [b*(u-v+3),b*(u+v+3)], w),
    return ctx.hypercomb(h, [u,v], **kwargs)

@defun
def lommels2(ctx, u, v, z, **kwargs):
    u = ctx._convert_param(u)[0]
    v = ctx._convert_param(v)[0]
    z = ctx.convert(z)
    # Asymptotic expansion (GR p. 947) -- need to be careful
    # not to use for small arguments
    # def h(u,v):
    #    b = ctx.mpq_1_2
    #    w = -(z/2)**(-2)
    #    return ([z], [u-1], [], [], [b*(1-u+v)], [b*(1-u-v)], w),
    def h(u,v):
        b = ctx.mpq_1_2
        w = ctx.square_exp_arg(z, mult=-0.25)
        T1 = [u-v+1, u+v+1, z], [-1, -1, u+1], [], [], [1], [b*(u-v+3),b*(u+v+3)], w
        T2 = [2, z], [u+v-1, -v], [v, b*(u+v+1)], [b*(v-u+1)], [], [1-v], w
        T3 = [2, z], [u-v-1, v], [-v, b*(u-v+1)], [b*(1-u-v)], [], [1+v], w
        #c1 = ctx.cospi((u-v)*b)
        #c2 = ctx.cospi((u+v)*b)
        #s = ctx.sinpi(v)
        #r1 = (u-v+1)*b
        #r2 = (u+v+1)*b
        #T2 = [c1, s, z, 2], [1, -1, -v, v], [], [-v+1], [], [-v+1], w
        #T3 = [-c2, s, z, 2], [1, -1, v, -v], [], [v+1], [], [v+1], w
        #T2 = [c1, s, z, 2], [1, -1, -v, v+u-1], [r1, r2], [-v+1], [], [-v+1], w
        #T3 = [-c2, s, z, 2], [1, -1, v, -v+u-1], [r1, r2], [v+1], [], [v+1], w
        return T1, T2, T3
    return ctx.hypercomb(h, [u,v], **kwargs)

@defun
def ber(ctx, n, z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/KelvinBer2/26/01/02/0001/
    def h(n):
        r = -(z/4)**4
        cos, sin = ctx.cospi_sinpi(-0.75*n)
        T1 = [cos, z/2], [1, n], [], [n+1], [], [0.5, 0.5*(n+1), 0.5*n+1], r
        T2 = [sin, z/2], [1, n+2], [], [n+2], [], [1.5, 0.5*(n+3), 0.5*n+1], r
        return T1, T2
    return ctx.hypercomb(h, [n], **kwargs)

@defun
def bei(ctx, n, z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/KelvinBei2/26/01/02/0001/
    def h(n):
        r = -(z/4)**4
        cos, sin = ctx.cospi_sinpi(0.75*n)
        T1 = [cos, z/2], [1, n+2], [], [n+2], [], [1.5, 0.5*(n+3), 0.5*n+1], r
        T2 = [sin, z/2], [1, n], [], [n+1], [], [0.5, 0.5*(n+1), 0.5*n+1], r
        return T1, T2
    return ctx.hypercomb(h, [n], **kwargs)

@defun
def ker(ctx, n, z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/KelvinKer2/26/01/02/0001/
    def h(n):
        r = -(z/4)**4
        cos1, sin1 = ctx.cospi_sinpi(0.25*n)
        cos2, sin2 = ctx.cospi_sinpi(0.75*n)
        T1 = [2, z, 4*cos1], [-n-3, n, 1], [-n], [], [], [0.5, 0.5*(1+n), 0.5*(n+2)], r
        T2 = [2, z, -sin1], [-n-3, 2+n, 1], [-n-1], [], [], [1.5, 0.5*(3+n), 0.5*(n+2)], r
        T3 = [2, z, 4*cos2], [n-3, -n, 1], [n], [], [], [0.5, 0.5*(1-n), 1-0.5*n], r
        T4 = [2, z, -sin2], [n-3, 2-n, 1], [n-1], [], [], [1.5, 0.5*(3-n), 1-0.5*n], r
        return T1, T2, T3, T4
    return ctx.hypercomb(h, [n], **kwargs)

@defun
def kei(ctx, n, z, **kwargs):
    n = ctx.convert(n)
    z = ctx.convert(z)
    # http://functions.wolfram.com/Bessel-TypeFunctions/KelvinKei2/26/01/02/0001/
    def h(n):
        r = -(z/4)**4
        cos1, sin1 = ctx.cospi_sinpi(0.75*n)
        cos2, sin2 = ctx.cospi_sinpi(0.25*n)
        T1 = [-cos1, 2, z], [1, n-3, 2-n], [n-1], [], [], [1.5, 0.5*(3-n), 1-0.5*n], r
        T2 = [-sin1, 2, z], [1, n-1, -n], [n], [], [], [0.5, 0.5*(1-n), 1-0.5*n], r
        T3 = [-sin2, 2, z], [1, -n-1, n], [-n], [], [], [0.5, 0.5*(n+1), 0.5*(n+2)], r
        T4 = [-cos2, 2, z], [1, -n-3, n+2], [-n-1], [], [], [1.5, 0.5*(n+3), 0.5*(n+2)], r
        return T1, T2, T3, T4
    return ctx.hypercomb(h, [n], **kwargs)

# TODO: do this more generically?
def c_memo(f):
    name = f.__name__
    def f_wrapped(ctx):
        cache = ctx._misc_const_cache
        prec = ctx.prec
        p,v = cache.get(name, (-1,0))
        if p >= prec:
            return +v
        else:
            cache[name] = (prec, f(ctx))
            return cache[name][1]
    return f_wrapped

@c_memo
def _airyai_C1(ctx):
    return 1 / (ctx.cbrt(9) * ctx.gamma(ctx.mpf(2)/3))

@c_memo
def _airyai_C2(ctx):
    return -1 / (ctx.cbrt(3) * ctx.gamma(ctx.mpf(1)/3))

@c_memo
def _airybi_C1(ctx):
    return 1 / (ctx.nthroot(3,6) * ctx.gamma(ctx.mpf(2)/3))

@c_memo
def _airybi_C2(ctx):
    return ctx.nthroot(3,6) / ctx.gamma(ctx.mpf(1)/3)

def _airybi_n2_inf(ctx):
    prec = ctx.prec
    try:
        v = ctx.power(3,'2/3')*ctx.gamma('2/3')/(2*ctx.pi)
    finally:
        ctx.prec = prec
    return +v

# Derivatives at z = 0
# TODO: could be expressed more elegantly using triple factorials
def _airyderiv_0(ctx, z, n, ntype, which):
    if ntype == 'Z':
        if n < 0:
            return z
        r = ctx.mpq_1_3
        prec = ctx.prec
        try:
            ctx.prec += 10
            v = ctx.gamma((n+1)*r) * ctx.power(3,n*r) / ctx.pi
            if which == 0:
                v *= ctx.sinpi(2*(n+1)*r)
                v /= ctx.power(3,'2/3')
            else:
                v *= abs(ctx.sinpi(2*(n+1)*r))
                v /= ctx.power(3,'1/6')
        finally:
            ctx.prec = prec
        return +v + z
    else:
        # singular (does the limit exist?)
        raise NotImplementedError

@defun
def airyai(ctx, z, derivative=0, **kwargs):
    z = ctx.convert(z)
    if derivative:
        n, ntype = ctx._convert_param(derivative)
    else:
        n = 0
    # Values at infinities
    if not ctx.isnormal(z) and z:
        if n and ntype == 'Z':
            if n == -1:
                if z == ctx.inf:
                    return ctx.mpf(1)/3 + 1/z
                if z == ctx.ninf:
                    return ctx.mpf(-2)/3 + 1/z
            if n < -1:
                if z == ctx.inf:
                    return z
                if z == ctx.ninf:
                    return (-1)**n * (-z)
        if (not n) and z == ctx.inf or z == ctx.ninf:
            return 1/z
        # TODO: limits
        raise ValueError("essential singularity of Ai(z)")
    # Account for exponential scaling
    if z:
        extraprec = max(0, int(1.5*ctx.mag(z)))
    else:
        extraprec = 0
    if n:
        if n == 1:
            def h():
                # http://functions.wolfram.com/03.07.06.0005.01
                if ctx._re(z) > 4:
                    ctx.prec += extraprec
                    w = z**1.5; r = -0.75/w; u = -2*w/3
                    ctx.prec -= extraprec
                    C = -ctx.exp(u)/(2*ctx.sqrt(ctx.pi))*ctx.nthroot(z,4)
                    return ([C],[1],[],[],[(-1,6),(7,6)],[],r),
                # http://functions.wolfram.com/03.07.26.0001.01
                else:
                    ctx.prec += extraprec
                    w = z**3 / 9
                    ctx.prec -= extraprec
                    C1 = _airyai_C1(ctx) * 0.5
                    C2 = _airyai_C2(ctx)
                    T1 = [C1,z],[1,2],[],[],[],[ctx.mpq_5_3],w
                    T2 = [C2],[1],[],[],[],[ctx.mpq_1_3],w
                    return T1, T2
            return ctx.hypercomb(h, [], **kwargs)
        else:
            if z == 0:
                return _airyderiv_0(ctx, z, n, ntype, 0)
            # http://functions.wolfram.com/03.05.20.0004.01
            def h(n):
                ctx.prec += extraprec
                w = z**3/9
                ctx.prec -= extraprec
                q13,q23,q43 = ctx.mpq_1_3, ctx.mpq_2_3, ctx.mpq_4_3
                a1=q13; a2=1; b1=(1-n)*q13; b2=(2-n)*q13; b3=1-n*q13
                T1 = [3, z], [n-q23, -n], [a1], [b1,b2,b3], \
                    [a1,a2], [b1,b2,b3], w
                a1=q23; b1=(2-n)*q13; b2=1-n*q13; b3=(4-n)*q13
                T2 = [3, z, -z], [n-q43, -n, 1], [a1], [b1,b2,b3], \
                    [a1,a2], [b1,b2,b3], w
                return T1, T2
            v = ctx.hypercomb(h, [n], **kwargs)
            if ctx._is_real_type(z) and ctx.isint(n):
                v = ctx._re(v)
            return v
    else:
        def h():
            if ctx._re(z) > 4:
                # We could use 1F1, but it results in huge cancellation;
                # the following expansion is better.
                # TODO: asymptotic series for derivatives
                ctx.prec += extraprec
                w = z**1.5; r = -0.75/w; u = -2*w/3
                ctx.prec -= extraprec
                C = ctx.exp(u)/(2*ctx.sqrt(ctx.pi)*ctx.nthroot(z,4))
                return ([C],[1],[],[],[(1,6),(5,6)],[],r),
            else:
                ctx.prec += extraprec
                w = z**3 / 9
                ctx.prec -= extraprec
                C1 = _airyai_C1(ctx)
                C2 = _airyai_C2(ctx)
                T1 = [C1],[1],[],[],[],[ctx.mpq_2_3],w
                T2 = [z*C2],[1],[],[],[],[ctx.mpq_4_3],w
                return T1, T2
        return ctx.hypercomb(h, [], **kwargs)

@defun
def airybi(ctx, z, derivative=0, **kwargs):
    z = ctx.convert(z)
    if derivative:
        n, ntype = ctx._convert_param(derivative)
    else:
        n = 0
    # Values at infinities
    if not ctx.isnormal(z) and z:
        if n and ntype == 'Z':
            if z == ctx.inf:
                return z
            if z == ctx.ninf:
                if n == -1:
                    return 1/z
                if n == -2:
                    return _airybi_n2_inf(ctx)
                if n < -2:
                    return (-1)**n * (-z)
        if not n:
            if z == ctx.inf:
                return z
            if z == ctx.ninf:
                return 1/z
        # TODO: limits
        raise ValueError("essential singularity of Bi(z)")
    if z:
        extraprec = max(0, int(1.5*ctx.mag(z)))
    else:
        extraprec = 0
    if n:
        if n == 1:
            # http://functions.wolfram.com/03.08.26.0001.01
            def h():
                ctx.prec += extraprec
                w = z**3 / 9
                ctx.prec -= extraprec
                C1 = _airybi_C1(ctx)*0.5
                C2 = _airybi_C2(ctx)
                T1 = [C1,z],[1,2],[],[],[],[ctx.mpq_5_3],w
                T2 = [C2],[1],[],[],[],[ctx.mpq_1_3],w
                return T1, T2
            return ctx.hypercomb(h, [], **kwargs)
        else:
            if z == 0:
                return _airyderiv_0(ctx, z, n, ntype, 1)
            def h(n):
                ctx.prec += extraprec
                w = z**3/9
                ctx.prec -= extraprec
                q13,q23,q43 = ctx.mpq_1_3, ctx.mpq_2_3, ctx.mpq_4_3
                q16 = ctx.mpq_1_6
                q56 = ctx.mpq_5_6
                a1=q13; a2=1; b1=(1-n)*q13; b2=(2-n)*q13; b3=1-n*q13
                T1 = [3, z], [n-q16, -n], [a1], [b1,b2,b3], \
                    [a1,a2], [b1,b2,b3], w
                a1=q23; b1=(2-n)*q13; b2=1-n*q13; b3=(4-n)*q13
                T2 = [3, z], [n-q56, 1-n], [a1], [b1,b2,b3], \
                    [a1,a2], [b1,b2,b3], w
                return T1, T2
            v = ctx.hypercomb(h, [n], **kwargs)
            if ctx._is_real_type(z) and ctx.isint(n):
                v = ctx._re(v)
            return v
    else:
        def h():
            ctx.prec += extraprec
            w = z**3 / 9
            ctx.prec -= extraprec
            C1 = _airybi_C1(ctx)
            C2 = _airybi_C2(ctx)
            T1 = [C1],[1],[],[],[],[ctx.mpq_2_3],w
            T2 = [z*C2],[1],[],[],[],[ctx.mpq_4_3],w
            return T1, T2
        return ctx.hypercomb(h, [], **kwargs)

def _airy_zero(ctx, which, k, derivative, complex=False):
    # Asymptotic formulas are given in DLMF section 9.9
    def U(t): return t**(2/3.)*(1-7/(t**2*48))
    def T(t): return t**(2/3.)*(1+5/(t**2*48))
    k = int(k)
    assert k >= 1
    assert derivative in (0,1)
    if which == 0:
        if derivative:
            return ctx.findroot(lambda z: ctx.airyai(z,1),
                -U(3*ctx.pi*(4*k-3)/8))
        return ctx.findroot(ctx.airyai, -T(3*ctx.pi*(4*k-1)/8))
    if which == 1 and complex == False:
        if derivative:
            return ctx.findroot(lambda z: ctx.airybi(z,1),
                -U(3*ctx.pi*(4*k-1)/8))
        return ctx.findroot(ctx.airybi, -T(3*ctx.pi*(4*k-3)/8))
    if which == 1 and complex == True:
        if derivative:
            t = 3*ctx.pi*(4*k-3)/8 + 0.75j*ctx.ln2
            s = ctx.expjpi(ctx.mpf(1)/3) * T(t)
            return ctx.findroot(lambda z: ctx.airybi(z,1), s)
        t = 3*ctx.pi*(4*k-1)/8 + 0.75j*ctx.ln2
        s = ctx.expjpi(ctx.mpf(1)/3) * U(t)
        return ctx.findroot(ctx.airybi, s)

@defun
def airyaizero(ctx, k, derivative=0):
    return _airy_zero(ctx, 0, k, derivative, False)

@defun
def airybizero(ctx, k, derivative=0, complex=False):
    return _airy_zero(ctx, 1, k, derivative, complex)

def _scorer(ctx, z, which, kwargs):
    z = ctx.convert(z)
    if ctx.isinf(z):
        if z == ctx.inf:
            if which == 0: return 1/z
            if which == 1: return z
        if z == ctx.ninf:
            return 1/z
        raise ValueError("essential singularity")
    if z:
        extraprec = max(0, int(1.5*ctx.mag(z)))
    else:
        extraprec = 0
    if kwargs.get('derivative'):
        raise NotImplementedError
    # Direct asymptotic expansions, to avoid
    # exponentially large cancellation
    try:
        if ctx.mag(z) > 3:
            if which == 0 and abs(ctx.arg(z)) < ctx.pi/3 * 0.999:
                def h():
                    return (([ctx.pi,z],[-1,-1],[],[],[(1,3),(2,3),1],[],9/z**3),)
                return ctx.hypercomb(h, [], maxterms=ctx.prec, force_series=True)
            if which == 1 and abs(ctx.arg(-z)) < 2*ctx.pi/3 * 0.999:
                def h():
                    return (([-ctx.pi,z],[-1,-1],[],[],[(1,3),(2,3),1],[],9/z**3),)
                return ctx.hypercomb(h, [], maxterms=ctx.prec, force_series=True)
    except ctx.NoConvergence:
        pass
    def h():
        A = ctx.airybi(z, **kwargs)/3
        B = -2*ctx.pi
        if which == 1:
            A *= 2
            B *= -1
        ctx.prec += extraprec
        w = z**3/9
        ctx.prec -= extraprec
        T1 = [A], [1], [], [], [], [], 0
        T2 = [B,z], [-1,2], [], [], [1], [ctx.mpq_4_3,ctx.mpq_5_3], w
        return T1, T2
    return ctx.hypercomb(h, [], **kwargs)

@defun
def scorergi(ctx, z, **kwargs):
    return _scorer(ctx, z, 0, kwargs)

@defun
def scorerhi(ctx, z, **kwargs):
    return _scorer(ctx, z, 1, kwargs)

@defun_wrapped
def coulombc(ctx, l, eta, _cache={}):
    if (l, eta) in _cache and _cache[l,eta][0] >= ctx.prec:
        return +_cache[l,eta][1]
    G3 = ctx.loggamma(2*l+2)
    G1 = ctx.loggamma(1+l+ctx.j*eta)
    G2 = ctx.loggamma(1+l-ctx.j*eta)
    v = 2**l * ctx.exp((-ctx.pi*eta+G1+G2)/2 - G3)
    if not (ctx.im(l) or ctx.im(eta)):
        v = ctx.re(v)
    _cache[l,eta] = (ctx.prec, v)
    return v

@defun_wrapped
def coulombf(ctx, l, eta, z, w=1, chop=True, **kwargs):
    # Regular Coulomb wave function
    # Note: w can be either 1 or -1; the other may be better in some cases
    # TODO: check that chop=True chops when and only when it should
    #ctx.prec += 10
    def h(l, eta):
        try:
            jw = ctx.j*w
            jwz = ctx.fmul(jw, z, exact=True)
            jwz2 = ctx.fmul(jwz, -2, exact=True)
            C = ctx.coulombc(l, eta)
            T1 = [C, z, ctx.exp(jwz)], [1, l+1, 1], [], [], [1+l+jw*eta], \
                [2*l+2], jwz2
        except ValueError:
            T1 = [0], [-1], [], [], [], [], 0
        return (T1,)
    v = ctx.hypercomb(h, [l,eta], **kwargs)
    if chop and (not ctx.im(l)) and (not ctx.im(eta)) and (not ctx.im(z)) and \
        (ctx.re(z) >= 0):
        v = ctx.re(v)
    return v

@defun_wrapped
def _coulomb_chi(ctx, l, eta, _cache={}):
    if (l, eta) in _cache and _cache[l,eta][0] >= ctx.prec:
        return _cache[l,eta][1]
    def terms():
        l2 = -l-1
        jeta = ctx.j*eta
        return [ctx.loggamma(1+l+jeta) * (-0.5j),
            ctx.loggamma(1+l-jeta) * (0.5j),
            ctx.loggamma(1+l2+jeta) * (0.5j),
            ctx.loggamma(1+l2-jeta) * (-0.5j),
            -(l+0.5)*ctx.pi]
    v = ctx.sum_accurately(terms, 1)
    _cache[l,eta] = (ctx.prec, v)
    return v

@defun_wrapped
def coulombg(ctx, l, eta, z, w=1, chop=True, **kwargs):
    # Irregular Coulomb wave function
    # Note: w can be either 1 or -1; the other may be better in some cases
    # TODO: check that chop=True chops when and only when it should
    if not ctx._im(l):
        l = ctx._re(l)  # XXX: for isint
    def h(l, eta):
        # Force perturbation for integers and half-integers
        if ctx.isint(l*2):
            T1 = [0], [-1], [], [], [], [], 0
            return (T1,)
        l2 = -l-1
        try:
            chi = ctx._coulomb_chi(l, eta)
            jw = ctx.j*w
            s = ctx.sin(chi); c = ctx.cos(chi)
            C1 = ctx.coulombc(l,eta)
            C2 = ctx.coulombc(l2,eta)
            u = ctx.exp(jw*z)
            x = -2*jw*z
            T1 = [s, C1, z, u, c], [-1, 1, l+1, 1, 1], [], [], \
                [1+l+jw*eta], [2*l+2], x
            T2 = [-s, C2, z, u],   [-1, 1, l2+1, 1],    [], [], \
                [1+l2+jw*eta], [2*l2+2], x
            return T1, T2
        except ValueError:
            T1 = [0], [-1], [], [], [], [], 0
            return (T1,)
    v = ctx.hypercomb(h, [l,eta], **kwargs)
    if chop and (not ctx._im(l)) and (not ctx._im(eta)) and (not ctx._im(z)) and \
        (ctx._re(z) >= 0):
        v = ctx._re(v)
    return v

def mcmahon(ctx,kind,prime,v,m):
    """
    Computes an estimate for the location of the Bessel function zero
    j_{v,m}, y_{v,m}, j'_{v,m} or y'_{v,m} using McMahon's asymptotic
    expansion (Abramowitz & Stegun 9.5.12-13, DLMF 20.21(vi)).

    Returns (r,err) where r is the estimated location of the root
    and err is a positive number estimating the error of the
    asymptotic expansion.
    """
    u = 4*v**2
    if kind == 1 and not prime: b = (4*m+2*v-1)*ctx.pi/4
    if kind == 2 and not prime: b = (4*m+2*v-3)*ctx.pi/4
    if kind == 1 and prime: b = (4*m+2*v-3)*ctx.pi/4
    if kind == 2 and prime: b = (4*m+2*v-1)*ctx.pi/4
    if not prime:
        s1 = b
        s2 = -(u-1)/(8*b)
        s3 = -4*(u-1)*(7*u-31)/(3*(8*b)**3)
        s4 = -32*(u-1)*(83*u**2-982*u+3779)/(15*(8*b)**5)
        s5 = -64*(u-1)*(6949*u**3-153855*u**2+1585743*u-6277237)/(105*(8*b)**7)
    if prime:
        s1 = b
        s2 = -(u+3)/(8*b)
        s3 = -4*(7*u**2+82*u-9)/(3*(8*b)**3)
        s4 = -32*(83*u**3+2075*u**2-3039*u+3537)/(15*(8*b)**5)
        s5 = -64*(6949*u**4+296492*u**3-1248002*u**2+7414380*u-5853627)/(105*(8*b)**7)
    terms = [s1,s2,s3,s4,s5]
    s = s1
    err = 0.0
    for i in range(1,len(terms)):
        if abs(terms[i]) < abs(terms[i-1]):
            s += terms[i]
        else:
            err = abs(terms[i])
    if i == len(terms)-1:
        err = abs(terms[-1])
    return s, err

def generalized_bisection(ctx,f,a,b,n):
    """
    Given f known to have exactly n simple roots within [a,b],
    return a list of n intervals isolating the roots
    and having opposite signs at the endpoints.

    TODO: this can be optimized, e.g. by reusing evaluation points.
    """
    assert n >= 1
    N = n+1
    points = []
    signs = []
    while 1:
        points = ctx.linspace(a,b,N)
        signs = [ctx.sign(f(x)) for x in points]
        ok_intervals = [(points[i],points[i+1]) for i in range(N-1) \
            if signs[i]*signs[i+1] == -1]
        if len(ok_intervals) == n:
            return ok_intervals
        N = N*2

def find_in_interval(ctx, f, ab):
    return ctx.findroot(f, ab, solver='illinois', verify=False)

def bessel_zero(ctx, kind, prime, v, m, isoltol=0.01, _interval_cache={}):
    prec = ctx.prec
    workprec = max(prec, ctx.mag(v), ctx.mag(m))+10
    try:
        ctx.prec = workprec
        v = ctx.mpf(v)
        m = int(m)
        prime = int(prime)
        assert v >= 0
        assert m >= 1
        assert prime in (0,1)
        if kind == 1:
            if prime: f = lambda x: ctx.besselj(v,x,derivative=1)
            else:     f = lambda x: ctx.besselj(v,x)
        if kind == 2:
            if prime: f = lambda x: ctx.bessely(v,x,derivative=1)
            else:     f = lambda x: ctx.bessely(v,x)
        # The first root of J' is very close to 0 for small
        # orders, and this needs to be special-cased
        if kind == 1 and prime and m == 1:
            if v == 0:
                return ctx.zero
            if v <= 1:
                # TODO: use v <= j'_{v,1} < y_{v,1}?
                r = 2*ctx.sqrt(v*(1+v)/(v+2))
                return find_in_interval(ctx, f, (r/10, 2*r))
        if (kind,prime,v,m) in _interval_cache:
            return find_in_interval(ctx, f, _interval_cache[kind,prime,v,m])
        r, err = mcmahon(ctx, kind, prime, v, m)
        if err < isoltol:
            return find_in_interval(ctx, f, (r-isoltol, r+isoltol))
        # An x such that 0 < x < r_{v,1}
        if kind == 1 and not prime: low = 2.4
        if kind == 1 and prime: low = 1.8
        if kind == 2 and not prime: low = 0.8
        if kind == 2 and prime: low = 2.0
        n = m+1
        while 1:
            r1, err = mcmahon(ctx, kind, prime, v, n)
            if err < isoltol:
                r2, err2 = mcmahon(ctx, kind, prime, v, n+1)
                intervals = generalized_bisection(ctx, f, low, 0.5*(r1+r2), n)
                for k, ab in enumerate(intervals):
                    _interval_cache[kind,prime,v,k+1] = ab
                return find_in_interval(ctx, f, intervals[m-1])
            else:
                n = n*2
    finally:
        ctx.prec = prec

@defun
def besseljzero(ctx, v, m, derivative=0):
    r"""
    For a real order `\nu \ge 0` and a positive integer `m`, returns
    `j_{\nu,m}`, the `m`-th positive zero of the Bessel function of the
    first kind `J_{\nu}(z)` (see :func:`~mpmath.besselj`). Alternatively,
    with *derivative=1*, gives the first nonnegative simple zero
    `j'_{\nu,m}` of `J'_{\nu}(z)`.

    The indexing convention is that used by Abramowitz & Stegun
    and the DLMF. Note the special case `j'_{0,1} = 0`, while all other
    zeros are positive. In effect, only simple zeros are counted
    (all zeros of Bessel functions are simple except possibly `z = 0`)
    and `j_{\nu,m}` becomes a monotonic function of both `\nu`
    and `m`.

    The zeros are interlaced according to the inequalities

    .. math ::

        j'_{\nu,k} < j_{\nu,k} < j'_{\nu,k+1}

        j_{\nu,1} < j_{\nu+1,2} < j_{\nu,2} < j_{\nu+1,2} < j_{\nu,3} < \cdots

    **Examples**

    Initial zeros of the Bessel functions `J_0(z), J_1(z), J_2(z)`::

        >>> from mpmath import *
        >>> mp.dps = 25; mp.pretty = True
        >>> besseljzero(0,1); besseljzero(0,2); besseljzero(0,3)
        2.404825557695772768621632
        5.520078110286310649596604
        8.653727912911012216954199
        >>> besseljzero(1,1); besseljzero(1,2); besseljzero(1,3)
        3.831705970207512315614436
        7.01558666981561875353705
        10.17346813506272207718571
        >>> besseljzero(2,1); besseljzero(2,2); besseljzero(2,3)
        5.135622301840682556301402
        8.417244140399864857783614
        11.61984117214905942709415

    Initial zeros of `J'_0(z), J'_1(z), J'_2(z)`::

        0.0
        3.831705970207512315614436
        7.01558666981561875353705
        >>> besseljzero(1,1,1); besseljzero(1,2,1); besseljzero(1,3,1)
        1.84118378134065930264363
        5.331442773525032636884016
        8.536316366346285834358961
        >>> besseljzero(2,1,1); besseljzero(2,2,1); besseljzero(2,3,1)
        3.054236928227140322755932
        6.706133194158459146634394
        9.969467823087595793179143

    Zeros with large index::

        >>> besseljzero(0,100); besseljzero(0,1000); besseljzero(0,10000)
        313.3742660775278447196902
        3140.807295225078628895545
        31415.14114171350798533666
        >>> besseljzero(5,100); besseljzero(5,1000); besseljzero(5,10000)
        321.1893195676003157339222
        3148.657306813047523500494
        31422.9947255486291798943
        >>> besseljzero(0,100,1); besseljzero(0,1000,1); besseljzero(0,10000,1)
        311.8018681873704508125112
        3139.236339643802482833973
        31413.57032947022399485808

    Zeros of functions with large order::

        >>> besseljzero(50,1)
        57.11689916011917411936228
        >>> besseljzero(50,2)
        62.80769876483536093435393
        >>> besseljzero(50,100)
        388.6936600656058834640981
        >>> besseljzero(50,1,1)
        52.99764038731665010944037
        >>> besseljzero(50,2,1)
        60.02631933279942589882363
        >>> besseljzero(50,100,1)
        387.1083151608726181086283

    Zeros of functions with fractional order::

        >>> besseljzero(0.5,1); besseljzero(1.5,1); besseljzero(2.25,4)
        3.141592653589793238462643
        4.493409457909064175307881
        15.15657692957458622921634

    Both `J_{\nu}(z)` and `J'_{\nu}(z)` can be expressed as infinite
    products over their zeros::

        >>> v,z = 2, mpf(1)
        >>> (z/2)**v/gamma(v+1) * \
        ...     nprod(lambda k: 1-(z/besseljzero(v,k))**2, [1,inf])
        ...
        0.1149034849319004804696469
        >>> besselj(v,z)
        0.1149034849319004804696469
        >>> (z/2)**(v-1)/2/gamma(v) * \
        ...     nprod(lambda k: 1-(z/besseljzero(v,k,1))**2, [1,inf])
        ...
        0.2102436158811325550203884
        >>> besselj(v,z,1)
        0.2102436158811325550203884

    """
    return +bessel_zero(ctx, 1, derivative, v, m)

@defun
def besselyzero(ctx, v, m, derivative=0):
    r"""
    For a real order `\nu \ge 0` and a positive integer `m`, returns
    `y_{\nu,m}`, the `m`-th positive zero of the Bessel function of the
    second kind `Y_{\nu}(z)` (see :func:`~mpmath.bessely`). Alternatively,
    with *derivative=1*, gives the first positive zero `y'_{\nu,m}` of
    `Y'_{\nu}(z)`.

    The zeros are interlaced according to the inequalities

    .. math ::

        y_{\nu,k} < y'_{\nu,k} < y_{\nu,k+1}

        y_{\nu,1} < y_{\nu+1,2} < y_{\nu,2} < y_{\nu+1,2} < y_{\nu,3} < \cdots

    **Examples**

    Initial zeros of the Bessel functions `Y_0(z), Y_1(z), Y_2(z)`::

        >>> from mpmath import *
        >>> mp.dps = 25; mp.pretty = True
        >>> besselyzero(0,1); besselyzero(0,2); besselyzero(0,3)
        0.8935769662791675215848871
        3.957678419314857868375677
        7.086051060301772697623625
        >>> besselyzero(1,1); besselyzero(1,2); besselyzero(1,3)
        2.197141326031017035149034
        5.429681040794135132772005
        8.596005868331168926429606
        >>> besselyzero(2,1); besselyzero(2,2); besselyzero(2,3)
        3.384241767149593472701426
        6.793807513268267538291167
        10.02347797936003797850539

    Initial zeros of `Y'_0(z), Y'_1(z), Y'_2(z)`::

        >>> besselyzero(0,1,1); besselyzero(0,2,1); besselyzero(0,3,1)
        2.197141326031017035149034
        5.429681040794135132772005
        8.596005868331168926429606
        >>> besselyzero(1,1,1); besselyzero(1,2,1); besselyzero(1,3,1)
        3.683022856585177699898967
        6.941499953654175655751944
        10.12340465543661307978775
        >>> besselyzero(2,1,1); besselyzero(2,2,1); besselyzero(2,3,1)
        5.002582931446063945200176
        8.350724701413079526349714
        11.57419546521764654624265

    Zeros with large index::

        >>> besselyzero(0,100); besselyzero(0,1000); besselyzero(0,10000)
        311.8034717601871549333419
        3139.236498918198006794026
        31413.57034538691205229188
        >>> besselyzero(5,100); besselyzero(5,1000); besselyzero(5,10000)
        319.6183338562782156235062
        3147.086508524556404473186
        31421.42392920214673402828
        >>> besselyzero(0,100,1); besselyzero(0,1000,1); besselyzero(0,10000,1)
        313.3726705426359345050449
        3140.807136030340213610065
        31415.14112579761578220175

    Zeros of functions with large order::

        >>> besselyzero(50,1)
        53.50285882040036394680237
        >>> besselyzero(50,2)
        60.11244442774058114686022
        >>> besselyzero(50,100)
        387.1096509824943957706835
        >>> besselyzero(50,1,1)
        56.96290427516751320063605
        >>> besselyzero(50,2,1)
        62.74888166945933944036623
        >>> besselyzero(50,100,1)
        388.6923300548309258355475

    Zeros of functions with fractional order::

        >>> besselyzero(0.5,1); besselyzero(1.5,1); besselyzero(2.25,4)
        1.570796326794896619231322
        2.798386045783887136720249
        13.56721208770735123376018

    """
    return +bessel_zero(ctx, 2, derivative, v, m)

max_stars_repo_name: benjaminmcdonald/sympy

max_issues_repo_name: jegerjensen/sympy

max_forks_repo_name: jegerjensen/sympy

reason: from mpmath


text:
from __future__ import division
import os
import shutil
import copy
from tqdm import tqdm
import time
import numpy as np
from PIL import Image, ImageDraw

def MaxSubarrayFL(array, w, T):
    for i in range(1, len(array)): array[i] += array[i-1]
    j1 = -1
    Smax = -1
    wr = -1
    w1, w2 = min(w), max(w)
    if T > 0 and w1 > 0:
        for st in range(1, len(array)-w2+1):
            for offset in range(w2-w1+1):
                S0 = array[st+offset+w1-1] - (array[st-1] if st>=1 else 0)
                if S0 >= T and S0 > Smax:
                    j1 = st
                    wr = w1 + offset
                    Smax = S0
    return j1, wr, Smax

def FixedAspRatioRectangle(G, alpha, ratio):
    m, n = G.shape
    GcSum = copy.deepcopy(G)
    for i in range(1, len(G)):
        GcSum[i] += GcSum[i-1]
    GSum = copy.deepcopy(GcSum)
    for i in range(1, len(G[0])):
        GSum[:,i] += GSum[:,i-1]
    i, j, w, h = 0, 0, n, m
    i1, i2, T, Smin = 0, 0, alpha*GSum[m-1,n-1], -1

    if GSum[m-1,n-1]/(255*m*n) < 0.1:
        return i, j, w, h, 0.10

    while i1 < m-1 and i2 < m:
        h0 = i2 - i1 + 1
        w0 = int(h0 / ratio)
        if w0 > n:
            i1 += 1
        else:
            array = GcSum[i2, :] - (GcSum[i1-1, :] if i1 >=1 else 0)
            j1, wr, S0 = MaxSubarrayFL(array, (h0, w0), T)
            if j1 > -1:
                if wr*h0 < w*h or (wr*h0 == w*h and S0 < Smin):
                    i, j, w, h = i1, j1, wr, h0
                    Smin = S0
                i1 += 1
            else:
                i2 += 1
    return i, j, w, h, Smin if Smin == -1 else Smin/GSum[m-1,n-1]

if os.path.exists('./crops'):
    shutil.rmtree('./crops')
os.mkdir('./crops')

filenames = sorted(os.listdir('./masks'))
for i in tqdm(range(len(filenames))):
    filename = filenames[i]
    image = Image.open('../test_images/{:s}'.format(filename))
    ori_shape = image.size
    hw_ratio = ori_shape[1] / ori_shape[0]
    if ori_shape[0] < ori_shape[1]:
        mask_shape = (224, int(224*hw_ratio))
    else:
        mask_shape = (int(224/hw_ratio), 224)

    mask = Image.open('./masks/{:s}'.format(filename)).resize(mask_shape, Image.ANTIALIAS)
    npMask = np.array(mask.convert('L')).astype(np.float32)

    #npMask[npMask < 100] = 0
    #npMask[npMask > 100] = 255

    i, j, w, h, trueAlpha = FixedAspRatioRectangle(npMask, 0.9, 4/3)
    w_resize_ratio, h_resize_ratio = np.array(ori_shape) / np.array(mask_shape)
    i, j, w, h = i*h_resize_ratio, j*w_resize_ratio, w*w_resize_ratio, h*h_resize_ratio

    draw = ImageDraw.Draw(image)
    draw.rectangle((j, i, j+w, i+h), outline=(255,0,0), width=2)
    image.save('crops/{:s}_{:.2f}.jpg'.format(filename[:-4], trueAlpha))


text:
from keras import backend as K
from keras.engine.topology import Layer
import numpy as np
from utilities import gram_matrix_sum

class pretrained_layer(Layer):
    def __init__(self, output_dim, pretrained_layer_functors, **kwargs):
        self.output_dim = output_dim
        self.pretrained_layer_functors = pretrained_layer_functors
        super(pretrained_layer, self).__init__(**kwargs)

    def build(self, input_shape):
        super(pretrained_layer, self).build(input_shape)

    def call(self, x):
        rssults = np.array([np.array(func([x, 1.]))[0][0] for func in self.pretrained_layer_functors])
        preds = []
        for layer in results[0]:
            preds.append(gram_matrix_sum(layer))
        print(preds)
        return preds

    def compute_output_shape(self, input_shape):
        return (input_shape[0], self.output_dim)


text:
# Copyright (c) Facebook, Inc. and its affiliates.

"""
Renders mesh using OpenDr / Pytorch-3d for visualization.
Part of code is modified from https://github.com/akanazawa/hmr
"""

import sys
import numpy as np
import cv2
import pdb
from PIL import Image, ImageDraw
from opendr.camera import ProjectPoints
from opendr.renderer import ColoredRenderer
from opendr.lighting import LambertianPointLight

class OpendrRenderer(object):
    def __init__(self,
                 img_size=224,
                 mesh_color=np.array([0.5, 0.5, 0.5]),):

        self.w = img_size
        self.h = img_size
        self.color = mesh_color
        self.img_size = img_size
        self.flength = 500.

    def render(self, verts, faces, bg_img):
        verts = verts.copy()
        faces = faces.copy()

        input_size = 500

        f = 10

        verts[:, 0] = (verts[:, 0] - input_size) / input_size
        verts[:, 1] = (verts[:, 1] - input_size) / input_size

        verts[:, 2] /= (5 * 112)
        verts[:, 2] += f

        cam_for_render = np.array([f, 1, 1]) * input_size

        rend_img = self.__call__(
            img=bg_img, cam=cam_for_render,
            verts=verts, faces=faces, color=self.color)

        return rend_img

    def __call__(self,
                 verts,
                 faces,
                 cam=None,
                 img=None,
                 do_alpha=False,
                 far=None,
                 near=None,
                 color = np.array([0, 0, 255]),
                 img_size=None):
        """
        cam is 3D [f, px, py]
        """
        if img is not None:
            h, w = img.shape[:2]
        elif img_size is not None:
            h = img_size[0]
            w = img_size[1]
        else:
            h = self.h
            w = self.w

        if cam is None:
            cam = [self.flength, w / 2., h / 2.]

        use_cam = ProjectPoints(
            f=cam[0] * np.ones(2),
            rt=np.zeros(3),
            t=np.zeros(3),
            k=np.zeros(5),
            c=cam[1:3])

        if near is None:
            near = np.maximum(np.min(verts[:, 2]) - 25, 0.1)
        if far is None:
            far = np.maximum(np.max(verts[:, 2]) + 25, 25)

        return_value = render_model(
            verts,
            faces,
            w,
            h,
            use_cam,
            do_alpha=do_alpha,
            img=img,
            far=far,
            near=near,
            color=color)

        imtmp = return_value
        image = (imtmp * 255).astype('uint8')
        return image

def _create_renderer(w=640,
                     h=480,
                     rt=np.zeros(3),
                     t=np.zeros(3),
                     f=None,
                     c=None,
                     k=None,
                     near=.5,
                     far=10.):

    f = np.array([w, w]) / 2. if f is None else f
    c = np.array([w, h]) / 2. if c is None else c
    k = np.zeros(5) if k is None else k

    rn = ColoredRenderer()

    rn.camera = ProjectPoints(rt=rt, t=t, f=f, c=c, k=k)
    rn.frustum = {'near': near, 'far': far, 'height': h, 'width': w}
    return rn

def _rotateY(points, angle):
    """Rotate the points by a specified angle."""
    ry = np.array([[np.cos(angle), 0., np.sin(angle)], [0., 1., 0.],
                   [-np.sin(angle), 0., np.cos(angle)]])
    return np.dot(points, ry)

def simple_renderer(rn,
                    verts,
                    faces,
                    yrot=np.radians(70),
                    color=np.array([0, 0, 255])
                    ):

    # Rendered model color
    rn.set(v=verts, f=faces, vc=color, bgcolor=np.ones(3))
    albedo = rn.vc

    # Construct Back Light (on back right corner)
    rn.vc = LambertianPointLight(
        f=rn.f,
        v=rn.v,
        num_verts=len(rn.v),
        light_pos=_rotateY(np.array([-200, -100, -100]), yrot),
        vc=albedo,
        light_color=np.array([1, 1, 1]))

    # Construct Left Light
    rn.vc += LambertianPointLight(
        f=rn.f,
        v=rn.v,
        num_verts=len(rn.v),
        # light_pos=_rotateY(np.array([800, 10, 300]), yrot),
        light_pos=_rotateY(np.array([800, 10, 300]), yrot),
        vc=albedo,
        light_color=np.array([1, 1, 1]))

    # Construct Right Light
    rn.vc += LambertianPointLight(
        f=rn.f,
        v=rn.v,
        num_verts=len(rn.v),
        light_pos=_rotateY(np.array([-500, 500, 1000]), yrot),
        # light_pos=_rotateY(np.array([-500, 500, 1000]), yrot),
        vc=albedo,
        light_color=np.array([.7, .7, .7]))

    return rn.r

def get_alpha(imtmp, bgval=1.):
    h, w = imtmp.shape[:2]
    alpha = (~np.all(imtmp == bgval, axis=2)).astype(imtmp.dtype)

    b_channel, g_channel, r_channel = cv2.split(imtmp)

    im_RGBA = cv2.merge((b_channel, g_channel, r_channel, alpha.astype(
        imtmp.dtype)))
    return im_RGBA

def append_alpha(imtmp):
    alpha = np.ones_like(imtmp[:, :, 0]).astype(imtmp.dtype)
    if np.issubdtype(imtmp.dtype, np.uint8):
        alpha = alpha * 255
    b_channel, g_channel, r_channel = cv2.split(imtmp)
    im_RGBA = cv2.merge((b_channel, g_channel, r_channel, alpha))
    return im_RGBA

def render_model(verts,
                 faces,
                 w,
                 h,
                 cam,
                 near=0.5,
                 far=25,
                 img=None,
                 do_alpha=False,
                 color=None):
    rn = _create_renderer(
        w=w, h=h, near=near, far=far, rt=cam.rt, t=cam.t, f=cam.f, c=cam.c)

    # Uses img as background, otherwise white background.
    if img is not None:
        rn.background_image = img / 255. if img.max() > 1.1 else img

    imtmp = simple_renderer(rn, verts, faces, color=color)

    # If white bg, make transparent.
    if img is None and do_alpha:
        imtmp = get_alpha(imtmp)
    elif img is not None and do_alpha:
        imtmp = append_alpha(imtmp)

    return imtmp

max_stars_repo_name: archonic/frankmocap

max_issues_repo_name: archonic/frankmocap

max_forks_repo_name: archonic/frankmocap