text:
from picamera.array import PiRGBArray
from picamera import PiCamera

from networktables import NetworkTables as nt

import time
import cv2
import numpy as np
import math

import socket

# camera setup

IMG_WIDTH = 320
IMG_HEIGHT = 240

camera = PiCamera(resolution=str(IMG_WIDTH) + "x" + str(IMG_HEIGHT), framerate = 60) # use 60 fps to avoid issues with rolling shutter
#camera.exposure_compensation = -25 # not sure if necessary
#camera.exposure_mode = 'off' # avoid automatic exposure compensation
#camera.shutter_speed = 5555 # low shutter speed to lower lighting besides the retroreflective tape
rawCapture = PiRGBArray(camera)
nt.initialize(server='roborio-424-frc.local')

time.sleep(3) # camera warmup, position setup

UDP_IP = "0.0.0.0"
UDP_PORT = 5803

sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
sock.bind((UDP_IP, UDP_PORT))

startTime = time.time() # for timing purposes

# kernel for dilate/erode
kernel = np.ones((10,10), np.uint8);

#font for debug drawing
font = cv2.FONT_HERSHEY_SIMPLEX

# main execution

color = 0;

vi = nt.getTable('vision')

for image in camera.capture_continuous(rawCapture, format="bgr", use_video_port=True):
    startTime = time.time()
    # now grab frame and store it
    frame = image.array
    rawCapture.truncate(0) # important

    # convert images to HSV color spectrum
    hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)

    # hsv thresholds
    lower = np.array([21,90,83])
    upper = np.array([39,255,255])

    thresh = cv2.inRange(hsv, lower, upper)

    # dilate and erode
    dilation = cv2.dilate(thresh, kernel, iterations = 1)

    erode = cv2.erode(dilation, kernel, iterations = 1)

    thresh = erode

    # find contours
    im2, contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    centerx = -1
    centery = 0
    width = -1
    contour = -1

    # sock.sendto(str(160),   ("roborio-424-frc.local", 5803))
    debug = frame

    # filter contours
    for i in range(0, len(contours)):
         if cv2.contourArea(contours[i]) > 1500:
            x,y,w,h = cv2.boundingRect(contours[i])
            rRect = cv2.minAreaRect(contours[i])
            rBox = cv2.boxPoints(rRect)
            rBox = np.int0(rBox)
            cv2.drawContours(debug, [rBox], 0, (0,255,0), 2)
            cv2.drawContours(debug, contours, i, (0, 0, 255), 1)
            cbb = cv2.contourArea(contours[i])/(rRect[1][0]*rRect[1][1])
            cv2.putText(debug, str(cbb)[:5], (x+w/2, y+h/2), font, 0.5, (255, 0, 0), 1, cv2.LINE_AA)
            # closest to center ?
            if(abs(IMG_WIDTH/2 - (x + w/2)) < abs(IMG_WIDTH/2 - centerx)):
                width = w
                centerx = x + w/2 # info about the target
                centery = y + h/2 # info about the target
                contour = i

    print(centerx)
    vi.putNumber('centerx', centerx);
    vi.putNumber('width', width);
    sock.sendto(str(centerx) + "," + str(width),   ("roborio-424-frc.local", 5803))

    print("--- %s seconds processing ---" % (time.time() - startTime))
    # save debug images
    cv2.imwrite("cube_orig.jpg", frame)

    cv2.drawContours(frame, contours, -1, (0, 0, 255), 2) # draw all contours in red
    cv2.imwrite('cube_found.jpg', frame)
    cv2.imwrite('cube_debug.jpg', debug)
    # cv2.imwrite("cube_dilate.jpg", dilation)
    # cv2.imwrite("cube_erode.jpg", erode)

print("--- %s seconds capturing ---" % (time.time() - startTime))
processingStart = time.time();

# end main execution

print("--- %s seconds processing ---" % (time.time() - processingStart))
print("--- %s seconds total ---" % (time.time() - startTime))

camera.close(); # important
"""
while True:
    led.fill((255, 0, 0))
    led.update()
    time.sleep(DELAY)
    led.fill((0, 255, 0))
    led.update()
    time.sleep(DELAY)
    led.fill((0, 0, 255))
    led.update()
    time.sleep(DELAY)
"""


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.

"""Utilities for data loading that are shared across multiple datasets.

Some datasets are very similar, so to prevent code duplication, shared utilities
are put into this class.
"""

# pylint:skip-file
import sys
import time
from collections import defaultdict

import numpy as np
import tensorflow as tf

from smurf import smurf_utils
from smurf import smurf_plotting

def parse_data(proto,
               include_flow,
               height=None,
               width=None,
               include_occlusion=False,
               include_invalid=False,
               resize_gt_flow=True,
               include_image_path=False,
               gt_flow_shape=None,
               include_segments=False):
  """Parse a data proto with flow.

  Args:
    proto: path to data proto file
    include_flow: bool, whether or not to include flow in the output
    height: int or None height to resize image to
    width: int or None width to resize image to
    include_occlusion: bool, whether or not to also return occluded pixels (will
      throw error if occluded pixels are not present)
    include_invalid: bool, whether or not to also return invalid pixels (will
      throw error if invalid pixels are not present)
    resize_gt_flow: bool, wether or not to resize flow ground truth as the image
    include_image_path: bool, if True, return the string for the key
      "image1_path" alongside the data.
    gt_flow_shape: list, shape of the original ground truth flow (only required
      to set a fixed ground truth flow shape for tensorflow estimator in case of
      supervised training at full resolution resize_gt_flow=False)
    include_segments: bool, if True, include the Sintel segmentation data.

  Returns:
    images, flow: A tuple of (image1, image2), flow
  """

  # Parse context and image sequence from protobuffer.
  context_features = {
      'height': tf.io.FixedLenFeature([], tf.int64),
      'width': tf.io.FixedLenFeature([], tf.int64),
  }
  sequence_features = {
      'images': tf.io.FixedLenSequenceFeature([], tf.string),
  }

  if include_invalid:
    sequence_features['invalid_masks'] = tf.io.FixedLenSequenceFeature(
        [], tf.string)

  if include_segments:
    sequence_features['segments'] = tf.io.FixedLenSequenceFeature(
        [], tf.string)
    sequence_features['segments_invalid'] = tf.io.FixedLenSequenceFeature(
        [], tf.string)

  if include_image_path:
    context_features['image1_path'] = tf.io.FixedLenFeature((), tf.string)

  if include_flow:
    context_features['flow_uv'] = tf.io.FixedLenFeature([], tf.string)

  if include_occlusion:
    context_features['occlusion_mask'] = tf.io.FixedLenFeature([], tf.string)

  context_parsed, sequence_parsed = tf.io.parse_single_sequence_example(
      proto,
      context_features=context_features,
      sequence_features=sequence_features,
  )

  def deserialize(s, dtype, dims):
    return tf.reshape(
        tf.io.decode_raw(s, dtype),
        [context_parsed['height'], context_parsed['width'], dims])

  images = tf.map_fn(
      lambda s: deserialize(s, tf.uint8, 3),
      sequence_parsed['images'],
      dtype=tf.uint8)

  images = tf.image.convert_image_dtype(images, tf.float32)
  if height is not None and width is not None:
    images = smurf_utils.resize(images, height, width, is_flow=False)
  output = {'images': images}

  if include_flow:
    flow_uv = deserialize(context_parsed['flow_uv'], tf.float32, 2)
    flow_uv = flow_uv[Ellipsis, ::-1]
    # Flying things has some images with erroneously large flow.
    # Mask out any values above / below 1000.
    invalid_cond = tf.math.logical_or(tf.greater(flow_uv, 1000),
                                      tf.less(flow_uv, -1000))
    mask = tf.where(invalid_cond, tf.zeros_like(flow_uv), tf.ones_like(flow_uv))
    flow_valid = tf.reduce_min(mask, axis=-1, keepdims=True)
    if height is not None and width is not None and resize_gt_flow:
      flow_uv = smurf_utils.resize(flow_uv, height, width, is_flow=True)
      flow_valid = smurf_utils.resize(flow_valid, height, width, is_flow=False)
    else:
      if gt_flow_shape is not None:
        flow_uv.set_shape(gt_flow_shape)
        flow_valid.set_shape((gt_flow_shape[0], gt_flow_shape[1], 1))
    # To be consistent with SMURF internals, we flip the ordering of flow.
    # create valid mask
    flow_valid = tf.ones_like(flow_uv[Ellipsis, :1], dtype=tf.float32)
    output['flow_valid'] = flow_valid
    output['flow'] = flow_uv

  if include_occlusion:
    occlusion_mask = deserialize(context_parsed['occlusion_mask'], tf.uint8, 1)
    if height is not None and width is not None:
      occlusion_mask = smurf_utils.resize(
          occlusion_mask, height, width, is_flow=False)
    output['occlusions'] = occlusion_mask

  if include_invalid:
    invalid_masks = tf.map_fn(
        lambda s: deserialize(s, tf.uint8, 1),
        sequence_parsed['invalid_masks'],
        dtype=tf.uint8)
    if height is not None and width is not None:
      invalid_masks = smurf_utils.resize(
          invalid_masks, height, width, is_flow=False)
    output['flow_valid'] = 1. - invalid_masks

  if include_image_path:
    output['image1_path'] = context_parsed['image1_path']

  if include_segments:
    segments = tf.map_fn(
        lambda s: deserialize(s, tf.uint8, 3),
        sequence_parsed['segments'],
        dtype=tf.uint8)
    segments = tf.image.convert_image_dtype(segments, tf.float32)
    segments_invalid = tf.map_fn(
        lambda s: deserialize(s, tf.uint8, 1),
        sequence_parsed['segments_invalid'],
        dtype=tf.uint8)
    segments_invalid = tf.image.convert_image_dtype(segments_invalid,
                                                    tf.float32)
    segments = tf.image.resize(segments, (height, width), method='nearest')
    segments_invalid = tf.image.resize(
        segments_invalid, (height, width), method='nearest')
    output['segments'] = segments
    output['segments_invalid'] = segments_invalid

  return output

def evaluate(
    inference_fn,
    dataset,
    height,
    width,
    progress_bar=False,
    plot_dir='',
    num_plots=0,
    max_num_evals=10000,
    prefix='',
    has_occlusion=True,
    weights=None,
):
  """Evaluate an inference function for flow.

  Args:
    inference_fn: An inference function that produces a flow_field from two
      images, e.g. the infer method of SMURF.
    dataset: A dataset produced by the method above with for_eval=True.
    height: int, the height to which the images should be resized for inference.
    width: int, the width to which the images should be resized for inference.
    progress_bar: boolean, flag to indicate whether the function should print a
      progress_bar during evaluaton.
    plot_dir: string, optional path to a directory in which plots are saved (if
      num_plots > 0).
    num_plots: int, maximum number of qualitative results to plot for the
      evaluation.
    max_num_evals: int, maxmim number of evaluations.
    prefix: str, prefix to prepend to filenames for saved plots and for keys in
      results dictionary.
    has_occlusion: bool indicating whether or not the dataset includes ground
      truth occlusion.
    weights: unsupervised loss weights

  Returns:
    A dictionary of floats that represent different evaluation metrics. The keys
    of this dictionary are returned by the method list_eval_keys (see below).
  """

  eval_start_in_s = time.time()

  it = tf.compat.v1.data.make_one_shot_iterator(dataset)
  epe_occ = []  # End point errors.
  errors_occ = []
  inference_times = []
  unsuper_losses = []
  all_occlusion_results = defaultdict(lambda: defaultdict(int))

  plot_count = 0
  eval_count = -1
  for test_batch in it:

    image_batch = test_batch['images']
    flow_gt = test_batch['flow']
    flow_valid = test_batch['flow_valid']
    if has_occlusion:
      occ_mask_gt = test_batch['occlusions']
    else:
      occ_mask_gt = tf.ones_like(flow_valid)

    if eval_count >= max_num_evals:
      break

    eval_count += 1
    if eval_count >= max_num_evals:
      break

    if progress_bar:
      sys.stdout.write(':')
      sys.stdout.flush()

    f = lambda: inference_fn(
        image_batch[0],
        image_batch[1],
        input_height=height,
        input_width=width,
        infer_occlusion=True,
        infer_bw=True)

    inference_time_in_ms, (flow, soft_occlusion_mask,
                           flow_bw) = smurf_utils.time_it(
                               f, execute_once_before=eval_count == 1)
    inference_times.append(inference_time_in_ms)

    if not has_occlusion:
      best_thresh = .5
    else:
      f_dict = compute_f_metrics(soft_occlusion_mask, occ_mask_gt)
      best_thresh = -1.
      best_f_score = -1.
      for thresh, metrics in f_dict.items():
        precision = metrics['tp'] / (metrics['tp'] + metrics['fp'] + 1e-6)
        recall = metrics['tp'] / (metrics['tp'] + metrics['fn'] + 1e-6)
        f1 = 2 * precision * recall / (precision + recall + 1e-6)
        if f1 > best_f_score:
          best_thresh = thresh
          best_f_score = f1
        all_occlusion_results[thresh]['tp'] += metrics['tp']
        all_occlusion_results[thresh]['fp'] += metrics['fp']
        all_occlusion_results[thresh]['tn'] += metrics['tn']
        all_occlusion_results[thresh]['fn'] += metrics['fn']

    final_flow = flow
    endpoint_error_occ = epe_elementwise(final_flow, flow_gt) * flow_valid
    outliers_occ = outliers_elementwise(final_flow, flow_gt) * flow_valid
    # NOTE(austinstone): The unsupervised loss function expects occluded areas
    # to be zeros, whereas the occlusion mask above returns 1s in areas
    # of occlusions. The unsupervised loss function below assumes some
    # parameters (smoothness edge info) is left at the default values.
    if weights is not None:
      inv_mask = tf.expand_dims(1. - soft_occlusion_mask, axis=0)
      occlusion_estimation_fn = lambda forward_flow, backward_flow: inv_mask
      loss_dict = smurf_utils.unsupervised_loss(
          images=tf.expand_dims(image_batch, axis=0),
          flows={
              (0, 1, 'augmented-student'): [tf.expand_dims(flow, axis=0)],
              (1, 0, 'augmented-student'): [tf.expand_dims(flow_bw, axis=0)]
          },
          weights=weights,
          occlusion_estimation_fn=occlusion_estimation_fn)
      loss_dict['total_loss'] = sum(loss_dict.values())
    else:
      loss_dict = {}
    unsuper_losses.append(loss_dict)
    epe_occ.append(tf.reduce_mean(input_tensor=endpoint_error_occ))
    errors_occ.append(tf.reduce_mean(input_tensor=outliers_occ))

    if plot_dir and plot_count < num_plots:
      plot_count += 1
      mask_thresh = tf.cast(
          tf.math.greater(soft_occlusion_mask, best_thresh), tf.float32)
      smurf_plotting.complete_paper_plot(
          plot_dir,
          plot_count,
          image_batch[0].numpy(),
          image_batch[1].numpy(),
          final_flow.numpy(),
          flow_gt.numpy(),
          np.ones_like(mask_thresh.numpy()),
          1. - mask_thresh.numpy(),
          1. - occ_mask_gt.numpy().astype('float32'),
          frame_skip=None)
  if progress_bar:
    sys.stdout.write('\n')
    sys.stdout.flush()

  fmax, best_thresh = get_fmax_and_best_thresh(all_occlusion_results)
  eval_stop_in_s = time.time()

  results = {
      'occl-f-max': fmax,
      'best-occl-thresh': best_thresh,
      'EPE': np.mean(np.array(epe_occ)),
      'ER': np.mean(np.array(errors_occ)),
      'inf-time(ms)': np.mean(inference_times),
      'eval-time(s)': eval_stop_in_s - eval_start_in_s,
  }
  for k in unsuper_losses[0].keys():
    results.update({k: np.mean([l[k] for l in unsuper_losses])})

  if prefix:
    return {prefix + '-' + k: v for k, v in results.items()}
  return results

def compute_f_metrics(mask_prediction, mask_gt, num_thresholds=40):
  """Return a dictionary of the true positives, etc. for two binary masks."""
  results = defaultdict(dict)
  mask_prediction = tf.cast(mask_prediction, tf.float32)
  mask_gt = tf.cast(mask_gt, tf.float32)
  for threshold in np.linspace(0, 1, num_thresholds):
    mask_thresh = tf.cast(
        tf.math.greater(mask_prediction, threshold), tf.float32)
    true_pos = tf.cast(tf.math.count_nonzero(mask_thresh * mask_gt), tf.float32)
    true_neg = tf.math.count_nonzero((mask_thresh - 1) * (mask_gt - 1))
    false_pos = tf.cast(
        tf.math.count_nonzero(mask_thresh * (mask_gt - 1)), tf.float32)
    false_neg = tf.cast(
        tf.math.count_nonzero((mask_thresh - 1) * mask_gt), tf.float32)
    results[threshold]['tp'] = true_pos
    results[threshold]['fp'] = false_pos
    results[threshold]['fn'] = false_neg
    results[threshold]['tn'] = true_neg
  return results

def get_fmax_and_best_thresh(results):
  """Select which threshold produces the best f1 score."""
  fmax = -1.
  best_thresh = -1.
  for thresh, metrics in results.items():
    precision = metrics['tp'] / (metrics['tp'] + metrics['fp'] + 1e-6)
    recall = metrics['tp'] / (metrics['tp'] + metrics['fn'] + 1e-6)
    f1 = 2 * precision * recall / (precision + recall + 1e-6)
    if f1 > fmax:
      fmax = f1
      best_thresh = thresh
  return fmax, best_thresh

def list_eval_keys(prefix=''):
  """List the keys of the dictionary returned by the evaluate function."""
  keys = [
      'EPE', 'ER', 'inf-time(ms)', 'eval-time(s)', 'occl-f-max',
      'best-occl-thresh'
  ]
  if prefix:
    return [prefix + '-' + k for k in keys]
  return keys

def epe_elementwise(estimation, groundtruth):
  """Computes the endpoint-error of each flow vector."""
  return tf.reduce_sum(
      (estimation - groundtruth)**2, axis=-1, keepdims=True)**0.5

def length_elementwise(flow):
  """Computes the length of each flow vector."""
  return tf.reduce_sum(flow**2, axis=-1, keepdims=True)**0.5

def outliers_elementwise(estimation, groundtruth, epe_threshold=3.0):
  """Computes the outlier criteria for the error rate per flow vector."""
  epe = epe_elementwise(estimation, groundtruth)
  length = length_elementwise(groundtruth)
  return tf.cast(
      tf.logical_and(epe > epe_threshold, epe > 0.05 * length), 'float32')

def angular_error_elementwise(estimation, groundtruth):
  """Computes the anuglar-error of each flow vector."""
  h, w, _ = tf.unstack(tf.shape(estimation))
  time_dim = tf.ones([h, w, 1])

  e_t = tf.concat([estimation, time_dim], -1)
  g_t = tf.concat([groundtruth, time_dim], -1)
  e_t = e_t / tf.reduce_sum(e_t**2, axis=-1, keepdims=True)**0.5
  g_t = g_t / tf.reduce_sum(g_t**2, axis=-1, keepdims=True)**0.5
  return tf.math.acos(
      tf.clip_by_value(tf.reduce_sum(e_t * g_t, -1, keepdims=True), -1., 1.))

ext: py

max_stars_repo_name: DionysisChristopoulos/google-research

max_issues_repo_name: DionysisChristopoulos/google-research

max_forks_repo_name: admariner/google-research

reason: import numpy


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

import os
import sys
import glob
import h5py
import numpy as np
from scipy.spatial.transform import Rotation
from torch.utils.data import Dataset

# Part of the code is referred from: https://github.com/charlesq34/pointnet

def download():
    BASE_DIR = os.path.dirname(os.path.abspath(__file__))
    DATA_DIR = os.path.join(BASE_DIR, 'data')
    if not os.path.exists(DATA_DIR):
        os.mkdir(DATA_DIR)

    datasets = [
        ('modelnet40_ply_hdf5_2048', 'https://shapenet.cs.stanford.edu/media/modelnet40_ply_hdf5_2048.zip'),
        ('3dmatch', 'http://web.tecnico.ulisboa.pt/sergio.agostinho/share/just-a-spoonful/3dmatch.zip'),
    ]
    for folder, www in datasets:
        if not os.path.exists(os.path.join(DATA_DIR, folder)):
            zipfile = os.path.basename(www)
            os.system('wget %s --no-check-certificate; unzip %s' % (www, zipfile))
            os.system('mv %s %s' % (zipfile[:-4], DATA_DIR))
            os.system('rm %s' % (zipfile))
    return DATA_DIR

def load_data(partition, prefix=None):
    if prefix:
        DATA_DIR = prefix
    else:
        download()
        BASE_DIR = os.path.dirname(os.path.abspath(__file__))
        DATA_DIR = os.path.join(BASE_DIR, 'data', 'modelnet40_ply_hdf5_2048')

    all_data = []
    all_label = []
    for h5_name in sorted(glob.glob(os.path.join(DATA_DIR, 'ply_data_%s*.h5' % partition))):
        f = h5py.File(h5_name, mode="r", swmr=True)
        data = f['data'][:].astype('float32')
        label = f['label'][:].astype('int64')
        f.close()
        all_data.append(data)
        all_label.append(label)
    all_data = np.concatenate(all_data, axis=0)
    all_label = np.concatenate(all_label, axis=0)
    return all_data, all_label

def translate_pointcloud(pointcloud):
    xyz1 = np.random.uniform(low=2. / 3., high=3. / 2., size=[3])
    xyz2 = np.random.uniform(low=-0.2, high=0.2, size=[3])

    translated_pointcloud = np.add(np.multiply(pointcloud, xyz1), xyz2).astype('float32')
    return translated_pointcloud

def jitter_pointcloud(pointcloud, sigma=0.01, clip=0.05):
    N, C = pointcloud.shape
    pointcloud += np.clip(sigma * np.random.randn(N, C), -1 * clip, clip)
    return pointcloud

class ModelNet40(Dataset):
    def __init__(self, num_points, partition='train', gaussian_noise=False, unseen=False, factor=4, prefix=None):
        self.data, self.label = load_data(partition, prefix=prefix)
        self.num_points = num_points
        self.partition = partition
        self.gaussian_noise = gaussian_noise
        self.unseen = unseen
        self.label = self.label.squeeze()
        self.factor = factor
        if self.unseen:
            ######## simulate testing on first 20 categories while training on last 20 categories
            if self.partition == 'test':
                self.data = self.data[self.label>=20]
                self.label = self.label[self.label>=20]
            elif self.partition == 'train':
                self.data = self.data[self.label<20]
                self.label = self.label[self.label<20]

    def __getitem__(self, item):
        pointcloud = self.data[item][:self.num_points]
        if self.gaussian_noise:
            pointcloud = jitter_pointcloud(pointcloud)
        if self.partition != 'train':
            np.random.seed(item)
        anglex = np.random.uniform() * np.pi / self.factor
        angley = np.random.uniform() * np.pi / self.factor
        anglez = np.random.uniform() * np.pi / self.factor

        cosx = np.cos(anglex)
        cosy = np.cos(angley)
        cosz = np.cos(anglez)
        sinx = np.sin(anglex)
        siny = np.sin(angley)
        sinz = np.sin(anglez)
        Rx = np.array([[1, 0, 0],
                        [0, cosx, -sinx],
                        [0, sinx, cosx]])
        Ry = np.array([[cosy, 0, siny],
                        [0, 1, 0],
                        [-siny, 0, cosy]])
        Rz = np.array([[cosz, -sinz, 0],
                        [sinz, cosz, 0],
                        [0, 0, 1]])
        R_ab = Rx.dot(Ry).dot(Rz)
        R_ba = R_ab.T
        translation_ab = np.array([np.random.uniform(-0.5, 0.5), np.random.uniform(-0.5, 0.5),
                                   np.random.uniform(-0.5, 0.5)])
        translation_ba = -R_ba.dot(translation_ab)

        pointcloud1 = pointcloud.T

        rotation_ab = Rotation.from_euler('zyx', [anglez, angley, anglex])
        pointcloud2 = rotation_ab.apply(pointcloud1.T).T + np.expand_dims(translation_ab, axis=1)

        euler_ab = np.asarray([anglez, angley, anglex])
        euler_ba = -euler_ab[::-1]

        pointcloud1 = np.random.permutation(pointcloud1.T).T
        pointcloud2 = np.random.permutation(pointcloud2.T).T

        return pointcloud1.astype('float32'), pointcloud2.astype('float32'), R_ab.astype('float32'), \
               translation_ab.astype('float32'), R_ba.astype('float32'), translation_ba.astype('float32'), \
               euler_ab.astype('float32'), euler_ba.astype('float32')

    def __len__(self):
        return self.data.shape[0]

def generate_random_poses(N, factor):
    euler_ab = np.random.rand(N, 3) * np.pi / factor
    euler_ba = -euler_ab[:,::-1]

    rot = Rotation.from_euler("zyx", euler_ab)
    R_ab = rot.as_matrix()
    R_ba = R_ab.transpose(0, 2, 1)

    translation_ab = np.random.rand(N, 3) - 0.5
    translation_ba = - np.squeeze(translation_ab[:, None] @ R_ab, axis=1)
    return dict(
        R_ab=R_ab,
        R_ba=R_ba,
        translation_ab=translation_ab,
        translation_ba=translation_ba,
        euler_ab=euler_ab,
        euler_ba=euler_ba,
    )

class ThreeDMatch(Dataset):

    def __init__(self, prefix, partition, minimum_overlap=0.3, factor=4):
        super().__init__()

        self.overlap_options = set([0.3, 0.5, 0.7])

        if minimum_overlap is not None and minimum_overlap not in self.overlap_options:
            msg = f"Accepted minimum_overlap values are the following: {self.overlap_options}"
            raise ValueError(msg)

        # Check and download data if needed
        if not prefix:
            data_dir = download()
            self.prefix = os.path.join(data_dir, "3dmatch")
        else:
            self.prefix = prefix

        # use stage information to populate list of files belonging to the split
        scenes = self._parse_scenes(partition)

        # retrieve all valid point cloud pairs
        self.pairs = self._parse_sequences(scenes, minimum_overlap)

        # generate random poses
        self.poses = generate_random_poses(len(self.pairs), factor)

    def _parse_scenes(self, partition):
        file = os.path.join(self.prefix, "splits", f"{partition}_3dmatch.txt")
        with open(file) as f:
            scenes = [line[:-1] for line in f.readlines()]
        return scenes

    def _parse_sequences(self, scenes, minimum_overlap):

        path = os.path.join(self.prefix, "preprocessed")
        pattern = "@seq-[0-9][0-9].txt" if minimum_overlap is None else f"@seq-[0-9][0-9]-{minimum_overlap:0.2f}.txt"

        pairs = []
        for scene in scenes:
            # sorting to ensure reproduceability across OSes
            sequences = sorted(list(glob.glob(os.path.join(path, scene + pattern))))
            for seq in sequences:
                with open(seq) as f:
                    pairs += [tuple(line.split()) for line in f.readlines()]
        return pairs

    def __len__(self) -> int:
        return len(self.pairs)

    def __getitem__(self, index: int):

        # load both point clouds
        file0, file1 = self.pairs[index][:2]
        pointcloud0 = np.load(os.path.join(self.prefix, "preprocessed", file0))["pcd"]
        pointcloud1 = np.load(os.path.join(self.prefix, "preprocessed", file1))["pcd"]

        # Rescale both point clouds, to lie inside a unit sphere
        scale = np.max(np.linalg.norm(np.stack([pointcloud0, pointcloud1]), axis=-1))
        pointcloud0 /= scale
        pointcloud1 /= scale

        # pose data
        R_ab = self.poses["R_ab"][index]
        R_ba = self.poses["R_ba"][index]
        euler_ab = self.poses["euler_ab"][index]
        euler_ba = self.poses["euler_ba"][index]
        translation_ab = self.poses["translation_ab"][index]
        translation_ba = self.poses["translation_ba"][index]

        # Point cloud data in 3DMatch is perfectly superimposed
        # apply transformation to target point cloud
        pointcloud1 = pointcloud1 @ R_ba + translation_ab
        out = (
            pointcloud0.T.astype(np.float32),
            pointcloud1.T.astype(np.float32),
            R_ab.astype(np.float32),
            translation_ab.astype(np.float32),
            R_ba.astype(np.float32),
            translation_ba.astype(np.float32),
            euler_ab.astype(np.float32),
            euler_ba.astype(np.float32),
        )
        return out

if __name__ == '__main__':
    train = ModelNet40(1024)
    test = ModelNet40(1024, 'test')
    for data in train:
        print(len(data))
        break

max_stars_repo_name: SergioRAgostinho/dcp

max_issues_repo_name: SergioRAgostinho/dcp

max_forks_repo_name: SergioRAgostinho/dcp

reason: import numpy,from scipy


text:
import numpy as np
import pandas
import sklearn.metrics
from skweak import utils
from spacy.tokens import Span # type: ignore

def evaluate(docs, all_labels, target_sources):
    """Extracts the evaluation results for one or more sources, and add them to a pandas DataFrame."""

    if isinstance(target_sources, str):
        target_sources = [target_sources]

    records = []
    for source in target_sources:
        results = get_results(docs, all_labels, source)
        labels = set(results["label_weights"].keys())
        for name in sorted(labels) + ["micro", "weighted", "macro"]:
            if name in results:
                record = results[name]
                record["label"] = name
                record["model"] = source
                if name in labels:
                    record["proportion"] = results["label_weights"][name]
                records.append(record)

    df = pandas.DataFrame.from_records(records)
    df["proportion"] = df.proportion.apply(lambda x: "%.1f %%"%(x*100) if not np.isnan(x) else "")
    df["tok_cee"] = df.tok_cee.apply(lambda x: str(x) if not np.isnan(x) else "")
    df["tok_acc"] = df.tok_acc.apply(lambda x: str(x) if not np.isnan(x) else "")
    df["coverage"] = df.coverage.apply(lambda x: str(x) if not np.isnan(x) else "")
    df = df.set_index(["label", "proportion", "model"]).sort_index()
    df = df[["tok_precision", "tok_recall", "tok_f1", "tok_cee", "tok_acc", "coverage",
             "ent_precision", "ent_recall", "ent_f1"]]
    return df

def get_results(docs, all_labels, target_source, conf_threshold=0.5):
    """Computes the usual metrics (precision, recall, F1, cross-entropy) on the dataset, using the spacy entities
    in each document as gold standard, and the annotations of a given source as the predicted values"""

    all_numbers = compute_raw_numbers(docs, all_labels, target_source, conf_threshold)
    tok_tp, tok_fp, tok_fn, tok_logloss, tok_nb, tok_tp_tn, ent_tp, ent_fp, ent_fn, ent_support, tok_support = all_numbers

    # We then compute the metrics themselves
    results = {}
    for label in ent_support:
        ent_pred = ent_tp[label]+ent_fp[label] + 1E-10
        ent_true = ent_tp[label]+ent_fn[label] + 1E-10
        tok_pred = tok_tp[label]+tok_fp[label] + 1E-10
        tok_true = tok_tp[label]+tok_fn[label] + 1E-10
        results[label] = {}
        results[label]["ent_precision"] = round(ent_tp[label]/ent_pred, 3)
        results[label]["ent_recall"] = round(ent_tp[label]/ent_true, 3)
        results[label]["tok_precision"] = round(tok_tp[label]/tok_pred, 3)
        results[label]["tok_recall"] = round(tok_tp[label]/tok_true, 3)

        ent_f1_numerator = (results[label]["ent_precision"] * results[label]["ent_recall"])
        ent_f1_denominator = (results[label]["ent_precision"] +results[label]["ent_recall"]) + 1E-10
        results[label]["ent_f1"] = 2*round(ent_f1_numerator / ent_f1_denominator, 3)

        tok_f1_numerator = (results[label]["tok_precision"] * results[label]["tok_recall"])
        tok_f1_denominator = (results[label]["tok_precision"] +results[label]["tok_recall"]) + 1E-10
        results[label]["tok_f1"] = 2*round(tok_f1_numerator / tok_f1_denominator, 3)

    results["macro"] = {"ent_precision":round(np.mean([results[l]["ent_precision"] for l in ent_support]), 3),
                       "ent_recall":round(np.mean([results[l]["ent_recall"] for l in ent_support]), 3),
                       "tok_precision":round(np.mean([results[l]["tok_precision"] for l in ent_support]), 3),
                       "tok_recall":round(np.mean([results[l]["tok_recall"] for l in ent_support]), 3)}

    label_weights = {l:ent_support[l]/sum(ent_support.values()) for l in ent_support}
    results["label_weights"] = label_weights
    results["weighted"] = {"ent_precision":round(np.sum([results[l]["ent_precision"]*label_weights[l]
                                                            for l in ent_support]), 3),
                           "ent_recall":round(np.sum([results[l]["ent_recall"]*label_weights[l]
                                                         for l in ent_support]), 3),
                           "tok_precision":round(np.sum([results[l]["tok_precision"]*label_weights[l]
                                                           for l in ent_support]), 3),
                           "tok_recall":round(np.sum([results[l]["tok_recall"]*label_weights[l]
                                                        for l in ent_support]), 3)}

    ent_pred = sum([ent_tp[l] for l in ent_support]) + sum([ent_fp[l] for l in ent_support]) + 1E-10
    ent_true = sum([ent_tp[l] for l in ent_support]) + sum([ent_fn[l] for l in ent_support]) + 1E-10
    tok_pred = sum([tok_tp[l] for l in ent_support]) + sum([tok_fp[l] for l in ent_support]) + 1E-10
    tok_true = sum([tok_tp[l] for l in ent_support])  + sum([tok_fn[l] for l in ent_support]) + 1E-10
    results["micro"] = {"ent_precision":round(sum([ent_tp[l] for l in ent_support]) / ent_pred, 3),
                        "ent_recall":round(sum([ent_tp[l] for l in ent_support]) / ent_true, 3),
                        "tok_precision":round(sum([tok_tp[l] for l in ent_support]) /tok_pred, 3),
                        "tok_recall":round(sum([tok_tp[l] for l in ent_support]) / tok_true, 3),
                        "tok_cee":round(tok_logloss/tok_nb, 3),
                        "tok_acc": round(tok_tp_tn/tok_nb, 3),
                        "coverage":round((sum(tok_tp.values()) +sum(tok_fp.values())) / sum(tok_support.values()), 3)}

    for metric in ["macro", "weighted", "micro"]:
        ent_f1_numerator = (results[metric]["ent_precision"] * results[metric]["ent_recall"])
        ent_f1_denominator = (results[metric]["ent_precision"] +results[metric]["ent_recall"]) + 1E-10
        results[metric]["ent_f1"] = 2*round(ent_f1_numerator / ent_f1_denominator, 3)

        tok_f1_numerator = (results[metric]["tok_precision"] * results[metric]["tok_recall"])
        tok_f1_denominator = (results[metric]["tok_precision"] +results[metric]["tok_recall"]) + 1E-10
        results[metric]["tok_f1"] = 2*round(tok_f1_numerator / tok_f1_denominator, 3)

    return results

def compute_raw_numbers(docs, all_labels, target_source, conf_threshold=0.5):
    """Computes the raw metrics (true positives, true negatives, ...) on the dataset, using the spacy entities
    in each document as gold standard, and the annotations of a given source as the predicted values"""

    # We start by computing the TP, FP and FN values
    tok_tp = {}
    tok_fp = {}
    tok_fn ={}

    tok_logloss = 0
    tok_nb = 0
    tok_tp_tn = 0

    ent_tp ={}
    ent_fp = {}
    ent_fn = {}
    ent_support = {}
    tok_support = {}

    for doc in docs:
        if target_source in doc.spans:
            spans = utils.get_spans_with_probs(doc, target_source)
        else:
            spans = []
        spans = [span for (span, prob) in spans if prob >= conf_threshold]

        for label in all_labels:
            true_spans = {(ent.start, ent.end) for ent in doc.ents if ent.label_==label}
            pred_spans = {(span.start,span.end) for span in spans if span.label_==label}

            ent_tp[label] = ent_tp.get(label,0) + len(true_spans.intersection(pred_spans))
            ent_fp[label] = ent_fp.get(label,0) + len(pred_spans - true_spans)
            ent_fn[label] = ent_fn.get(label,0) +  len(true_spans - pred_spans)
            ent_support[label] = ent_support.get(label, 0) + len(true_spans)

            true_tok_labels = {i for start,end in true_spans for i in range(start, end)}
            pred_tok_labels = {i for start,end in pred_spans for i in range(start, end)}
            tok_tp[label] = tok_tp.get(label, 0) + len(true_tok_labels.intersection(pred_tok_labels))
            tok_fp[label] = tok_fp.get(label, 0) + len(pred_tok_labels - true_tok_labels)
            tok_fn[label] = tok_fn.get(label,0) + len(true_tok_labels - pred_tok_labels)
            tok_support[label] = tok_support.get(label, 0) + len(true_tok_labels)

        gold_probs, pred_probs = _get_probs(doc, all_labels, target_source)
        tok_logloss += sklearn.metrics.log_loss(gold_probs, pred_probs, normalize=False)
        tok_tp_tn += sum(gold_probs.argmax(axis=1) == pred_probs.argmax(axis=1))
        tok_nb += len(doc)

    return (tok_tp, tok_fp, tok_fn, tok_logloss, tok_nb, tok_tp_tn, ent_tp,
            ent_fp, ent_fn, ent_support, tok_support)

def _get_probs(doc, all_labels, target_source):
    """Retrieves the gold and predicted probabilities (as matrices)"""

    out_label_indices = {"O":0}
    for label in all_labels:
        for prefix in "BI":
            out_label_indices["%s-%s" % (prefix, label)] = len(out_label_indices)

    gold_probs = np.zeros((len(doc), len(out_label_indices)), dtype=np.int16)
    for ent in doc.ents:
        gold_probs[ent.start, out_label_indices.get("B-%s" % ent.label_, 0)] = 1
        for i in range(ent.start+1, ent.end):
            gold_probs[i, out_label_indices.get("I-%s" % ent.label_, 0)] = 1

    pred_probs = np.zeros(gold_probs.shape)
    if target_source in doc.spans and "probs" in doc.spans[target_source].attrs:
        for tok_pos, labels in doc.spans[target_source].attrs["probs"].items():
            for label, label_prob in labels.items():
                pred_probs[tok_pos, out_label_indices[label]] = label_prob
        pred_probs[:,0] = np.clip(1-pred_probs[:,1:].sum(axis=1), 0.0, 1.0)
    else:
        vector = utils.spans_to_array(doc, all_labels, [target_source])[:,0]
        pred_probs[np.arange(vector.size), vector] = True

    return gold_probs, pred_probs

def show_errors(docs, all_labels, target_source, conf_threshold=0.5):
    """Utilities to display the errors/omissions of a given source"""

    for i, doc in enumerate(docs):

        spans = utils.get_spans_with_probs(doc, target_source)

        print("Doc %i:"%i, doc)
        true_spans = {(ent.start, ent.end):ent.label_ for ent in doc.ents}
        pred_spans = {(span.start,span.end):span.label_ for span, prob in spans if prob >=conf_threshold}

        for start,end in true_spans:
            if (start,end) not in pred_spans:
                print("Not found: %s [%i:%i] -> %s"%(doc[start:end], start, end, true_spans[(start,end)]))
            elif true_spans[(start,end)]!=pred_spans[(start,end)]:
                print("Wrong label: %s [%i:%i] -> %s but predicted as %s"%(doc[start:end], start, end,
                true_spans[(start,end)], pred_spans[(start,end)]))

        for start,end in pred_spans:
            if (start,end) not in true_spans:
                print("Spurious: %s [%i:%i] -> %s"%(doc[start:end], start, end, pred_spans[(start,end)]))


text:
#!/usr/bin/env python
# ------------------------------------------------------------------------------------------------------%
# Created by "Thieu Nguyen" at 09:48, 16/03/2020                                                        %
#                                                                                %
#       Email:      nguyenthieu2102@gmail.com                                                           %
#       Homepage:   https://www.researchgate.net/profile/Thieu_Nguyen6                                  %
#       Github:     https://github.com/thieu1995                                                        %
#-------------------------------------------------------------------------------------------------------%

import concurrent.futures as parallel
from functools import partial
import numpy as np
from mealpy.optimizer import Optimizer
from scipy.stats import cauchy
"""
BaseDE: - the very first DE algorithm (Novel mutation strategy for enhancing SHADE and LSHADE algorithms for global numerical optimization)
    strategy = 0: DE/current-to-rand/1/bin
             = 1: DE/best/1/bin
             = 2: DE/best/2/bin
             = 3: DE/rand/2/bin
             = 4: DE/current-to-best/1/bin
             = 5: DE/current-to-rand/1/bin
"""

class BaseDE(Optimizer):
    """
        The original version of: Differential Evolution (DE)
    """

    def __init__(self, problem, epoch=10000, pop_size=100, wf=0.8, cr=0.9, strategy=0, **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
            wf (float): weighting factor, default = 0.8
            cr (float): crossover rate, default = 0.9
            strategy (int): There are lots of variant version of DE algorithm, default is DE/current-to-rand/1/bin
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size
        self.weighting_factor = wf
        self.crossover_rate = cr
        self.strategy = strategy

    def _mutation__(self, current_pos, new_pos):
        pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < self.crossover_rate, current_pos, new_pos)
        return self.amend_position_faster(pos_new)

    def create_child(self, idx, pop_copy, g_best):
        if self.strategy == 0:
            # Choose 3 random element and different to i
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 3, replace=False)
            pos_new = pop_copy[idx_list[0]][self.ID_POS] + self.weighting_factor * \
                      (pop_copy[idx_list[1]][self.ID_POS] - pop_copy[idx_list[2]][self.ID_POS])
        elif self.strategy == 1:
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 2, replace=False)
            pos_new = g_best[self.ID_POS] + self.weighting_factor * (pop_copy[idx_list[0]][self.ID_POS] - pop_copy[idx_list[1]][self.ID_POS])
        elif self.strategy == 2:
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 4, replace=False)
            pos_new = g_best[self.ID_POS] + self.weighting_factor * (pop_copy[idx_list[0]][self.ID_POS] - pop_copy[idx_list[1]][self.ID_POS]) + \
                      self.weighting_factor * (pop_copy[idx_list[2]][self.ID_POS] - pop_copy[idx_list[3]][self.ID_POS])
        elif self.strategy == 3:
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 5, replace=False)
            pos_new = pop_copy[idx_list[0]][self.ID_POS] + self.weighting_factor * \
                      (pop_copy[idx_list[1]][self.ID_POS] - pop_copy[idx_list[2]][self.ID_POS]) + \
                      self.weighting_factor * (pop_copy[idx_list[3]][self.ID_POS] - pop_copy[idx_list[4]][self.ID_POS])
        elif self.strategy == 4:
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 2, replace=False)
            pos_new = pop_copy[idx][self.ID_POS] + self.weighting_factor * (g_best[self.ID_POS] - pop_copy[idx][self.ID_POS]) + \
                      self.weighting_factor * (pop_copy[idx_list[0]][self.ID_POS] - pop_copy[idx_list[1]][self.ID_POS])
        else:
            idx_list = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 3, replace=False)
            pos_new = pop_copy[idx][self.ID_POS] + self.weighting_factor * (pop_copy[idx_list[0]][self.ID_POS] - pop_copy[idx][self.ID_POS]) + \
                      self.weighting_factor * (pop_copy[idx_list[1]][self.ID_POS] - pop_copy[idx_list[2]][self.ID_POS])
        pos_new = self._mutation__(pop_copy[idx][self.ID_POS], pos_new)
        fit_new = self.get_fitness_position(pos_new)
        return [pos_new, fit_new]

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()
        pop_idx = np.array(range(0, self.pop_size))

        ## Reproduction
        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, g_best=g_best), pop_idx)
            child = [x for x in pop_child]
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, g_best=g_best), pop_idx)
            child = [x for x in pop_child]
        else:
            child = [self.create_child(idx, pop_copy, g_best=g_best) for idx in pop_idx]

        # create new pop by comparing fitness of corresponding each member in pop and children
        pop = self.greedy_selection_population(pop, child)
        return pop

class JADE(Optimizer):
    """
        The original version of: Differential Evolution (JADE)
        Link:
            JADE: Adaptive Differential Evolution with Optional External Archive
    """

    def __init__(self, problem, epoch=10000, pop_size=100, miu_f=0.5, miu_cr=0.5, pt=0.1, ap=0.1, **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
            miu_f (float): cross-over probability, default = 0.5
            miu_cr (float): mutation probability, default = 0.5
            pt (float): The percent of top best agents (p in the paper), default = 0.1
            ap (float): The Adaptation Parameter control value of f and cr (c in the paper), default=0.1
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size
        self.miu_f = miu_f  # the initial f, location is changed then that f is good
        self.miu_cr = miu_cr  # the initial cr,
        self.pt = pt  # np.random.uniform(0.05, 0.2) # the x_best is select from the top 100p % solutions
        self.ap = ap  # np.random.uniform(1/20, 1/5) # the adaptation parameter control value of f and cr

        ## Dynamic variable, changing in run time
        self.dyn_miu_cr = self.miu_cr
        self.dyn_miu_f = self.miu_f
        self.dyn_pop_archive = list()

    ### Survivor Selection
    def lehmer_mean(self, list_objects):
        return sum(list_objects**2) / sum(list_objects)

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()
        pop_idx = np.array(range(0, self.pop_size))

        list_f = list()
        list_cr = list()
        pop_sorted = sorted(pop, key=lambda x: x[self.ID_FIT][self.ID_TAR])

        ## Reproduction
        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted, list_f=list_f, list_cr=list_cr), pop_idx)
            child = [x for x in pop_child]
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted, list_f=list_f, list_cr=list_cr), pop_idx)
            child = [x for x in pop_child]
        else:
            child = [self.create_child(idx, pop_copy, pop_sorted, list_f, list_cr) for idx in pop_idx]
        pop = child.copy()
        # Randomly remove solution
        temp = len(self.dyn_pop_archive) - self.pop_size
        if temp > 0:
            idx_list = np.random.choice(range(0, len(self.dyn_pop_archive)), temp, replace=False)
            archive_pop_new = []
            for idx, solution in enumerate(self.dyn_pop_archive):
                if idx not in idx_list:
                    archive_pop_new.append(solution)
            self.dyn_pop_archive = archive_pop_new.copy()

        # Update miu_cr and miu_f
        self.dyn_miu_cr = (1 - self.ap) * self.dyn_miu_cr + self.ap * np.mean(np.array(list_cr))
        self.dyn_miu_f = (1 - self.ap) * self.dyn_miu_f + self.ap * self.lehmer_mean(np.array(list_f))
        return pop

    def create_child(self, idx, pop_copy, pop_sorted, list_f, list_cr):
        ## Calculate adaptive parameter cr and f
        cr = np.random.normal(self.dyn_miu_cr, 0.1)
        cr = np.clip(cr, 0, 1)
        while True:
            f = cauchy.rvs(self.dyn_miu_f, 0.1)
            if f < 0:
                continue
            elif f > 1:
                f = 1
            break
        top = int(self.pop_size * self.pt)
        x_best = pop_sorted[np.random.randint(0, top)]
        x_r1 = pop_copy[np.random.choice(list(set(range(0, self.pop_size)) - {idx}))]
        new_pop = pop_copy + self.dyn_pop_archive
        while True:
            x_r2 = new_pop[np.random.randint(0, len(new_pop))]
            if np.any(x_r2[self.ID_POS] - x_r1[self.ID_POS]) and np.any(x_r2[self.ID_POS] - pop_copy[idx][self.ID_POS]):
                break
        x_new = pop_copy[idx][self.ID_POS] + f * (x_best[self.ID_POS] - pop_copy[idx][self.ID_POS]) + f * (x_r1[self.ID_POS] - x_r2[self.ID_POS])
        pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < cr, x_new, pop_copy[idx][self.ID_POS])
        j_rand = np.random.randint(0, self.problem.n_dims)
        pos_new[j_rand] = x_new[j_rand]
        pos_new = self.amend_position_faster(pos_new)
        fit_new = self.get_fitness_position(pos_new)

        if self.compare_agent([pos_new, fit_new], pop_copy[idx]):
            self.dyn_pop_archive.append(pop_copy[idx].copy())
            list_cr.append(cr)
            list_f.append(f)
            return [pos_new, fit_new]
        return pop_copy[idx]

class SADE(Optimizer):
    """
        The original version of: Self-Adaptive Differential Evolution(SADE)
        Link:
            Self-adaptive differential evolution algorithm for numerical optimization
    """

    def __init__(self, problem, epoch=10000, pop_size=100, **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size

        self.loop_probability = 50
        self.loop_cr = 5
        self.ns1 = self.ns2 = self.nf1 = self.nf2 = 0
        self.crm = 0.5
        self.p1 = 0.5

        # Dynamic variable
        self.dyn_list_cr = list()

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()
        pop_idx = np.array(range(0, self.pop_size))

        ## Reproduction
        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, g_best=g_best), pop_idx)
            child = [x for x in pop_child]
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, g_best=g_best), pop_idx)
            child = [x for x in pop_child]
        else:
            child = [self.create_child(idx, pop_copy, g_best) for idx in pop_idx]

        # Update cr and p1
        if (epoch + 1) / self.loop_cr == 0:
            self.crm = np.mean(self.dyn_list_cr)
            self.dyn_list_cr = list()

        if (epoch + 1) / self.loop_probability == 0:
            self.p1 = self.ns1 * (self.ns2 + self.nf2) / (self.ns2 * (self.ns1 + self.nf1) + self.ns1 * (self.ns2 + self.nf2))
            self.ns1 = self.ns2 = self.nf1 = self.nf2 = 0

        return child.copy()

    def create_child(self, idx, pop_copy, g_best):
        ## Calculate adaptive parameter cr and f
        cr = np.random.normal(self.crm, 0.1)
        cr = np.clip(cr, 0, 1)
        while True:
            f = np.random.normal(0.5, 0.3)
            if f < 0:
                continue
            elif f > 1:
                f = 1
            break

        id1, id2, id3 = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 3, replace=False)
        if np.random.rand() < self.p1:
            x_new = pop_copy[id1][self.ID_POS] + f * (pop_copy[id2][self.ID_POS] - pop_copy[id3][self.ID_POS])
            pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < cr, x_new, pop_copy[idx][self.ID_POS])
            j_rand = np.random.randint(0, self.problem.n_dims)
            pos_new[j_rand] = x_new[j_rand]
            pos_new = self.amend_position_faster(pos_new)
            fit_new = self.get_fitness_position(pos_new)
            if self.compare_agent([pos_new, fit_new], pop_copy[idx]):
                self.ns1 += 1
                return [pos_new, fit_new]
            else:
                self.nf1 += 1
                return pop_copy[idx]
        else:
            x_new = pop_copy[idx][self.ID_POS] + f * (g_best[self.ID_POS] - pop_copy[idx][self.ID_POS]) + \
                    f * (pop_copy[id1][self.ID_POS] - pop_copy[id2][self.ID_POS])
            pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < cr, x_new, pop_copy[idx][self.ID_POS])
            j_rand = np.random.randint(0, self.problem.n_dims)
            pos_new[j_rand] = x_new[j_rand]
            pos_new = self.amend_position_faster(pos_new)
            fit_new = self.get_fitness_position(pos_new)
            if self.compare_agent([pos_new, fit_new], pop_copy[idx]):
                self.ns2 += 1
                self.dyn_list_cr.append(cr)
                return [pos_new, fit_new]
            else:
                self.nf2 += 1
                return pop_copy[idx]

class SHADE(Optimizer):
    """
        The original version of: Success-History Adaptation Differential Evolution (SHADE)
        Link:
            Success-History Based Parameter Adaptation for Differential Evolution
    """

    def __init__(self, problem, epoch=750, pop_size=100, miu_f=0.5, miu_cr=0.5, **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
            miu_f (float): cross-over probability, default = 0.5
            miu_cr (float): mutation probability, default = 0.5
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size

        # Dynamic variable
        self.dyn_miu_f = miu_f * np.ones(self.pop_size)     # list the initial f,
        self.dyn_miu_cr = miu_cr * np.ones(self.pop_size)   # list the initial cr,
        self.dyn_pop_archive = list()
        self.k_counter = 0

    ### Survivor Selection
    def weighted_lehmer_mean(self, list_objects, list_weights):
        up = list_weights * list_objects ** 2
        down = list_weights * list_objects
        return sum(up) / sum(down)

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()  # Save all old elements --> Use to update cr value
        pop_idx = np.array(range(0, self.pop_size))

        list_f = list()
        list_cr = list()
        list_f_index = list()
        list_cr_index = list()

        list_f_new = np.ones(self.pop_size)
        list_cr_new = np.ones(self.pop_size)
        pop_old = pop.copy()
        pop_sorted = sorted(pop, key=lambda x: x[self.ID_FIT][self.ID_TAR])

        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted,
                                                 list_cr_new=list_cr_new, list_f_new=list_f_new), pop_idx)
            pop_new = [x for x in pop_child]  # Save all new elements --> Use to update the list_cr and list_f
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted,
                                                 list_cr_new=list_cr_new, list_f_new=list_f_new), pop_idx)
            pop_new = [x for x in pop_child]
        else:
            pop_new = [self.create_child(idx, pop_copy, pop_sorted, list_cr_new, list_f_new) for idx in pop_idx]
        for i in range(0, self.pop_size):
            if self.compare_agent(pop_new[i], pop[i]):
                list_cr.append(list_cr_new[i])
                list_f.append(list_f_new[i])
                list_f_index.append(i)
                list_cr_index.append(i)
                pop[i] = pop_new[i].copy()
                self.dyn_pop_archive.append(pop[i].copy())

        # Randomly remove solution
        temp = len(self.dyn_pop_archive) - self.pop_size
        if temp > 0:
            idx_list = np.random.choice(range(0, len(self.dyn_pop_archive)), temp, replace=False)
            archive_pop_new = []
            for idx, solution in enumerate(self.dyn_pop_archive):
                if idx not in idx_list:
                    archive_pop_new.append(solution)
            self.dyn_pop_archive = archive_pop_new.copy()

        # Update miu_cr and miu_f
        if len(list_f) != 0 and len(list_cr) != 0:
            # Eq.13, 14, 10
            list_fit_old = np.ones(len(list_cr_index))
            list_fit_new = np.ones(len(list_cr_index))
            idx_increase = 0
            for i in range(0, self.pop_size):
                if i in list_cr_index:
                    list_fit_old[idx_increase] = pop_old[i][self.ID_FIT][self.ID_TAR]
                    list_fit_new[idx_increase] = pop_new[i][self.ID_FIT][self.ID_TAR]
                    idx_increase += 1
            list_weights = abs(list_fit_new - list_fit_old) / sum(abs(list_fit_new - list_fit_old))
            self.dyn_miu_cr[self.k_counter] = sum(list_weights * np.array(list_cr))
            self.dyn_miu_f[self.k_counter] = self.weighted_lehmer_mean(np.array(list_f), list_weights)
            self.k_counter += 1
            if self.k_counter >= self.pop_size:
                self.k_counter = 0
        return pop

    def create_child(self, idx, pop_copy, pop_sorted, list_cr_new, list_f_new):
        ## Calculate adaptive parameter cr and f
        idx_rand = np.random.randint(0, self.pop_size)
        cr = np.random.normal(self.dyn_miu_cr[idx_rand], 0.1)
        cr = np.clip(cr, 0, 1)
        while True:
            f = cauchy.rvs(self.dyn_miu_f[idx_rand], 0.1)
            if f < 0:
                continue
            elif f > 1:
                f = 1
            break
        list_cr_new[idx] = cr
        list_f_new[idx] = f
        p = np.random.uniform(2 / self.pop_size, 0.2)
        top = int(self.pop_size * p)
        x_best = pop_sorted[np.random.randint(0, top)]
        x_r1 = pop_copy[np.random.choice(list(set(range(0, self.pop_size)) - {idx}))]
        new_pop = pop_copy + self.dyn_pop_archive
        while True:
            x_r2 = new_pop[np.random.randint(0, len(new_pop))]
            if np.any(x_r2[self.ID_POS] - x_r1[self.ID_POS]) and np.any(x_r2[self.ID_POS] - pop_copy[idx][self.ID_POS]):
                break
        x_new = pop_copy[idx][self.ID_POS] + f * (x_best[self.ID_POS] - pop_copy[idx][self.ID_POS]) + f * (x_r1[self.ID_POS] - x_r2[self.ID_POS])
        pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < cr, x_new, pop_copy[idx][self.ID_POS])
        j_rand = np.random.randint(0, self.problem.n_dims)
        pos_new[j_rand] = x_new[j_rand]
        pos_new = self.amend_position_faster(pos_new)
        fit_new = self.get_fitness_position(pos_new)
        return [pos_new, fit_new]

class L_SHADE(Optimizer):
    """
        The original version of: Linear Population Size Reduction Success-History Adaptation Differential Evolution (LSHADE)
        Link:
            Improving the Search Performance of SHADE Using Linear Population Size Reduction
    """

    def __init__(self, problem, epoch=750, pop_size=100, miu_f=0.5, miu_cr=0.5, **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
            miu_f (float): cross-over probability, default = 0.5
            miu_cr (float): mutation probability, default = 0.5
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size

        # Dynamic variable
        self.dyn_miu_f = miu_f * np.ones(self.pop_size)  # list the initial f,
        self.dyn_miu_cr = miu_cr * np.ones(self.pop_size)  # list the initial cr,
        self.dyn_pop_archive = list()
        self.dyn_pop_size = self.pop_size
        self.k_counter = 0
        self.n_min = int(self.pop_size / 5)

    ### Survivor Selection
    def weighted_lehmer_mean(self, list_objects, list_weights):
        up = list_weights * list_objects ** 2
        down = list_weights * list_objects
        return sum(up) / sum(down)

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()  # Save all old elements --> Use to update cr value
        pop_idx = np.array(range(0, self.pop_size))

        list_f = list()
        list_cr = list()
        list_f_index = list()
        list_cr_index = list()

        list_f_new = np.ones(self.pop_size)
        list_cr_new = np.ones(self.pop_size)
        pop_old = pop.copy()
        pop_sorted = sorted(pop, key=lambda x: x[self.ID_FIT][self.ID_TAR])

        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted,
                                                 list_cr_new=list_cr_new, list_f_new=list_f_new), pop_idx)
            pop_new = [x for x in pop_child]  # Save all new elements --> Use to update the list_cr and list_f
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy, pop_sorted=pop_sorted,
                                                 list_cr_new=list_cr_new, list_f_new=list_f_new), pop_idx)
            pop_new = [x for x in pop_child]
        else:
            pop_new = [self.create_child(idx, pop_copy, pop_sorted, list_cr_new, list_f_new) for idx in pop_idx]

        for i in range(0, self.pop_size):
            if self.compare_agent(pop_new[i], pop[i]):
                list_cr.append(list_cr_new[i])
                list_f.append(list_f_new[i])
                list_f_index.append(i)
                list_cr_index.append(i)
                pop[i] = pop_new[i].copy()
                self.dyn_pop_archive.append(pop[i].copy())

        # Randomly remove solution
        temp = len(self.dyn_pop_archive) - self.pop_size
        if temp > 0:
            idx_list = np.random.choice(range(0, len(self.dyn_pop_archive)), temp, replace=False)
            archive_pop_new = []
            for idx, solution in enumerate(self.dyn_pop_archive):
                if idx not in idx_list:
                    archive_pop_new.append(solution)
            self.dyn_pop_archive = archive_pop_new.copy()

        # Update miu_cr and miu_f
        if len(list_f) != 0 and len(list_cr) != 0:
            # Eq.13, 14, 10
            list_fit_old = np.ones(len(list_cr_index))
            list_fit_new = np.ones(len(list_cr_index))
            idx_increase = 0
            for i in range(0, self.dyn_pop_size):
                if i in list_cr_index:
                    list_fit_old[idx_increase] = pop_old[i][self.ID_FIT][self.ID_TAR]
                    list_fit_new[idx_increase] = pop_new[i][self.ID_FIT][self.ID_TAR]
                    idx_increase += 1
            list_weights = np.abs(list_fit_new - list_fit_old) / sum(np.abs(list_fit_new - list_fit_old))
            self.dyn_miu_cr[self.k_counter] = sum(list_weights * np.array(list_cr))
            self.dyn_miu_f[self.k_counter] = self.weighted_lehmer_mean(np.array(list_f), list_weights)
            self.k_counter += 1
            if self.k_counter >= self.dyn_pop_size:
                self.k_counter = 0

        # Linear Population Size Reduction
        self.dyn_pop_size = round(self.pop_size + epoch * ((self.n_min - self.pop_size) / self.epoch))
        return pop

    def create_child(self, idx, pop_copy, pop_sorted, list_cr_new, list_f_new):
        ## Calculate adaptive parameter cr and f
        idx_rand = np.random.randint(0, self.pop_size)
        cr = np.random.normal(self.dyn_miu_cr[idx_rand], 0.1)
        cr = np.clip(cr, 0, 1)
        while True:
            f = cauchy.rvs(self.dyn_miu_f[idx_rand], 0.1)
            if f < 0:
                continue
            elif f > 1:
                f = 1
            break
        list_cr_new[idx] = cr
        list_f_new[idx] = f
        p = np.random.uniform(0.15, 0.2)
        top = int(self.dyn_pop_size * p)
        x_best = pop_sorted[np.random.randint(0, top)]
        x_r1 = pop_copy[np.random.choice(list(set(range(0, self.dyn_pop_size)) - {idx}))]
        new_pop = pop_copy + self.dyn_pop_archive
        while True:
            x_r2 = new_pop[np.random.randint(0, len(new_pop))]
            if np.any(x_r2[self.ID_POS] - x_r1[self.ID_POS]) and np.any(x_r2[self.ID_POS] - pop_copy[idx][self.ID_POS]):
                break
        x_new = pop_copy[idx][self.ID_POS] + f * (x_best[self.ID_POS] - pop_copy[idx][self.ID_POS]) + f * (x_r1[self.ID_POS] - x_r2[self.ID_POS])
        pos_new = np.where(np.random.uniform(0, 1, self.problem.n_dims) < cr, x_new, pop_copy[idx][self.ID_POS])
        j_rand = np.random.randint(0, self.problem.n_dims)
        pos_new[j_rand] = x_new[j_rand]
        pos_new = self.amend_position_faster(pos_new)
        fit_new = self.get_fitness_position(pos_new)
        return [pos_new, fit_new]

class SAP_DE(Optimizer):
    """
        The original version of: Differential Evolution with Self-Adaptive Populations
        Link:
            Exploring dynamic self-adaptive populations in differential evolution
    """
    ID_CR = 2
    ID_MR = 3
    ID_PS = 4

    def __init__(self, problem, epoch=750, pop_size=100, wf=0.8, cr=0.9, F=1, branch="ABS", **kwargs):
        """
        Args:
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
        """
        super().__init__(problem, kwargs)
        self.nfe_per_epoch = pop_size
        self.sort_flag = False

        self.epoch = epoch
        self.pop_size = pop_size
        self.weighting_factor = wf
        self.crossover_rate = cr
        self.F = F
        self.M = pop_size
        self.branch = branch            # np.absolute (ABS) or relative (REL)

    def create_solution(self):
        """
        Returns:
            The position position with 2 element: index of position/location and index of fitness wrapper
            The general format: [position, [target, [obj1, obj2, ...]], strategy, times_win]

        ## To get the position, fitness wrapper, target and obj list
        ##      A[self.ID_POS]                  --> Return: position
        ##      A[self.ID_FIT]                  --> Return: [target, [obj1, obj2, ...]]
        ##      A[self.ID_FIT][self.ID_TAR]     --> Return: target
        ##      A[self.ID_FIT][self.ID_OBJ]     --> Return: [obj1, obj2, ...]
        """
        position = np.random.uniform(self.problem.lb, self.problem.ub)
        fitness = self.get_fitness_position(position=position)
        crossover_rate = np.random.uniform(0, 1)
        mutation_rate = np.random.uniform(0, 1)
        if self.branch == "ABS":
            pop_size = int(10 * self.problem.n_dims + np.random.normal(0, 1))
        else:    #elif self.branch == "REL":
            pop_size = int(10 * self.problem.n_dims + np.random.uniform(-0.5, 0.5))
        return [position, fitness, crossover_rate, mutation_rate, pop_size]

    def edit_to_range(self, var=None, lower=0, upper=1, func_value=None):
        while var <= lower or var >= upper:
            if var <= lower:
                var += func_value()
            if var >= upper:
                var -= func_value()
        return var

    def evolve(self, mode='sequential', epoch=None, pop=None, g_best=None):
        """
        Args:
            mode (str): 'sequential', 'thread', 'process'
                + 'sequential': recommended for simple and small task (< 10 seconds for calculating objective)
                + 'thread': recommended for IO bound task, or small computing task (< 2 minutes for calculating objective)
                + 'process': recommended for hard and big task (> 2 minutes for calculating objective)

        Returns:
            [position, fitness value]
        """
        pop_copy = pop.copy()  # Save all old elements --> Use to update cr value
        pop_idx = np.array(range(0, self.pop_size))

        if mode == "thread":
            with parallel.ThreadPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy), pop_idx)
            pop_new = [x for x in pop_child]  # Save all new elements --> Use to update the list_cr and list_f
        elif mode == "process":
            with parallel.ProcessPoolExecutor() as executor:
                pop_child = executor.map(partial(self.create_child, pop_copy=pop_copy), pop_idx)
            pop_new = [x for x in pop_child]
        else:
            pop_new = [self.create_child(idx, pop_copy) for idx in pop_idx]
        pop = pop_new.copy()

        # Calculate new population size
        total = sum([pop[i][self.ID_PS] for i in range(0, self.pop_size)])
        if self.branch == "ABS":
            m_new = int(total / self.pop_size)
        else: # elif self.branch == "REL":
            m_new = int(self.pop_size + total)
        if m_new <= 4:
            m_new = self.M + int(np.random.uniform(0, 4))
        elif m_new > 4 * self.M:
            m_new = self.M - int(np.random.uniform(0, 4))

        ## Change population by population size
        if m_new <= self.pop_size:
            pop = pop[:m_new]
        else:
            pop_sorted = sorted(pop, key=lambda x: x[self.ID_FIT][self.ID_TAR])
            best = pop_sorted[0].copy()
            pop_best = [best for _ in range(0, m_new - self.pop_size)]
            pop = pop + pop_best
        self.pop_size = m_new
        return pop

    def create_child(self, idx, pop_copy):
        # Choose 3 random element and different to idx
        idxs = np.random.choice(list(set(range(0, self.pop_size)) - {idx}), 3, replace=False)
        j = np.random.randint(0, self.pop_size)
        self.F = np.random.uniform(0, 1)

        sol_new = pop_copy[idxs[0]].copy()
        ## Crossover
        if np.random.uniform(0, 1) < pop_copy[idx][self.ID_CR] or idx == j:
            pos_new = pop_copy[idxs[0]][self.ID_POS] + self.F * (pop_copy[idxs[1]][self.ID_POS] - pop_copy[idxs[2]][self.ID_POS])
            cr_new = pop_copy[idxs[0]][self.ID_CR] + self.F * (pop_copy[idxs[1]][self.ID_CR] - pop_copy[idxs[2]][self.ID_CR])
            mr_new = pop_copy[idxs[0]][self.ID_MR] + self.F * (pop_copy[idxs[1]][self.ID_MR] - pop_copy[idxs[2]][self.ID_MR])
            if self.branch == "ABS":
                ps_new = pop_copy[idxs[0]][self.ID_PS] + int(self.F * (pop_copy[idxs[1]][self.ID_PS] - pop_copy[idxs[2]][self.ID_PS]))
            else:       # elif self.branch == "REL":
                ps_new = pop_copy[idxs[0]][self.ID_PS] + self.F * (pop_copy[idxs[1]][self.ID_PS] - pop_copy[idxs[2]][self.ID_PS])
            pos_new = self.amend_position_faster(pos_new)
            fit_new = self.get_fitness_position(pos_new)
            cr_new = self.edit_to_range(cr_new, 0, 1, np.random.random)
            mr_new = self.edit_to_range(mr_new, 0, 1, np.random.random)
            sol_new = [pos_new, fit_new, cr_new, mr_new, ps_new]
        ## Mutation
        if np.random.uniform(0, 1) < pop_copy[idxs[0]][self.ID_MR]:
            pos_new = pop_copy[idx][self.ID_POS] + np.random.normal(0, pop_copy[idxs[0]][self.ID_MR])
            cr_new = np.random.normal(0, 1)
            mr_new = np.random.normal(0, 1)
            if self.branch == "ABS":
                ps_new = pop_copy[idx][self.ID_PS] + int(np.random.normal(0.5, 1))
            else:       # elif self.branch == "REL":
                ps_new = pop_copy[idx][self.ID_PS] + np.random.normal(0, pop_copy[idxs[0]][self.ID_MR])
            pos_new = self.amend_position_faster(pos_new)
            fit_new = self.get_fitness_position(pos_new)
            sol_new = [pos_new, fit_new, cr_new, mr_new, ps_new]
        return sol_new


text:
# Copyright 2021 DeepMind Technologies Limited
#
# 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.
# ==============================================================================

"""Utilities for parsing actions and observations produced by the environment."""

import collections
import enum
from typing import Optional, Sequence, Tuple

import numpy as np

# Indices of bits representing global state.
NUM_POWERS = 7
# Max number of units that might need an order in total, across all players.
MAX_ACTIONS = 34

# Bits representing interesting things in per-province state vectors
OBSERVATION_UNIT_ARMY = 0
OBSERVATION_UNIT_FLEET = 1
OBSERVATION_UNIT_ABSENT = 2
OBSERVATION_UNIT_POWER_START = 3
OBSERVATION_BUILDABLE = 11
OBSERVATION_REMOVABLE = 12
OBSERVATION_DISLODGED_ARMY = 13
OBSERVATION_DISLODGED_FLEET = 14
OBSERVATION_DISLODGED_START = 16
OBSERVATION_SC_POWER_START = 27
PROVINCE_VECTOR_LENGTH = 35

# Areas information. The observation has NUM_AREAS area vectors:
# - First, a vector for each of the SINGLE_COASTED_PROVINCES provinces.
# - Then, three vectors for each of the BICOASTAL_PROVINCES. The first one is
#   the land area, and the other two are the coasts.
SINGLE_COASTED_PROVINCES = 72
BICOASTAL_PROVINCES = 3
NUM_AREAS = SINGLE_COASTED_PROVINCES + 3 * BICOASTAL_PROVINCES
NUM_PROVINCES = SINGLE_COASTED_PROVINCES + BICOASTAL_PROVINCES
NUM_SUPPLY_CENTERS = 34

OBSERVATION_BOARD_SHAPE = [NUM_AREAS, PROVINCE_VECTOR_LENGTH]

# Agents may want a list of areas that they can take actions in. This
# doesn't make sense in the builds phase, so we use the value -2 to indicate
# that, and will use lists like [BUILD_PHASE_AREA]*num_builds as 'area lists'
BUILD_PHASE_AREA_FLAG = -2
INVALID_AREA_FLAG = -1

# Main areas_ids of bicoastal provinces
BICOASTAL_PROVINCES_MAIN_AREAS = [72, 75, 78]

ProvinceID = int  # in [0, 1, ...74]
AreaID = int  # in [0, 1, ...80]
# AreaIndex is 0 for the "land" area of bicoastal provinces and for any province
# that has either 0 or 1 coasts. AreaIndex is 1 for the first coast and 2 for
# the second coast of bicoastal provinces.
AreaIndex = int  # in [0, 1, 2]
# CoastFlag is 0 for the first and 1 for the second coast of bicoastal
# provinces. CoastFlag is 0 for all provinces that have 0 or 1 coast.
CoastFlag = int  # in [0, 1]
ProvinceWithFlag = Tuple[ProvinceID, CoastFlag]

class Season(enum.Enum):
  """Diplomacy season."""
  SPRING_MOVES = 0
  SPRING_RETREATS = 1
  AUTUMN_MOVES = 2
  AUTUMN_RETREATS = 3
  BUILDS = 4

  def is_moves(self):
    return self == Season.SPRING_MOVES or self == Season.AUTUMN_MOVES

  def is_retreats(self):
    return self == Season.SPRING_RETREATS or self == Season.AUTUMN_RETREATS

  def is_builds(self):
    return self == Season.BUILDS

NUM_SEASONS = len(Season)

class UnitType(enum.Enum):
  ARMY = 0
  FLEET = 1

Observation = collections.namedtuple(
    'Observation', ['season', 'board', 'build_numbers', 'last_actions'])

class ProvinceType(enum.Enum):
  LAND = 0
  SEA = 1
  COASTAL = 2
  BICOASTAL = 3

def province_type_from_id(province_id: ProvinceID) -> ProvinceType:
  """Returns the ProvinceType for the province."""
  if province_id < 14:
    return ProvinceType.LAND
  elif province_id < 33:
    return ProvinceType.SEA
  elif province_id < 72:
    return ProvinceType.COASTAL
  elif province_id < 75:
    return ProvinceType.BICOASTAL
  else:
    raise ValueError('Invalid ProvinceID (too large)')

def province_id_and_area_index(area: AreaID) -> Tuple[ProvinceID, AreaIndex]:
  """Returns the province_id and the area index within the province.

  Args:
    area: the ID of the area in the observation vector, an integer from 0 to 80

  Returns:
    province_id: This is between 0 and
      SINGLE_COASTED_PROVINCES + BICOASTAL_PROVINCES - 1, and corresponds to the
      representation of this area used in orders.
    area_index: this is 0 for the main area of a province, and 1 or 2 for a
      coast in a bicoastal province.
  """
  if area < SINGLE_COASTED_PROVINCES:
    return area, 0
  province_id = (
      SINGLE_COASTED_PROVINCES + (area - SINGLE_COASTED_PROVINCES) // 3)
  area_index = (area - SINGLE_COASTED_PROVINCES) % 3
  return province_id, area_index

_prov_and_area_id_to_area = {province_id_and_area_index(area): area
                             for area in range(NUM_AREAS)}

def area_from_province_id_and_area_index(province_id: ProvinceID,
                                         area_index: AreaIndex) -> AreaID:
  """The inverse of province_id_and_area_index.

  Args:
    province_id: This is between 0 and
      SINGLE_COASTED_PROVINCES + BICOASTAL_PROVINCES - 1, and corresponds to the
      representation of this area used in orders.
    area_index: this is 0 for the main area of a province, and 1 or 2 for a
      coast in a bicoastal province.

  Returns:
    area: the id of the area in the observation vector

  Raises:
    KeyError: If the province_id and area_index are invalid
  """
  return _prov_and_area_id_to_area[(province_id, area_index)]

def area_index_for_fleet(
    province_tuple: ProvinceWithFlag) -> AreaIndex:
  if province_type_from_id(province_tuple[0]) == ProvinceType.BICOASTAL:
    return province_tuple[1] + 1
  else:
    return 0

def obs_index_start_and_num_areas(
    province_id: ProvinceID) -> Tuple[AreaID, int]:
  """Returns the area_id of the province's main area, and the number of areas.

  Args:
    province_id: the id of the province.
  """
  if province_id < SINGLE_COASTED_PROVINCES:
    return province_id, 1
  area_start = (
      SINGLE_COASTED_PROVINCES + (province_id - SINGLE_COASTED_PROVINCES) * 3)
  return area_start, 3

def moves_phase_areas(country_index: int, board_state: np.ndarray,
                      retreats: bool) -> Sequence[AreaID]:
  """Returns the areas with country_index's units active for this phase."""
  offset = (OBSERVATION_DISLODGED_START if retreats else
            OBSERVATION_UNIT_POWER_START)
  our_areas = np.where(board_state[:, country_index + offset])[0]
  filtered_areas = []
  provinces = set()
  for area in our_areas:
    province_id, area_index = province_id_and_area_index(area)
    if retreats:
      u_type = dislodged_unit_type(province_id, board_state)
    else:
      u_type = unit_type(province_id, board_state)
    # This area is valid, unless the unit is a fleet, this is the first area of
    # its province, and the province is bicoastal.
    if (u_type == UnitType.FLEET and
        area_index == 0 and
        obs_index_start_and_num_areas(province_id)[1] > 1):
      continue
    filtered_areas.append(area)
    if province_id in provinces:
      raise ValueError('Duplicate province in move phase areas')
    provinces.add(province_id)
  return sorted(list(filtered_areas))

def order_relevant_areas(observation: Observation, player: int,
                         topological_index=None) -> Sequence[AreaID]:
  """Areas with moves sorted according to topological_index."""
  season = observation.season
  if season.is_moves():
    areas = moves_phase_areas(player, observation.board, False)
  elif season.is_retreats():
    areas = moves_phase_areas(player, observation.board, True)
  else:
    areas = [BUILD_PHASE_AREA_FLAG] * abs(
        observation.build_numbers[player])
    return areas

  provinces_to_areas = dict()
  for area in areas:
    province, _ = province_id_and_area_index(area)
    if (province not in provinces_to_areas or
        area > provinces_to_areas[province]):  # This selects coasts over land
      provinces_to_areas[province] = area

  areas_without_repeats = list(provinces_to_areas.values())
  if topological_index:
    areas_without_repeats.sort(key=topological_index.index)

  return areas_without_repeats

def unit_type(province_id: ProvinceID,
              board_state: np.ndarray) -> Optional[UnitType]:
  """Returns the unit type in the province."""
  main_area, _ = obs_index_start_and_num_areas(province_id)
  return unit_type_from_area(main_area, board_state)

def unit_type_from_area(area_id: AreaID,
                        board_state: np.ndarray) -> Optional[UnitType]:
  if board_state[area_id, OBSERVATION_UNIT_ARMY] > 0:
    return UnitType(UnitType.ARMY)
  elif board_state[area_id, OBSERVATION_UNIT_FLEET] > 0:
    return UnitType(UnitType.FLEET)
  return None

def dislodged_unit_type(province_id: ProvinceID,
                        board_state: np.ndarray) -> Optional[UnitType]:
  """Returns the type of any dislodged unit in the province."""
  main_area, _ = obs_index_start_and_num_areas(province_id)
  return dislodged_unit_type_from_area(main_area, board_state)

def dislodged_unit_type_from_area(
    area_id: AreaID, board_state: np.ndarray) -> Optional[UnitType]:
  """Returns the type of any dislodged unit in the province."""
  if board_state[area_id, OBSERVATION_DISLODGED_ARMY] > 0:
    return UnitType(UnitType.ARMY)
  elif board_state[area_id, OBSERVATION_DISLODGED_FLEET] > 0:
    return UnitType(UnitType.FLEET)
  return None

def unit_power(province_id: ProvinceID,
               board_state: np.ndarray) -> Optional[int]:
  """Returns which power controls the unit province (None if no unit there)."""
  main_area, _ = obs_index_start_and_num_areas(province_id)
  return unit_power_from_area(main_area, board_state)

def unit_power_from_area(area_id: AreaID,
                         board_state: np.ndarray) -> Optional[int]:
  if unit_type_from_area(area_id, board_state) is None:
    return None

  for power_id in range(NUM_POWERS):
    if board_state[area_id, OBSERVATION_UNIT_POWER_START + power_id]:
      return power_id
  raise ValueError('Expected a unit there, but none of the powers indicated')

def dislodged_unit_power(province_id: ProvinceID,
                         board_state: np.ndarray) -> Optional[int]:
  """Returns which power controls the unit province (None if no unit there)."""
  main_area, _ = obs_index_start_and_num_areas(province_id)
  return dislodged_unit_power_from_area(main_area, board_state)

def dislodged_unit_power_from_area(area_id: AreaID,
                                   board_state: np.ndarray) -> Optional[int]:
  if unit_type_from_area(area_id, board_state) is None:
    return None

  for power_id in range(NUM_POWERS):
    if board_state[area_id, OBSERVATION_DISLODGED_START + power_id]:
      return power_id
  raise ValueError('Expected a unit there, but none of the powers indicated')

def build_areas(country_index: int,
                board_state: np.ndarray) -> Sequence[AreaID]:
  """Returns all areas where it is legal for a power to build.

  Args:
    country_index: The power to get provinces for.
    board_state: Board from observation.
  """
  return np.where(
      np.logical_and(
          board_state[:, country_index + OBSERVATION_SC_POWER_START] > 0,
          board_state[:, OBSERVATION_BUILDABLE] > 0))[0]

def build_provinces(country_index: int,
                    board_state: np.ndarray) -> Sequence[ProvinceID]:
  """Returns all provinces where it is legal for a power to build.

  This returns province IDs, not area numbers.

  Args:
    country_index: The power to get provinces for.
    board_state: Board from observation.
  """
  buildable_provinces = []
  for a in build_areas(country_index, board_state):
    province_id, area_index = province_id_and_area_index(a)
    if area_index != 0:
      # We get only the main province.
      continue
    buildable_provinces.append(province_id)
  return buildable_provinces

def sc_provinces(country_index: int,
                 board_state: np.ndarray) -> Sequence[ProvinceID]:
  """Returns all supply centres the power owns.

  This returns province IDs, not area IDs.

  Args:
    country_index: The power to get provinces for.
    board_state: Board from observation.
  """
  sc_areas = np.where(
      board_state[:, country_index + OBSERVATION_SC_POWER_START] > 0)[0]
  provinces = []
  for a in sc_areas:
    province_id, area_index = province_id_and_area_index(a)
    if area_index != 0:
      # We get only the main province.
      continue
    provinces.append(province_id)
  return provinces

def removable_areas(country_index: int,
                    board_state: np.ndarray) -> Sequence[AreaID]:
  """Get all areas where it is legal for a power to remove."""
  return np.where(
      np.logical_and(
          board_state[:, country_index + OBSERVATION_UNIT_POWER_START] > 0,
          board_state[:, OBSERVATION_REMOVABLE] > 0))[0]

def removable_provinces(country_index: int,
                        board_state: np.ndarray) -> Sequence[ProvinceID]:
  """Get all provinces where it is legal for a power to remove."""
  remove_provinces = []
  for a in removable_areas(country_index, board_state):
    province_id, area_index = province_id_and_area_index(a)
    if area_index != 0:
      # We get only the main province.
      continue
    remove_provinces.append(province_id)
  return remove_provinces

def area_id_for_unit_in_province_id(province_id: ProvinceID,
                                    board_state: np.ndarray) -> AreaID:
  """AreaID from [0..80] of the unit in board_state."""
  if unit_type(province_id, board_state) is None:
    raise ValueError('No unit in province {}'.format(province_id))

  if (obs_index_start_and_num_areas(province_id)[1] == 3 and
      unit_type(province_id, board_state) == UnitType.FLEET):
    first_coast = area_from_province_id_and_area_index(province_id, 1)
    if unit_type_from_area(first_coast, board_state) is not None:
      return area_from_province_id_and_area_index(province_id, 1)
    else:
      return area_from_province_id_and_area_index(province_id, 2)
  else:
    return area_from_province_id_and_area_index(province_id, 0)

max_stars_repo_name: akhti/diplomacy

max_issues_repo_name: akhti/diplomacy

max_forks_repo_name: akhti/diplomacy


text:
import os
import shutil
from PIL import Image
import PIL
import numpy as np
import cv2
import torch
import numpy as np
import torch.nn.functional as F

def _fspecial_gauss_1d(size, sigma):
    r"""Create 1-D gauss kernel
    Args:
        size (int): the size of gauss kernel
        sigma (float): sigma of normal distribution
    Returns:
        torch.Tensor: 1D kernel (1 x 1 x size)
    """
    coords = torch.arange(size).to(dtype=torch.float)
    coords -= size//2

    g = torch.exp(-(coords**2) / (2*sigma**2))
    g /= g.sum()

    return g.unsqueeze(0).unsqueeze(0)

def gaussian_filter(input, win):
    r""" Blur input with 1-D kernel
    Args:
        input (torch.Tensor): a batch of tensors to be blured
        window (torch.Tensor): 1-D gauss kernel
    Returns:
        torch.Tensor: blured tensors
    """
    N, C, H, W = input.shape
    out = F.conv2d(input, win, stride=1, padding=0, groups=C)
    out = F.conv2d(out, win.transpose(2, 3), stride=1, padding=0, groups=C)
    return out

def _ssim(X, Y,
          data_range,
          win,
          size_average=True,
          K=(0.01,0.03)):

    r""" Calculate ssim index for X and Y
    Args:
        X (torch.Tensor): images
        Y (torch.Tensor): images
        win (torch.Tensor): 1-D gauss kernel
        data_range (float or int, optional): value range of input images. (usually 1.0 or 255)
        size_average (bool, optional): if size_average=True, ssim of all images will be averaged as a scalar
    Returns:
        torch.Tensor: ssim results.
    """
    K1, K2 = K
    batch, channel, height, width = X.shape
    compensation = 1.0

    C1 = (K1 * data_range)**2
    C2 = (K2 * data_range)**2

    win = win.to(X.device, dtype=X.dtype)

    mu1 = gaussian_filter(X, win)
    mu2 = gaussian_filter(Y, win)

    mu1_sq = mu1.pow(2)
    mu2_sq = mu2.pow(2)
    mu1_mu2 = mu1 * mu2

    sigma1_sq = compensation * ( gaussian_filter(X * X, win) - mu1_sq )
    sigma2_sq = compensation * ( gaussian_filter(Y * Y, win) - mu2_sq )
    sigma12   = compensation * ( gaussian_filter(X * Y, win) - mu1_mu2 )

    cs_map = (2 * sigma12 + C2) / (sigma1_sq + sigma2_sq + C2) # set alpha=beta=gamma=1
    ssim_map = ((2 * mu1_mu2 + C1) / (mu1_sq + mu2_sq + C1)) * cs_map

    ssim_per_channel = torch.flatten(ssim_map, 2).mean(-1)
    cs = torch.flatten( cs_map, 2 ).mean(-1)
    return ssim_per_channel, cs

def ssim(X, Y,
         data_range=255,
         size_average=True,
         win_size=11,
         win_sigma=1.5,
         win=None,
         K=(0.01, 0.03),
         nonnegative_ssim=False):
    r""" interface of ssim
    Args:
        X (torch.Tensor): a batch of images, (N,C,H,W)
        Y (torch.Tensor): a batch of images, (N,C,H,W)
        data_range (float or int, optional): value range of input images. (usually 1.0 or 255)
        size_average (bool, optional): if size_average=True, ssim of all images will be averaged as a scalar
        win_size: (int, optional): the size of gauss kernel
        win_sigma: (float, optional): sigma of normal distribution
        win (torch.Tensor, optional): 1-D gauss kernel. if None, a new kernel will be created according to win_size and win_sigma
        K (list or tuple, optional): scalar constants (K1, K2). Try a larger K2 constant (e.g. 0.4) if you get a negative or NaN results.
        nonnegative_ssim (bool, optional): force the ssim response to be nonnegative with relu
    Returns:
        torch.Tensor: ssim results
    """

    if len(X.shape) != 4:
        raise ValueError('Input images should be 4-d tensors.')

    if not X.type() == Y.type():
        raise ValueError('Input images should have the same dtype.')

    if not X.shape == Y.shape:
        raise ValueError('Input images should have the same shape.')

    if win is not None: # set win_size
        win_size = win.shape[-1]

    if not (win_size % 2 == 1):
        raise ValueError('Window size should be odd.')

    if win is None:
        win = _fspecial_gauss_1d(win_size, win_sigma)
        win = win.repeat(X.shape[1], 1, 1, 1)

    ssim_per_channel, cs = _ssim(X, Y,
                                data_range=data_range,
                                win=win,
                                size_average=False,
                                K=K)
    if nonnegative_ssim:
        ssim_per_channel = torch.relu(ssim_per_channel)

    if size_average:
        return ssim_per_channel.mean()
    else:
        return ssim_per_channel.mean(1)

def calc_psnr(sr, hr, scale=2, rgb_range=255, dataset=None):
    if hr.nelement() == 1: return 0

    diff = (sr - hr) / rgb_range
    if dataset and dataset.dataset.benchmark:
        shave = scale
        if diff.size(1) > 1:
            gray_coeffs = [65.738, 129.057, 25.064]
            convert = diff.new_tensor(gray_coeffs).view(1, 3, 1, 1) / 256
            diff = diff.mul(convert).sum(dim=1)
    else:
        shave = scale + 6

    valid = diff[..., shave:-shave, shave:-shave]
    mse = valid.pow(2).mean()

    return -10 * math.log10(mse)

def calc_ssim(sr, hr, scale, rgb_range):
    if hr.nelement() == 1: return 0

    shave = scale + 6

    X = sr[..., shave:-shave, shave:-shave]
    Y = hr[..., shave:-shave, shave:-shave]

    return ssim(X, Y,
         data_range=rgb_range,
         size_average=True,
         win_size=11,
         win_sigma=1.5,
         win=None,
         K=(0.01, 0.03),
         nonnegative_ssim=False)

predicted_img_dir = '/nfsdata1/AIM/X2/valid/HR'
ground_truth_dir = ''

idx = 0
score = 0.
for root, dirs, files in os.walk(predicted_img_dir):
    for name in files:
        if name[-3:] == 'png':
            print('Processing : ', name)
            idx += 1
            in_sr_path = os.path.join(root, name)
            sr = Image.open(in_sr_path)
            sr = np.array(sr, dtype=np.float32)
            sr = torch.from_numpy(sr).float().cuda()

            in_hr_path = os.path.join(ground_truth_dir, in_sr_path.split('/')[-1].replace('.', '_HR.'))
            hr = Image.open(in_hr_path)
            hr = np.array(hr, dtype=np.float32)
            hr = torch.from_numpy(hr).float().cuda()

            score += 0.5*calc_psnr(sr, hr, scale=2, rgb_range=255)/50. + 0.5*(calc_ssim(sr, hr, scale=2, rgb_range=255)-0.4)/0.6

print('Score :', score/idx)