text:
import os
import time
import numpy as np
import itertools
import scipy.misc
import matplotlib.pyplot as plt
import imageio

from zhiqiang.envs import AbstractEnv

def int_round(value):
    return int(round(value, 0))

class Pong(AbstractEnv):
    """
    """
    def __init__(self, settings):
        """
        """
        self.settings = settings
        self.sizeH = self.settings.env_settings["size_h"]
        self.sizeW = self.settings.env_settings["size_w"]
        self.ball_v = self.settings.env_settings["ball_v"]
        self.max_v = self.settings.env_settings["max_v"]
        #
        self.num_players = 1
        self.reward_hit = 1
        self.reward_loss = -1
        #
        self.reset()
        #

    def reset(self):
        """
        """
        self._reset_env()
        return self._reset_game()

    #
    def _reset_env(self):
        """
        """
        self.list_actions = list(range(-self.max_v, self.max_v+1))
        self.num_actions = len(self.list_actions)
        #
        self.racket_h = 3
        self.racket_w = 1
        self.ball_size = 1
        #
        # board
        self.size_all_h = self.sizeH + 2
        self.size_all_w = self.sizeW + 2
        self.pa_h = self.size_all_h - 1   # all
        self.pa_w = self.size_all_w - 1
        self.pd_h = self.size_all_h - 2   # board
        self.pd_w = self.size_all_w - 2
        #
        self.size_center_h = self.size_all_h // 2
        self.size_center_w = self.size_all_w // 2
        #
        self.flag_reset_game = 1
        #
        self.pic_h = self.size_all_h
        self.pic_w = self.size_all_w
        #

    def _reset_game(self):
        """
        """
        # start, position, player, ball
        self.pp_0 = self.size_center_h    # player
        self.pp_1 = self.size_center_h
        self.pb_h = self.size_center_h    # ball
        self.pb_w = self.size_center_w
        #
        self.mb_h = np.random.randint(-10, 10)  # movement, ball
        self.mb_w = np.random.randint(2, 10)
        self._normalize_mb()
        #
        self.action_player_0 = 0
        self.action_player_1 = 0
        #
        self.flag_reset_game = 0
        #
        # state
        self.state_prev = self.render()
        self.state_curr = self.render()
        self.state = (self.state_prev, self.state_curr)
        return self.state
        #

    def render(self):
        """
        """
        a = np.ones([self.size_all_h, self.size_all_w, 3])
        a[1:-1,1:-1, 0] = 0
        a[1:-1,1:-1, 1] = 1
        a[1:-1,1:-1, 2] = 1
        #
        if self.num_players == 1:
            pass
        else:
            a[self.pp_0-1:self.pp_0+2, 0, 0] = 1
            a[self.pp_0-1:self.pp_0+2, 0, 1] = 0
            a[self.pp_0-1:self.pp_0+2, 0, 2] = 0
        #
        a[self.pp_1-1:self.pp_1+2, self.pa_w, 0] = 0
        a[self.pp_1-1:self.pp_1+2, self.pa_w, 1] = 0
        a[self.pp_1-1:self.pp_1+2, self.pa_w, 2] = 1
        #
        a[self.pb_h, self.pb_w, :] = 0
        #
        return a
        #

    def step(self, action, action_0=0):
        """
        """
        self.action_player_0 = self.list_actions[ action_0 ]
        self.action_player_1 = self.list_actions[ action ]
        #
        if self.flag_reset_game:
            return self._reset_game(), self._score(0.0, 0.0), False, {"reset": 1}
        #
        # ball
        d_h = self.mb_h
        d_w = self.mb_w
        #
        # check, reward
        #
        ## first
        flag_reset_game, reward_0, reward_1 = self._check_bound_w(d_h, d_w)
        #
        # reset
        if flag_reset_game:
            #
            self.flag_reset_game = 1
            #
            self.state_prev = self.state_curr
            self.state_curr = self.render()
            self.state = (self.state_prev, self.state_curr)
            #
            return self.state, self._score(reward_0, reward_1), False, None
            #
        #
        self._check_bound_h(d_h, d_w)
        #

        #
        # player
        self._execute_player_action(0)
        self._execute_player_action(1)
        #
        # ball_v
        self._normalize_mb()
        #

        #
        self.state_prev = self.state_curr
        self.state_curr = self.render()
        self.state = (self.state_prev, self.state_curr)
        #
        return self.state, self._score(reward_0, reward_1), False, None
        #

    def close(self):
        pass
    #

    #
    def map_to_pic(self, state_all):
        """
        """
        state = state_all[1]
        #
        h = self.pic_h
        w = self.pic_w
        b = scipy.misc.imresize(state[:,:,0],[h,w,1],interp='nearest')
        c = scipy.misc.imresize(state[:,:,1],[h,w,1],interp='nearest')
        d = scipy.misc.imresize(state[:,:,2],[h,w,1],interp='nearest')
        pic = np.stack([b,c,d], axis=2)
        #
        return pic

    def display(self, state=None, show=True, step=None, score=None):
        """
        """
        if state is None:
            state = self.state
        #
        # grid = plt.GridSpec(1, 5, wspace=0.5, hspace=0.5)
        # plt.subplot(grid[0, 0:3])
        pic = self.map_to_pic(state)
        plt.imshow(pic, interpolation="nearest")
        # plt.subplot(grid[0, 4])
        if step is not None:
            if score is not None:
                plt.title("step: %d, score: %f" % (step, score))
            else:
                plt.title("step: %d" % (step, ))
        else:
            if score is not None:
                plt.title("score: %f" % (score, ))
            #
        #
        if show:
            plt.show()
        #
        filename = "pong_temp.eps"  # eps < png < jpg
        plt.savefig(filename)
        pic_with_score = imageio.imread(filename)
        os.remove(filename)
        #
        return pic, pic_with_score
        #
    #

    #
    def _execute_player_action(self, player_idx):
        """
        """
        if player_idx == 0:
            self.pp_0 += self.action_player_0
            #
            if self.pp_0 < 1:
                self.pp_0 = 1
            elif self.pp_0 > self.pd_h:   # board
                self.pp_0 = self.pd_h
            #
        elif player_idx == 1:
            self.pp_1 += self.action_player_1
            #
            if self.pp_1 < 1:
                self.pp_1 = 1
            elif self.pp_1 > self.pd_h:
                self.pp_1 = self.pd_h
            #
    #
    def _normalize_mb(self):
        """
        """
        norm = np.sqrt(self.mb_h ** 2 + self.mb_w ** 2)
        ratio = self.ball_v / norm
        self.mb_h = int_round(self.mb_h * ratio)
        self.mb_w = int_round(self.mb_w * ratio)
    #
    def _check_bound_w(self, d_h, d_w):
        """
        """
        # check, reward
        flag_reset_game = 0
        reward_0 = 0
        reward_1 = 0
        #
        if self.pb_w + d_w <= 0:
            #
            if self.num_players == 1:
                # reflection
                self.pb_w += d_w
                self.pb_w = 1 - self.pb_w
                self.mb_w = -self.mb_w
                #
            else:
                #
                # when pb_w = 1,
                d_pbw = 1 - self.pb_w
                d_pbh = int_round(d_h * d_pbw / d_w)
                d_p0 = int_round(self.action_player_0 * d_pbw / d_w)
                #
                pb_h = self.pb_h + d_pbh
                p0 = self.pp_0 + d_p0
                #
                if pb_h <= 0:
                    pb_h = 1 - pb_h
                elif pb_h >= self.pa_h:
                    pb_h = 2 * self.pd_h - pb_h - 1
                #
                if p0 <= 0:
                    p0 = 1
                elif p0 >= self.pa_h:
                    p0 = self.pd_h
                #
                if pb_h >= p0-1 and pb_h <= p0+1:
                    # reflection
                    self.pb_w += d_w
                    self.pb_w = 1 - self.pb_w
                    self.mb_w = -self.mb_w
                    #
                    reward_0 = self.reward_hit
                    #
                else:
                    reward_0 = self.reward_loss
                    #
                    flag_reset_game = 1
                    self.pb_w = 0
                    self.pb_h = pb_h
                    self.pp_0 = p0
                    self._execute_player_action(1)
                    #
            #
        elif self.pb_w + d_w >= self.pa_w:
            #
            # when pb_w = self.pd_w,
            d_pbw = self.pd_w - self.pb_w
            d_pbh = int_round(d_h * d_pbw / d_w)
            d_p1 = int_round(self.action_player_1 * d_pbw / d_w)
            #
            pb_h = self.pb_h + d_pbh
            p1 = self.pp_1 + d_p1
            #
            if pb_h <= 0:
                pb_h = 1 - pb_h
            elif pb_h >= self.pa_h:
                pb_h = 2 * self.pd_h - pb_h - 1
            #
            if p1 <= 0:
                p1 = 1
            elif p1 >= self.pa_h:
                p1 = self.pd_h
            #
            if pb_h >= p1-1 and pb_h <= p1+1:
                # reflection
                self.pb_w += d_w
                self.pb_w = 2 * self.pd_w - self.pb_w + 1
                self.mb_w = -self.mb_w
                #
                reward_1 = self.reward_hit
                #
            else:
                reward_1 = self.reward_loss
                #
                flag_reset_game = 1
                self.pb_w = self.pa_w
                self.pb_h = pb_h
                self.pp_1 = p1
                self._execute_player_action(0)
                #
            #
        else:
            self.pb_w += d_w
            #
        #
        return flag_reset_game, reward_0, reward_1
        #
    #
    def _check_bound_h(self, d_h, d_w):
        """
        """
        # check, reflection, h
        self.pb_h += d_h
        #
        if self.pb_h <= 0:
            self.pb_h = 1 - self.pb_h
            self.mb_h = -self.mb_h -1
        elif self.pb_h >= self.pa_h:
            self.pb_h = 2 * self.pd_h - self.pb_h - 1
            self.mb_h = -self.mb_h +1
        #
    #
    def _score(self, reward_0, reward_1):
        """
        """
        if self.num_players == 1:
            return reward_1
        else:
            return np.array( [reward_0, reward_1] )
        #

#

max_stars_repo_name: Li-Ming-Fan/zhiqiang

max_issues_repo_name: Li-Ming-Fan/zhiqiang

max_forks_repo_name: Li-Ming-Fan/zhiqiang


text:
import numpy
import warnings
from skimage.transform import AffineTransform, EuclideanTransform, SimilarityTransform

from . import defaults
from .exceptions import CoordIOUserWarning

from . import libcoordio

# low level, functional conversions

def fieldAngle2Cart(xField, yField):
    """Convert (ZEMAX-style) field angles in degrees to a cartesian point on
    the unit sphere

    Zemax defines field angles like this:
    Positive field angles imply positive slope for the ray in that direction,
    and thus refer to **negative** coordinates on distant objects. ZEMAX converts
    x field angles ( αx ) and y field angles ( αy ) to ray direction cosines
    using the following formulas:

    tanαx = l/n
    tanαy = m/n
    l^2 + m^2 + n^2 = 1

    where l, m, and n are the x, y, and z direction cosines.

    Parameters
    -----------
    xField: scalar or 1D array
        zemax style x field
    yField: scalar or 1D array
        zemax style y field

    Returns
    ---------
    result: list
        [x,y,z] coordinates on unit sphere
    """
    # xField, yField = fieldXY

    # invert field degrees, fields represent slope
    # of ray, so a positivly sloped ray is coming from
    # below its respective plane, the resulting vector
    # from this transform should thus point downward, not up
    tanx = numpy.tan(numpy.radians(-1*xField))
    tany = numpy.tan(numpy.radians(-1*yField))

    # calculate direction cosines
    z = numpy.sqrt(1/(tanx**2 + tany**2 + 1))
    x = tanx * z
    y = tany * z
    # if numpy.isnan(n) or numpy.isnan(l) or numpy.isnan(m):
    #     raise RuntimeError("NaN output [%.2f, %.2f, %.2f] for input [%.2f, %.2f]"%(n, l, m, xField, yField))

    return [x,y,z]

def cart2FieldAngle(x, y, z):
    """Convert cartesian point on unit sphere of sky to (ZEMAX-style) field
    angles in degrees.

    Zemax defines field angles like this:
    Positive field angles imply positive slope for the ray in that direction,
    and thus refer to **negative** coordinates on distant objects. ZEMAX converts
    x field angles ( αx ) and y field angles ( αy ) to ray direction cosines
    using the following formulas:

    tanαx = l/n
    tanαy = m/n
    l^2 + m^2 + n^2 = 1

    where l, m, and n are the x, y, and z direction cosines.

    Parameters
    -----------
    x: scalar or 1D array
    y: scalar or 1D array
    z: scalar or 1D array

    Returns
    ---------
    result: list
        [xField, yField] ZEMAX style field coordinates (degrees)
    """
    # if numpy.abs(numpy.linalg.norm(cartXYZ) - 1) > SMALL_NUM:
    #     raise RuntimeError("cartXYZ must be a vector on unit sphere with L2 norm = 1")
    # l, m, n = x, y, z

    # invert field degrees, fields represent slope
    # of ray, so a positivly sloped ray is coming from
    # below the optical axis, the resulting vector
    # from this transform should thus point downward, not up

    xField = -1 * numpy.degrees(numpy.arctan2(x, z))
    yField = -1 * numpy.degrees(numpy.arctan2(y, z))
    # if numpy.isnan(xField) or numpy.isnan(yField):
    #     raise RuntimeError("NaN output [%.2f, %.2f] for input [%.2f, %.2f, %.2f]"%(xField, yField, cartXYZ[0], cartXYZ[1], cartXYZ[2]))
    return [xField, yField]

def cart2Sph(x, y, z):
    """Convert cartesian coordinates to spherical
    coordinates theta, phi in degrees

    phi is polar angle measure from z axis
    theta is azimuthal angle measured from x axis

    Parameters
    -----------
    x: scalar or 1D array
    y: scalar or 1D array
    z: scalar or 1D array

    Returns
    ---------
    result: list
        [theta, phi] degrees
    """

    # if numpy.abs(numpy.linalg.norm(cartXYZ) - 1) > SMALL_NUM:
    #     raise RuntimeError("cartXYZ must be a vector on unit sphere with L2 norm = 1")
    # if cartXYZ[2] < 0:
    #     raise RuntimeError("z direction must be positive for cartesian field coord")
    # x, y, z = cartXYZ
    theta = numpy.degrees(numpy.arctan2(y, x))
    # wrap theta to be between 0 and 360 degrees
    theta = theta % 360

    phi = numpy.degrees(numpy.arccos(z))
    # if numpy.isnan(theta) or numpy.isnan(phi):
    #     raise RuntimeError("NaN output [%.2f, %.2f] from input [%.2f, %.2f, %.2f]"%(theta, phi, cartXYZ[0], cartXYZ[1], cartXYZ[2]))
    return [theta, phi]

def sph2Cart(theta, phi, r=1):
    """Convert spherical coordinates theta, phi in degrees
    to cartesian coordinates on unit sphere.

    phi is polar angle measure from z axis
    theta is azimuthal angle measured from x axis

    Parameters
    -----------
    theta: scalar or 1D array
        degrees, azimuthal angle
    phi: scalar or 1D array
        degrees, polar angle
    r: scalar
        radius of curvature. Default to 1 for unit sphere

    Returns
    ---------
    result: list
        [x,y,z] coordinates on sphere

    """

    theta, phi = numpy.radians(theta), numpy.radians(phi)
    x = r*numpy.cos(theta) * numpy.sin(phi)
    y = r*numpy.sin(theta) * numpy.sin(phi)
    z = r*numpy.cos(phi)

    # if numpy.isnan(x) or numpy.isnan(y) or numpy.isnan(z):
    #     raise RuntimeError("NaN output [%.2f, %.2f, %.2f] for input [%.2f, %.2f]"%(x, y, z, numpy.degrees(theta), numpy.degrees(phi)))

    return [x, y, z]

def fieldToObserved(x, y, z, altCenter, azCenter, pa):
    """Convert xyz unit-spherical field coordinates to Alt/Az

    the inverse of observetToField

    Az=0 is north, Az=90 is east

    +x points along +RA, +y points along +Dec, +z points along boresight
    from telescope toward sky

    Parameters
    ------------
    x: scalar or 1D array
        unit-spherical x field coord (aligned with +RA)
    y: scalar or 1D array
        unit-spherical y field coord (aligned with +Dec)
    z: scalar or 1D array
        unit-spherical z coord
    altCenter: float
        altitude in degrees, field center
    azCenter: float
        azimuth in degrees, field center
    pa: float
        position angle of field center

    Returns
    --------
    alt : float or 1D array
        altitude of coordinate in degrees
    az : float or 1D array

    """
    cosQ = numpy.cos(numpy.radians(-1*pa))
    sinQ = numpy.sin(numpy.radians(-1*pa))
    rotQ = numpy.array([
        [ cosQ, sinQ, 0],
        [-sinQ, cosQ, 0],
        [    0,    0, 1]
    ])

    coords = numpy.array(
        [x, y, z]
    ).T

    coords = rotQ.dot(coords.T).T

    sinPhi = numpy.sin(-1*numpy.radians(90-altCenter))
    cosPhi = numpy.cos(-1*numpy.radians(90-altCenter))
    rotPhi = numpy.array([
        [1,       0,      0],
        [0,  cosPhi, sinPhi],
        [0, -sinPhi, cosPhi]
    ])
    coords = rotPhi.dot(coords.T).T

    sinTheta = numpy.sin(-1*numpy.radians(90-azCenter))
    cosTheta = numpy.cos(-1*numpy.radians(90-azCenter))
    rotTheta = numpy.array([
        [ cosTheta, sinTheta, 0],
        [-sinTheta, cosTheta, 0],
        [        0,        0, 1]
    ])

    coords = rotTheta.dot(coords.T).T

    x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]
    theta, phi = cart2Sph(x, y, z)

    # convert sph theta, phi to az, alt
    az = -1 * theta
    alt = 90 - phi

    az = az % 360  # wrap to 0,360

    return alt, az

def observedToField(alt, az, altCenter, azCenter, pa):
    """Convert Az/Alt coordinates to cartesian coords on the
    unit sphere with the z axis aligned with field center (boresight)

    Az=0 is north, Az=90 is east

    Resulting Cartesian coordinates have the following
    convention +x points along +RA, +y points along +Dec

    Parameters
    ------------
    alt: float or 1D array
        altitude in degrees, positive above horizon, negative below
    az: float or 1D array
        azimuth in degrees az=0 is north az=90 is east
    altCenter: float
        altitude in degrees, field center
    azCenter: float
        azimuth in degrees, field center
    pa: float
        position angle of field center

    Returns
    --------
    theta : float or 1D array
        azimuthal angle of field coordinate (deg)
    phi : float or 1D array
        polar angle of field coordinat (deg off axis)
    """
    # convert from alt/az to spherical
    phis = 90 - alt  # alt
    thetas = -1 * az  # az

    # work in cartesian frame
    coords = sph2Cart(thetas, phis)
    coords = numpy.array(coords).T

    # rotate the xyz coordinate system about z axis
    # such that -y axis is aligned with the azimuthal angle
    # of the field center

    sinTheta = numpy.sin(numpy.radians(90 - azCenter))
    cosTheta = numpy.cos(numpy.radians(90 - azCenter))
    rotTheta = numpy.array([
        [ cosTheta, sinTheta, 0],
        [-sinTheta, cosTheta, 0],
        [        0,        0, 1]
    ])

    coords = rotTheta.dot(coords.T).T

    # rotate the xyz coordinate system about the x axis
    # such that +z points to the field center.

    sinPhi = numpy.sin(numpy.radians(90 - altCenter))
    cosPhi = numpy.cos(numpy.radians(90 - altCenter))
    rotPhi = numpy.array([
        [1,       0,      0],
        [0,  cosPhi, sinPhi],
        [0, -sinPhi, cosPhi]
    ])
    coords = rotPhi.dot(coords.T).T
    # return coords

    # finally rotate about z by the parallactic angle
    # this puts +RA along +X and +DEC along +Y
    cosQ = numpy.cos(numpy.radians(pa))
    sinQ = numpy.sin(numpy.radians(pa))
    rotQ = numpy.array([
        [ cosQ, sinQ, 0],
        [-sinQ, cosQ, 0],
        [    0,    0, 1]
    ])

    coords = rotQ.dot(coords.T).T

    x = coords[:, 0]
    y = coords[:, 1]
    z = coords[:, 2]

    # finally convert back from cartesian to spherical (Field)
    thetaPhi = cart2Sph(x, y, z)
    thetaPhi = numpy.array(thetaPhi).T
    theta, phi = thetaPhi[:,0], thetaPhi[:,1]
    return theta, phi

def fieldToFocal(thetaField, phiField, site, waveCat):
    """Convert spherical field coordinates to the modeled spherical
    position on the focal plane.  Origin is the M1 vertex.

    Raises warnings for coordinates at large field angles.

    Parameters
    -----------
    thetaField: float or numpy.ndarray
        azimuthal angle of field coordinate (deg)
    phiField: float or numpy.ndarray
        polar angle of field coordinate (deg)
    site: str
        "APO" or "LCO"
    waveCat: str
        "Boss", "Apogee", or "GFA"

    Result
    -------
    xFocal: float or numpy.ndarray
        spherical x focal coord mm
    yFocal: float or numpy.ndarray
       spherical y focal coord mm
    zFocal: float or numpy.ndarray
        spherical z focal coord mm
    R: float
        radius of curvature for sphere
    b: float
        z location of center of sphere
    fieldWarn: bool or boolean array
        True if angle off-axis is large enough to be suspicious
    """

    if site == "APO":
        fieldWarn = numpy.array(phiField) > defaults.APO_MAX_FIELD_R
    if site == "LCO":
        fieldWarn = numpy.array(phiField) > defaults.LCO_MAX_FIELD_R

    if hasattr(fieldWarn, "__len__"):
        if True in fieldWarn:
            warnings.warn(
                "Warning! Coordinate far off telescope optical axis "\
                "conversion may be bogus",
                CoordIOUserWarning
            )
    elif fieldWarn is True:
        warnings.warn(
            "Warning! Coordinate far off telescope optical axis, "\
            "conversion may be bogus",
            CoordIOUserWarning
        )

    direction = "focal"
    R, b, c0, c1, c2, c3, c4 = defaults.getFPModelParams(site, direction, waveCat)

    phiFocal = c0 * phiField + c1 * phiField**3 + c2 * phiField**5 \
        + c3 * phiField**7 + c4 * phiField**9
    # thetaField = thetaFocal by definition
    # convert back to xyz coords
    rTheta = numpy.radians(thetaField)
    rPhi = numpy.radians(180 - phiFocal)
    xFocal = R * numpy.cos(rTheta) * numpy.sin(rPhi)
    yFocal = R * numpy.sin(rTheta) * numpy.sin(rPhi)
    zFocal = R * numpy.cos(rPhi) + b

    return xFocal, yFocal, zFocal, R, b, fieldWarn

def focalToField(xFocal, yFocal, zFocal, site, waveCat):
    """Convert xyz focal position to unit-spherical field position.
    xyzFocal do not need to lie on the modeled sphere, and the spherical
    radius of curvature is not needed, only the origin of the sphere
    is needed (b).

    Raises warnings for coordinates at large field angles.

    Parameters
    -----------
    xFocal: float or 1D array
        x focal coord mm
    yFocal: float or 1D array
        y focal coord mm
    zFocal: float or 1D array
        z focal coord mm, +z points along boresight toward sky.
    site: str
        "APO" or "LCO"
    waveCat: str
        "Boss", "Apogee", or "GFA"

    Result
    -------
    thetaField: float or 1D array
        azimuthal field coordinate (deg)
    phiField: float or 1D array
       polar field coordinate (deg)
    fieldWarn: bool or boolean array
        True if angle off-axis is large enough to be suspicious
    """
    direction = "field"
    R, b, c0, c1, c2, c3, c4 = defaults.getFPModelParams(site, direction, waveCat)
    # note by definition thetaField==thetaFocal
    thetaField = numpy.degrees(numpy.arctan2(yFocal, xFocal))

    # generate focal phis (degree off boresight)
    # angle off-axis from optical axis
    rFocal = numpy.sqrt(xFocal**2 + yFocal**2)
    v = numpy.array([rFocal, zFocal]).T
    v[:, 1] = v[:, 1] - b

    # unit vector pointing toward object on focal plane from circle center
    # arc from vertex to off axis
    v = v / numpy.vstack([numpy.linalg.norm(v, axis=1)] * 2).T
    # FP lands at -Z so, Angle from sphere center towards ground
    downwardZaxis = numpy.array([0, -1])
    # phi angle between ray and optical axis measured from sphere center
    phiFocal = numpy.degrees(numpy.arccos(v.dot(downwardZaxis)))
    phiField = c0 * phiFocal + c1 * phiFocal**3 + c2 * phiFocal**5 \
        + c3 * phiFocal**7 + c4 * phiFocal**9

    if site == "APO":
        fieldWarn = numpy.array(phiField) > defaults.APO_MAX_FIELD_R
    if site == "LCO":
        fieldWarn = numpy.array(phiField) > defaults.LCO_MAX_FIELD_R

    if hasattr(fieldWarn, "__len__"):
        if True in fieldWarn:
            warnings.warn(
                "Warning! Coordinate far off telescope optical axis, conversion may be bogus",
                CoordIOUserWarning
            )
    elif fieldWarn is True:
        warnings.warn(
            "Warning! Coordinate far off telescope optical axis, conversion may be bogus",
            CoordIOUserWarning
        )

    return thetaField, phiField, fieldWarn

def focalToWok(
    xFocal, yFocal, zFocal, positionAngle=0,
    xOffset=0, yOffset=0, zOffset=0, tiltX=0, tiltY=0
):
    """Convert xyz focal coordinates in mm to xyz wok coordinates in mm.

    The origin of the focal coordinate system is the
    M1 vertex. focal +y points toward North, +x points toward E.
    The origin of the wok coordinate system is the wok vertex.  -x points
    toward the boss slithead.  +z points from telescope to sky.

    Tilt is applied about x axis then y axis.

    Parameters
    -------------
    xFocal: scalar or 1D array
        x position of object on focal plane mm
        (+x aligned with +RA on image)
    yFocal: scalar or 1D array
        y position of object on focal plane mm
        (+y aligned with +Dec on image)
    zFocal: scalar or 1D array
        z position of object on focal plane mm
        (+z aligned boresight and increases from the telescope to the sky)
    positionAngle: scalar
        position angle deg.  Angle measured from (image) North through East to wok +y.
        So position angle of 45 deg, wok +y points NE
    site : str
        either "APO" or "LCO"
    xOffset: scalar or None
        x position (mm) of wok origin (vertex) in focal coords
        calibrated
    yOffset: scalar
        y position (mm) of wok origin (vertex) in focal coords
        calibratated
    zOffset: scalar
        z position (mm) of wok origin (vertex) in focal coords
        calibratated
    tiltX: scalar
        tilt (deg) of wok about focal x axis at PA=0
        calibrated
    tiltY: scalar
        tilt (deg) of wok about focal y axis at PA=0
        calibrated

    Returns
    ---------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    """

    # apply calibrated tilts and translation (at PA=0)
    # where they should be fit
    # tilts are defined https://mathworld.wolfram.com/RotationMatrix.html
    # as coordinate system rotations counter clockwise when looking
    # down the positive axis toward the origin
    coords = numpy.array([xFocal, yFocal, zFocal])
    transXYZ = numpy.array([xOffset, yOffset, zOffset])

    rotX = numpy.radians(tiltX)
    rotY = numpy.radians(tiltY)
    # rotation about z axis is position angle
    # position angle is clockwise positive for rotation measured from
    # north to wok +y (when looking from above the wok)
    # however rotation matrices are positinve for counter-clockwise rotation
    # hence the sign flip that's coming
    rotZ = -1*numpy.radians(positionAngle)

    rotMatX = numpy.array([
        [1, 0, 0],
        [0, numpy.cos(rotX), numpy.sin(rotX)],
        [0, -1*numpy.sin(rotX), numpy.cos(rotX)]
    ])

    rotMatY = numpy.array([
        [numpy.cos(rotY), 0, -1*numpy.sin(rotY)],
        [0, 1, 0],
        [numpy.sin(rotY), 0, numpy.cos(rotY)]
    ])

    # rotates coordinate system
    rotMatZ = numpy.array([
        [numpy.cos(rotZ), numpy.sin(rotZ), 0],
        [-numpy.sin(rotZ), numpy.cos(rotZ), 0],
        [0, 0, 1]
    ])

    # first apply rotation about x axis (x tilt)
    coords = rotMatX.dot(coords)
    # next apply rotation about y axis (y tilt)
    coords = rotMatY.dot(coords)
    # apply rotation about z axis (angle of observation)
    coords = rotMatZ.dot(coords)

    # apply translation
    if hasattr(xFocal, "__len__"):
        # list of coords fed in
        transXYZ = numpy.array([transXYZ]*len(xFocal)).T
        xWok, yWok, zWok = coords - transXYZ
    else:
        # single set of xyz coords fed in
        xWok, yWok, zWok = coords - transXYZ

    return xWok, yWok, zWok

def wokToFocal(
    xWok, yWok, zWok, positionAngle=0,
    xOffset=0, yOffset=0, zOffset=0, tiltX=0, tiltY=0
):
    """Convert xyz wok coordinates in mm to xyz focal coordinates in mm.

    The origin of the focal coordinate system is the
    M1 vertex. focal +y points toward North, +x points toward E.
    The origin of the wok coordinate system is the wok vertex.  -x points
    toward the boss slithead.  +z points from telescope to sky.

    Tilt is applied first about x axis then y axis.

    Parameters
    -------------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image at PA=0)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image at PA=0)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    positionAngle: scalar
        position angle deg.  Angle measured from (image) North through East to wok +y.
        So position angle of 45 deg, wok +y points NE
    xOffset: scalar or None
        x position (mm) of wok origin (vertex) in focal coords
        calibrated
    yOffset: scalar
        y position (mm) of wok origin (vertex) in focal coords
        calibratated
    zOffset: scalar
        z position (mm) of wok origin (vertex) in focal coords
        calibratated
    tiltX: scalar
        tilt (deg) of wok about focal x axis at PA=0
        calibrated
    tiltY: scalar
        tilt (deg) of wok about focal y axis at PA=0
        calibrated

    Returns
    ---------
    xFocal: scalar or 1D array
        x position of object in focal coord sys mm
        (+x aligned with +RA on image)
    yFocal: scalar or 1D array
        y position of object in focal coord sys mm
        (+y aligned with +Dec on image)
    zFocal: scalar or 1D array
        z position of object in focal coord sys mm
        (+z aligned boresight and increases from the telescope to the sky)
    """
    # this routine is a reversal of the steps
    # in the function focalToWok, with rotational
    # angles inverted and translations applied in reverse
    coords = numpy.array([xWok, yWok, zWok])

    rotX = numpy.radians(-1*tiltX)
    rotY = numpy.radians(-1*tiltY)
    rotZ = numpy.radians(positionAngle)

    rotMatX = numpy.array([
        [1, 0, 0],
        [0, numpy.cos(rotX), numpy.sin(rotX)],
        [0, -numpy.sin(rotX), numpy.cos(rotX)]
    ])

    rotMatY = numpy.array([
        [numpy.cos(rotY), 0, -numpy.sin(rotY)],
        [0, 1, 0],
        [numpy.sin(rotY), 0, numpy.cos(rotY)]
    ])

    rotMatZ = numpy.array([
        [numpy.cos(rotZ), numpy.sin(rotZ), 0],
        [-numpy.sin(rotZ), numpy.cos(rotZ), 0],
        [0, 0, 1]
    ])

    transXYZ = numpy.array([xOffset, yOffset, zOffset])

    # add offsets for reverse transform
    if hasattr(xWok, "__len__"):
        # list of coords fed in
        transXYZ = numpy.array([transXYZ]*len(xWok)).T
        coords = coords + transXYZ
    else:
        # single set of xyz coords fed in
        coords = coords + transXYZ

    coords = rotMatZ.dot(coords)
    coords = rotMatY.dot(coords)
    xFocal, yFocal, zFocal = rotMatX.dot(coords)
    return xFocal, yFocal, zFocal

def _verify3Vector(checkMe, label):
    if not hasattr(checkMe, "__len__"):
        raise RuntimeError("%s must be a 3-vector"%label)
    checkMe = numpy.array(checkMe).squeeze()
    if len(checkMe) != 3:
        raise RuntimeError("%s must be a 3-vector"%label)
    if checkMe.ndim != 1:
        raise RuntimeError("%s must be a 3-vector"%label)
    return checkMe

def wokToTangent(xWok, yWok, zWok, b, iHat, jHat, kHat,
                 elementHeight=defaults.POSITIONER_HEIGHT, scaleFac=1,
                 dx=0, dy=0, dz=0):
    """
    Convert from wok coordinates to tangent coordinates.

    xyz Wok coords are mm with orgin at wok vertex. +z points from wok toward M2.
    -x points toward the boss slithead

    In the tangent coordinate frame the xy plane is tangent to the wok
    surface at xyz position b.  The origin is set to be elementHeight above
    the wok surface.

    Parameters
    -------------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image at PA = 0)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image at PA = 0)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    b: 3-vector
        x,y,z position (mm) of element hole on wok surface measured in wok coords
    iHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate x axis
    jHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate y axis
    kHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate z axis
    elementHeight: scalar
        height (mm) of positioner/GFA chip above wok surface
    scaleFac: scalar
        scale factor to apply to b to account for thermal expansion of wok.
        scale is applied to b radially
    dx: scalar
        x offset (mm), calibration to capture small displacements of
        tangent x
    dy: scalar
        y offset (mm), calibration to capture small displacements of
        tangent y
    dz: scalar
        z offset (mm), calibration to capture small displacements of
        tangent x

    Returns
    ---------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    zTangent: scalar or 1D array
        z position (mm) in tangent coordinates
    """
    b = _verify3Vector(b, "b")
    iHat = _verify3Vector(iHat, "iHat")
    jHat = _verify3Vector(jHat, "jHat")
    kHat = _verify3Vector(kHat, "kHat")
    # format into list of lists
    if hasattr(xWok, "__len__"):
        xyzWok = numpy.array([xWok, yWok, zWok]).T.tolist()

        xyzTangent = libcoordio.wokToTangentArr(
            xyzWok, list(b), list(iHat), list(jHat), list(kHat),
            elementHeight, scaleFac, dx, dy, dz)

        xyzTangent = numpy.array(xyzTangent)
        xTangent = xyzTangent[:, 0]
        yTangent = xyzTangent[:, 1]
        zTangent = xyzTangent[:, 2]
    else:
        xTangent, yTangent, zTangent = libcoordio.wokToTangent(
            [xWok, yWok, zWok], list(b), list(iHat), list(jHat), list(kHat),
            elementHeight, scaleFac, dx, dy, dz
        )
    return xTangent, yTangent, zTangent

def _wokToTangent(xWok, yWok, zWok, b, iHat, jHat, kHat,
                 elementHeight=defaults.POSITIONER_HEIGHT, scaleFac=1,
                 dx=0, dy=0, dz=0):
    """
    Convert from wok coordinates to tangent coordinates.

    xyz Wok coords are mm with orgin at wok vertex. +z points from wok toward M2.
    -x points toward the boss slithead

    In the tangent coordinate frame the xy plane is tangent to the wok
    surface at xyz position b.  The origin is set to be elementHeight above
    the wok surface.

    Parameters
    -------------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image at PA = 0)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image at PA = 0)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    b: 3-vector
        x,y,z position (mm) of element hole on wok surface measured in wok coords
    iHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate x axis
    jHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate y axis
    kHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate z axis
    elementHeight: scalar
        height (mm) of positioner/GFA chip above wok surface
    scaleFac: scalar
        scale factor to apply to b to account for thermal expansion of wok.
        scale is applied to b radially
    dx: scalar
        x offset (mm), calibration to capture small displacements of
        tangent x
    dy: scalar
        y offset (mm), calibration to capture small displacements of
        tangent y
    dz: scalar
        z offset (mm), calibration to capture small displacements of
        tangent x

    Returns
    ---------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    zTangent: scalar or 1D array
        z position (mm) in tangent coordinates
    """
    # check for problems
    b = _verify3Vector(b, "b")
    iHat = _verify3Vector(iHat, "iHat")
    jHat = _verify3Vector(jHat, "jHat")
    kHat = _verify3Vector(kHat, "kHat")

    coords = numpy.array([xWok,yWok,zWok])

    # apply radial scale factor to b
    # assume that z is unaffected by
    # thermal expansion (probably reasonable)
    if scaleFac != 1:
        r = numpy.sqrt(b[0]**2+b[1]**2)
        theta = numpy.arctan2(b[1], b[0])
        r = r*scaleFac
        b[0] = r*numpy.cos(theta)
        b[1] = r*numpy.sin(theta)

    isArr = hasattr(xWok, "__len__")
    calibOff = numpy.array([dx,dy,dz])
    if isArr:
        b = numpy.array([b]*len(xWok)).T
        calibOff = numpy.array([calibOff]*len(xWok)).T

    # offset to wok base position
    coords = coords - b

    # rotate normal to wok surface at point b
    rotNorm = numpy.array([iHat, jHat, kHat], dtype="float64")
    coords = rotNorm.dot(coords)

    # offset xy plane to focal surface
    if isArr:
        coords[2,:] = coords[2,:] - elementHeight
    else:
        coords[2] = coords[2] - elementHeight

    # apply rotational calibration
    # if dRot != 0:
    #     rotZ = numpy.radians(dRot)
    #     rotMatZ = numpy.array([
    #         [numpy.cos(rotZ), numpy.sin(rotZ), 0],
    #         [-numpy.sin(rotZ), numpy.cos(rotZ), 0],
    #         [0, 0, 1]
    #     ])

    #     coords = rotMatZ.dot(coords)

    coords = coords - calibOff

    return coords[0], coords[1], coords[2]

def tangentToWok(xTangent, yTangent, zTangent, b, iHat, jHat, kHat,
                 elementHeight=defaults.POSITIONER_HEIGHT, scaleFac=1,
                 dx=0, dy=0, dz=0):
    """
    Convert from tangent coordinates at b to wok coordinates.

    xyz Wok coords are mm with orgin at wok vertex. +z points from wok toward M2.
    -x points toward the boss slithead

    In the tangent coordinate frame the xy plane is tangent to the wok
    surface at xyz position b.  The origin is set to be elementHeight above
    the wok surface.

    Parameters
    -------------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    zTangent: scalar or 1D array
        z position (mm) in tangent coordinates
    b: 3-vector
        x,y,z position (mm) of element hole on wok surface measured in wok coords
    iHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate x axis
    jHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate y axis
    kHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate z axis
    elementHeight: scalar
        height (mm) of positioner/GFA chip above wok surface
    scaleFac: scalar
        scale factor to apply to b to account for thermal expansion of wok.
        scale is applied to b radially
    dx: scalar
        x offset (mm), calibration to capture small displacements of
        tangent x
    dy: scalar
        y offset (mm), calibration to capture small displacements of
        tangent y
    dz: scalar
        z offset (mm), calibration to capture small displacements of
        tangent x

    Returns
    ---------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image at PA = 0)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image at PA = 0)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    """
    b = _verify3Vector(b, "b")
    iHat = _verify3Vector(iHat, "iHat")
    jHat = _verify3Vector(jHat, "jHat")
    kHat = _verify3Vector(kHat, "kHat")
    # format into list of lists
    if hasattr(xTangent, "__len__"):
        xyzTangent = numpy.array([xTangent,yTangent,zTangent]).T.tolist()
        xyzWok = libcoordio.tangentToWokArr(
            xyzTangent, list(b), list(iHat), list(jHat), list(kHat),
            elementHeight, scaleFac, dx, dy, dz)

        xyzWok = numpy.array(xyzWok)
        xWok = xyzWok[:,0]
        yWok = xyzWok[:,1]
        zWok = xyzWok[:,2]
    else:
        xWok, yWok, zWok = libcoordio.tangentToWok(
            [xTangent, yTangent, zTangent], list(b), list(iHat), list(jHat),
            list(kHat), elementHeight, scaleFac, dx, dy, dz
        )
    return xWok, yWok, zWok

def _tangentToWok(xTangent, yTangent, zTangent, b, iHat, jHat, kHat,
                 elementHeight=defaults.POSITIONER_HEIGHT, scaleFac=1,
                 dx=0, dy=0, dz=0):
    """
    Convert from tangent coordinates at b to wok coordinates.

    xyz Wok coords are mm with orgin at wok vertex. +z points from wok toward M2.
    -x points toward the boss slithead

    In the tangent coordinate frame the xy plane is tangent to the wok
    surface at xyz position b.  The origin is set to be elementHeight above
    the wok surface.

    Parameters
    -------------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    zTangent: scalar or 1D array
        z position (mm) in tangent coordinates
    b: 3-vector
        x,y,z position (mm) of element hole on wok surface measured in wok coords
    iHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate x axis
    jHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate y axis
    kHat: 3-vector
        x,y,z unit vector in wok coords that indicate the direction
        of the tangent coordinate z axis
    elementHeight: scalar
        height (mm) of positioner/GFA chip above wok surface
    scaleFac: scalar
        scale factor to apply to b to account for thermal expansion of wok.
        scale is applied to b radially
    dx: scalar
        x offset (mm), calibration to capture small displacements of
        tangent x
    dy: scalar
        y offset (mm), calibration to capture small displacements of
        tangent y
    dz: scalar
        z offset (mm), calibration to capture small displacements of
        tangent x

    Returns
    ---------
    xWok: scalar or 1D array
        x position of object in wok space mm
        (+x aligned with +RA on image at PA = 0)
    yWok: scalar or 1D array
        y position of object in wok space mm
        (+y aligned with +Dec on image at PA = 0)
    zWok: scalar or 1D array
        z position of object in wok space mm
        (+z aligned boresight and increases from the telescope to the sky)
    """
    # check for problems
    b = _verify3Vector(b, "b")
    iHat = _verify3Vector(iHat, "iHat")
    jHat = _verify3Vector(jHat, "jHat")
    kHat = _verify3Vector(kHat, "kHat")

    calibOff = numpy.array([dx,dy,dz])
    isArr = hasattr(xTangent, "__len__")
    if isArr:
        calibOff = numpy.array([calibOff]*len(xTangent)).T

    coords = numpy.array([xTangent,yTangent,zTangent])

    # apply calibration offsets
    coords = coords + calibOff

    # apply rotational calibration
    # if dRot != 0:
    #     rotZ = -1*numpy.radians(dRot)
    #     rotMatZ = numpy.array([
    #         [numpy.cos(rotZ), numpy.sin(rotZ), 0],
    #         [-numpy.sin(rotZ), numpy.cos(rotZ), 0],
    #         [0, 0, 1]
    #     ])

    #     coords = rotMatZ.dot(coords)

    # offset xy plane by element height
    if isArr:
        coords[2, :] = coords[2, :] + elementHeight
    else:
        coords[2] = coords[2] + elementHeight

    # rotate normal to wok surface at point b
    # invRotNorm = numpy.linalg.inv(numpy.array([iHat, jHat, kHat], dtype="float64"))
    # transpose is inverse! ortho-normal rows or unit vectors!
    invRotNorm = numpy.array([iHat, jHat, kHat], dtype="float64").T

    coords = invRotNorm.dot(coords)

    # offset to wok base position
    # apply radial scale factor to b
    # assume that z is unaffected by
    # thermal expansion (probably reasonable)
    if scaleFac != 1:
        r = numpy.sqrt(b[0]**2+b[1]**2)
        theta = numpy.arctan2(b[1], b[0])
        r = r * scaleFac
        b[0] = r * numpy.cos(theta)
        b[1] = r * numpy.sin(theta)

    if isArr:
        b = numpy.array([b]*len(xTangent)).T
    coords = coords + b

    return coords[0], coords[1], coords[2]

def proj2XYplane(x, y, z, rayOrigin):
    """Given a point x, y, z orginating from rayOrigin, project
    this point onto the xy plane.

    Parameters
    ------------
    x: scalar or 1D array
        x position of point to be projected
    y: scalar or 1D array
        y position of point to be projected
    z: scalar or 1D array
        z position of point to be projected
    rayOrigin: 3-vector
        [x,y,z] origin of ray(s)

    Returns
    --------
    xProj: scalar or 1D array
        projected x value
    yProj: scalar or 1D array
        projected y value
    zProj: scalar or 1D array
        projected z value (should be zero always!)
    projDist: scalar or 1D array
        projected distance (a proxy for something like focus offset)
    """

    # intersections of lines to planes...
    # http://geomalgorithms.com/a05-_intersect-1.html
    # intersect3D_SegmentPlane
    rayOrigin = _verify3Vector(rayOrigin, "rayOrigin")
    if hasattr(x, "__len__"):
        rayOrigin = numpy.array([rayOrigin] * len(z), dtype="float64").T
        x = numpy.array(x, dtype="float64")
        y = numpy.array(y, dtype="float64")
        z = numpy.array(z, dtype="float64")

    xyz = numpy.array([x, y, z], dtype="float64")

    u = xyz - rayOrigin
    w = rayOrigin
    zHat = numpy.array([0, 0, 1])  # normal to xy plane
    D = zHat.dot(u)
    N = -1 * zHat.dot(w)

    sI = N / D
    xyzProj = rayOrigin + sI * u
    projDist = numpy.linalg.norm(xyz - xyzProj, axis=0)
    # negative projDist for points below xy plane
    projDist *= numpy.sign(z)

    return xyzProj[0], xyzProj[1], xyzProj[2], projDist

def tangentToPositioner(
    xTangent, yTangent, xBeta, yBeta, la=7.4, alphaOffDeg=0, betaOffDeg=0,
    lefthand=False,
):
    """
    Determine alpha/beta positioner angles that place xBeta, yBeta coords in mm
    at xTangent, yTangent.  CPP Version.

    todo: include hooks for positioner non-linearity

    Parameters
    -------------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    xBeta: scalar or 1D array
        x position (mm) in beta arm frame
    yBeta: scalar or 1D array
        y position (mm) in beta arm frame
    la: scalar or 1D array
        length (mm) of alpha arm
    alphaOffDeg: scalar
        alpha zeropoint offset (deg)
    betaOffDeg: scalar
        beta zeropoint offset (deg)
    lefthand: boolean
        whether to return a lefthand or righthand parity robot

    Returns
    ---------
    alphaDeg: scalar or 1D array
        alpha angle in degrees
    betaDeg: scalar or 1D array
        beta angle in degrees
    isOK: boolean or 1D boolean array
        True if point physically accessible, False otherwise
    """
    # C++ wrapped stuff wants lists, not numpy arrays
    isArr = hasattr(xTangent, "__len__")
    if isArr:
        xyTangent = numpy.array([xTangent, yTangent]).T.tolist()

        if hasattr(xBeta, "__len__"):
            xyBeta = numpy.array([xBeta, yBeta]).T.tolist()
        else:
            xyBeta = [[xBeta, yBeta]]*len(xyTangent)

        alphaBeta = libcoordio.tangentToPositionerArr(
            xyTangent, xyBeta, la, alphaOffDeg, betaOffDeg, lefthand
        )

        alphaBeta = numpy.array(alphaBeta)
        alpha = alphaBeta[:, 0]
        beta = alphaBeta[:, 1]
    else:
        alpha, beta = libcoordio.tangentToPositioner(
            [xTangent, yTangent], [xBeta, yBeta],
            la, alphaOffDeg, betaOffDeg, lefthand
        )

    isOKAlpha = numpy.isfinite(alpha)
    isOKBeta = numpy.isfinite(beta)
    isOK = isOKAlpha & isOKBeta

    if not isArr:
        isOK = bool(isOK)

    return alpha, beta, isOK

def _tangentToPositioner(
    xTangent, yTangent, xBeta, yBeta, la=7.4, alphaOffDeg=0, betaOffDeg=0
):
    """
    Determine alpha/beta positioner angles that place xBeta, yBeta coords in mm
    at xTangent, yTangent.  Python Version.

    todo: include hooks for positioner non-linearity

    note, didn't add left hand option here...use C++

    Parameters
    -------------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    xBeta: scalar or 1D array
        x position (mm) in beta arm frame
    yBeta: scalar or 1D array
        y position (mm) in beta arm frame
    la: scalar or 1D array
        length (mm) of alpha arm
    alphaOffDeg: scalar
        alpha zeropoint offset (deg)
    betaOffDeg: scalar
        beta zeropoint offset (deg)

    Returns
    ---------
    alphaDeg: scalar or 1D array
        alpha angle in degrees
    betaDeg: scalar or 1D array
        beta angle in degrees
    isOK: boolean or 1D boolean array
        True if point physically accessible, False otherwise
    """

    isArr = hasattr(xTangent, "__len__")
    if isArr:
        xTangent = numpy.array(xTangent, dtype="float64")
        yTangent = numpy.array(yTangent, dtype="float64")

    if hasattr(xBeta, "__len__"):
        xBeta = numpy.array(xBeta, dtype="float64")
        yBeta = numpy.array(yBeta, dtype="float64")

    # polar coords jive better for this calculation
    thetaTangent = numpy.arctan2(yTangent, xTangent)
    rTangentSq = xTangent**2 + yTangent**2

    # convert xy Beta to radial coords
    # the origin of the beta coord system is the
    # beta axis of rotation
    thetaBAC = numpy.arctan2(yBeta, xBeta) # radians!
    rBacSq = xBeta**2+yBeta**2

    gamma = numpy.arccos(
        (la**2 + rBacSq - rTangentSq) / (2 * la * numpy.sqrt(rBacSq))
    )
    xi = numpy.arccos(
        (la**2 + rTangentSq - rBacSq) / (2 * la * numpy.sqrt(rTangentSq))
    )

    thetaTangent = numpy.degrees(thetaTangent)
    thetaBAC = numpy.degrees(thetaBAC)
    gamma = numpy.degrees(gamma)
    xi = numpy.degrees(xi)

    alphaDeg = thetaTangent - xi - alphaOffDeg
    betaDeg = 180 - gamma - thetaBAC - betaOffDeg

    # look for nans
    isOKAlpha = numpy.isfinite(alphaDeg)
    isOKBeta = numpy.isfinite(betaDeg)
    isOK = isOKAlpha & isOKBeta

    if not isArr:
        isOK = bool(isOK)

    # handle wrapping? not sure it's necessary
    alphaDeg = alphaDeg % 360
    betaDeg = betaDeg % 360  # should already be there, but...
    # if hasattr(betaDeg, "__len__"):
    #     if True in betaDeg > 180:
    #         raise RuntimeError("problem in alpha/beta conversion! beta > 180")
    # elif betaDeg > 180:
    #     raise RuntimeError("problem in alpha/beta conversion! beta > 180")

    return alphaDeg, betaDeg, isOK

def positionerToTangent(
    alphaDeg, betaDeg, xBeta, yBeta, la=7.4, alphaOffDeg=0, betaOffDeg=0
):
    """
    Determine tangent coordinates (mm) of xBeta, yBeta coords in mm
    from alpha/beta angle.  CPP version.

    todo: include hooks for positioner non-linearity

    Parameters
    -------------
    alphaDeg: scalar or 1D array
        alpha angle in degrees
    betaDeg: scalar or 1D array
        beta angle in degrees
    xBeta: scalar or 1D array
        x position (mm) in beta arm frame
    yBeta: scalar or 1D array
        y position (mm) in beta arm frame
    la: scalar or 1D array
        length (mm) of alpha arm
    alphaOffDeg: scalar
        alpha zeropoint offset (deg)
    betaOffDeg: scalar
        beta zeropoint offset (deg)

    Returns
    ---------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    """
    if hasattr(alphaDeg, "__len__"):
        alphaBeta = numpy.array([alphaDeg, betaDeg]).T.tolist()

        if hasattr(xBeta, "__len__"):
            xyBeta = numpy.array([xBeta, yBeta]).T.tolist()
        else:
            xyBeta = [[xBeta, yBeta]] * len(alphaBeta)

        xyTangent = libcoordio.positionerToTangentArr(
            alphaBeta, xyBeta, la, alphaOffDeg, betaOffDeg
        )

        xyTangent = numpy.array(xyTangent)
        xTangent = xyTangent[:, 0]
        yTangent = xyTangent[:, 1]
    else:
        xTangent, yTangent = libcoordio.positionerToTangent(
            [alphaDeg, betaDeg], [xBeta, yBeta],
            la, alphaOffDeg, betaOffDeg
        )
    return xTangent, yTangent

def _positionerToTangent(
    alphaDeg, betaDeg, xBeta, yBeta, la=7.4, alphaOffDeg=0, betaOffDeg=0
):
    """
    Determine tangent coordinates (mm) of xBeta, yBeta coords in mm
    from alpha/beta angle.  Python version.

    todo: include hooks for positioner non-linearity

    Parameters
    -------------
    alphaDeg: scalar or 1D array
        alpha angle in degrees
    betaDeg: scalar or 1D array
        beta angle in degrees
    xBeta: scalar or 1D array
        x position (mm) in beta arm frame
    yBeta: scalar or 1D array
        y position (mm) in beta arm frame
    la: scalar or 1D array
        length (mm) of alpha arm
    alphaOffDeg: scalar
        alpha zeropoint offset (deg)
    betaOffDeg: scalar
        beta zeropoint offset (deg)

    Returns
    ---------
    xTangent: scalar or 1D array
        x position (mm) in tangent coordinates
    yTangent: scalar or 1D array
        y position (mm) in tangent coordinates
    """
    # convert xy Beta to radial coords
    # the origin of the beta coord system is the
    # beta axis of rotation
    if hasattr(xBeta, "__len__"):
        xBeta = numpy.array(xBeta, dtype="float64")
        yBeta = numpy.array(yBeta, dtype="float64")

    thetaBAC = numpy.arctan2(yBeta, xBeta)  # radians!
    rBAC = numpy.sqrt(xBeta**2 + yBeta**2)
    alphaRad = numpy.radians(alphaDeg+alphaOffDeg)
    betaRad = numpy.radians(betaDeg+betaOffDeg)

    cosAlpha = numpy.cos(alphaRad)
    sinAlpha = numpy.sin(alphaRad)
    cosAlphaBeta = numpy.cos(alphaRad + betaRad + thetaBAC)
    sinAlphaBeta = numpy.sin(alphaRad + betaRad + thetaBAC)

    xTangent = la * cosAlpha + rBAC * cosAlphaBeta
    yTangent = la * sinAlpha + rBAC * sinAlphaBeta

    return xTangent, yTangent

class FVCUW(object):
    # right now this is for UWs test bench, fix image scale and use
    # euclidean transforms
    defaultScale = 0.026926930620282834
    defaultRotation = -0.005353542811111436
    defaultTranslation = numpy.array([-75.69606047, -48.68561274])
    # transform from CCD pixels to mm
    def __init__(self):

        self.tform = SimilarityTransform(
            scale=self.defaultScale,
            rotation=self.defaultRotation,
            translation=self.defaultTranslation
         )

    def fit(self, xpix, ypix, xWok, yWok):
        xypixels = numpy.array([xpix, ypix]).T
        xyWok = numpy.array([xWok,yWok]).T
        self.tform.estimate(xypixels, xyWok)

    def fvcToWok(self, xpix, ypix):
        xypixels = numpy.array([xpix, ypix]).T
        xyWok = self.tform(xypixels)
        return xyWok[:,0], xyWok[:,1]

    def wokToFVC(self, xWok, yWok):
        xyWok = numpy.array([xWok, yWok]).T
        xypix = self.tform.inverse(xyWok)
        return xypix[:,0], xypix[:,1]

ext: py

max_stars_repo_name: sdss/coordio

max_issues_repo_name: sdss/coordio

max_forks_repo_name: sdss/coordio

reason: import numpy


text:
import networkx as nx
import numpy as np
from autode.log import logger
from autode.mol_graphs import is_isomorphic

def get_sum_energy_mep(saddle_point_r1r2, pes_2d):
    """
    Calculate the sum of the minimum energy path that traverses reactants (r)
    to products (p) via the saddle point (s)::

                /          p
               /     s
          r2  /
             /r
             ------------
                  r1

    Arguments:
        saddle_point_r1r2 (tuple(float)):

        pes_2d (autode.pes_2d.PES2d):

    Returns:
        (float): Path energy (Ha)
    """
    logger.info('Finding the total energy along the minimum energy pathway')

    reactant_point = (0, 0)
    product_point, product_energy = None, 9999

    # The saddle point indexes are those that are closest tp the saddle points
    # r1 and r2 distances
    saddle_point = (np.argmin(np.abs(pes_2d.r1s - saddle_point_r1r2[0])),
                    np.argmin(np.abs(pes_2d.r2s - saddle_point_r1r2[1])))

    # Generate a grid graph (i.e. all nodes are corrected
    energy_graph = nx.grid_2d_graph(pes_2d.n_points_r1, pes_2d.n_points_r2)

    min_energy = min([species.energy for species in pes_2d.species.flatten()])

    # For energy point on the 2D surface
    for i in range(pes_2d.n_points_r1):
        for j in range(pes_2d.n_points_r2):
            point_rel_energy = pes_2d.species[i, j].energy - min_energy

            # Populate the relative energy of each node in the graph
            energy_graph.nodes[i, j]['energy'] = point_rel_energy

            # Find the point where products are made
            if is_isomorphic(graph1=pes_2d.species[i, j].graph,
                             graph2=pes_2d.product_graph):

                # If products have not yet found, or they have and the energy
                # are lower but are still isomorphic
                if product_point is None or point_rel_energy < product_energy:
                    product_point = (i, j)
                    product_energy = point_rel_energy

    logger.info(f'Reactants at r1={pes_2d.r1s[0]:.4f} , '
                f'r2={pes_2d.r2s[0]:.4f} Å and '
                f'products r1={pes_2d.rs[product_point][0]:.4f}, '
                f'r2={pes_2d.rs[product_point][1]:.4f} Å')

    def energy_diff(curr_node, final_node, d):
        """Energy difference between the twp points on the graph. d is required
         to satisfy nx. Must only increase in energy to a saddle point so take
          the magnitude to prevent traversing s mistakenly"""
        return (np.abs(energy_graph.nodes[final_node]['energy']
                       - energy_graph.nodes[curr_node]['energy']))

    # Calculate the energy along the MEP up to the saddle point from reactants
    # and products
    path_energy = 0.0

    for point in (reactant_point, product_point):
        path_energy += nx.dijkstra_path_length(energy_graph,
                                               source=point,
                                               target=saddle_point,
                                               weight=energy_diff)

    logger.info(f'Path energy to {saddle_point} is {path_energy:.4f} Hd')
    return path_energy

max_stars_repo_name: Crossfoot/autodE

max_issues_repo_name: skphy/autodE

reason: import numpy,import networkx


text:
import astropy.units as u
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from plyades.config import config

# Values taken from:
# Archinal, et al. "Report of the IAU working group on cartographic
# coordinates and rotational elements: 2009." Celestial Mechanics
# and Dynamical Astronomy 109.2 (2011): 101-135.
#
# Archinal, et al. "Erratum to: Reports of the IAU Working Group on
# Cartographic Coordinates and Rotational Elements: 2006 & 2009."
# Celestial Mechanics and Dynamical Astronomy 110.4 (2011): 401-403.
#
# Vallado, David A., and Wayne D. McClain.
# Fundamentals of astrodynamics and applications.
# Fourth Edition. Springer Science & Business Media, 2013.

# symbols = {
#     "Sun": u"\u2609", "Mercury": u"\u263F", "Venus": u"\u2640",
#     "Earth": u"\u2641", "Mars": u"\u2642", "Jupiter": u"\u2643",
#     "Saturn": u"\u2644", "Uranus": u"\u26E2", "Neptune": u"\u2646",
#     "Moon": u"\u263E", "Pluto": u"\u2647"}

class Planet:
    def __init__(
        self,
        name,
        mu,
        mean_radius,
        equatorial_radius,
        polar_radius,
        j2,
        jpl_id,
        symbol,
        # rotational_elements,
    ):
        self.name = name
        self.mu = mu
        self.mean_radius = mean_radius
        self.equatorial_radius = equatorial_radius
        self.polar_radius = polar_radius
        self.j2 = j2
        self.jpl_id = jpl_id
        self.symbol = symbol
        # self._rotational_elements = rotational_elements

    def _repr_latex_(self):
        strs = [
            "{} {}".format(self.name, self.symbol),
            "Gravitational parameter:",
            "$\mu={}$ {}".format(self.mu.value, self.mu.unit),
            "Mean radius:",
            "$r_m={}$: {}".format(self.mean_radius.value, self.mean_radius.unit),
            "Equatorial radius:",
            "$r_e={}$ {}".format(self.equatorial_radius.value, self.equatorial_radius.unit),
            "Polar radius:",
            "$r_p={}$ {}".format(self.polar_radius.value, self.polar_radius.unit),
            "J2:",
            "$J_2={}$".format(self.j2),
        ]
        return "<br />".join(strs)

    def rv(self, jd, jd2=0.0):
        return config['ephemeris'].rv(self.jpl_id, jd, jd2)

    def wrt(self, body, jd, jd2=0.0):
        r_target, v_target = self.rv(jd, jd2)
        r_origin, v_origin = body.rv(jd, jd2)
        return r_target-r_origin, v_target-v_origin

    def plot3d(self, ax=None):
        if not ax:
            fig = plt.figure("Plyades Plot")
            ax = fig.add_subplot(111, projection='3d')
        re = self.equatorial_radius
        rp = self.polar_radius

        u = np.linspace(0, 2 * np.pi, 100)
        v = np.linspace(0, np.pi, 100)

        x = re * np.outer(np.cos(u), np.sin(v))
        y = re * np.outer(np.sin(u), np.sin(v))
        z = rp * np.outer(np.ones(np.size(u)), np.cos(v))
        ax.plot_surface(x, y, z,
            rstride=4,
            cstride=4,
            color='b',
            alpha=.3,
            linewidth=1,
            edgecolor="b"
        )

class Moon:
    pass

class SmallBody:
    pass

MERCURY = Planet(
    'Mercury',
    2.2032e4*u.km**3/u.s**2,
    2439.7*u.km,
    2439.7*u.km,
    2439.7*u.km,
    0.00006,
    199,
    u"\u263F",
)
VENUS = Planet(
    'Venus',
    3.257e5*u.km**3/u.s**2,
    6051.8*u.km,
    6051.8*u.km,
    6051.8*u.km,
    0.000027,
    299,
    u"\u2640",
)
EARTH = Planet(
    'Earth',
    398600.4418*u.km**3/u.s**2,
    6371.0084*u.km,
    6378.1366*u.km,
    6356.7519*u.km,
    0.0010826269,
    399,
    u"\u2641",
)
MARS = Planet(
    'Mars',
    4.305e4*u.km**3/u.s**2,
    3389.50*u.km,
    3396.19*u.km,
    3376.20*u.km,
    0.001964,
    4,
    u"\u2642",
)
JUPITER = Planet(
    'Jupiter',
    1.268e8*u.km**3/u.s**2,
    69911.0*u.km,
    71492.0*u.km,
    66854.0*u.km,
    0.01475,
    5,
    u"\u2643",
)
SATURN = Planet(
    'Saturn',
    3.794e7*u.km**3/u.s**2,
    58232.0*u.km,
    60268.0*u.km,
    54364.0*u.km,
    0.01645,
    6,
    u"\u2644",
)
URANUS = Planet(
    'Uranus',
    5.794e6*u.km**3/u.s**2,
    25362.0*u.km,
    25559.0*u.km,
    24973.0*u.km,
    0.012,
    7,
    u"\u26E2",
)
NEPTUNE = Planet(
    'Neptune',
    6.809e6*u.km**3/u.s**2,
    24622.0*u.km,
    24764.0*u.km,
    24341.0*u.km,
    0.004,
    8,
    u"\u2646",
)

reason: import numpy,import astropy


text:
from neuron import h, rxd
h.load_file("stdrun.hoc")

sec = h.Section(name="sec")
sec.L = 10
sec.nseg = 11
sec.diam = 5
rxd.set_solve_type(dimension=3)

cyt = rxd.Region(h.allsec(), name="cyt", nrn_region="i")
ip3 = rxd.Species(cyt, name="ip3", initial=lambda nd: 1000 if nd.segment == sec(0.3) else 0 )

def callbackfun():
    return 1000

for nd in ip3.nodes(sec(0.1)):
    nd.include_flux(1000)

for nd in ip3.nodes(sec(0.5)):
    nd.include_flux(callbackfun)

for nd in ip3.nodes(sec(0.9)):
    nd.include_flux(sec(0.3)._ref_ip3i)

h.finitialize(-70)
h.continuerun(10)
import numpy
ip33d = ip3._intracellular_instances[cyt]
con = []
for seg in sec:
    con.append(numpy.mean([nd.concentration*nd.volume for nd in ip3.nodes(seg)])/seg.volume())
print(con)

con1D = [0.0, 0.0009302760807252256, 0.0, 1000.0, 0.0, 0.0009302760807252256, 0.0, 0.0, 0.0, 0.000930276080725227, 0.0]

diff = numpy.array(con) - numpy.array(con1D)

max_stars_repo_name: tommorse/nrn

max_issues_repo_name: tommorse/nrn

max_forks_repo_name: tommorse/nrn


text:
# 二次罚函数法（不等式约束）
def penalty_quadratic(funcs, args, cons, x_0, draw=True, output_f=False, method="gradient_descent", sigma=10, p=0.4, epsilon=1e-10, k=0):
    '''
    Parameters
    ----------
    funcs : sympy.matrices.dense.MutableDenseMatrix
        当前目标方程

    args : sympy.matrices.dense.MutableDenseMatrix
        参数列表

    cons : sympy.matrices.dense.MutableDenseMatrix
        不等式参数约束列表

    x_0 : list
        初始迭代点列表

    draw : bool
        绘图接口参数

    output_f : bool
        输出迭代函数值列表

    method : string
        无约束优化方法内核

    sigma : double
        罚函数因子

    p : double
        修正参数

    epsilon : double
        迭代停机准则

    k : int
        迭代次数

    Returns
    -------
    tuple
        最终收敛点, 迭代次数, (迭代函数值列表)

    '''
    import numpy as np
    import sympy as sp
    from optimtool.functions.tools import function_f_x_k, function_plot_iteration, function_data_convert
    from optimtool.unconstrain.gradient_descent import barzilar_borwein
    from optimtool.unconstrain.newton import CG
    from optimtool.unconstrain.newton_quasi import L_BFGS
    from optimtool.unconstrain.trust_region import steihaug_CG
    assert sigma > 0
    assert p > 0
    assert p < 1
    funcs, args, _, cons = function_data_convert(funcs, args, None, cons)
    point = []
    f = []
    while 1:
        point.append(np.array(x_0))
        f.append(function_f_x_k(funcs, args, x_0))
        reps = dict(zip(args, x_0))
        consv = np.array(cons.subs(reps)).astype(np.float64)
        consv = np.where(consv <= 0, consv, 1)
        consv = np.where(consv > 0, consv, 0)
        pe = sp.Matrix([funcs + (sigma / 2) * cons.T * consv])
        if method == "gradient_descent":
            x_0, _ = barzilar_borwein(pe, args, tuple(x_0), draw=False)
        elif method == "newton":
            x_0, _ = CG(pe, args, tuple(x_0), draw=False)
        elif method == "newton_quasi":
            x_0, _ = L_BFGS(pe, args, tuple(x_0), draw=False)
        elif method == "trust_region":
            x_0, _ = steihaug_CG(pe, args, tuple(x_0), draw=False)
        k = k + 1
        if np.linalg.norm(x_0 - point[k - 1]) < epsilon:
            point.append(np.array(x_0))
            f.append(function_f_x_k(funcs, args, x_0))
            break
        sigma = p * sigma
    function_plot_iteration(f, draw, "penalty_quadratic_unequal")
    if output_f is True:
        return x_0, k, f
    else:
        return x_0, k

# 内点罚函数法（不等式约束）
'''
保证点在定义域内
'''
def penalty_interior_log(funcs, args, cons, x_0, draw=True, output_f=False, sigma=12, p=0.6, epsilon=1e-10, k=0):
    '''
    Parameters
    ----------
    funcs : sympy.matrices.dense.MutableDenseMatrix
        当前目标方程

    args : sympy.matrices.dense.MutableDenseMatrix
        参数列表

    cons : sympy.matrices.dense.MutableDenseMatrix
        不等式参数约束列表

    x_0 : list
        初始迭代点列表

    draw : bool
        绘图接口参数

    output_f : bool
        输出迭代函数值列表

    sigma : double
        罚函数因子

    p : double
        修正参数

    epsilon : double
        迭代停机准则

    k : int
        迭代次数

    Returns
    -------
    tuple
        最终收敛点, 迭代次数, (迭代函数值列表)

    '''
    import numpy as np
    import sympy as sp
    from optimtool.functions.tools import function_f_x_k, function_plot_iteration, function_data_convert
    from optimtool.hybrid.approximate_point_gradient import neg_log
    assert sigma > 0
    assert p > 0
    assert p < 1
    funcs, args, _, cons = function_data_convert(funcs, args, None, cons)
    point = []
    f = []
    while 1:
        point.append(np.array(x_0))
        f.append(function_f_x_k(funcs, args, x_0))
        x_0, _ = neg_log(funcs, sigma, -cons, args, tuple(x_0), draw=False)
        k = k + 1
        sigma = p * sigma
        print(np.linalg.norm(x_0 - point[k - 1]))
        if np.linalg.norm(x_0 - point[k - 1]) < epsilon:
            point.append(np.array(x_0))
            f.append(function_f_x_k(funcs, args, x_0))
            break
    function_plot_iteration(f, draw, "penalty_interior_fraction")
    if output_f is True:
        return x_0, k, f
    else:
        return x_0, k

# 分式
def penalty_interior_fraction(funcs, args, cons, x_0, draw=True, output_f=False, method="gradient_descent", sigma=12, p=0.6, epsilon=1e-6, k=0):
    '''
    Parameters
    ----------
    funcs : sympy.matrices.dense.MutableDenseMatrix
        当前目标方程

    args : sympy.matrices.dense.MutableDenseMatrix
        参数列表

    cons : sympy.matrices.dense.MutableDenseMatrix
        不等式参数约束列表

    x_0 : list
        初始迭代点列表

    draw : bool
        绘图接口参数

    output_f : bool
        输出迭代函数值列表

    method : string
        无约束优化方法内核

    sigma : double
        罚函数因子

    p : double
        修正参数

    epsilon : double
        迭代停机准则

    k : int
        迭代次数

    Returns
    -------
    tuple
        最终收敛点, 迭代次数, (迭代函数值列表)

    '''
    import numpy as np
    import sympy as sp
    from optimtool.functions.tools import function_f_x_k, function_plot_iteration, function_data_convert
    from optimtool.unconstrain.gradient_descent import barzilar_borwein
    from optimtool.unconstrain.newton import CG
    from optimtool.unconstrain.newton_quasi import L_BFGS
    from optimtool.unconstrain.trust_region import steihaug_CG
    assert sigma > 0
    assert p > 0
    assert p < 1
    funcs, args, _, cons = function_data_convert(funcs, args, None, cons)
    point = []
    f = []
    sub_pe = 0
    for i in cons:
        sub_pe += 1 / i
    sub_pe = sp.Matrix([sub_pe])
    while 1:
        point.append(np.array(x_0))
        f.append(function_f_x_k(funcs, args, x_0))
        pe = sp.Matrix([funcs - sigma * sub_pe])
        if method == "gradient_descent":
            x_0, _ = barzilar_borwein(pe, args, tuple(x_0), draw=False)
        elif method == "newton":
            x_0, _ = CG(pe, args, tuple(x_0), draw=False)
        elif method == "newton_quasi":
            x_0, _ = L_BFGS(pe, args, tuple(x_0), draw=False)
        elif method == "trust_region":
            x_0, _ = steihaug_CG(pe, args, tuple(x_0), draw=False)
        k = k + 1
        sigma = p * sigma
        if np.linalg.norm(x_0 - point[k - 1]) < epsilon:
            point.append(np.array(x_0))
            f.append(function_f_x_k(funcs, args, x_0))
            break
    function_plot_iteration(f, draw, "penalty_interior_fraction")
    if output_f is True:
        return x_0, k, f
    else:
        return x_0, k

# 增广拉格朗日函数法（不等式约束）
def lagrange_augmented(funcs, args, cons, x_0, draw=True, output_f=False, method="gradient_descent", muk=10, sigma=8, alpha=0.2, beta=0.7, p=2, eta=1e-1, epsilon=1e-4, k=0):
    '''
    Parameters
    ----------
    funcs : sympy.matrices.dense.MutableDenseMatrix
        当前目标方程

    args : sympy.matrices.dense.MutableDenseMatrix
        参数列表

    cons : sympy.matrices.dense.MutableDenseMatrix
        不等式参数约束列表

    x_0 : list
        初始迭代点列表

    draw : bool
        绘图接口参数

    output_f : bool
        输出迭代函数值列表

    method : string
        无约束优化方法内核

    muk : float
        因子

    sigma : float
        罚函数因子

    alpha : float
        初始步长

    beta : float
        修正参数

    p : float
        修正参数

    eta : float
        常数

    epsilon : float
        迭代停机准则

    k : int
        迭代次数

    Returns
    -------
    tuple
        最终收敛点, 迭代次数, (迭代函数值列表)

    '''
    import sympy as sp
    import numpy as np
    from optimtool.functions.tools import function_f_x_k, function_plot_iteration, function_cons_unequal_L, function_renew_mu_k, function_v_k, function_data_convert
    from optimtool.unconstrain.gradient_descent import barzilar_borwein
    from optimtool.unconstrain.newton import CG
    from optimtool.unconstrain.newton_quasi import L_BFGS
    from optimtool.unconstrain.trust_region import steihaug_CG
    assert sigma > 0
    assert p > 1
    assert alpha > 0
    assert alpha <= beta
    assert beta < 1
    funcs, args, _, cons = function_data_convert(funcs, args, None, cons)
    f = []
    muk = np.array([muk for i in range(cons.shape[0])]).reshape(cons.shape[0], 1)
    while 1:
        etak = 1 / sigma
        epsilonk = 1 / sigma**alpha
        cons_uneuqal_modifyed = function_cons_unequal_L(cons, args, muk, sigma, x_0)
        L = sp.Matrix([funcs + (sigma / 2) * cons_uneuqal_modifyed])
        f.append(function_f_x_k(funcs, args, x_0))
        if method == "gradient_descent":
            x_0, _ = barzilar_borwein(L, args, x_0, draw=False, epsilon=etak)
        elif method == "newton":
            x_0, _ = CG(L, args, x_0, draw=False, epsilon=etak)
        elif method == "newton_quasi":
            x_0, _ = L_BFGS(L, args, x_0, draw=False, epsilon=etak)
        elif method == "trust_region":
            x_0, _ = steihaug_CG(L, args, x_0, draw=False, epsilon=etak)
        k = k + 1
        vkx = function_v_k(None, cons, args, muk, sigma, x_0)
        if vkx <= epsilonk:
            res = L.jacobian(args)
            if (vkx <= epsilon) and (np.linalg.norm(np.array(res.subs(dict(zip(args, x_0)))).astype(np.float64)) <= eta):
                f.append(function_f_x_k(funcs, args, x_0))
                break
            else:
                muk = function_renew_mu_k(cons, args, muk, sigma, x_0)
                sigma = sigma
                etak = etak / sigma
                epsilonk = epsilonk / sigma**beta
        else:
            sigma = p * sigma
            etak = 1 / sigma
            epsilonk  = 1 / sigma**alpha
    function_plot_iteration(f, draw, "lagrange_augmented_unequal")
    if output_f is True:
        return x_0, k, f
    else:
        return x_0, k

max_stars_repo_name: linjing-lab/optimtool

max_issues_repo_name: linjing-lab/optimtool

max_forks_repo_name: linjing-lab/optimtool

reason: import numpy,import sympy


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

"""
utils.py
~~~~~~~~

A selection of useful helper functions used throughout the Tesseract library.

"""

import logging
from datetime import datetime, date
from functools import wraps
from timeit import default_timer as timer

import numpy as np

def resolve_date(d):
    """Convert a str or date to an appropriate datetime.

    Strings should be of the format '%Y', '%Y-%m or '%Y-%m-%d', for example:
    '2012', '1994-02' or '1991-12-11'. Date objects with no time information
    will be rounded down to the midnight beginning that date.

    Args:
        d (Union[str, date]): The string or date to convert.

    Returns:
        datetime: The parsed datetime equivalent of d.
    """
    if isinstance(d, datetime):
        return d

    if isinstance(d, date):
        return datetime.combine(d, datetime.min.time())

    for fmt in ('%Y', '%Y-%m', '%Y-%m-%d'):
        try:
            return datetime.strptime(d, fmt)
        except ValueError:
            pass

    raise ValueError('date string format not recognized.')

def check_for_raw_scores(y_pred):
    # Heuristic to check if input are raw scores
    if y_pred.ndim > 1:
        for v in y_pred:
            if ((np.linalg.norm(v, 0),
                 np.linalg.norm(v), 2) != (1, 1)):
                return True
    return False

def select_prediction_function(clf, scores_only=False, labels_only=False):
    if hasattr(clf, 'predict_proba') and not labels_only:
        prediction_function = clf.predict_proba
    elif hasattr(clf, 'decision_function') and not labels_only:
        prediction_function = clf.decision_function
    elif hasattr(clf, 'predict') and not scores_only:
        prediction_function = clf.predict
    else:
        raise TypeError(
            'Unsure how to handle predictions with '
            'classifier of type {}.'.format(clf.__class__))
    return prediction_function

def resolve_categorical(y):
    return np.argmax(y, 1) if y.ndim > 1 else y

def binary_labels(array, positive='malicious', negative='benign'):
    return [positive if x else negative for x in array]

def parse_percentage(n):
    return float(n[:-1]) / 100

def resolve_percentage(n):
    return parse_percentage(n) if isinstance(n, str) else n

def seconds_to_time(seconds):
    """Return a nicely formatted time given the number of seconds."""
    m, s = divmod(seconds, 60)
    h, m = divmod(m, 60)
    d, h = divmod(h, 24)
    return "%d days, %02d hours, %02d minutes, %02d seconds" % (d, h, m, s)

def timing(f):
    @wraps(f)
    def wrap(*args, **kwargs):
        start = timer()
        result = f(*args, **kwargs)
        elapsed = seconds_to_time(timer() - start)
        logging.debug('{} took: {}'.format(f.__name__, elapsed))
        return result

    return wrap


text:
import numpy as np
import pandas as pd

class Data:
    def __init__(self, metadata, pixels, fc6):
        self.metadata = metadata
        self.pixels = pixels
        self.fc6 = fc6
    def filter_out(self, kinds = {'invalid' : True,
                                  'outlier' : True}):
        '''
        input:
            frame = pandas dataframe
            remove = boolean flag, if True will remove flagged data
            kind = str, can be either 'invalid' or 'outlier'
        output:
            pandas dataframe (subset of frame)
            pixel-level features (subsetted)
            fc6 features (subsetted)
        '''
        kind2col = {'invalid': 'isInvalid', 'outlier': 'isOutlier'} # mapping between keyword arg and column name in frame
        for kind in kinds:
            self.metadata = self.metadata[self.metadata[kind2col[kind]] != True] if kinds[kind] else self.metadata
        self.pixels = self.pixels[self.metadata.index]
        self.fc6 = self.fc6[self.metadata.index]
        self.metadata = self.metadata.reset_index(drop=True)
    def preprocess(self, mean=True, std=False):
        '''
        input:
            mean = boolean, if true, will normalize features to mean = 0
            std  = boolean, if true, will normalize features to std = 1
        output:
            updates feature parameters to be normalized/standardized
        '''
        stds = self.pixels.std(axis = 0)
        stds[stds == 0.] = 1
        means = self.pixels.mean(axis = 0)
        self.pixels = self.pixels - means if mean else self.pixels
        self.pixels = self.pixels / stds if std else self.pixels

        stds = self.fc6.std(axis = 0)
        stds[stds == 0.] = 1
        means = self.fc6.mean(axis = 0)
        self.fc6 = self.fc6 - means if mean else self.fc6
        self.fc6 = self.fc6 / stds if std else self.fc6

        if std == True:
            assert round(np.sum(self.pixels.std(axis=0)) + np.sum((self.pixels.std(axis=0) == 0.))) == len(self.pixels[0])
            assert round(np.sum(self.fc6.std(axis=0)) + np.sum((self.fc6.std(axis=0) == 0.))) == len(self.fc6[0])
        if mean == True:
            assert round(np.sum(self.pixels.mean(axis=0))) == 0
            assert round(np.sum(self.fc6.mean(axis=0))) == 0


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

#/*##########################################################################
# Copyright (C) 2016 K. Kummer, A. Tamborino, European Synchrotron Radiation
# Facility
#
# This file is part of the ID32 RIXSToolBox developed at the ESRF by the ID32
# staff and the ESRF Software group.
#
# 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.
#
#############################################################################*/

from __future__ import division, print_function

__author__ = "K. Kummer - ESRF ID32"
__contact__ = "kurt.kummer@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
___doc__ = """
    ...
"""

import numpy as np

class RTB_Math(object):
    def gaussian_kernel_1d(self, sigma, order, radius):
        """
        Computes a 1D Gaussian convolution kernel.
        """
        if order < 0:
            raise ValueError('order must be non-negative')
        p = np.polynomial.Polynomial([0, 0, -0.5 / (sigma * sigma)])
        x = np.arange(-radius, radius + 1)
        phi_x = np.exp(p(x), dtype=np.double)
        phi_x /= phi_x.sum()
        if order > 0:
            q = np.polynomial.Polynomial([1])
            p_deriv = p.deriv()
            for _ in range(order):
                # f(x) = q(x) * phi(x) = q(x) * exp(p(x))
                # f'(x) = (q'(x) + q(x) * p'(x)) * phi(x)
                q = q.deriv() + q * p_deriv
            phi_x *= q(x)
        return phi_x

    def gaussian_filter(self, y, sigma, order=0, truncate=4.0):
        """
        Applies a Gaussian filter to a 1D input array.
        Adapted from scipy.ndimage
        """
        if abs(sigma) < 1e-15:
            retval = y
        else:
            sd = float(sigma)
            lw = int(truncate * sd + 0.5)
            weights = self.gaussian_kernel_1d(sigma, order, lw)[::-1]
            ytemp = np.hstack((len(weights)//len(y)+1)*[y, y[::-1],y, y[::-1]]
                + [y] + (len(weights)//len(y)+1)*[y[::-1], y, y[::-1], y])
            retval = np.correlate(ytemp, weights, mode='same')
            retval = retval[(len(weights)//len(y)+1)*4*len(y):
                ((len(weights)//len(y)+1)*4+1)*len(y)]
            retval[:len(weights)//2] = y[:len(weights)//2]
            retval[-len(weights)//2:] = y[-len(weights)//2:]
        return retval

    def gaussian(self, x, x0, amp, fwhm):
        return amp*np.exp(-(x-x0)**2/2/(fwhm/2.3548)**2)

    def minimize(self, func, x0, args=(), xatol=1e-4, fatol=1e-4, maxiter=None,
                    maxfev=None, disp=1):
        """
        Minimize a function using the downhill simplex algorithm.
        Adapted from scipy.optimize
        """

        def wrap_function(function, args):
            ncalls = [0]
            if function is None:
                return ncalls, None

            def function_wrapper(*wrapper_args):
                ncalls[0] += 1
                return function(*(wrapper_args + args))

            return ncalls, function_wrapper

        maxfun = maxfev
        fcalls, func = wrap_function(func, args)

        rho = 1
        chi = 2
        psi = 0.5
        sigma = 0.5

        nonzdelt = 0.05
        zdelt = 0.00025

        x0 = np.asfarray(x0).flatten()

        N = len(x0)

        sim = np.zeros((N + 1, N), dtype=x0.dtype)
        sim[0] = x0
        for k in range(N):
            y = np.array(x0, copy=True)
            if y[k] != 0:
                y[k] = (1 + nonzdelt)*y[k]
            else:
                y[k] = zdelt
            sim[k + 1] = y

        # If neither are set, then set both to default
        if maxiter is None and maxfun is None:
            maxiter = N * 200
            maxfun = N * 200
        elif maxiter is None:
            # Convert remaining Nones, to np.inf, unless the other is np.inf, in
            # which case use the default to avoid unbounded iteration
            if maxfun == np.inf:
                maxiter = N * 200
            else:
                maxiter = np.inf
        elif maxfun is None:
            if maxiter == np.inf:
                maxfun = N * 200
            else:
                maxfun = np.inf

        one2np1 = list(range(1, N + 1))
        fsim = np.zeros((N + 1,), float)

        for k in range(N + 1):
            fsim[k] = func(sim[k])

        ind = np.argsort(fsim)
        fsim = np.take(fsim, ind, 0)
        # sort so sim[0,:] has the lowest function value
        sim = np.take(sim, ind, 0)

        iterations = 1

        while (fcalls[0] < maxfun and iterations < maxiter):
            if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
                    np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
                break

            xbar = np.add.reduce(sim[:-1], 0) / N
            xr = (1 + rho) * xbar - rho * sim[-1]
            fxr = func(xr)
            doshrink = 0

            if fxr < fsim[0]:
                xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
                fxe = func(xe)

                if fxe < fxr:
                    sim[-1] = xe
                    fsim[-1] = fxe
                else:
                    sim[-1] = xr
                    fsim[-1] = fxr
            else:  # fsim[0] <= fxr
                if fxr < fsim[-2]:
                    sim[-1] = xr
                    fsim[-1] = fxr
                else:  # fxr >= fsim[-2]
                    # Perform contraction
                    if fxr < fsim[-1]:
                        xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
                        fxc = func(xc)

                        if fxc <= fxr:
                            sim[-1] = xc
                            fsim[-1] = fxc
                        else:
                            doshrink = 1
                    else:
                        # Perform an inside contraction
                        xcc = (1 - psi) * xbar + psi * sim[-1]
                        fxcc = func(xcc)

                        if fxcc < fsim[-1]:
                            sim[-1] = xcc
                            fsim[-1] = fxcc
                        else:
                            doshrink = 1

                    if doshrink:
                        for j in one2np1:
                            sim[j] = sim[0] + sigma * (sim[j] - sim[0])
                            fsim[j] = func(sim[j])

            ind = np.argsort(fsim)
            sim = np.take(sim, ind, 0)
            fsim = np.take(fsim, ind, 0)
            iterations += 1

        x = sim[0]
        fval = np.min(fsim)
        warnflag = 0

        if fcalls[0] >= maxfun:
            warnflag = 1
            msg = 'Maximum number of function evaluations has been exceeded.'
            if disp:
                print('Warning: ' + msg)
        elif iterations >= maxiter:
            warnflag = 2
            msg = 'Maximum number of iterations has been exceeded.'
            if disp:
                print('Warning: ' + msg)
        else:
            msg = 'Optimization terminated successfully.'
            if disp:
                print(msg)
                print("         Current function value: %f" % fval)
                print("         Iterations: %d" % iterations)
                print("         Function evaluations: %d" % fcalls[0])

        result = dict(fun=fval, nit=iterations, nfev=fcalls[0],
                                status=warnflag, success=(warnflag == 0),
                                message=msg, x=x, final_simplex=(sim, fsim))

        return x, result

    def curve_fit(self, theory, x, y, p0=(1, 1, 1)):
        """
        Covinience function for curve fitting.
        """
        def func2minimize(p0, theory, data):
            return ((data[:,1]-theory(data[:,0], *p0))**2).sum()

        popt, msg = self.minimize(func2minimize, x0=p0,
                    args=(theory, np.vstack([x, y]).T), disp=False)
        return popt, msg

    def interpolate_on_grid(self, qvals, xvals, yvals, grid, method='nearest',
            fill_value=None):
        """
        Interpolate one-dimensional irregular data on a regular, two-dimensional
        grid. Mimics the functionality of scipy.interpolate.griddata.
        """

        gridx, gridy = grid
        qvals = np.array(qvals)

        gridz = np.vstack([np.interp(gridx[:,0], xvals[i], yvals[i],
            left=fill_value, right=fill_value) for i in range(len(qvals))]).T
        if method == 'nearest':
            gridz = np.vstack([gridz[:,np.argmin(np.abs(qvals-q))]
                                for q in gridy[0,:]]).T
        else:
            idx = np.argsort(qvals)
            qvals = qvals[idx]
            gridz = gridz[:,idx]
            gridz = np.vstack([np.interp(gridy[0,:], qvals, gridz[i,:],
                                left=fill_value, right=fill_value)
                                for i in range(len(gridz[:,0]))])
        return gridz

if __name__ == "__main__":
    try:
        import matplotlib.pyplot as plt
        PLOT = True
    except ImportError:
        print('Failed to import matplotlib. Skip plotting.')
        PLOT = False

    a = RTB_Math()
    x = np.linspace(-1, 1, 101)
    y = 3*(np.random.random(len(x))-0.5) + a.gaussian(x, 0, 25, .2)
    ys = a.gaussian_filter(y, 5)

    def func2minimize(p0, theory, data):
        return ((data[:,1]-theory(data[:,0], *p0))**2).sum()

    popt, msg = a.minimize(func2minimize, x0=[1,1,1],
                    args=(a.gaussian, np.vstack([x, y]).T), disp=True)

    if PLOT:
        plt.figure().add_subplot(111)
        plt.errorbar(x, y, yerr=3, marker='.', ls='none', label='Raw data')
        plt.plot(x, ys, label='Smoothed')
        plt.plot(x, a.gaussian(x, *popt), label='fit\n(%.3f, %.3f, %.3f)' % tuple(popt))
        plt.legend()

        plt.show()

max_stars_repo_name: kkummer/RixsToolBox

max_issues_repo_name: kkummer/RixsToolBox

max_forks_repo_name: kkummer/RixsToolBox

reason: import numpy,from scipy


text:
import numpy as np
import string
import os,re
import random

to1letter = {
    "ALA":'A', "ARG":'R', "ASN":'N', "ASP":'D', "CYS":'C',
    "GLN":'Q', "GLU":'E', "GLY":'G', "HIS":'H', "ILE":'I',
    "LEU":'L', "LYS":'K', "MET":'M', "PHE":'F', "PRO":'P',
    "SER":'S', "THR":'T', "TRP":'W', "TYR":'Y', "VAL":'V' }

# read A3M and convert letters into
# integers in the 0..20 range,
def parse_a3m(filename):

    msa = []

    table = str.maketrans(dict.fromkeys(string.ascii_lowercase))

    # read file line by line
    for line in open(filename,"r"):

        # skip labels
        if line[0] == '>':
            continue

        # remove right whitespaces
        line = line.rstrip()

        # remove lowercase letters and append to MSA
        msa.append(line.translate(table))

    # convert letters into numbers
    alphabet = np.array(list("ARNDCQEGHILKMFPSTWYV-"), dtype='|S1').view(np.uint8)
    msa = np.array([list(s) for s in msa], dtype='|S1').view(np.uint8)
    for i in range(alphabet.shape[0]):
        msa[msa == alphabet[i]] = i

    # treat all unknown characters as gaps
    msa[msa > 20] = 20

    return msa


text:
import numpy as np
import sophus as sp

class Frame(object):
    """docstring for Frame"""
    factory_id = -1

    def __init__(self, nid, time_stamp=0, Tcw=None, camera=None, color=None, depth=None):
        self.id = nid
        self.time_stamp = time_stamp
        self.Tcw = Tcw if Tcw is not None else sp.SE3()
        self.camera = camera
        self.color = color
        self.depth = depth

    @staticmethod
    def createFrame():
        Frame.factory_id += 1
        return Frame(Frame.factory_id)

    def getCamCenter(self):
        return self.Tcw.inverse().translation()

    def isInFrame(self, pt_world):
        p_cam = self.camera.world2camera(pt_world, self.Tcw)
        if p_cam[2] < 0:
            return False
        pixel = self.camera.world2camera(pt_world, self.Tcw)

        lower_left = np.array([0, 0])
        uper_right = self.color.shape
        return np.logical_and(lower_left < pixel, pixel < uper_right).all()

    def findDepth(self, kp):
        x, y = np.around(kp.pt).astype(int)
        d = self.depth[y, x]
        if d != 0:
            return d / self.camera.depth_scale
        else:
            # check the nearby points
            dx = [-1, 0, 1, 0]
            dy = [0, -1, 0, 1]

            for i in range(4):
                d = self.depth[y + dy[i], x + dx[i]]
                if d != 0: return d / self.camera.depth_scale
        return -1


text:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from sklearn.mixture import BayesianGaussianMixture

def main():
    cluster_nums = [2,4,8,12]
    samples_nums = [100,500,900]
    dimensions = [2,3,4]
    k_range = list(range(1,16))

    result_aic = []
    result_bic = []
    result_vbem = []
    cond_cluster = []
    cond_dimension = []
    cond_sample = []

    np.random.seed(0)
    cluster_data = []
    for dimension in dimensions:
        for samples_num in samples_nums:
            for cluster_num in cluster_nums:
                data = []
                for i in range(cluster_num):
                    cluster_data += [np.random.randn(np.random.randint(100,100+samples_num), dimension) + np.random.randint(10,50,size=(1,dimension))]
                    data.extend(cluster_data[i])
                data = np.array(data)

                lowest_aic = np.infty
                best_aic_k = 0
                aic = []

                lowest_bic = np.infty
                best_bic_k = 0
                bic = []

                lowest_bound = np.infty
                best_vbem_k = 0
                vbem = []

                for k in k_range:
                    gmm = GaussianMixture(n_components=k, covariance_type='full')
                    gmm.fit(data)
                    aic.append(gmm.aic(data))
                    bic.append(gmm.bic(data))
                    if aic[-1] < lowest_aic:
                        lowest_aic = aic[-1]
                        best_aic_k = k
                    if bic[-1] < lowest_bic:
                        lowest_bic = bic[-1]
                        best_bic_k = k

                    vbgmm = BayesianGaussianMixture( weight_concentration_prior_type="dirichlet_distribution", n_components=k,covariance_type='full')
                    vbgmm.fit(data)
                    vbem.append(vbgmm.lower_bound_ * (-1))
                    if vbem[-1] < lowest_bound:
                        lowest_bound = vbem[-1]
                        best_vbem_k = k

                cond_cluster += [cluster_num]
                cond_dimension += [dimension]
                cond_sample += ['<' + str(samples_num)]
                result_aic += [best_aic_k]
                result_bic += [best_bic_k]
                result_vbem += [best_vbem_k]

                print(str(cluster_num) + '\t\t' + str(dimension) + '\t\t<' + str(100 + samples_num) + '\t\t' + str(best_aic_k) + '\t\t' + str(best_bic_k) + '\t\t' + str(best_vbem_k))
                plt.figure(str(cluster_num) + '-cluster ' + str(dimension) + 'D ' + 'with random <' + str(100 + samples_num) + ' samples')
                plt.title(str(cluster_num) + ' clusters ' + str(dimension) + 'D ' + ' <' + str(100 + samples_num) + ' samples')
                plt.plot(k_range, aic, 'b-o', label='AIC Optimal k = ' + str(best_aic_k))
                plt.plot(k_range, bic, 'r-o', label='BIC Optimal k = ' + str(best_bic_k))
                plt.plot(k_range, vbem, 'g-o', label='VBEM\'s (minus) lower bound Optimal k = ' + str(best_vbem_k))
                plt.legend(loc='upper right')
                plt.savefig('./'+ str(cluster_num) + 'cluster_' + str(dimension) + 'D_' + 'with_random_' + str(100 + samples_num) + '_samples')
                # plt.show()
    cond_cluster = np.array(cond_cluster)
    result_aic = np.array(result_aic)
    result_bic = np.array(result_bic)
    result_vbem = np.array(result_vbem)
    error_aic = result_aic - cond_cluster
    error_bic = result_bic - cond_cluster
    error_vbem = result_vbem - cond_cluster
    error_ratio_aic = sum(list(map(lambda x: 1 if x != 0 else 0,(result_aic - cond_cluster)))) * 1.0 / len(cond_cluster)
    error_ratio_bic = sum(list(map(lambda x: 1 if x != 0 else 0,(result_bic - cond_cluster)))) * 1.0 / len(cond_cluster)
    error_ratio_vbem = sum(list(map(lambda x: 1 if x != 0 else 0,(result_vbem - cond_cluster)))) * 1.0 / len(cond_cluster)

    print(cond_cluster)
    print(cond_dimension)
    print(cond_sample)
    print(result_aic)
    print(result_bic)
    print(result_vbem)
    print(error_aic)
    print(error_bic)
    print(error_vbem)
    print(error_ratio_aic)
    print(error_ratio_bic)
    print(error_ratio_vbem)

main()

max_stars_repo_name: AlexKnightDu/Lessons

max_issues_repo_name: AlexKnightDu/Lessons

max_forks_repo_name: AlexKnightDu/Lessons