text:
# -*- coding: utf-8 -*-
"""Routines for constructing estimates of dynamical quantities on trajectory
data using Galerkin expansion.

@author: Erik

"""
from __future__ import absolute_import, division, print_function, unicode_literals

import numpy as np
import scipy.linalg as spl

from .data_manipulation import get_initial_final_split, tlist_to_flat, flat_to_tlist

def compute_mfpt(basis, stateA, lag=1, dt=1.):
    """Calculates the mean first passage time into state A as a function of
    each point.

    Parameters
    ----------
    basis : list of trajectories
        Basis for the Galerkin expansion. Must be zero in state A.
    state_A : list of trajectories
        List of trajectories where each element is 1 or 0, corresponding to whether or not the datapoint is in state A.
    lag : int, optional
        Number of timepoints in the future to use for the finite difference in the discrete-time generator.  If not provided, uses value in the generator.
    timestep : scalar, optional
        Time between timepoints in the trajectory data.

    Returns
    -------
    mfpt : list of trajectories
        List of trajectories containing the values of the mean first passage time at each timepoint.

    """
    complement = [(A_i - 1.) for A_i in stateA]
    soln = compute_FK(basis, complement, lag=lag, dt=dt)

    post_processed_soln = []
    for ht_i in soln:
        ht_i[ht_i < 0.] = 0.
        post_processed_soln.append(ht_i)
    return post_processed_soln

def compute_committor(basis, guess_committor, lag=1):
    """Calculates the forward committor into state A as a function of
    each point.

    Parameters
    ----------
    basis : list of trajectories
        Basis for the Galerkin expansion. Must be zero in state A and B
    guess_committor : list of trajectories, optional
        The value of the guess function obeying the inhomogenous boundary conditions.
    lag : int, optional
        Number of timepoints in the future to use for the finite difference in the discrete-time generator.

    Returns
    -------
    committor: dynamical basis object
        List of trajectories containing the values of the forward committor at each point.

    """
    h = [np.zeros(gi.shape) for gi in guess_committor]
    soln = compute_FK(basis, h, r=guess_committor, lag=lag, dt=1.)

    post_processed_soln = []
    for qi in soln:
        qi[qi > 1.] = 1.
        qi[qi < 0.] = 0.
        post_processed_soln.append(qi)
    return post_processed_soln

def compute_bwd_committor(basis, guess_committor, stationary_com, lag=1):
    """Calculates the backward into state A as a function of
    each point.

    Parameters
    ----------
    basis : list of trajectories
        Basis for the Galerkin expansion. Must be zero in state A and B
    guess_committor : list of trajectories, optional
        The value of the guess function obeying the inhomogenous boundary conditions.
    stationary_com : list of trajectories
        Values of the change of measure to the stationary distribution.
    lag : int, optional
        Number of timepoints in the future to use for the finite difference in the discrete-time generator.

    Returns
    -------
    bwd_committor: dynamical basis object
        List of trajectories containing the values of the backward_committor at each point.
    """
    h = [np.zeros(gi.shape) for gi in guess_committor]
    soln = compute_adj_FK(basis, h, com=stationary_com, r=guess_committor, lag=lag, dt=1.)

    post_processed_soln = []
    for qi in soln:
        qi[qi > 1.] = 1.
        qi[qi < 0.] = 0.
        post_processed_soln.append(qi)
    return post_processed_soln

def compute_change_of_measure(basis, lag=1):
    """Calculates the value of the change of measure to the stationary distribution for each datapoint.

    Parameters
    ----------
    basis : list of trajectories
        Basis for the Galerkin expansion. Must be zero in state A and B
    lag : int, optional
        Number of timepoints in the future to use for the finite difference in the discrete-time generator.

    Returns
    -------
    change_of_measure : dynamical basis object
        List of trajectories containing the values of the change of measure to the stationary distribution at each point.

    """
    evals, evecs = compute_esystem(basis, lag, left=True, right=False)
    com = [np.real(ev_i[:, 0]) for ev_i in evecs]

    # Ensure positivity
    com_sign = np.sign(np.sum(com))
    com = [ci * com_sign for ci in com]
    post_processed_soln = []
    for ci in com:
        ci[ci < 0.] = 0.
        post_processed_soln.append(ci)
    return post_processed_soln

def compute_esystem(basis, lag=1, dt=1., left=False, right=True):
    """Calculates the eigenvectors and eigenvalues of the generator through
    Galerkin expansion.

    Parameters
    ----------
    basis : list of trajectories
        List of trajectories containing the basis for the Galerkin expansion.
        This method works much better if the basis set is zero on states A and B, however this is not a necessity.
    lag : int, optional
        Number of timepoints in the future to use for the finite difference in the discrete-time generator.
    left : bool, optional
        Whether or not to calculate the left eigenvectors  of the system.
    right : bool, optional
        Whether or not to calculate the right eigenvectors  of the system.

    Returns
    -------
    eigenvalues : numpy array
        Numpy array containing the eigenvalues of the generator.
    left_eigenvectors : list of trajectories, optional
        If left was set to true, the left eigenvectors are returned as a list of trajectories.
    right_eigenvectors : list of trajectories, optional
        If right was set to true, the right eigenvectors are returned as a list of trajectories.

    """
    if lag is None:
        lag = basis.lag
    basis_list = basis
    basis_flat_traj, traj_edges = tlist_to_flat(basis_list)
    K = compute_correlation_mat(basis_list, lag=lag)
    S = compute_stiffness_mat(basis_list, lag=lag)
    L = (K - S) / (lag * dt)

    # Calculate, sort eigensystem
    if (left and right):
        evals, evecs_l, evecs_r = spl.eig(L, b=S, left=True, right=True)
        evals, [evecs_l, evecs_r] = _sort_esystem(evals, [evecs_l, evecs_r])
        evecs_l = flat_to_tlist(np.dot(basis_flat_traj, evecs_l), traj_edges)
        evecs_r = flat_to_tlist(np.dot(basis_flat_traj, evecs_r), traj_edges)
        return evals, evecs_l, evecs_r
    elif (left or right):
        evals, evecs = spl.eig(L, b=S, left=left, right=right)
        evals, [evecs] = _sort_esystem(evals, [evecs])
        evecs = np.dot(basis_flat_traj, evecs)
        evecs = flat_to_tlist(evecs, traj_edges)
        return evals, evecs
    else:
        evals = spl.eig(L, b=S, left=False, right=False)
        return np.sort(evals)[::-1]

def _sort_esystem(evals, evec_collection):
    """Utility function that sorts a collection of eigenvetors in desceding
    order, and then sorts the eigenvectors accordingly.

    Parameters
    ----------
    evals : numpy array
        The unsorted eigenvalues.
    evec_collection : list of arrays
        List where each element in the list is a collection of eigenvectors.
        Each collection is sorted according to the eigenvalues.

    Returns
    -------
    sorted_evals : numpy array
        The sorted eigenvalues
    sorted_evecs : of arrays
        list where each element is a collection of sorted eigenvectors.

    """
    idx = evals.argsort()[::-1]
    sorted_evals = evals[idx]
    sorted_evecs = [evecs[:, idx] for evecs in evec_collection]
    return sorted_evals, sorted_evecs

def compute_FK(basis, h, r=None, lag=1, dt=1., return_coeffs=False):
    """
    Solves the forward Feynman-Kac problem Lg=h on a domain D, with boundary
    conditions g=b on the complement of D. To account for the boundary
    conditions, we solve the homogeneous problem Lg = h - Lr, where r is the
    provided guess.

    Parameters
    ----------
    traj_data : list of arrays OR single numpy array
        Value of the basis functions at every time point.  Should only be nonzero
        for points on the domain.
    h : list of 1d arrays or single 1d array
        Value of the RHS of the FK formula.  This should only be nonzero at
        points on the domain, Domain.
    r : list of 1d arrays or single 1d array, optional
        Value of the guess function.  Should be equal to b every point off of
        the domain.  IF not provided, the boundary conditions are assumed to be
        homogeneous.
    lag : int
        Number of timepoints in the future to use for the finite difference in
        the discrete-time generator.  If not provided, defaults to 1.
    timestep : scalar, optional
        Time between timepoints in the trajectory data.  Defaults to 1.

    Returns
    -------
    g : list of arrays
        Estimated solution to the Feynman-Kac problem.
    coeffs : ndarray
        Coefficients for the solution, only returned if return_coeffs is True.
    """
    h = [h_i.reshape(-1, 1).astype('float') for h_i in h]
    L_basis = compute_generator(basis, lag=lag, dt=dt)
    h_i = compute_stiffness_mat(basis, h, lag=lag)
    if r is not None:
        r = [r_i.reshape(-1, 1) for r_i in r]
        L_guess = compute_generator(basis, r, lag=lag, dt=dt)
        h_i -= L_guess
    coeffs = spl.solve(L_basis, h_i.ravel())

    # Construct solution vector
    N_traj = len(basis)
    soln = [np.dot(basis[i], coeffs) for i in range(N_traj)]
    if r is not None:
        soln = [soln[i] + r[i].ravel() for i in range(N_traj)]

    # Return calculated values
    if return_coeffs:
        return soln, coeffs
    else:
        return soln

def compute_adj_FK(basis, h, com=None, r=None, lag=1, dt=1., return_coeffs=False):
    """
    Solves the Feynman-Kac problem L^t dagger g=h on a domain D, with
    boundary conditions g=b on the complement of D. Here L^t is the adjoint of
    the generator with respect to the provided change of measure.  To account
    for the boundary conditions, we solve the homogeneous problem
    L^t g = h - L^t r, where r is the provided guess.

    Parameters
    ----------
    traj_data : list of arrays OR single numpy array
        Value of the basis functions at every time point.  Should only be nonzero
        for points on the domain.
    h : list of 1d arrays or single 1d array
        Value of the RHS of the FK formula.  This should only be nonzero at
        points on the domain, Domain.
    com : list of 1d arrays or single 1d array, optional
        Values of the change of measure against which to take the desired
        adjoint.  If not provided, takes the adjoint against the sampled meaure
    r : list of 1d arrays or single 1d array, optional
        Value of the guess function.  Should be equal to b every point off of
        the domain.  If not provided, the boundary conditions are assumed to be
        homogeneous.
    lag : int
        Number of timepoints in the future to use for the finite difference in
        the discrete-time generator.  If not provided, defaults to 1.
    timestep : scalar, optional
        Time between timepoints in the trajectory data.  Defaults to 1.

    Returns
    -------
    g : list of arrays
        Estimated solution to the Feynman-Kac problem.
    coeffs : ndarray
        Coefficients for the solution, only returned if return_coeffs is True.
    """
    h = [h_i.reshape(-1, 1).astype('float') for h_i in h]

    L_basis = compute_generator(basis, lag=lag, dt=dt, com=com)
    h_i = compute_stiffness_mat(basis, h, lag=lag, com=com)
    if r is not None:
        r = [r_i.reshape(-1, 1) for r_i in r]
        L_guess = compute_generator(r, basis, lag=lag, dt=dt, com=com).T
        h_i -= L_guess
    coeffs = spl.solve(L_basis.T, h_i.ravel())

    # Construct solution vector
    N_traj = len(basis)
    soln = [np.dot(basis[i], coeffs) for i in range(N_traj)]
    if r is not None:
        soln = [soln[i] + r[i].ravel() for i in range(N_traj)]

    # Return calculated values
    if return_coeffs:
        return soln, coeffs
    else:
        return soln

def compute_correlation_mat(Xs, Ys=None, lag=1, com=None):
    """Computes the time-lagged correlation matrix between two sets of observables.

    Parameters
    ----------
    Xs : list of trajectories
        List of trajectories for the first set of observables.
    Ys : list of trajectories, optional
        List of trajectories for the second set of observables.  If None, set to be X.
    lag : int, optional
        Lag to use in the correlation matrix.  Default is one step.
    com : list of trajectories
        Values of the change of measure against which to compute the average

    Returns
    -------
    K : numpy array
        The time-lagged correlation matrix between X and Y.
    """
    if Ys is None:
        Ys = Xs
    # Move to flat convention for easier processing.
    X_flat, traj_edges = tlist_to_flat(Xs)
    Y_flat = tlist_to_flat(Ys)[0]
    if com is not None:
        com_flat = tlist_to_flat(com)[0]
        X_flat *= com_flat

    # Split into initial, final points
    initial_indices, final_indices = get_initial_final_split(traj_edges, lag=lag)
    Y_t_lag = Y_flat[final_indices]
    X_t_0 = X_flat[initial_indices]

    # Calculate correlation, stiffness matrices.
    N = X_t_0.shape[0]
    K = X_t_0.T.dot(Y_t_lag) / float(N)
    return K

def compute_stiffness_mat(Xs, Ys=None, lag=1, com=None):
    """Computes the stiffness matrix between two sets of observables.

    Parameters
    ----------
    Xs : list of trajectories
        List of trajectories for the first set of observables.
    Ys : list of trajectories, optional
        List of trajectories for the second set of observables.  If None, set to be X.
    lag : int, optional
        Lag to use in the correlation matrix.  Default is one step.
        This is required as the stiffness is only evaluated over the initial points.
    com : list of trajectories
        Values of the change of measure against which to compute the average

    Returns
    -------
    S : numpy array
        The time-lagged stiffness matrix between X and Y.
    """
    if Ys is None:
        Ys = Xs
    # Move to flat convention for easier processing.
    X_flat, traj_edges = tlist_to_flat(Xs)
    Y_flat = tlist_to_flat(Ys)[0]
    if com is not None:
        com_flat = tlist_to_flat(com)[0]
        X_flat *= com_flat

    # Split into initial, final points
    initial_indices, final_indices = get_initial_final_split(traj_edges, lag=lag)
    Y_t_0 = Y_flat[initial_indices]
    X_t_0 = X_flat[initial_indices]
    N = X_t_0.shape[0]
    S = X_t_0.T.dot(Y_t_0) / float(N)
    return S

def compute_generator(Xs, Ys=None, lag=1, dt=1., com=None):
    """
    Computes the matrix of inner product elements against the generator.

    Parameters
    ----------
    Xs : list of trajectories
        List of trajectories for the first set of observables.
    Ys : list of trajectories, optional
        List of trajectories for the second set of observables.  If None, set to be X.
    lag : int, optional
        Lag to use in the correlation matrix.  Default is one step.
    dt: float, optional
        time per step of dynamics. Default is one time unit.
    com : list of trajectories
        Values of the change of measure against which to compute the average.

    Returns
    -------
    L : numpy array
        The approximation to the inner product <X, L Y>.
    """
    if Ys is None:
        Ys = Xs
    K = compute_correlation_mat(Xs, Ys, lag=lag, com=com)
    S = compute_stiffness_mat(Xs, Ys, lag=lag, com=com)
    return (K - S) / (lag * dt)

max_stars_repo_name: ehthiede/PyEDGAR

max_issues_repo_name: ehthiede/PyEDGAR

max_forks_repo_name: ehthiede/PyEDGAR

reason: import numpy,import scipy


text:
#!/usr/bin/env python3

import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

def fit_line(x, y):
    return 0, 0

def main():

if __name__ == "__main__":
    main()

ext: py

reason: import numpy


text:
import re
import numpy as np
import pytinydiffsim as dp
import meshcat_utils_dp
import meshcat

import time
np.set_printoptions(precision=3, suppress=True, linewidth=200)

def vecx_to_np(a):
  lst = []
  for i in range (a.size()):
    lst.append(a[i])
  return np.array(lst)

def vec3_to_np(a):
  lst = []
  for i in range(3):
    lst.append(a[i])
  return np.array(lst)
def vec4_to_np(a):
  lst = []
  for i in range(4):
    lst.append(a[i])
  return np.array(lst)

URDF_NAME = "laikago/laikago_toes_zup.urdf"
START_POS = [1.5, -2.0, 0.50]
START_ORN = [0,0,0,1]
MPC_BODY_MASS = 215/9.8
MPC_BODY_INERTIA = (0.07335, 0, 0, 0, 0.25068, 0, 0, 0, 0.25447)
MPC_BODY_HEIGHT = 0.42

MPC_VELOCITY_MULTIPLIER = 1.0

ACTION_REPEAT = 10

_IDENTITY_ORIENTATION=[0,0,0,1]
CHASSIS_NAME_PATTERN = re.compile(r"\w+_chassis_\w+")
MOTOR_NAME_PATTERN = re.compile(r"\w+_hip_motor_\w+")
KNEE_NAME_PATTERN = re.compile(r"\w+_lower_leg_\w+")
TOE_NAME_PATTERN = re.compile(r"jtoe\d*")
_DEFAULT_HIP_POSITIONS = (
    (0.21, -0.1157, 0),
    (0.21, 0.1157, 0),
    (-0.21, -0.1157, 0),
    (-0.21, 0.1157, 0),
)
_BODY_B_FIELD_NUMBER = 2
_LINK_A_FIELD_NUMBER = 3

HIP_JOINT_OFFSET = 0.0
UPPER_LEG_JOINT_OFFSET = 0#-0.6
KNEE_JOINT_OFFSET = 0#0.66

LAIKAGO_DEFAULT_ABDUCTION_ANGLE = 0.0#0.2
LAIKAGO_DEFAULT_HIP_ANGLE = 0.67-0.6
LAIKAGO_DEFAULT_KNEE_ANGLE = -1.25 +0.66
NUM_LEGS = 4
NUM_MOTORS = 12
# Bases on the readings from Laikago's default pose.
INIT_MOTOR_ANGLES = np.array([
    LAIKAGO_DEFAULT_ABDUCTION_ANGLE,
    LAIKAGO_DEFAULT_HIP_ANGLE,
    LAIKAGO_DEFAULT_KNEE_ANGLE
] * NUM_LEGS)
MOTOR_NAMES = [
    "FR_hip_motor_2_chassis_joint",
    "FR_upper_leg_2_hip_motor_joint",
    "FR_lower_leg_2_upper_leg_joint",
    "FL_hip_motor_2_chassis_joint",
    "FL_upper_leg_2_hip_motor_joint",
    "FL_lower_leg_2_upper_leg_joint",
    "RR_hip_motor_2_chassis_joint",
    "RR_upper_leg_2_hip_motor_joint",
    "RR_lower_leg_2_upper_leg_joint",
    "RL_hip_motor_2_chassis_joint",
    "RL_upper_leg_2_hip_motor_joint",
    "RL_lower_leg_2_upper_leg_joint",
]

#Use a PD controller
MOTOR_CONTROL_POSITION = 1
# Apply motor torques directly.
MOTOR_CONTROL_TORQUE = 2
# Apply a tuple (q, qdot, kp, kd, tau) for each motor. Here q, qdot are motor
# position and velocities. kp and kd are PD gains. tau is the additional
# motor torque. This is the most flexible control mode.
MOTOR_CONTROL_HYBRID = 3
MOTOR_CONTROL_PWM = 4 #only for Minitaur

MOTOR_COMMAND_DIMENSION = 5

# These values represent the indices of each field in the motor command tuple
POSITION_INDEX = 0
POSITION_GAIN_INDEX = 1
VELOCITY_INDEX = 2
VELOCITY_GAIN_INDEX = 3
TORQUE_INDEX = 4

class LaikagoMotorModel(object):
  """A simple motor model for Laikago.

    When in POSITION mode, the torque is calculated according to the difference
    between current and desired joint angle, as well as the joint velocity.
    For more information about PD control, please refer to:
    https://en.wikipedia.org/wiki/PID_controller.

    The model supports a HYBRID mode in which each motor command can be a tuple
    (desired_motor_angle, position_gain, desired_motor_velocity, velocity_gain,
    torque).

  """

  def __init__(self,
               kp,
               kd,
               torque_limits=None,
               motor_control_mode=MOTOR_CONTROL_POSITION):
    self._kp = kp
    self._kd = kd
    self._torque_limits = torque_limits
    if torque_limits is not None:
      if isinstance(torque_limits, (collections.Sequence, np.ndarray)):
        self._torque_limits = np.asarray(torque_limits)
      else:
        self._torque_limits = np.full(NUM_MOTORS, torque_limits)
    self._motor_control_mode = motor_control_mode
    self._strength_ratios = np.full(NUM_MOTORS, 1)

  def set_strength_ratios(self, ratios):
    """Set the strength of each motors relative to the default value.

    Args:
      ratios: The relative strength of motor output. A numpy array ranging from
        0.0 to 1.0.
    """
    self._strength_ratios = ratios

  def set_motor_gains(self, kp, kd):
    """Set the gains of all motors.

    These gains are PD gains for motor positional control. kp is the
    proportional gain and kd is the derivative gain.

    Args:
      kp: proportional gain of the motors.
      kd: derivative gain of the motors.
    """
    self._kp = kp
    self._kd = kd

  def set_voltage(self, voltage):
    pass

  def get_voltage(self):
    return 0.0

  def set_viscous_damping(self, viscous_damping):
    pass

  def get_viscous_dampling(self):
    return 0.0

  def convert_to_torque(self,
                        motor_commands,
                        motor_angle,
                        motor_velocity,
                        true_motor_velocity,
                        motor_control_mode=None):
    """Convert the commands (position control or torque control) to torque.

    Args:
      motor_commands: The desired motor angle if the motor is in position
        control mode. The pwm signal if the motor is in torque control mode.
      motor_angle: The motor angle observed at the current time step. It is
        actually the true motor angle observed a few milliseconds ago (pd
        latency).
      motor_velocity: The motor velocity observed at the current time step, it
        is actually the true motor velocity a few milliseconds ago (pd latency).
      true_motor_velocity: The true motor velocity. The true velocity is used to
        compute back EMF voltage and viscous damping.
      motor_control_mode: A MotorControlMode enum.

    Returns:
      actual_torque: The torque that needs to be applied to the motor.
      observed_torque: The torque observed by the sensor.
    """
    del true_motor_velocity
    if not motor_control_mode:
      motor_control_mode = self._motor_control_mode

    # No processing for motor torques
    if motor_control_mode is MOTOR_CONTROL_TORQUE:
      assert len(motor_commands) == NUM_MOTORS
      motor_torques = self._strength_ratios * motor_commands
      return motor_torques, motor_torques

    desired_motor_angles = None
    desired_motor_velocities = None
    kp = None
    kd = None
    additional_torques = np.full(NUM_MOTORS, 0)
    if motor_control_mode is MOTOR_CONTROL_POSITION:
      assert len(motor_commands) == NUM_MOTORS
      kp = self._kp
      kd = self._kd
      desired_motor_angles = motor_commands
      desired_motor_velocities = np.full(NUM_MOTORS, 0)
    elif motor_control_mode is MOTOR_CONTROL_HYBRID:
      # The input should be a 60 dimension vector
      assert len(motor_commands) == MOTOR_COMMAND_DIMENSION * NUM_MOTORS
      kp = motor_commands[POSITION_GAIN_INDEX::MOTOR_COMMAND_DIMENSION]
      kd = motor_commands[VELOCITY_GAIN_INDEX::MOTOR_COMMAND_DIMENSION]
      desired_motor_angles = motor_commands[
          POSITION_INDEX::MOTOR_COMMAND_DIMENSION]
      desired_motor_velocities = motor_commands[
          VELOCITY_INDEX::MOTOR_COMMAND_DIMENSION]
      additional_torques = motor_commands[TORQUE_INDEX::MOTOR_COMMAND_DIMENSION]
    motor_torques = -1 * (kp * (motor_angle - desired_motor_angles)) - kd * (
        motor_velocity - desired_motor_velocities) + additional_torques
    motor_torques = self._strength_ratios * motor_torques
    if self._torque_limits is not None:
      if len(self._torque_limits) != len(motor_torques):
        raise ValueError(
            "Torque limits dimension does not match the number of motors.")
      motor_torques = np.clip(motor_torques, -1.0 * self._torque_limits,
                              self._torque_limits)

    return motor_torques, motor_torques

class SimpleRobot(object):
  def __init__(self, simulation_time_step):
    self.MPC_BODY_MASS  = MPC_BODY_MASS
    self.MPC_BODY_INERTIA = MPC_BODY_INERTIA
    self.MPC_BODY_HEIGHT = MPC_BODY_HEIGHT

    self.time_step = simulation_time_step

    self.num_legs = NUM_LEGS
    self.num_motors = NUM_MOTORS

    self.skip_sync=0
    self.world = dp.TinyWorld()
    self.world.friction = 10.0

    self.mb_solver = dp.TinyMultiBodyConstraintSolver()
    self.vis = meshcat.Visualizer(zmq_url='tcp://127.0.0.1:6000')
    self.vis.delete()

    urdf_parser = dp.TinyUrdfParser()

    plane_urdf_data = urdf_parser.load_urdf("../../../data/plane_implicit.urdf")
    plane2vis = meshcat_utils_dp.convert_visuals(plane_urdf_data, "../../../data/checker_purple.png", self.vis, "../../../data/")
    self.plane_mb = dp.TinyMultiBody(False)
    plane2mb = dp.UrdfToMultiBody2()
    res = plane2mb.convert2(plane_urdf_data, self.world, self.plane_mb)

    self.urdf_data = urdf_parser.load_urdf("../../../data/laikago/laikago_toes_zup_joint_order.urdf")
    print("robot_name=",self.urdf_data.robot_name)
    self.b2vis = meshcat_utils_dp.convert_visuals(self.urdf_data, "../../../data/laikago/laikago_tex.jpg", self.vis, "../../../data/laikago/")
    is_floating=True
    self.mb = dp.TinyMultiBody(is_floating)
    urdf2mb = dp.UrdfToMultiBody2()
    self.res = urdf2mb.convert2(self.urdf_data, self.world, self.mb)
    self.mb.set_base_position(dp.TinyVector3(0,0,0.6))

    knee_angle = -0.5
    abduction_angle = 0.2

    self.initial_poses = [abduction_angle, 0., knee_angle, abduction_angle, 0., knee_angle,
                     abduction_angle, 0., knee_angle, abduction_angle, 0., knee_angle]

    qcopy = self.mb.q
    print("mb.q=",self.mb.q)
    print("qcopy=",qcopy)
    for q_index in range (12):
      qcopy[q_index+7]=self.initial_poses[q_index]
    print("qcopy=",qcopy)
    self.mb.set_q(qcopy)
    print("2 mb.q=",self.mb.q)
    print("self.mb.links=", self.mb.links)
    print("TDS:")
    for i in range (len(self.mb.links)):
      print("i=",i)
      l = self.mb.links[i]
      #print("l.link_name=", l.link_name)
      print("l.joint_name=", l.joint_name)
    dp.forward_kinematics(self.mb, self.mb.q, self.mb.qd)

    self._BuildJointNameToIdDict()
    self._BuildUrdfIds()
    self._BuildMotorIdList()
    self.ResetPose()
    self._motor_enabled_list = [True] * self.num_motors
    self._step_counter = 0
    self._state_action_counter = 0
    self._motor_offset= np.array([0.]*12)
    self._motor_direction= np.array([1.]*12)
    self.ReceiveObservation()
    self._kp = self.GetMotorPositionGains()
    self._kd = self.GetMotorVelocityGains()
    self._motor_model = LaikagoMotorModel(kp=self._kp, kd=self._kd, motor_control_mode=MOTOR_CONTROL_HYBRID)

    self.ReceiveObservation()
    self.mb.set_base_position(dp.TinyVector3(START_POS[0],START_POS[1],START_POS[2]))
    self.mb.set_base_orientation(dp.TinyQuaternion(START_ORN[0],START_ORN[1],START_ORN[2],START_ORN[3]))

    dp.forward_kinematics(self.mb, self.mb.q, self.mb.qd)
    meshcat_utils_dp.sync_visual_transforms(self.mb, self.b2vis, self.vis)

    self._SettleDownForReset(reset_time=1.0)

  def ResetPose(self):
    pass

  def _SettleDownForReset(self, reset_time):
    self.ReceiveObservation()
    if reset_time <= 0:
      return
    for _ in range(1500):
      self._StepInternal(
          INIT_MOTOR_ANGLES,
          motor_control_mode=MOTOR_CONTROL_POSITION)

  def _GetMotorNames(self):
    return MOTOR_NAMES

  def _BuildMotorIdList(self):
    self._motor_id_list = [
        self._joint_name_to_id[motor_name]
        for motor_name in self._GetMotorNames()
    ]

  def GetMotorPositionGains(self):
    return [220.]*self.num_motors

  def GetMotorVelocityGains(self):
    return np.array([2.]*self.num_motors)#np.array([1., 2., 2., 1., 2., 2., 1., 2., 2., 1., 2., 2.])

  def compute_jacobian(self, robot, link_id):
    """Computes the Jacobian matrix for the given link.

    Args:
      robot: A robot instance.
      link_id: The link id as returned from loadURDF.

    Returns:
      The 3 x N transposed Jacobian matrix. where N is the total DoFs of the
      robot. For a quadruped, the first 6 columns of the matrix corresponds to
      the CoM translation and rotation. The columns corresponds to a leg can be
      extracted with indices [6 + leg_id * 3: 6 + leg_id * 3 + 3].
    """

    world_point = self.mb.links[link_id].world_transform.translation
    local_point = self.mb.get_world_transform(-1).apply_inverse(world_point)
    jacdp = dp.point_jacobian(self.mb, link_id, local_point, True)
    #jacdp.print("jacdp")
    #print("number of rows=", jacdp.num_rows)
    #print("number of columns=", jacdp.num_columns)
    jakkes=[]
    diff=[]
    for r in range (jacdp.num_rows):
      for c in range (jacdp.num_columns):
        v = jacdp.get_at(r,c)
        jakkes.append(v)

    jacobian2 = np.array(jakkes).reshape(jacdp.num_rows,jacdp.num_columns)
    #print("jacobian2=", jacobian2)
    #print("diffjac[",link_id,"]=",diffjac)
    assert jacobian2.shape[0] == 3
    return jacobian2

  def ComputeJacobian(self, leg_id):
    """Compute the Jacobian for a given leg."""
    # Does not work for Minitaur which has the four bar mechanism for now.
    assert len(self._foot_link_ids) == self.num_legs
    return self.compute_jacobian(
        robot=self,
        link_id=self._foot_link_ids[leg_id],
    )

  def MapContactForceToJointTorques(self, leg_id, contact_force):
    """Maps the foot contact force to the leg joint torques."""
    jv = self.ComputeJacobian(leg_id)
    all_motor_torques = np.matmul(contact_force, jv)
    motor_torques = {}
    motors_per_leg = self.num_motors // self.num_legs
    com_dof = 6
    for joint_id in range(leg_id * motors_per_leg,
                          (leg_id + 1) * motors_per_leg):
      motor_torques[joint_id] = all_motor_torques[
          com_dof + joint_id] * self._motor_direction[joint_id]

    return motor_torques

  def GetBaseRollPitchYaw(self):
    """Get minitaur's base orientation in euler angle in the world frame.

    Returns:
      A tuple (roll, pitch, yaw) of the base in world frame.
    """
    orn = self.mb.get_base_orientation()

    roll_pitch_yaw2 = vec3_to_np(orn.get_euler_rpy())
    #print("roll_pitch_yaw=",roll_pitch_yaw)
    #print("roll_pitch_yaw2=",roll_pitch_yaw2)

    return np.asarray(roll_pitch_yaw2)

  def joint_angles_from_link_position(  self,    robot,    link_position,    link_id,    joint_ids,    position_in_world_frame,
              base_translation = (0, 0, 0),    base_rotation = (0, 0, 0, 1)):
    """Uses Inverse Kinematics to calculate joint angles.

    Args:
      robot: A robot instance.
      link_position: The (x, y, z) of the link in the body or the world frame,
        depending on whether the argument position_in_world_frame is true.
      link_id: The link id as returned from loadURDF.
      joint_ids: The positional index of the joints. This can be different from
        the joint unique ids.
      position_in_world_frame: Whether the input link_position is specified
        in the world frame or the robot's base frame.
      base_translation: Additional base translation.
      base_rotation: Additional base rotation.

    Returns:
      A list of joint angles.
    """

    if not position_in_world_frame:
      #tds
      base_world_tr = self.mb.get_world_transform(-1)
      local_base_tr = dp.TinySpatialTransform()
      local_base_tr.translation = dp.TinyVector3(base_translation[0],base_translation[1],base_translation[2])
      local_base_tr.rotation.setRotation(dp.TinyQuaternion(base_rotation[0],base_rotation[1],base_rotation[2],base_rotation[3]))
      world_tr = base_world_tr*local_base_tr
      local_link_tr = dp.TinySpatialTransform()
      local_link_tr.rotation.set_identity()
      local_link_tr.translation = dp.TinyVector3(link_position[0],link_position[1],link_position[2])
      link_tr = world_tr* local_link_tr
      world_link_pos = [link_tr.translation[0],link_tr.translation[1],link_tr.translation[2]]

    else:
      world_link_pos = link_position

    ik_solver = 0

    target_world_pos = dp.TinyVector3(world_link_pos[0],world_link_pos[1],world_link_pos[2])
    #target_local_pos = self.mb.get_world_transform(-1).apply_inverse(target_world_pos)

    all_joint_angles2 = dp.inverse_kinematics(self.mb, link_id, target_world_pos)
    #print("--")
    #print("len(all_joint_angles)=",len(all_joint_angles))
    #print("all_joint_angles=",all_joint_angles)

    vec = []
    for i in range (all_joint_angles2.size()-7):
      vec.append(all_joint_angles2[i+7])
    #print("len(all_joint_angles2)=",len(vec))
    #print("all_joint_angles2=",vec)
    vec = np.array(vec)
    #diff = vec - all_joint_angles
    #print("diff=",diff)
    # Extract the relevant joint angles.
    joint_angles = [vec[i] for i in joint_ids]
    return joint_angles

  def ComputeMotorAnglesFromFootLocalPosition(self, leg_id,
                                              foot_local_position):
    """Use IK to compute the motor angles, given the foot link's local position.

    Args:
      leg_id: The leg index.
      foot_local_position: The foot link's position in the base frame.

    Returns:
      A tuple. The position indices and the angles for all joints along the
      leg. The position indices is consistent with the joint orders as returned
      by GetMotorAngles API.
    """
    return self._EndEffectorIK(
        leg_id, foot_local_position, position_in_world_frame=False)
  def _EndEffectorIK(self, leg_id, position, position_in_world_frame):
    """Calculate the joint positions from the end effector position."""
    assert len(self._foot_link_ids) == self.num_legs
    toe_id = self._foot_link_ids[leg_id]
    motors_per_leg = self.num_motors // self.num_legs
    joint_position_idxs = [
        i for i in range(leg_id * motors_per_leg, leg_id * motors_per_leg +
                         motors_per_leg)
    ]
    joint_angles = self.joint_angles_from_link_position(
        robot=self,
        link_position=position,
        link_id=toe_id,
        joint_ids=joint_position_idxs,
        position_in_world_frame=position_in_world_frame)
    # Joint offset is necessary for Laikago.
    joint_angles = np.multiply(
        np.asarray(joint_angles) -
        np.asarray(self._motor_offset)[joint_position_idxs],
        self._motor_direction[joint_position_idxs])
    # Return the joing index (the same as when calling GetMotorAngles) as well
    # as the angles.
    return joint_position_idxs, joint_angles.tolist()

  def GetTimeSinceReset(self):
    return self._step_counter * self.time_step

  def GetHipPositionsInBaseFrame(self):
    return _DEFAULT_HIP_POSITIONS

  def GetBaseAngularVelocity(self):
    angvel = [self.mb.qd[0],self.mb.qd[1],self.mb.qd[2]]
    return angvel
  def GetBaseVelocity(self):
    """Get the linear velocity of minitaur's base.

    Returns:
      The velocity of minitaur's base.
    """
    #print("tds linvel")
    linvel = [self.mb.qd[3],self.mb.qd[4],self.mb.qd[5]]
    return linvel

  def GetTrueBaseOrientation(self):
    rn = self.mb.get_base_orientation()
    orn=[rn[0],rn[1],rn[2],rn[3]]
    return orn

  def TransformAngularVelocityToLocalFrame(self, angular_velocity, orientation):
    """Transform the angular velocity from world frame to robot's frame.

    Args:
      angular_velocity: Angular velocity of the robot in world frame.
      orientation: Orientation of the robot represented as a quaternion.

    Returns:
      angular velocity of based on the given orientation.
    """

    # Treat angular velocity as a position vector, then transform based on the
    # orientation given by dividing (or multiplying with inverse).
    # Get inverse quaternion assuming the vector is at 0,0,0 origin.
    tr = dp.TinySpatialTransform()
    tr.translation = dp.TinyVector3(1,2,3)
    tr.translation.set_zero()
    orn = dp.TinyQuaternion(orientation[0],orientation[1],orientation[2],orientation[3])
    tr.rotation.setRotation(orn)
    tr_inv = tr.get_inverse()
    rel_vel = tr.apply_inverse(dp.TinyVector3(angular_velocity[0],angular_velocity[1],angular_velocity[2]))
    #tr.print("tds trans")
    #print("tds rel_vel=",rel_vel)

    #print(tr)

    return vec3_to_np(rel_vel)

  def GetBaseRollPitchYawRate(self):
    """Get the rate of orientation change of the minitaur's base in euler angle.

    Returns:
      rate of (roll, pitch, yaw) change of the minitaur's base.
    """
    angular_velocity = self.GetBaseAngularVelocity()
    orientation = self.GetTrueBaseOrientation()

    return self.TransformAngularVelocityToLocalFrame(angular_velocity,orientation)

  def GetFootContacts(self):
    contacts = [False, False, False, False]
    return contacts

    #all_contacts = self.pybullet_client.getContactPoints(bodyA=self.quadruped)
    #
    #print("len(all_contacts)=", len(all_contacts))
    #for contact in all_contacts:
    #  # Ignore self contacts
    #  if contact[_BODY_B_FIELD_NUMBER] == self.quadruped:
    #    continue
    #  try:
    #    toe_link_index = self._foot_link_ids.index(
    #        contact[_LINK_A_FIELD_NUMBER])
    #    contacts[toe_link_index] = True
    #  except ValueError:
    #    continue
    #return contacts

  def GetTrueMotorAngles(self):
    """Gets the eight motor angles at the current moment, mapped to [-pi, pi].

    Returns:
      Motor angles, mapped to [-pi, pi].
    """
    self.ReceiveObservation()

    #motor_angles = [state[0] for state in self._joint_states]
    motor_angles = np.multiply(
        np.asarray(self.motor_angles) - np.asarray(self._motor_offset),
        self._motor_direction)
    return motor_angles

  def GetPDObservation(self):
    self.ReceiveObservation()
    observation = []
    observation.extend(self.GetTrueMotorAngles())
    observation.extend(self.GetTrueMotorVelocities())
    q = observation[0:self.num_motors]
    qdot = observation[self.num_motors:2 * self.num_motors]
    return (np.array(q), np.array(qdot))

  def GetTrueMotorVelocities(self):
    """Get the velocity of all eight motors.

    Returns:
      Velocities of all eight motors.
    """
    #motor_velocities = [state[1] for state in self._joint_states]

    motor_velocities = np.multiply(self.motor_velocities, self._motor_direction)
    return motor_velocities

  #def GetTrueObservation(self):
  #  self.ReceiveObservation()
  #  observation = []
  #  observation.extend(self.GetTrueMotorAngles())
  #  observation.extend(self.GetTrueMotorVelocities())
  #  observation.extend(self.GetTrueMotorTorques())
  #  observation.extend(self.GetTrueBaseOrientation())
  #  observation.extend(self.GetTrueBaseRollPitchYawRate())
  #  return observation

  def ApplyAction(self, motor_commands, motor_control_mode):
    """Apply the motor commands using the motor model.

    Args:
      motor_commands: np.array. Can be motor angles, torques, hybrid commands
      motor_control_mode: A MotorControlMode enum.
    """

    dp.forward_kinematics(self.mb, self.mb.q, self.mb.qd)
    dp.forward_dynamics(self.mb, dp.TinyVector3(0.,0.,-10.))

    mb_q = vecx_to_np(self.mb.q)[7:]
    mb_qd = vecx_to_np(self.mb.qd)[6:]

    actual_torque, observed_torque = self._motor_model.convert_to_torque(
                          motor_commands,
                          mb_q,
                          mb_qd,
                          mb_qd,
                          motor_control_mode=motor_control_mode)
    #print("actual_torque=",actual_torque)
    for i in range(len(actual_torque)):
      self.mb.tau[i] = actual_torque[i]

  def _SetMotorTorqueByIds(self, motor_ids, torques):
    pass

  def ReceiveObservation(self):
    self.motor_angles = vecx_to_np(self.mb.q)[7:]
    #print("tds motor_angles=",self.motor_angles)
    self.motor_velocities = vecx_to_np(self.mb.qd)[6:]
    #print("tds motor_velocities=",self.motor_velocities)

  def _StepInternal(self, action, motor_control_mode):
    self.ApplyAction(action, motor_control_mode)

    dp.integrate_euler_qdd(self.mb, self.time_step)

    multi_bodies = [self.plane_mb, self.mb]
    dispatcher = self.world.get_collision_dispatcher()
    contacts = self.world.compute_contacts_multi_body(multi_bodies,dispatcher)

    sync_to_pybullet = False

    #collision solver
    for cps in contacts:
      self.mb_solver.resolve_collision(cps, self.time_step)
    dp.integrate_euler(self.mb, self.time_step)

    joint_states = []
    for q_index in range (12):
      #qcopy[q_index+7]=self.initial_poses[q_index]
      joint_position = self.mb.q[q_index+7]
      joint_velocity = self.mb.qd[q_index+6]
      if sync_to_pybullet:
        self.pybullet_client.resetJointState(self.quadruped, self._motor_id_list[q_index], joint_position,joint_velocity)

    ps = self.mb.get_base_position()
    rn = self.mb.get_base_orientation()
    angvel = [self.mb.qd[0],self.mb.qd[1],self.mb.qd[2]]
    linvel = [self.mb.qd[3],self.mb.qd[4],self.mb.qd[5]]
    pos=[ps[0],ps[1],ps[2]]
    orn=[rn[0],rn[1],rn[2],rn[3]]

    if sync_to_pybullet:
      self.pybullet_client.resetBasePositionAndOrientation(self.quadruped, pos, orn)
      self.pybullet_client.resetBaseVelocity(self.quadruped, linvel, angvel)
    self.ReceiveObservation()

    self.skip_sync-=1
    if self.skip_sync<=0:
      #print("frame=",frame)
      self.skip_sync=16
      meshcat_utils_dp.sync_visual_transforms(self.mb, self.b2vis, self.vis)

    self._state_action_counter += 1

  def Step(self, action):
    """Steps simulation."""
    #if self._enable_action_filter:
    #  action = self._FilterAction(action)

    for i in range(ACTION_REPEAT):
      #proc_action = self.ProcessAction(action, i)
      proc_action = action
      self._StepInternal(proc_action, motor_control_mode=MOTOR_CONTROL_HYBRID)
      self._step_counter += 1

  def _BuildJointNameToIdDict(self):
    self._joint_name_to_id = {}
    print("TDS:")
    for joint_id in range(len(self.urdf_data.joints)):
      j = self.urdf_data.joints[joint_id]
      print("joint_id=",joint_id)
      print(j.joint_name)
      self._joint_name_to_id[j.joint_name] = joint_id

  def _BuildUrdfIds(self):
    """Build the link Ids from its name in the URDF file.

    Raises:
      ValueError: Unknown category of the joint name.
    """
    self._chassis_link_ids = [-1]
    self._leg_link_ids = []
    self._motor_link_ids = []
    self._knee_link_ids = []
    self._foot_link_ids = []

    print("TDS")
    num_joints = len(self.urdf_data.joints)
    for joint_id in range(num_joints):
      j = self.urdf_data.joints[joint_id]
      joint_name = j.joint_name
      print("link index=", joint_id)
      print("joint_name=",joint_name)
      joint_id = self._joint_name_to_id[joint_name]
      if CHASSIS_NAME_PATTERN.match(joint_name):
        self._chassis_link_ids.append(joint_id)
      elif MOTOR_NAME_PATTERN.match(joint_name):
        self._motor_link_ids.append(joint_id)
      # We either treat the lower leg or the toe as the foot link, depending on
      # the urdf version used.
      elif KNEE_NAME_PATTERN.match(joint_name):
        self._knee_link_ids.append(joint_id)
      elif TOE_NAME_PATTERN.match(joint_name):
        #assert self._urdf_filename == URDF_WITH_TOES
        self._foot_link_ids.append(joint_id)
      else:
        raise ValueError("Unknown category of joint %s" % joint_name)

    self._leg_link_ids.extend(self._knee_link_ids)
    self._leg_link_ids.extend(self._foot_link_ids)

    #assert len(self._foot_link_ids) == NUM_LEGS
    self._chassis_link_ids.sort()
    self._motor_link_ids.sort()
    self._knee_link_ids.sort()
    self._foot_link_ids.sort()
    self._leg_link_ids.sort()

    return

  def link_position_in_base_frame( self,   link_id ):
    """Computes the link's local position in the robot frame.

    Args:
      robot: A robot instance.
      link_id: The link to calculate its relative position.

    Returns:
      The relative position of the link.
    """
    link_local_frame = dp.link_transform_base_frame(self.mb, self.mb.q, link_id)
    link_local_position = vec3_to_np(link_local_frame.translation)

    return np.array(link_local_position)

  def GetFootLinkIDs(self):
    """Get list of IDs for all foot links."""
    return self._foot_link_ids

  def GetFootPositionsInBaseFrame(self):
    """Get the robot's foot position in the base frame."""
    assert len(self._foot_link_ids) == self.num_legs
    foot_positions = []
    for foot_id in self.GetFootLinkIDs():
      foot_positions.append(
          self.link_position_in_base_frame(link_id=foot_id)
          )
    return np.array(foot_positions)

max_stars_repo_name: eric-heiden/tds-merge

max_issues_repo_name: eric-heiden/tds-merge

max_forks_repo_name: eric-heiden/tds-merge


text:
import numpy as np
import time
import sys
from os.path import expanduser
home = expanduser("~")
sys.path.append(home+ "/pyrobot/src")
from pyrobot.core import Robot

def main():
    bot = Robot('tm700')
    bot.gripper.open()
    time.sleep(3)
    bot.gripper.close()

if __name__ == "__main__":
    main()


text:
import cv2
import tensorflow as tf
import numpy as np
import glob
import imageio
import os
import configparser
import re
import time

from auto_pose.meshrenderer import meshrenderer_phong
from auto_pose.ae.pysixd_stuff import misc

from auto_pose.ae import factory, utils

import argparse

import copy
import yaml
import tarfile
import six.moves.urllib as urllib
from tensorflow.core.framework import graph_pb2

# Protobuf Compilation (once necessary)
#os.system('protoc object_detection/protos/*.proto --python_out=.')

from auto_pose.test.googledet_utils import label_map_util
from auto_pose.test.googledet_utils.helper import FPS2, WebcamVideoStream, SessionWorker

import time

# helper function for split model
def _node_name(n):
    if n.startswith("^"):
        return n[1:]
    else:
        return n.split(":")[0]

# Load a (frozen) Tensorflow model into memory.
def load_frozenmodel():
    print('> Loading frozen model into memory')
    config = tf.ConfigProto(allow_soft_placement=True, log_device_placement=log_device)
    config.gpu_options.allow_growth=allow_memory_growth
    config.gpu_options.per_process_gpu_memory_fraction = 0.5 ###Jetson only
    if not split_model:
        detection_graph = tf.Graph()
        with detection_graph.as_default():
          od_graph_def = tf.GraphDef()
          with tf.gfile.GFile(model_path, 'rb') as fid:
            serialized_graph = fid.read()
            od_graph_def.ParseFromString(serialized_graph)
            tf.import_graph_def(od_graph_def, name='')
        return detection_graph, None, None

    else:
        # load a frozen Model and split it into GPU and CPU graphs
        # Hardcoded for ssd_mobilenet
        input_graph = tf.Graph()
        with tf.Session(graph=input_graph,config=config):
            if ssd_shape == 600:
                shape = 7326
                print 'ssd_shape = 600 :('
                exit()
            else:
                shape = 1917
            score = tf.placeholder(tf.float32, shape=(None, shape, num_classes), name="Postprocessor/convert_scores")
            expand = tf.placeholder(tf.float32, shape=(None, shape, 1, 4), name="Postprocessor/ExpandDims_1")
            for node in input_graph.as_graph_def().node:
                if node.name == "Postprocessor/convert_scores":
                    score_def = node
                if node.name == "Postprocessor/ExpandDims_1":
                    expand_def = node

        detection_graph = tf.Graph()
        with detection_graph.as_default():
            od_graph_def = tf.GraphDef()
            with tf.gfile.GFile(model_path, 'rb') as fid:
                serialized_graph = fid.read()
                od_graph_def.ParseFromString(serialized_graph)
                dest_nodes = ['Postprocessor/convert_scores','Postprocessor/ExpandDims_1']

                edges = {}
                name_to_node_map = {}
                node_seq = {}
                seq = 0
                for node in od_graph_def.node:
                    n = _node_name(node.name)
                    name_to_node_map[n] = node
                    edges[n] = [_node_name(x) for x in node.input]
                    node_seq[n] = seq
                    seq += 1
                for d in dest_nodes:
                    assert d in name_to_node_map, "%s is not in graph" % d

                nodes_to_keep = set()
                next_to_visit = dest_nodes[:]

                while next_to_visit:
                    n = next_to_visit[0]
                    del next_to_visit[0]
                    if n in nodes_to_keep: continue
                    nodes_to_keep.add(n)
                    next_to_visit += edges[n]

                nodes_to_keep_list = sorted(list(nodes_to_keep), key=lambda n: node_seq[n])
                nodes_to_remove = set()

                for n in node_seq:
                    if n in nodes_to_keep_list: continue
                    nodes_to_remove.add(n)
                nodes_to_remove_list = sorted(list(nodes_to_remove), key=lambda n: node_seq[n])

                keep = graph_pb2.GraphDef()
                for n in nodes_to_keep_list:
                    keep.node.extend([copy.deepcopy(name_to_node_map[n])])

                remove = graph_pb2.GraphDef()
                remove.node.extend([score_def])
                remove.node.extend([expand_def])
                for n in nodes_to_remove_list:
                    remove.node.extend([copy.deepcopy(name_to_node_map[n])])

                with tf.device('/gpu:0'):
                    tf.import_graph_def(keep, name='')
                with tf.device('/cpu:0'):
                    tf.import_graph_def(remove, name='')

        return detection_graph, score, expand

def load_labelmap():
    print('> Loading label map')
    label_map = label_map_util.load_labelmap(label_path)
    print label_map
    categories = label_map_util.convert_label_map_to_categories(label_map, max_num_classes=num_classes, use_display_name=True)
    print categories
    category_index = label_map_util.create_category_index(categories)
    print category_index
    return category_index

def detection(detection_graph, category_index, score, expand):
    print("> Building Graph")
    print category_index
    # Session Config: allow seperate GPU/CPU adressing and limit memory allocation
    config = tf.ConfigProto(allow_soft_placement=True, log_device_placement=log_device)
    config.gpu_options.allow_growth=allow_memory_growth
    config.gpu_options.per_process_gpu_memory_fraction = 0.4  ###Jetson only
    cur_frames = 0
    with detection_graph.as_default():
        #run_meta = tf.RunMetadata()
        with tf.Session(graph=detection_graph,config=config) as sess:
            # Define Input and Ouput tensors
            image_tensor = detection_graph.get_tensor_by_name('image_tensor:0')
            detection_boxes = detection_graph.get_tensor_by_name('detection_boxes:0')
            detection_scores = detection_graph.get_tensor_by_name('detection_scores:0')
            detection_classes = detection_graph.get_tensor_by_name('detection_classes:0')
            num_detections = detection_graph.get_tensor_by_name('num_detections:0')
            if split_model:
                score_out = detection_graph.get_tensor_by_name('Postprocessor/convert_scores:0')
                expand_out = detection_graph.get_tensor_by_name('Postprocessor/ExpandDims_1:0')
                score_in = detection_graph.get_tensor_by_name('Postprocessor/convert_scores_1:0')
                expand_in = detection_graph.get_tensor_by_name('Postprocessor/ExpandDims_1_1:0')
                # Threading
                gpu_worker = SessionWorker("GPU",detection_graph,config)
                cpu_worker = SessionWorker("CPU",detection_graph,config)
                gpu_opts = [score_out, expand_out]
                cpu_opts = [detection_boxes, detection_scores, detection_classes, num_detections]
                gpu_counter = 0
                cpu_counter = 0

            for i,experiment_name in enumerate(arguments.experiment_names):

                full_name = experiment_name.split('/')
                experiment_name = full_name.pop()
                experiment_group = full_name.pop() if len(full_name) > 0 else ''

                train_cfg_file_path = utils.get_config_file_path(workspace_path, experiment_name, experiment_group)
                train_args = configparser.ConfigParser()
                train_args.read(train_cfg_file_path)
                h_train, w_train, c = train_args.getint('Dataset','H'),train_args.getint('Dataset','W'), train_args.getint('Dataset','C')
                model_paths.append(train_args.get('Paths','MODEL_PATH'))
                all_train_args.append(train_args)

                log_dir = utils.get_log_dir(workspace_path,experiment_name,experiment_group)
                ckpt_dir = utils.get_checkpoint_dir(log_dir)

                all_codebooks.append(factory.build_codebook_from_name(experiment_name, experiment_group, return_dataset=False))
                factory.restore_checkpoint(sess, tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope=experiment_name)), ckpt_dir)

            #opts = tf.profiler.ProfileOptionBuilder.float_operation()
            #flops = tf.profiler.profile(sess.graph, run_meta=run_meta, cmd='op', options=opts)
            #exit()
            # i_class_mapping = {v: k for k, v in class_i_mapping.iteritems()}
            renderer = meshrenderer_phong.Renderer(
                model_paths,
                1
            )

            # Start Video Stream and FPS calculation
            fps = FPS2(fps_interval).start()
            video_stream = WebcamVideoStream(video_input,width,height).start()
            cur_frames = 0
            print("> Press 'q' to Exit, 'a' to start auto_pose")
            print('> Starting Detection')
            while video_stream.isActive():
                # actual Detection
                if split_model:
                    # split model in seperate gpu and cpu session threads
                    if gpu_worker.is_sess_empty():
                        # read video frame, expand dimensions and convert to rgb
                        image = video_stream.read()
                        image_expanded = np.expand_dims(cv2.cvtColor(image, cv2.COLOR_BGR2RGB), axis=0)
                        # put new queue
                        gpu_feeds = {image_tensor: image_expanded}
                        if visualize:
                            gpu_extras = image # for visualization frame
                        else:
                            gpu_extras = None
                        gpu_worker.put_sess_queue(gpu_opts,gpu_feeds,gpu_extras)

                    g = gpu_worker.get_result_queue()
                    if g is None:
                        # gpu thread has no output queue. ok skip, let's check cpu thread.
                        gpu_counter += 1
                    else:
                        # gpu thread has output queue.
                        gpu_counter = 0
                        score,expand,image = g["results"][0],g["results"][1],g["extras"]

                        if cpu_worker.is_sess_empty():
                            # When cpu thread has no next queue, put new queue.
                            # else, drop gpu queue.
                            cpu_feeds = {score_in: score, expand_in: expand}
                            cpu_extras = image
                            cpu_worker.put_sess_queue(cpu_opts,cpu_feeds,cpu_extras)

                    c = cpu_worker.get_result_queue()
                    if c is None:
                        # cpu thread has no output queue. ok, nothing to do. continue
                        cpu_counter += 1
                        time.sleep(0.005)
                        continue # If CPU RESULT has not been set yet, no fps update
                    else:
                        cpu_counter = 0
                        boxes, scores, classes, num, image = c["results"][0],c["results"][1],c["results"][2],c["results"][3],c["extras"]
                else:
                    # default session
                    image = video_stream.read()
                    image_expanded = np.expand_dims(cv2.cvtColor(image, cv2.COLOR_BGR2RGB), axis=0)
                    boxes, scores, classes, num = sess.run(
                            [detection_boxes, detection_scores, detection_classes, num_detections],
                            feed_dict={image_tensor: image_expanded})

                # Visualization of the results of a detection.

                H, W = image.shape[:2]

                img_crops = []
                det_bbs = []
                det_classes = []
                det_scores = []

                det_aae_bbs = []
                det_aae_objects_k = []
                #print vis_img.shape
                boxes = np.squeeze(boxes)
                scores = np.squeeze(scores)
                classes = np.squeeze(classes).astype(np.int32)

                highest_class_score = {clas:0.0 for clas in classes}
                for box,score, clas in zip(boxes,scores, classes):
                    if score > det_th and score > highest_class_score[clas]:

                        highest_class_score[clas] = score
                        ymin, xmin, ymax, xmax = (np.array(box)*np.array([height,width,height,width])).astype(np.int32)

                        h, w = (ymax-ymin,xmax-xmin)
                        det_bbs.append([xmin,ymin,w,h])
                        det_classes.append(clas)
                        det_scores.append(score)
                        if clas in clas_k_map:
                            det_aae_bbs.append([xmin,ymin,w,h])

                            det_aae_objects_k.append(clas_k_map[clas])

                            size = int(np.maximum(h, w) * train_args.getfloat('Dataset','PAD_FACTOR'))
                            cx = xmin + (xmax - xmin)/2
                            cy = ymin + (ymax - ymin)/2

                            left = np.maximum(cx-size/2, 0)
                            top = np.maximum(cy-size/2, 0)

                            img_crop = image[top:cy+size/2,left:cx+size/2]
                            img_crop = cv2.resize(img_crop, (h_train,w_train))

                            img_crop = img_crop / 255.
                            img_crops.append(img_crop)

                if len(det_aae_bbs) > 0:

                    Rs = []
                    ts = []
                    for k,bb,img_crop in zip(det_aae_objects_k,det_aae_bbs,img_crops):
                        R, t = all_codebooks[k].auto_pose6d(sess, img_crop, bb, K_test, 1, all_train_args[k], upright=False)
                        Rs.append(R.squeeze())
                        ts.append(t.squeeze())

                    Rs = np.array(Rs)
                    ts = np.array(ts)

                    bgr_y,_,_ = renderer.render_many(
                        obj_ids=np.array(det_aae_objects_k).astype(np.int32),
                        W=width/arguments.down,
                        H=height/arguments.down,
                        K=K_down,
                        Rs=Rs,
                        ts=ts,
                        near=1.,
                        far=10000.,
                        random_light=False,
                        # calc_bbs=False,
                        # depth=False
                    )

                    bgr_y = cv2.resize(bgr_y,(width,height))

                    g_y = np.zeros_like(bgr_y)
                    g_y[:,:,1]= bgr_y[:,:,1]
                    im_bg = cv2.bitwise_and(image,image,mask=(g_y[:,:,1]==0).astype(np.uint8))
                    image = cv2.addWeighted(im_bg,1,g_y,1,0)

                for bb,score,clas in zip(det_bbs,det_scores,det_classes):
                    xmin,ymin,xmax,ymax = bb[0],bb[1],bb[2]+bb[0],bb[1]+bb[3]
                    cv2.putText(image, '%s : %1.3f' % (category_index[clas]['name'],score), (xmin, ymax+20), cv2.FONT_ITALIC, .5, color_dict[clas-1], 2)
                    cv2.rectangle(image, (xmin,ymin),(xmax,ymax), color_dict[clas-1], 2)

                if vis_text:
                    cv2.putText(image,"fps: {}".format(fps.fps_local()), (10,30),
                                cv2.FONT_HERSHEY_SIMPLEX, 0.75, (77, 255, 9), 2)
                cv2.imshow('object_detection', image)
                # Exit Option
                key = cv2.waitKey(1)
                if key == ord('q'):
                    break

                fps.update()

    # End everything
    if split_model:
        gpu_worker.stop()
        cpu_worker.stop()
    fps.stop()
    video_stream.stop()
    cv2.destroyAllWindows()
    print('> [INFO] elapsed time (total): {:.2f}'.format(fps.elapsed()))
    print('> [INFO] approx. FPS: {:.2f}'.format(fps.fps()))

## LOAD CONFIG PARAMS ##
try:
    with open("googledet_utils/googledet_config.yml", 'r') as ymlfile:
        cfg = yaml.load(ymlfile)
except:
    print 'no config file found'
    exit()

video_input         = cfg['video_input']
visualize           = cfg['visualize']
vis_text            = cfg['vis_text']
width               = cfg['width']
height              = cfg['height']
fps_interval        = cfg['fps_interval']
allow_memory_growth = cfg['allow_memory_growth']
det_th              = cfg['det_th']
model_path          = cfg['model_path']
label_path          = cfg['label_path']
num_classes         = cfg['num_classes']
split_model         = cfg['split_model']
log_device          = cfg['log_device']
ssd_shape           = cfg['ssd_shape']
K_test              = cfg['K_test']

parser = argparse.ArgumentParser()
parser.add_argument("experiment_names", nargs='+',type=str)
parser.add_argument('-down', default=1, type=int)
parser.add_argument("-s", action='store_true', default=False)

# parser.add_argument("-gt_bb", action='store_true', default=False)
arguments = parser.parse_args()

workspace_path = os.environ.get('AE_WORKSPACE_PATH')
if workspace_path == None:
    print 'Please define a workspace path:\n'
    print 'export AE_WORKSPACE_PATH=/path/to/workspace\n'
    exit(-1)

all_codebooks = []
all_train_args = []
model_paths = []

K_test =  np.array(K_test).reshape(3,3)
K_down = K_test.copy()
if arguments.down > 1:
    K_down[:2,:] = K_down[:2,:] / arguments.down

aae_list = [full_name.pop() for full_name in [experiment_name.split('/') for experiment_name in arguments.experiment_names]]
color_dict = [(0,255,0),(0,0,255),(255,0,0),(255,255,0)] * 10

graph, score, expand = load_frozenmodel()
category = load_labelmap()

clas_k_map = {}
for _,val in category.iteritems():
    if val['name'] in aae_list:
        clas_k_map[int(val['id'])] = aae_list.index(val['name'])

detection(graph, category, score, expand)


text:
from sage_utils.amqp.workers import BaseRegisterWorker

class MicroserviceRegisterWorker(BaseRegisterWorker):

    def get_microservice_data(self, app):
        return {
            'name': self.app.config['SERVICE_NAME'],
            'version': self.app.config['SERVICE_VERSION'],
            'permissions': []
        }

max_stars_repo_name: OpenMatchmaking/microservice-strategist

max_issues_repo_name: OpenMatchmaking/microservice-strategist

max_forks_repo_name: OpenMatchmaking/microservice-strategist

reason: from sage


text:
import keras
from keras.models import load_model
import imageio
import os
import numpy as np
from scipy.misc import imresize
import sys, getopt, re, argparse
import tensorflow as tf
from time import time

# Keras' definition converted to numpy...
def softmax(x, axis=-1):
	ndim = np.ndim(x)
	if ndim >= 2:
		e = np.exp(x - np.max(x, axis=axis, keepdims=True))
		s = np.sum(e, axis=axis, keepdims=True)
		return e / s
	else:
		raise ValueError('Cannot apply softmax to a tensor that is 1D')

def main(argv):
	parser = argparse.ArgumentParser(description='Inference 3DConv Models')
	parser.add_argument('--mov_name', help='Name of movie to process')

	args = parser.parse_args()
	arg_dict = args.__dict__

	# Helix
	net = load_model('3Dconv_Keras.h5')
	# Remove the softmax layer...
	net.pop()
	net.compile('adam','categorical_crossentropy') # rebuild to set the output correctly
	print(net.summary())

	if 'mov_name' in arg_dict.keys() and arg_dict['mov_name'] is not None:
		video_pattern = os.path.splitext(args.mov_name)[0]

	reader = imageio.get_reader(video_pattern+'.avi')
	writer = imageio.get_writer(video_pattern+'_out_keras_BN.avi', fps=30, codec='mpeg4', quality=10)
	im_iter = reader.iter_data()

	input_size = 112
	time_depth = 16
	frames = [np.zeros([input_size, input_size, 1]) for x in range(time_depth)]

	framenum = 0
	while True:
		try:
			start_time = time()
			frames[0:time_depth-1] = np.copy(frames[1:time_depth])
			frame = np.uint8(next(im_iter))
			frame = imresize(frame, (input_size, input_size, 3))
			frame = frame[:,:,0]
			frame = np.reshape(frame, [input_size, input_size, 1])
			frames[time_depth-1] = frame
			batch_input_single = np.reshape(frames,[time_depth, input_size, input_size, 1])
			batch_input = np.reshape([batch_input_single,
									np.flipud(batch_input_single),
									np.fliplr(batch_input_single),
									np.fliplr(np.flipud(batch_input_single)),
									np.transpose(batch_input_single,(0,2,1,3)),
									np.transpose(np.flipud(batch_input_single),(0,2,1,3)),
									np.transpose(np.fliplr(batch_input_single),(0,2,1,3)),
									np.transpose(np.fliplr(np.flipud(batch_input_single)),(0,2,1,3))], [8, time_depth, input_size, input_size, 1])
			if framenum % 1000 == 0:
				print('Batch ' + str(framenum))
				print('Batch Assembled in: ' + str(time()-start_time))
			start_time = time()
			results_nosoftmax = net.predict(batch_input, batch_size=8)
			if framenum % 1000 == 0:
				print('Batch Computed in: ' + str(time()-start_time))
			start_time = time()
			results = softmax(results_nosoftmax)
			#print(results)
			predictions = np.argmax(results, 1)
			#prob = np.exp(results)[:,1] / np.sum(np.exp(results), axis=1)
			new_frame = np.zeros([input_size, input_size, 3])
			new_frame[:,:,0] = np.squeeze(frame)
			new_frame[:,:,1] = np.squeeze(frame)
			new_frame[:,:,2] = np.squeeze(frame)
			new_frame[:,0:15,:] = 255
			for i in range(np.shape(results)[0]):
				if results[i][1] < 0.50000001: # Negative
					level = 112-int(results[i][1]*112)
					new_frame[level:111,(0+2*i):(1+2*i),0] = 200
					new_frame[level:111,(0+2*i):(1+2*i),1] = 0
					new_frame[level:111,(0+2*i):(1+2*i),2] = 0
				else: # Positive
					level = 112-int(results[i][1]*112)
					new_frame[level:111,(0+2*i):(1+2*i),1] = 200
					new_frame[level:111,(0+2*i):(1+2*i),0] = 0
					new_frame[level:111,(0+2*i):(1+2*i),2] = 0
			# Tick marks at 25, 50, 75%
			new_frame[int(112/2),0:5,:] = 0
			new_frame[int(3*112/4),0:5,:] = 0
			new_frame[int(112/4),0:5,:] = 0
			# Consensus predictions
			mean_pred = np.mean(results[:,1]) > 0.5
			vote_pred = np.sum(predictions) > 4 # only majority, not half
			maxpool_pred = np.argmax(np.diag(results_nosoftmax[np.argmax(results_nosoftmax,0)])) == 1
			if mean_pred:
				new_frame[1:6,20:25,1] = 200
				new_frame[1:6,20:25,0] = 0
				new_frame[1:6,20:25,2] = 0
			else:
				new_frame[1:6,20:25,0] = 200
				new_frame[1:6,20:25,1] = 0
				new_frame[1:6,20:25,2] = 0
			if vote_pred:
				new_frame[1:6,30:35,1] = 200
				new_frame[1:6,30:35,0] = 0
				new_frame[1:6,30:35,2] = 0
			else:
				new_frame[1:6,30:35,0] = 200
				new_frame[1:6,30:35,1] = 0
				new_frame[1:6,30:35,2] = 0
			if maxpool_pred:
				new_frame[1:6,40:45,1] = 200
				new_frame[1:6,40:45,0] = 0
				new_frame[1:6,40:45,2] = 0
			else:
				new_frame[1:6,40:45,0] = 200
				new_frame[1:6,40:45,1] = 0
				new_frame[1:6,40:45,2] = 0
			writer.append_data(new_frame.astype('u1'))
			if framenum % 1000 == 0:
				print('Batch Saved in: ' + str(time()-start_time))
			framenum = framenum + 1
		except StopIteration:
			break
	reader.close()
	writer.close()

if __name__ == '__main__':
	main(sys.argv[1:])

max_stars_repo_name: KumarLabJax/MouseGrooming

max_issues_repo_name: KumarLabJax/MouseGrooming

max_forks_repo_name: KumarLabJax/MouseGrooming

reason: import numpy,from scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov  6 15:42:56 2021

@author: mohsenr
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import loadmat
from scipy.stats import zscore
import torch
from torch import nn
from torch.utils.data import DataLoader, random_split
from eelbrain import *
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from mne.decoding import GeneralizingEstimator, LinearModel, SlidingEstimator
import matplotlib.pyplot as plt
import os
import sklearn
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import roc_auc_score
#%%
r_model = []
r_null = []
data = []
scores = np.zeros((21,1394,1394))
for subj in range(21):
    tmp = []
    dataSub = loadmat(f"datasetCND_musicImagery/musicImagery/dataCND/dataSub{subj+1}.mat")
    for n in range(88):
        tmp.append(dataSub['eeg']['data'][0,0][0,n].T)
    data.append(np.asarray(tmp))

data = np.asarray(data)

chorals = np.asarray([i//11 for i in range(44)])
confusion = np.zeros((21, len(set(chorals)), len(set(chorals))))
met = loadmat('datasetCND_musicImagery/musicImagery/dataCND/metronome.mat')
met = met['met']
idx_met = np.where(met[:,0] != 0)[0]

for imet in idx_met:
    data[:,:,:,imet:imet+34] = 0

a_ma = np.ma.masked_equal(data,0)
data = np.array([i.compressed() for i in a_ma]).reshape(21,88,64,1803 - 34*len(idx_met))

#%%

device = "cpu"#"cuda" if torch.cuda.is_available() else "cpu"
print("Using {} device".format(device))

class Lis2Img(nn.Module):
    def __init__(self):
        super().__init__()
        # self.subject_embedding = nn.Embedding(21, 16)
        self.encoder = nn.Sequential(
            ########nn.Conv1d(64+16, 64, 4, stride=2, padding=2),
            nn.Conv1d(64, 32, 4, stride=2, padding=2),
            nn.ReLU(),
            nn.Conv1d(32,32, 4, stride=2, padding=2),
            nn.ReLU(),

        )
        # self.lstm = nn.LSTM(input_size=64, hidden_size=64, num_layers=2, batch_first=True)
        self.decoder = nn.Sequential(
            nn.ConvTranspose1d(32, 32, 4, stride=2, padding=2),
            nn.ReLU(),
            nn.ConvTranspose1d(32, 64, 4, stride=2, padding=2)
        )

    def forward(self, x):
        # x has shape [4, 64, 1803]
        # subj_id has shape [4]

        # subj = self.subject_embedding(subj_id)              # subj has shape [4, 16]
        # # Copy embedding vector into all time points
        # subj = subj.view(4, -1, 1).expand(-1, -1, 1802)     # subj has shape [4, 16, 1803]
        # Now, subj[:, :, 0] == subj[:, :, 1] == .... == subj[:, :, 1801]

        # Append the 16 subject embedding channels onto the 64 EEG channels
        ########x = torch.cat((x, subj), dim=1)
        # x now has shape [4, 64+16, 1803]

        # Conv1d input shape: [batch size, # channels, signal length]
        x = self.encoder(x)

        # LSTM input shape: [batch size, signal length, # features] -- when batch_first=True!
        x = x.transpose(1, 2)
        # x, _ = self.lstm(x)
        x = x.transpose(2, 1)

        x = self.decoder(x)
        return x

# model = Lis2Img().to(device)
# model.to(dtype=torch.double)
# print(model)

# loss_fn = nn.MSELoss()
# optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

def reset_weights(model):
    for layer in model.children():
        if hasattr(layer, 'reset_parameters'):
            layer.reset_parameters()

def R_value(a, b):                      # Inputs shape [4, 64, 1803]
    dim = -1                            # Correlation over time points dimension
    a = a - a.mean(axis=dim, keepdim=True)
    b = b - b.mean(axis=dim, keepdim=True)
    cov = (a * b).mean(dim)
    na, nb = [(i**2).mean(dim)**0.5 for i in [a, b]]
    norms = na * nb
    R_matrix = cov / norms
    return R_matrix                     # [4, 64] - R values per channel for each trial

def train(x_dataloader, y_dataloader, model, loss_fn, optimizer):
    assert (len(x_dataloader.dataset) == len(y_dataloader.dataset)), \
                "Input, output dataloaders have different lengths! :O"

    size = len(x_dataloader.dataset)
    for batch, (X, Y) in enumerate(zip(x_dataloader, y_dataloader)):
        X, Y = X.to(device), Y.to(device)
        # Compute prediction error
        pred = model(X)
        loss = loss_fn(pred, Y)

        # Backpropagation
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        if batch % 50 == 0:
            loss, current = loss.item(), batch * len(X)
            print(f"[{current:>4d}/{size:>4d}]   loss: {loss}")

def test(x_dataloader, y_dataloader, model, loss_fn):
    assert (len(x_dataloader) == len(y_dataloader)), \
                "Input, output dataloaders have different lengths! :O"

    num_batches = len(x_dataloader)
    model.eval()
    avg_loss = 0
    R_values = []

    with torch.no_grad():
        for (X, Y) in zip(x_dataloader, y_dataloader):
            X, Y = X.to(device), Y.to(device)
            pred = model(X)
            avg_loss += loss_fn(pred, Y).item()
            R = R_value(pred, Y)
            R_values.extend(R.flatten().tolist())

    avg_loss /= num_batches
    avg_R = np.asarray(R_values).mean()
    print("Test:")
    print(f"\tAvg loss: {avg_loss}")
    print(f"\tAvg R value: {avg_R}\n")

acc = np.zeros((21,))
for test_sbj_idx in range(0,21):
    # model.apply(reset_weights)

    # loss_fn = nn.MSELoss()
    # optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
    ##
    model = Lis2Img().to(device)
    model.to(dtype=torch.double)
    loss_fn = nn.MSELoss()
    optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
    ##

    train_sbjs_idx = list(range(21))
    train_sbjs_idx.remove(test_sbj_idx)

    train_sbjs = data[train_sbjs_idx,:,:,:1394]#:1802]
    test_sbj = data[test_sbj_idx,:,:,:1394]#:1802]

    # (1) make data to pre-train the model as and AE for 20 subjects (all data from 20 sbjs)
    # (2) fix the encoder train the decoder part for lis_to_img (devide into train and test)
    # (3) test the network for the left-out sbj (all data from left_out sbjs)

    data_tr = train_sbjs.reshape(20*88,64,1394)                     #(1)
    data_lis = train_sbjs[:,:44,:,:].reshape(20*44,64,1394)         #(2a)
    data_img = train_sbjs[:,-44:,:,:].reshape(20*44,64,1394)        #(2b)
    tst_sbj_lis = test_sbj[:44,:,:]                                 #(3a)
    tst_sbj_img = test_sbj[-44:,:,:]                                #(3b)

    [all_train, all_test] = random_split(data_tr, [1700, 60],        #(1)
                                         generator=torch.Generator().manual_seed(42))

    [lis_train, lis_test] = random_split(data_lis, [800, 80],       #(2)
                                 generator=torch.Generator().manual_seed(42))
    [img_train, img_test] = random_split(data_img, [800, 80],        #(2))
                                 generator=torch.Generator().manual_seed(42))

    # Create dataloaders
    batch_size = 4
    all_train_dataloader = DataLoader(all_train, batch_size=batch_size)
    all_test_dataloader = DataLoader(all_test, batch_size=batch_size)

    lis_train_dataloader = DataLoader(lis_train, batch_size=batch_size)
    lis_test_dataloader = DataLoader(lis_test, batch_size=batch_size)
    img_train_dataloader = DataLoader(img_train, batch_size=batch_size)
    img_test_dataloader = DataLoader(img_test, batch_size=batch_size)

    epochs = 30
    print("[[Pre-training autoencoder]]")
    for t in range(epochs):
        print(f"Epoch {t+1}\n-------------------------------")
        train(all_train_dataloader, all_train_dataloader, model, loss_fn, optimizer)
        test(all_test_dataloader, all_test_dataloader, model, loss_fn)

    print("[[Fix decoder weights]]")
    model.encoder.requires_grad_(False)
    epochs = 30
    print("[[Pre-training autoencoder]]")
    for t in range(epochs):
        print(f"Epoch {t+1}\n-------------------------------")
        train(lis_train_dataloader, img_train_dataloader, model, loss_fn, optimizer)
        test(lis_test_dataloader, img_test_dataloader, model, loss_fn)
    print("Done!\n")

    ###### lis2img for the unseen subject
    [test_sbj_lis, _] = random_split(tst_sbj_lis, [44, 0], generator=torch.Generator().manual_seed(42))
    [test_sbj_img, _] = random_split(tst_sbj_img, [44, 0], generator=torch.Generator().manual_seed(42))
    test_lis = DataLoader(test_sbj_lis, batch_size=batch_size)
    test_img = DataLoader(test_sbj_img, batch_size=batch_size)
    print(f'Test the network on Subject: {test_sbj_idx + 1}')
    test(test_lis, test_img, model, loss_fn)

    X = torch.tensor(test_sbj_lis)
    Y = torch.tensor(test_sbj_img)
    Pred = model(X).detach()
    # tmp = R_value(Y,Pred)
    # r_model.append(np.asarray(tmp))
    # save.pickle(np.asarray(Pred), f'../Pred_img/sbj{test_sbj_idx}_Pred.pkl')
    ############ Linear Classification
    clf = make_pipeline(StandardScaler(),  # z-score normalization
                    LinearModel(LogisticRegression(penalty = 'l2', C=1,
                                solver = 'liblinear', class_weight='balanced')))
    decoder = SlidingEstimator(clf,scoring=None,n_jobs=-1,verbose=True)
    X = np.asarray(Pred)
    y = chorals[test_sbj_img.indices]
    classes = set(y)
    cv = StratifiedKFold(n_splits=5, random_state=0, shuffle=True)
    ## Compute confusion matrix for each cross-validation fold
    y_pred = np.zeros((len(y), 1394 ,len(classes)))
    # decoder.fit(X = np.asarray(Y), y = chorals[test_sbj_img.indices])
    # scores[test_sbj_idx,:,:] = decoder.score(X = np.asarray(Pred), y = chorals[test_sbj_img.indices])
    for train_t, test_t in cv.split(X, y):
        # Fit
        decoder.fit(X[train_t], y[train_t])
        # Probabilistic prediction (necessary for ROC-AUC scoring metric)
        y_pred[test_t] = decoder.predict_proba(X[test_t])

    tmp = y_pred.mean(axis = 1)
    acc[test_sbj_idx] = sum(tmp.argmax(1) == y)/44
    for ii, train_class in enumerate(classes):
        for jj in range(ii, len(classes)):
            confusion[test_sbj_idx, ii, jj] = roc_auc_score(y == train_class,tmp[:, jj])
            confusion[test_sbj_idx, jj, ii] = confusion[test_sbj_idx, ii, jj]

    # fig, ax = plt.subplots(1)
    # im = ax.matshow(scores[test_sbj_idx,:,:], vmin=.0, vmax=.5,
    #                 cmap='RdBu_r', origin='lower',)# extent = [0, 100, 0, 100])
    # # ax.contour(p_val, levels = [-0.1, 0.05], colors = 'k', linestyles = '--',
    # #              linewidths =.35)
    # ax.axhline(0., color='k')
    # ax.axvline(0., color='k')
    # ax.xaxis.set_ticks_position('bottom')
    # ax.set_xlabel('Testing Time (s)')
    # ax.set_ylabel('Training Time (s)')
    # ax.set_title(f'train on True, test on Pred: sbject{test_sbj_idx+1}')
    # clb = plt.colorbar(im)
    # clb.set_label('Accuracy  (chance = 0.25)')
    # plt.tight_layout()

    labels = [f'Choral {i}' for i in range(1,5)]
    fig, ax = plt.subplots(1)
    im = ax.matshow(confusion[test_sbj_idx,:,:], cmap='RdBu_r', clim=[0., 1])
    ax.set_yticks(range(len(classes)))
    ax.set_yticklabels(labels)
    ax.set_xticks(range(len(classes)))
    ax.set_xticklabels(labels, rotation=0, ha='center')
    ax.set_title(f'predicted Imagined Subject {test_sbj_idx+1}_{round(acc[test_sbj_idx],2)}')
    # ax.axhline(4, color='k')
    # ax.axvline(4, color='k')
    clt = plt.colorbar(im)
    clt.set_label('AUC')
    plt.tight_layout()
    plt.show()

    ###### lis2img for control
    [test_sbj_lis, _] = random_split(tst_sbj_lis, [44, 0], generator=torch.Generator().manual_seed(4))
    [test_sbj_img, _] = random_split(tst_sbj_img, [44, 0], generator=torch.Generator().manual_seed(42))
    test_lis = DataLoader(test_sbj_lis, batch_size=batch_size)
    test_img = DataLoader(test_sbj_img, batch_size=batch_size)
    print(f'Test the network on CONTROL Subject: {test_sbj_idx + 1}')
    test(test_lis, test_img, model, loss_fn)

    X = torch.tensor(test_sbj_lis)
    Y = torch.tensor(test_sbj_img)
    Pred = model(X).detach()
    tmp = R_value(Y,Pred)
    r_null.append(np.asarray(tmp))

#%%

idx = np.argsort(r_model.mean(1))
idx = idx[::-1]
plt.scatter(np.arange(1,23,1),r_model[idx,:].mean(axis = 1))
plt.errorbar(np.arange(1,23,1),r_model[idx,:].mean(axis = 1),
             yerr=r_model[idx,:].std(axis = 1)/np.sqrt(21))
plt.scatter(np.arange(1,23,1),r_null[idx,:].mean(axis = 1))
plt.errorbar(np.arange(1,23,1),r_null[idx,:].mean(axis = 1),
             yerr=r_null[idx,:].std(axis = 1)/np.sqrt(21))
plt.xticks(np.arange(1, 23, step=1))
plt.tight_layout()
# plt.savefig('../Figures/Participants_N-1.pdf')


text:
# Copyright 2022 Huawei Technologies Co., Ltd
#
# 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.
# ============================================================================
"""ONNX evaluation"""
import pickle
import time

import numpy as np
import onnxruntime as ort
from src.dataset import load_dataset
from src.gnmt_model.bleu_calculate import bleu_calculate
from src.dataset.tokenizer import Tokenizer
from src.utils.get_config import get_config

from model_utils.config import config as default_config

def create_session(checkpoint_path, target_device):
    """Create ONNX runtime session"""
    if target_device == 'GPU':
        providers = ['CUDAExecutionProvider']
    elif target_device in ('CPU', 'Ascend'):
        providers = ['CPUExecutionProvider']
    else:
        raise ValueError(f"Unsupported target device '{target_device}'. Expected one of: 'CPU', 'GPU', 'Ascend'")
    session = ort.InferenceSession(checkpoint_path, providers=providers)
    input_names = [x.name for x in session.get_inputs()]
    return session, input_names

def infer(config):
    """Run inference"""
    session, [ids_name, mask_name] = create_session(config.file_name, config.device_target)
    eval_dataset = load_dataset(data_files=config.test_dataset,
                                batch_size=config.batch_size,
                                sink_mode=config.dataset_sink_mode,
                                drop_remainder=False,
                                is_translate=True,
                                shuffle=False)

    predictions = []
    source_sentences = []

    batch_index = 1
    pad_idx = 0
    sos_idx = 2
    eos_idx = 3
    source_ids_pad = np.tile(np.array([[sos_idx, eos_idx] + [pad_idx] * (config.seq_length - 2)], np.int32),
                             [config.batch_size, 1])
    source_mask_pad = np.tile(np.array([[1, 1] + [0] * (config.seq_length - 2)], np.int32),
                              [config.batch_size, 1])
    for batch in eval_dataset.create_dict_iterator(output_numpy=True):
        source_sentences.append(batch["source_eos_ids"])
        source_ids = batch["source_eos_ids"]
        source_mask = batch["source_eos_mask"]

        active_num = source_ids.shape[0]
        if active_num < config.batch_size:
            source_ids = np.concatenate((source_ids, source_ids_pad[active_num:, :]))
            source_mask = np.concatenate((source_mask, source_mask_pad[active_num:, :]))

        start_time = time.time()
        [predicted_ids] = session.run(None, {ids_name: source_ids, mask_name: source_mask})

        print(f" | BatchIndex = {batch_index}, Batch size: {config.batch_size}, active_num={active_num}, "
              f"Time cost: {time.time() - start_time}.")
        if active_num < config.batch_size:
            predicted_ids = predicted_ids[:active_num, :]
        batch_index = batch_index + 1
        predictions.append(predicted_ids)

    output = []
    for inputs, batch_out in zip(source_sentences, predictions):
        for i, _ in enumerate(batch_out):
            if batch_out.ndim == 3:
                batch_out = batch_out[:, 0]

            example = {
                "source": inputs[i].tolist(),
                "prediction": batch_out[i].tolist()
            }
            output.append(example)

    return output

def run_onnx_eval():
    """ONNX eval"""
    config = get_config(default_config)
    result = infer(config)

    with open(config.output, "wb") as f:
        pickle.dump(result, f, 1)

    result_npy_addr = config.output
    vocab = config.vocab
    bpe_codes = config.bpe_codes
    test_tgt = config.test_tgt
    tokenizer = Tokenizer(vocab, bpe_codes, 'en', 'de')
    scores = bleu_calculate(tokenizer, result_npy_addr, test_tgt)
    print(f"BLEU scores is :{scores}")

if __name__ == '__main__':
    run_onnx_eval()

max_stars_repo_name: mindspore-ai/models

max_issues_repo_name: mindspore-ai/models

max_forks_repo_name: mindspore-ai/models


text:
#!/usr/bin/env python3
# coding: utf-8
"""
PanelResolver class

Copyright 2017 MicaSense, Inc.

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in the
Software without restriction, including without limitation the rights to use,
copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the
Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""

import math
import numpy as np
import cv2
import re
import pyzbar.pyzbar as pyzbar

from skimage import measure
import matplotlib.pyplot as plt

import micasense.image as image

class Panel(object):
    def __init__(
        self, img: image.Image, panel_corners=None, ignore_autocalibration=False
    ):
        # if we have panel images with QR metadata, panel detection is not called,
        # so this can be forced here
        if img is None:
            raise IOError("Must provide an image")

        self.image = img
        bias = img.radiance().min()
        scale = img.radiance().max() - bias
        self.gray8b = np.zeros(img.radiance().shape, dtype="uint8")
        cv2.convertScaleAbs(
            img.undistorted(img.radiance()),
            self.gray8b,
            256.0 / scale,
            -1.0 * scale * bias,
        )

        if (self.image.auto_calibration_image) and ~ignore_autocalibration:
            self.__panel_type = "auto"  # panels the camera found we call auto
            if panel_corners is not None:
                self.__panel_bounds = np.array(panel_corners)
            else:
                self.__panel_bounds = np.array(self.image.panel_region)
            self.panel_albedo = self.image.panel_albedo
            self.serial = self.image.panel_serial
            self.qr_area = None
            self.qr_bounds = None
            self.panel_std = None
            self.saturated_panel_pixels_pct = None
            self.panel_pixels_mean = None
            self.panel_version = None
            if re.search(r"RP\d{2}-(\d{7})-\D{2}", self.image.panel_serial):
                self.serial = self.image.panel_serial
                self.panel_version = int(self.image.panel_serial[2:4])
        else:
            self.__panel_type = "search"  # panels we search for we call search
            self.serial = None
            self.qr_area = None
            self.qr_bounds = None
            self.panel_std = None
            self.saturated_panel_pixels_pct = None
            self.panel_pixels_mean = None
            self.panel_version = None
            if panel_corners is not None:
                self.__panel_bounds = np.array(panel_corners)
            else:
                self.__panel_bounds = None

    def __expect_panel(self):
        return self.image.band_name.upper() != "LWIR"

    def __find_qr(self):
        decoded = pyzbar.decode(self.gray8b, symbols=[pyzbar.ZBarSymbol.QRCODE])
        for symbol in decoded:
            serial_str = symbol.data.decode("UTF-8")
            m = re.search(r"RP\d{2}-(\d{7})-\D{2}", serial_str)
            if m:
                self.serial = serial_str
                self.panel_version = int(self.serial[2:4])
                self.qr_bounds = []
                for point in symbol.polygon:
                    self.qr_bounds.append([point.x, point.y])
                self.qr_bounds = np.asarray(self.qr_bounds, np.int32)
                self.qr_area = cv2.contourArea(self.qr_bounds)
                # print (symbol.polygon)
                # print (self.qr_bounds)
                break

    def __pt_in_image_bounds(self, pt):
        width, height = self.image.size()
        if pt[0] >= width or pt[0] < 0:
            return False
        if pt[1] >= height or pt[1] < 0:
            return False
        return True

    def reflectance_from_panel_serial(self):
        if self.__panel_type == "auto":
            return self.panel_albedo

        if self.serial is None:
            self.__find_qr()
        if self.serial is None:
            raise ValueError("Panel serial number not found")
        if self.panel_version >= 4:
            min_wl = float(self.serial[-14:-10])
            min_rf = float(self.serial[-10:-7]) / 1000.0
            max_wl = float(self.serial[-7:-3])
            max_rf = float(self.serial[-3:]) / 1000.0
            c = np.polyfit([min_wl, max_wl], [min_rf, max_rf], 1)
            p = np.poly1d(c)
            return p(self.image.center_wavelength)
        else:
            return None

    def qr_corners(self):
        if self.__panel_type == "auto":
            return None

        if self.qr_bounds is None:
            self.__find_qr()
        return self.qr_bounds

    def panel_detected(self):
        if self.__expect_panel() is False:
            return False

        if self.__panel_type == "auto":
            return True

        if self.serial is None:
            self.__find_qr()
        return self.qr_bounds is not None

    def panel_corners(self):
        """get the corners of a panel region based on the qr code location
        Our algorithm to do this uses a 'reference' qr code location and
        it's associate panel region.  We find the affine transform
        between the reference qr and our qr, and apply that same transform to the
        reference panel region to find our panel region. Because of a limitation
        of the pyzbar library, the rotation of the absolute QR code isn't known,
        so we then try all 4 rotations and test against a cost function which is the
        minimum of the standard devation divided by the mean value for the panel region"""
        if self.__panel_bounds is not None:
            return self.__panel_bounds
        if self.serial is None:
            self.__find_qr()
        if self.serial is None:  # didn't find a panel in this image
            return None

        if self.panel_version < 3:
            # reference_panel_pts = np.asarray(
            #     [[894, 469], [868, 232], [630, 258], [656, 496]],
            #     dtype=np.int32,
            # )
            # reference_qr_pts = np.asarray(
            #     [[898, 748], [880, 567], [701, 584], [718, 762]],
            #     dtype=np.int32
            # )
            # use the actual panel measures here - we use units of [mm]
            # the panel is 154.4 x 152.4 mm , vs. the 84 x 84 mm for the QR code
            # it is left 143.20 mm from the QR code
            # use the inner 50% square of the panel
            s = 76.2
            p = 42
            t_off = np.array([-143.2, 0])

        elif (self.panel_version >= 3) and (self.panel_version < 6):
            s = 50
            p = 45
            t_off = np.array([-145.8, 0])
            # reference_panel_pts = np.asarray(
            #     [[557, 350], [550, 480], [695, 480], [700, 350]], dtype=np.int32
            # )
            # reference_qr_pts = np.asarray(
            #     [[821, 324], [819, 506], [996, 509], [999, 330]], dtype=np.int32
            # )
        elif self.panel_version >= 6:
            # use the actual panel measures here - we use units of [mm]
            # the panel is 100 x 100 mm , vs. the 91 x 91 mm for the QR code
            # it is down 125.94 mm from the QR code
            # use the inner 50% square of the panel
            p = 41
            s = 50
            t_off = np.array([0, -130.84])

        reference_panel_pts = (
            np.asarray([[-s, s], [s, s], [s, -s], [-s, -s]], dtype=np.float32) * 0.5
            + t_off
        )
        reference_qr_pts = np.asarray(
            [[-p, p], [p, p], [p, -p], [-p, -p]], dtype=np.float32
        )
        bounds = []
        costs = []
        for rotation in range(0, 4):
            qr_points = np.roll(reference_qr_pts, rotation, axis=0)

            src = np.asarray([tuple(row) for row in qr_points[:]], np.float32)
            dst = np.asarray([tuple(row) for row in self.qr_corners()[:]], np.float32)

            # we determine the homography from the 4 corner points
            warp_matrix = cv2.getPerspectiveTransform(src, dst)

            # warp_matrix = cv2.getAffineTransform(src, dst)

            pts = np.asarray([reference_panel_pts], "float32")
            panel_bounds = cv2.convexHull(
                cv2.perspectiveTransform(pts, warp_matrix), clockwise=False
            )
            panel_bounds = np.squeeze(panel_bounds)  # remove nested lists

            bounds_in_image = True
            for i, point in enumerate(panel_bounds):
                if not self.__pt_in_image_bounds(point):
                    bounds_in_image = False
            if bounds_in_image:
                mean, std, _, _ = self.region_stats(
                    self.image.raw(), panel_bounds, sat_threshold=65000
                )
                bounds.append(panel_bounds.astype(np.int32))
                costs.append(std / mean)

        idx = costs.index(min(costs))
        self.__panel_bounds = bounds[idx]
        return self.__panel_bounds

    def ordered_panel_coordinates(self):
        """
        Return panel region coordinates in a predictable order. Panel region
        coordinates that are automatically detected by the camera are ordered
        differently than coordinates detected by Panel.panel_corners().

        :return: [ (ur), (ul), (ll), (lr) ] to mirror Image.panel_region
        attribute order
        """
        pc = self.panel_corners()
        pc = sorted(pc, key=lambda x: x[0])

        # get the coordinates on the "left" and "right" side of the bounding box
        left_coords = pc[:2]
        right_coords = pc[2:]

        # sort y values ascending for correct order
        left_coords = sorted(left_coords, key=lambda y: y[0])
        right_coords = sorted(right_coords, key=lambda y: y[0])

        return [
            tuple(right_coords[1]),
            tuple(left_coords[1]),
            tuple(left_coords[0]),
            tuple(right_coords[0]),
        ]

    def region_stats(self, img, region, sat_threshold=None):
        """Provide regional statistics for a image over a region
        Inputs: img is any image ndarray, region is a skimage shape
        Outputs: mean, std, count, and saturated count tuple for the region"""
        rev_panel_pts = np.fliplr(region)  # skimage and opencv coords are reversed
        w, h = img.shape
        mask = measure.grid_points_in_poly((w, h), rev_panel_pts)
        num_pixels = mask.sum()
        panel_pixels = img[mask]
        stdev = panel_pixels.std()
        mean_value = panel_pixels.mean()
        saturated_count = 0
        if sat_threshold is not None:
            saturated_count = (panel_pixels > sat_threshold).sum()
            # set saturated pixels here
            if num_pixels > 0:
                self.saturated_panel_pixels_pct = (100.0 * saturated_count) / num_pixels
        return mean_value, stdev, num_pixels, saturated_count

    def raw(self):
        raw_img = self.image.undistorted(self.image.raw())
        return self.region_stats(raw_img, self.panel_corners(), sat_threshold=65000)

    def intensity(self):
        intensity_img = self.image.undistorted(self.image.intensity())
        return self.region_stats(intensity_img, self.panel_corners(), sat_threshold=65000)

    def radiance(self):
        radiance_img = self.image.undistorted(self.image.radiance())
        return self.region_stats(radiance_img, self.panel_corners())

    def reflectance_mean(self):
        reflectance_image = self.image.reflectance()
        if reflectance_image is None:
            print(
                "First calculate the reflectance image by providing a\n"
                " band specific irradiance to the calling image.reflectance(irradiance)"
            )
        mean, _, _, _ = self.region_stats(reflectance_image, self.panel_corners())
        return mean

    def irradiance_mean(self, reflectance):
        radiance_mean, _, _, _ = self.radiance()
        return radiance_mean * math.pi / reflectance

    def plot_image(self):
        display_img = cv2.cvtColor(self.gray8b, cv2.COLOR_GRAY2RGB)
        if self.panel_detected():
            if self.qr_corners() is not None:
                cv2.drawContours(display_img, [self.qr_corners()], 0, (255, 0, 0), 3)
            cv2.drawContours(display_img, [self.panel_corners()], 0, (0, 0, 255), 3)

        font = cv2.FONT_HERSHEY_DUPLEX
        if self.panel_detected():
            if self.qr_corners() is not None:
                xloc = self.qr_corners()[0][0] - 100
                yloc = self.qr_corners()[0][1] + 100
            else:
                xloc = self.panel_corners()[0][0] - 100
                yloc = self.panel_corners()[0][1] + 100
            cv2.putText(
                display_img,
                str(self.serial).split("_")[0],
                (xloc, yloc),
                font,
                1,
                255,
                2,
            )
        return display_img

    def plot(self, figsize=(14, 14)):
        display_img = self.plot_image()
        fig, ax = plt.subplots(figsize=figsize)
        ax.imshow(display_img)
        plt.tight_layout()
        plt.show()
        return fig, ax

max_stars_repo_name: Oceancolour-RG/imageprocessing

max_issues_repo_name: Oceancolour-RG/imageprocessing

max_forks_repo_name: Oceancolour-RG/imageprocessing