text:
"""
This file covers relevant cases for non-extremal spin SC/Kerr holes (a < 1).

This file is translated from ConstantsOfMotion.m from the Black Hole
Perturbation Toolkit by Niels Warburton et al.
"""

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

from mpmath import sqrt, mp
from sys import exit

# ------------------------------------------------------------------------------
#  SC orbits (a = 0)
# ------------------------------------------------------------------------------

def calc_sc_constants(slr, ecc, x):
    """
    Energy, angular momentum, and carter constant calculation.

    Schwarzschild case (spin parameter a = 0)

    Parameters:
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc) (0, 1]
                   negative x -> retrograde
                   positive x -> prograde

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    try:
        prec = mp.prec
        mp.prec += 10

        ecc2 = ecc * ecc
        slr2 = slr * slr
        x2 = x * x

        mp.prec += 10
        En = sqrt(
            (-4 * ecc2 + (-2 + slr)**2) /
            (slr * (-3 - ecc2 + slr))
        )
        Lz = (slr * x) / sqrt(-3 - ecc2 + slr)
        Q = (slr2 * (-1 + x2)) / (3 + ecc2 - slr)

        return En, Lz, Q
    finally:
        mp.prec = prec

# ------------------------------------------------------------------------------
#  Kerr equatorial orbits (x = +/- 1) -> sign is pro/retrograde motion
# ------------------------------------------------------------------------------

def eq_energy(aa, slr, ecc, x):
    """
    Energy for the Kerr equatorial case (inclination value x = 1).

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [~2, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x = -1 -> retrograde
                   x = +1 -> prograde

    Returns:
        En (mpf): energy
    """
    try:
        prec = mp.prec
        mp.prec += 10

        ecc2 = ecc * ecc
        eta = ecc2 - 1
        eta2 = eta * eta
        aa2 = aa * aa
        aa4 = aa2 * aa2
        aa6 = aa4 * aa2
        mp.prec += 10
        return (
            sqrt(
                1 - (-eta * (1 + (eta * (aa2 * (1 + 3 * ecc2 + slr) +
                     slr * (-3 - ecc2 + slr -
                     2 * sqrt((aa6 * eta2 +
                              aa2 * (-4 * ecc2 + (-2 + slr)**2) * slr**2 +
                              2 * aa4 * slr * (-2 + slr + ecc2 * (2 + slr))) /
                             (slr**3 * x**2)) * x))) /
                    (-4 * aa2 * eta2 + (3 + ecc2 - slr)**2 * slr))) / slr)
        )
    finally:
        # print('here')
        mp.prec = prec

def eq_ang_momentum(aa, slr, ecc, x):
    """
    Angular momentum for the Kerr equatorial case (inclination value x = 1).

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x = -1 -> retrograde
                   x = +1 -> prograde

    Returns:
        Lz (mpf): angular momentum
    """
    try:
        prec = mp.prec
        mp.prec += 10

        ecc2 = ecc * ecc
        eta = ecc2 - 1
        eta2 = eta * eta
        aa2 = aa * aa
        aa4 = aa2 * aa2
        aa6 = aa4 * aa2
        slr2 = slr * slr
        slr3 = slr2 * slr
        x2 = x * x
        mp.prec += 10
        num_root = slr * (-3 - ecc2 + slr - 2 *
                          sqrt((aa6 * eta2 + aa2 *
                               (-4 * ecc2 + (-2 + slr)**2) *
                               slr2 + 2 * aa4 * slr *
                               (-2 + slr + ecc2 * (2 + slr))) /
                               (slr3 * x2)) * x)
        denom = (-4 * aa2 * eta2 + (3 + ecc2 - slr)**2 * slr)
        return (
            slr * x *
            sqrt((aa2 * (1 + 3 * ecc2 + slr) + num_root) / (denom * x2)) +
            aa * sqrt(1 - (-eta * (1 + (eta * (aa2 * (1 + 3 * ecc2 + slr) +
                      num_root)) / denom)) / slr)
        )
    finally:
        mp.prec = prec

def calc_eq_constants(aa, slr, ecc, x):
    """
    Calculate equatorial constants in one function.

    Note that for the equatorial case, the Carter constant, Q = 0.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x = -1 -> retrograde
                   x = +1 -> prograde

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    En = eq_energy(aa, slr, ecc, x)
    Lz = eq_ang_momentum(aa, slr, ecc, x)
    Q = 0  # equatorial case
    return En, Lz, Q

# ------------------------------------------------------------------------------
#  Kerr spherical orbits (ecc = 0)
# ------------------------------------------------------------------------------

def spherical_energy(aa, slr, x):
    """
    Compute energy for the Kerr spherical case (ecc = 0).

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        x (mpf): inclination value given by cos(theta_inc) (0, 1]
                   negative x -> retrograde
                   positive x -> prograde

    Returns:
        En (mpf): energy
    """
    try:
        prec = mp.prec
        mp.prec += 10

        slrm2 = (-2 + slr)**2
        slr2 = slr * slr
        slr3 = slr * slr2
        slr4 = slr * slr3
        slr5 = slr3 * slr2
        slr7 = slr5 * slr2
        x2 = x * x
        x3 = x2 * x
        aa2 = aa * aa
        aa3 = aa2 * aa
        aa4 = aa2 * aa2
        aa5 = aa4 * aa
        aa6 = aa4 * aa2
        x2m1 = x2 - 1
        x2m12 = x2m1 * x2m1

        mp.prec += 10
        return (
            sqrt(((-3 + slr) * slrm2 * slr5 -
                 2 * aa5 * x * x2m1 * sqrt(slr3 + aa2 * slr * x2m1) +
                 aa4 * slr2 * x2m1 * (4 - 5 * slr * (-1 + x2) +
                 3 * slr2 * x2m1) - aa6 * x2m12 *
                 (x2 + slr2 * x2m1 - slr * (1 + 2 * x2)) +
                 aa2 * slr3 * (4 - 4 * x2 + slr * (12 - 7 * x2) -
                               3 * slr3 * (-1 + x2) +
                               slr2 * (-13 + 10 * x2)) +
                 aa * (-2 * slr**4.5 * x * sqrt(slr2 + aa2 * x2m1) +
                       4 * slr3 * x * sqrt(slr3 + aa2 * slr * x2m1)) +
                 2 * aa3 * (2 * slr * x * x2m1 *
                            sqrt(slr3 + aa2 * slr * x2m1) -
                            x3 * sqrt(slr7 + aa2 * slr5 * x2m1))) /
                 ((slr2 - aa2 * x2m1) *
                 ((-3 + slr)**2 * slr4 - 2 * aa2 * slr2 *
                  (3 + 2 * slr - 3 * x2 + slr2 * x2m1) +
                  aa4 * x2m1 * (-1 + x2 + slr2 * x2m1 -
                  2 * slr * (1 + x2)))))
        )
    finally:
        mp.prec = prec

def spherical_ang_momentum(En, aa, slr, x):
    """
    Compute angular momentum for the Kerr spherical case (ecc = 0).

    Parameters:
        En (mpf): energy
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        x (mpf): inclination value given by cos(theta_inc) (0, 1]
                   negative x -> retrograde
                   positive x -> prograde

    Returns:
        Lz (mpf): angular momentum
    """
    try:
        prec = mp.prec
        mp.prec += 10

        aa2 = aa * aa
        slr2 = slr * slr
        x2 = x * x
        mp.prec += 10
        g = 2 * aa * slr
        d = (aa2 + (-2 + slr) * slr) * (slr2 - aa2 * (-1 + x2))
        h = ((-2 + slr) * slr - aa2 * (-1 + x2)) / x2
        f = (slr**4 + aa2 * (slr * (2 + slr) - (aa2 + (-2 + slr) * slr) *
             (-1 + x2)))

        return (-(En * g) +
                sqrt((-(d * h) + En**2 * (g**2 + f * h)) / x2) * x) / h
    finally:
        mp.prec = prec

# ------------------------------------------------------------------------------
#  Generic orbit constant calculation
# ------------------------------------------------------------------------------

def calc_delta(r, aa):
    """
    Calculate ubiquitous function on Kerr spacetimes.

    Parameters:
        r (mpf): radius
        aa (mpf): spin parameter (0, 1)

    Returns:
        delta (mpf)
    """
    try:
        prec = mp.prec
        mp.prec += 20
        # return r * r - 2 * r + aa * aa
        return r * (r - 2) + aa * aa
    finally:
        mp.prec = prec

def calc_f(r, zm, aa):
    """
    DESCRIPTION TBD (FUNCTION USED IN GENERAL CONSTANTS).

    Parameters:
        r (mpf): radius
        zm (mpf): 1 - x * x
        aa (mpf): spin parameter (0, 1)

    Returns:
        f (mpf)
    """
    try:
        prec = mp.prec
        mp.prec += 10
        r2 = r * r
        mp.prec += 10
        r4 = r2 * r2
        aa2 = aa * aa
        zm2 = zm * zm

        delta = calc_delta(r, aa)
        mp.prec += 10
        return 2 * aa2 * r + aa2 * r2 + r4 + aa2 * zm2 * delta
    finally:
        mp.prec = prec

def calc_g(r, aa):
    """
    DESCRIPTION TBD (FUNCTION USED IN GENERAL CONSTANTS).

    Parameters:
        r (mpf): radius
        aa (mpf): spin parameter (0, 1)

    Returns:
        g (mpf)
    """
    try:
        prec = mp.prec
        mp.prec += 10
        return 2 * aa * r
    finally:
        mp.prec = prec

def calc_h(r, zm, aa):
    """
    DESCRIPTION TBD (FUNCTION USED IN GENERAL CONSTANTS).

    Parameters:
        r (mpf): radius
        zm (mpf): 1 - x * x
        aa (mpf): spin parameter (0, 1)

    Returns:
        h (mpf)
    """
    try:
        prec = mp.prec
        mp.prec += 10
        zm2 = zm * zm
        delta = calc_delta(r, aa)
        mp.prec += 10
        return (-2 + r) * r + (zm2 * delta) / (1 - zm2)
    finally:
        mp.prec = prec

def calc_d(r, zm, aa):
    """
    DESCRIPTION TBD (FUNCTION USED IN GENERAL CONSTANTS).

    Parameters:
        r (mpf): radius
        zm (mpf): 1 - x * x
        aa (mpf): spin parameter (0, 1)

    Returns:
        d (mpf)
    """
    try:
        prec = mp.prec
        mp.prec += 10
        r2 = r * r
        aa2 = aa * aa
        zm2 = zm * zm

        delta = calc_delta(r, aa)
        mp.prec += 10
        return (r2 + aa2 * zm2) * delta
    finally:
        mp.prec = prec

def gen_energy(zm, aa, slr, ecc, x):
    """
    Compute energy for generic orbit case.

    Parameters:
        zm (mpf): 1 - x * x
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [~2, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        En (mpf): energy
    """
    try:
        prec = mp.prec
        mp.prec += 20
        r1 = slr / (1 - ecc)
        r2 = slr / (1 + ecc)

        dr1 = calc_d(r1, zm, aa)
        dr2 = calc_d(r2, zm, aa)
        gr1 = calc_g(r1, aa)
        gr2 = calc_g(r2, aa)
        hr1 = calc_h(r1, zm, aa)
        hr2 = calc_h(r2, zm, aa)
        fr1 = calc_f(r1, zm, aa)
        fr2 = calc_f(r2, zm, aa)

        mp.prec += 10
        kappa = dr1 * hr2 - hr1 * dr2
        epsilon = dr1 * gr2 - gr1 * dr2
        rho = fr1 * hr2 - hr1 * fr2
        eta = fr1 * gr2 - gr1 * fr2
        sigma = gr1 * hr2 - hr1 * gr2

        mp.prec += 10
        kappa2 = kappa * kappa
        epsilon2 = epsilon * epsilon
        rho2 = rho * rho
        x2 = x * x

        mp.prec += 10
        En = sqrt(
            (kappa * rho + 2 * epsilon * sigma -
             2 * sqrt((sigma * (-(eta * kappa2) + epsilon * kappa * rho +
                      epsilon2 * sigma)) / x2) * x) / (rho2 + 4 * eta * sigma))

        return En
    finally:
        mp.prec = prec

def gen_ang_momentum(En, aa, slr, ecc, x):
    """
    Compute angular momentum for generic orbit case.

    Parameters:
        En (mpf): energy
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        z_ang_mom (mpf): angular momentum
    """
    try:
        prec = mp.prec
        mp.prec += 10
        r1 = slr / (1 - ecc)
        zm = sqrt(1 - x * x)

        # // delta_r1 = calc_delta(r1, aa)
        mp.prec += 10
        fr1 = calc_f(r1, zm, aa)
        gr1 = calc_g(r1, aa)
        hr1 = calc_h(r1, zm, aa)
        dr1 = calc_d(r1, zm, aa)

        mp.prec += 10
        En2 = En * En
        x2 = x * x
        mp.prec += 10
        Lz = ((-(En * gr1) + x * sqrt((-(dr1 * hr1) + En2 *
                     (gr1 * gr1 + fr1 * hr1)) / x2)) / hr1)
        return Lz
    finally:
        mp.prec = prec

def gen_carter_const(En, Lz, aa, slr, ecc, x):
    """
    Compute Carter constant for generic orbit case.

    Parameters:
        En (mpf): energy
        Lz (mpf): angular momentum
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        Q (mpf): Carter constant
    """
    try:
        prec = mp.prec
        mp.prec += 10
        zm = sqrt(1 - x * x)
        mp.prec += 10
        zm2 = zm * zm
        mp.prec += 10
        return (zm2 * (aa * aa * (1 - En * En) + Lz * Lz /
                (1 - zm2)))
    finally:
        mp.prec = prec

def calc_gen_constants(aa, slr, ecc, x):
    """
    Call generic orbit constant calculating functions in one function.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    try:
        prec = mp.prec
        mp.prec += 30
        zm = sqrt(1 - x * x)

        En = gen_energy(zm, aa, slr, ecc, x)
        Lz = gen_ang_momentum(En, aa, slr, ecc, x)
        Q = gen_carter_const(En, Lz, aa, slr, ecc, x)
        return En, Lz, Q
    finally:
        mp.prec = prec

def pol_energy(aa, slr, ecc):
    """
    Calculate energy for polar Kerr case.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity

    Returns:
        En (mpf): energy
    """
    try:
        prec = mp.prec
        mp.prec += 10
        aa2 = aa * aa
        aa4 = aa2 * aa2
        ecc2 = ecc * ecc
        ecc4 = ecc2 * ecc2
        slr2 = slr * slr
        slr4 = slr2 * slr2
        mp.prec += 10
        return sqrt(-((slr * (aa4 * (-1 + ecc2)**2 +
                    (-4 * ecc2 + (-2 + slr)**2) * slr2 +
                    2 * aa2 * slr * (-2 + slr + ecc2 * (2 + slr)))) /
                    (aa4 * (-1 + ecc2)**2 * (-1 + ecc2 - slr) +
                     (3 + ecc2 - slr) * slr4 - 2 * aa2 * slr2 *
                     (-1 - ecc4 + slr + ecc2*(2 + slr)))))
    finally:
        mp.prec = prec

def pol_carter_const(aa, slr, ecc):
    """
    Calculate energy for polar Kerr case.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity

    Returns:
        Q (mpf): Carter constant
    """
    try:
        prec = mp.prec
        mp.prec += 10
        aa2 = aa * aa
        aa4 = aa2 * aa2
        ecc2 = ecc * ecc
        ecc4 = ecc2 * ecc2
        slr2 = slr * slr
        slr4 = slr2 * slr2
        mp.prec += 10
        return -((slr2 * (aa4 * (-1 + ecc2)**2 + slr4 +
                 2 * aa2 * slr * (-2 + slr + ecc2 * (2 + slr)))) /
                 (aa4 * (-1 + ecc2)**2 *
                  (-1 + ecc2 - slr) + (3 + ecc2 - slr) * slr4 -
                 2 * aa2 * slr2 * (-1 - ecc4 + slr + ecc2 * (2 + slr))))
    finally:
        mp.prec = prec

def calc_sph_constants(aa, slr, x):
    """
    Call spherical orbit (ecc = 0) functions.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    ecc = 0

    En = spherical_energy(aa, slr, x)
    Lz = spherical_ang_momentum(En, aa, slr, x)
    Q = gen_carter_const(En, Lz, aa, slr, ecc, x)

    return En, Lz, Q

def calc_pol_constants(aa, slr, ecc):
    """
    Call Kerr polar orbit (x = 0) functions.

    Parameters:
        aa (mpf): spin parameter [0, 1)
        slr (mpf): semi-latus rectum
        ecc (mpf): eccentricity [0, 1)

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    Lz = 0

    En = pol_energy(aa, slr, ecc)
    Q = pol_carter_const(aa, slr, ecc)

    return En, Lz, Q

def mp_calc_constants(aa, slr, ecc, x):
    """
    Choose which function to call based on input parameters.

    This version uses mpmath for extended precision.

    Parameters:
        aa (mpf): spin parameter (0, 1)
        slr (mpf): semi-latus rectum [6, inf)
        ecc (mpf): eccentricity [0, 1)
        x (mpf): inclination value given by cos(theta_inc)
                   x < 0 -> retrograde
                   x > 0 -> prograde

    Returns:
        En (mpf): energy
        Lz (mpf): angular momentum
        Q (mpf): Carter constant
    """
    if aa == 0:
        if slr < 6:
            print("slr value is too small for the SC case.")
            exit("Error in input values.")
        return calc_sc_constants(slr, ecc, x)
    elif x == 0:
        return calc_pol_constants(aa, slr, ecc)
    elif x**2 == 1:
        return calc_eq_constants(aa, slr, ecc, x)
    elif ecc == 0:
        return calc_sph_constants(aa, slr, x)
    else:
        return calc_gen_constants(aa, slr, ecc, x)

max_stars_repo_name: AaronDJohnson/renormnu

max_issues_repo_name: AaronDJohnson/fbtpoint

max_forks_repo_name: AaronDJohnson/fbtpoint

reason: from mpmath


text:
#!/usr/bin/env python

'''Take a photo.

Take a photo using a USB or Raspberry Pi camera.
'''

from __future__ import print_function
import os
import sys
from time import time, sleep
import subprocess

WIDTH = os.getenv('take_photo_width', '640')
HEIGHT = os.getenv('take_photo_height', '480')
CAMERA_DISABLED_MSG = 'No camera selected. Choose a camera on the device page.'

def _log(text):
    try:
        import json, socket, struct
        s = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM)
        r = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM)
        s.connect(os.environ['FARMWARE_API_V2_REQUEST_PIPE'])
        r.connect(os.environ['FARMWARE_API_V2_RESPONSE_PIPE'])
        message = bytes(json.dumps({
            'kind': 'rpc_request', 'args': {'label': ''},
            'body': [{
                'kind': 'send_message',
                'args': {'message': text, 'message_type': 'error'}}]}), 'utf-8')
        s.sendall(struct.pack('!Hii', 0xFBFB, 0, len(message)) + message)
        r.recv(10)
    except (KeyError, TypeError):
        std_print(text)
    finally:
        s.close()
        r.close()

try:
    MissingError = FileNotFoundError
except NameError:
    MissingError = OSError

def get_camera_selection():
    'Fetch camera type selected.'
    return os.getenv('camera', 'USB').upper()

def rotation_disabled():
    'Check if rotation is disabled via environment variable.'
    return '1' in os.getenv('take_photo_disable_rotation_adjustment', '1')

def std_print(text):
    'Print.'
    if not 'quiet' in os.getenv('take_photo_logging', '').lower():
        try:
            print(text, flush=True)
        except TypeError:
            print(text)

def get_video_port_list():
    'Get available video ports from /dev.'
    return [d for d in os.listdir('/dev') if d.startswith('video')]

def usb_camera_call(savepath):
    'Call fswebcam.'
    size = '{}x{}'.format(WIDTH, HEIGHT)
    args = ['fswebcam', '-r', size, '-S', '10', '--no-banner', savepath]
    std_print('Calling `{}`...'.format(' '.join(args)))
    try:
        return subprocess.call(args)
    except MissingError:
        return 1

def rpi_photo_call(savepath):
    'Call raspistill.'
    width = min(int(WIDTH), 4056)
    height = min(int(HEIGHT), 3040)
    size = ['-w', str(width), '-h', str(height)]
    if height > 1500:
        size = ['-md', '3']
    args = ['raspistill'] + size + ['-o', savepath]
    std_print('Calling `{}`...'.format(' '.join(args)))
    try:
        return subprocess.call(args)
    except MissingError:
        return 1

# Takes photo and exits if rotation was disabled via environment variable.
# Without imports, logs, or processing, this is a much quicker path.
if rotation_disabled():
    SAVEPATH = '/tmp/images/{}.jpg'.format(int(time()))
    SELECTED_CAMERA = get_camera_selection()
    return_code = 0
    if 'NONE' in SELECTED_CAMERA:
        _log(CAMERA_DISABLED_MSG)
    elif 'RPI' in SELECTED_CAMERA:
        return_code = rpi_photo_call(SAVEPATH)
    else:
        ports = get_video_port_list()
        if len(ports) < 1:
            _log('USB Camera not detected.')
            sys.exit(0)
        return_code = usb_camera_call(SAVEPATH)
    if return_code == 0:
        sys.exit(0)
    else:
        std_print('command not found. Trying OpenCV...')

# start timer
START_TIME = time()

import requests
import numpy as np

FIRST_IMPORTS_COMPLETE_TIME = time()

def verbose_log(text, time_override=None):
    'Print text with time elapsed since start.'
    now = time_override if time_override is not None else time()
    elapsed = round(now - START_TIME, 4)
    timed_log = '[{:>8}] {}'.format(elapsed, text)
    log_level = os.getenv('take_photo_logging', '').lower()
    if 'quiet' in log_level:
        return
    if 'verbose' not in log_level:
        std_print(timed_log)
        return
    log_content = timed_log if 'timed' in log_level else text
    try:
        log(log_content, 'debug')
    except NameError:
        pass

def _farmware_api_url():
    major_version = int(os.getenv('FARMBOT_OS_VERSION', '0.0.0')[0])
    base_url = os.environ['FARMWARE_URL']
    return base_url + 'api/v1/' if major_version > 5 else base_url

def legacy_log(message, message_type):
    'Send a message to the log.'
    try:
        os.environ['FARMWARE_URL']
    except KeyError:
        print(message)
    else:
        log_message = '[take-photo] ' + str(message)
        headers = {
            'Authorization': 'bearer {}'.format(os.environ['FARMWARE_TOKEN']),
            'content-type': 'application/json'}
        payload = {
            'kind': 'send_message',
            'args': {'message': log_message, 'message_type': message_type}}
        requests.post(_farmware_api_url() + 'celery_script',
                      json=payload, headers=headers)

try:
    ft_import_start_msg = 'Importing Farmware Tools...'
    FT_IMPORT_START_TIME = time()
    from farmware_tools import env, device
except ImportError:
    ft_import_result_msg = 'farmware_tools import error. Using legacy logger.'
    log = legacy_log
    IMAGES_DIR = os.getenv('IMAGES_DIR')
else:
    ft_import_result_msg = 'Farmware Tools import complete.'
    IMAGES_DIR = env.Env().images_dir

    def log(message, message_type):
        'Send a log message.'
        device.log('[take-photo] {}'.format(message), message_type)

verbose_log('First imports complete.', FIRST_IMPORTS_COMPLETE_TIME)
verbose_log(ft_import_start_msg, FT_IMPORT_START_TIME)
verbose_log(ft_import_result_msg)

try:
    verbose_log('Importing OpenCV...')
    os.environ['OPENCV_VIDEOIO_DEBUG'] = '1'
    import cv2
except ImportError:
    log('OpenCV import error.', 'error')
    sys.exit(0)
else:
    verbose_log('OpenCV import complete.')

def rotate(image):
    'Rotate image if calibration data exists.'
    if rotation_disabled():
        raise KeyError('Rotation disabled.')
    angle = float(os.environ['CAMERA_CALIBRATION_total_rotation_angle'])
    sign = -1 if angle < 0 else 1
    turns, remainder = -int(angle / 90.), abs(angle) % 90  # 165 --> -1, 75
    if remainder > 45: turns -= 1 * sign  # 75 --> -1 more turn (-2 turns total)
    angle += 90 * turns                   #        -15 degrees
    image = np.rot90(image, k=turns)
    height, width, _ = image.shape
    matrix = cv2.getRotationMatrix2D((int(width / 2), int(height / 2)), angle, 1)
    return cv2.warpAffine(image, matrix, (width, height))

def image_filename():
    'Prepare filename with timestamp.'
    epoch = int(time())
    filename = '{timestamp}.jpg'.format(timestamp=epoch)
    return filename

def upload_path(filename):
    'Filename with path for uploading an image.'
    images_dir = IMAGES_DIR or '/tmp/images'
    if not os.path.isdir(images_dir):
        log('{} directory does not exist.'.format(images_dir), 'error')
    path = images_dir + os.sep + filename
    return path

def save_image(image):
    'Save an image to file after attempting rotation.'
    filename = image_filename()
    # Try to rotate the image
    try:
        verbose_log('Considering rotation...')
        final_image = rotate(image)
    except:
        verbose_log('Did not rotate image.')
        final_image = image
    else:
        verbose_log('Rotated image.')
        filename = 'rotated_' + filename
    # Save the image to file
    filename_path = upload_path(filename)
    cv2.imwrite(filename_path, final_image)
    verbose_log('Image saved: {}'.format(filename_path))

def _get_usb_device_list():
    try:
        raw_usb_results = subprocess.check_output(['lsusb'])
    except MissingError:
        verbose_log('Unable to check USB devices.')
        device_list_str = ''
    except subprocess.CalledProcessError:
        verbose_log('USB device check error.')
        device_list_str = ''
    else:
        usb_results = raw_usb_results.decode().strip().split('\n')
        usb_devices = [result.strip()[28:] for result in usb_results]
        usb_list_str = '|'.join(usb_devices)
        verbose_log('{} USB device entries detected: {}'.format(
            len(usb_results), usb_list_str))
        device_list_str = ' (Device list: {})'.format(usb_list_str)
    return device_list_str

def _open_camera(port):
    verbose_log('Opening camera...')
    try:
        camera = cv2.VideoCapture(port)
    except Exception as error:
        verbose_log(error)
    try:
        backend = camera.getBackendName()
    except:
        backend = 'not available'
    verbose_log('using backend: ' + backend)
    sleep(0.1)
    try:
        camera_open = camera.isOpened()
    except NameError:
        camera_open = False
    if not camera_open:
        verbose_log('Camera is not open.')
        log('Could not connect to camera.', 'error')
        return
    verbose_log('Camera opened successfully.')
    return camera

def _adjust_settings(camera, image_width, image_height):
    try:
        camera.set(cv2.CAP_PROP_FRAME_WIDTH, image_width)
        camera.set(cv2.CAP_PROP_FRAME_HEIGHT, image_height)
    except AttributeError:
        camera.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, image_width)
        camera.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, image_height)

def _check_camera_availability(camera_path):
    try:
        pids = subprocess.check_output(['fuser', camera_path])
    except MissingError:
        verbose_log('Unable to check if busy.')
    except subprocess.CalledProcessError:
        verbose_log('Camera not busy.')
    else:
        verbose_log('{} busy. Attempting to close...'.format(camera_path))
        for pid in pids.strip().split(b' '):
            subprocess.call(['kill', '-9', pid])

def _capture_usb_image(camera):
    try:
        return camera.read()
    except Exception as error:
        verbose_log(error)
        log('Image capture error.', 'error')
        return 0, None

def _log_no_image():
    verbose_log('No image.')
    log('Problem getting image.', 'error')

def usb_camera_photo():
    'Take a photo using a USB camera.'
    # Settings
    camera_port = 0      # default USB camera port
    max_port_num = 1     # highest port to try if not detected on port
    discard_frames = 10  # number of frames to discard for auto-adjust
    max_attempts = 5     # number of failed discard frames before quit
    image_width = int(WIDTH)
    image_height = int(HEIGHT)

    # Check USB devices for camera
    device_list_str = _get_usb_device_list()
    # Check video ports for camera
    video_ports = get_video_port_list()
    verbose_log('{} video ports detected: {}'.format(
        len(video_ports), ','.join(video_ports)))
    if len(video_ports) < 1:
        log('USB Camera not detected.{}'.format(device_list_str), 'error')
        return
    max_port_num = len(video_ports) - 1
    verbose_log('Adjusting max port number to {}.'.format(max_port_num))
    ret = False
    while camera_port <= max_port_num:
        camera_path = '/dev/video' + str(camera_port)
        verbose_log('Trying {}'.format(camera_path))
        if not os.path.exists(camera_path):
            verbose_log('{} missing'.format(camera_path))
            camera_port += 1
            continue

        # Close process using camera (if open)
        _check_camera_availability(camera_path)

        # Open the camera
        camera = _open_camera(camera_port)
        if camera is None:
            return

        verbose_log('Adjusting image with test captures...')
        # Set image size
        _adjust_settings(camera, image_width, image_height)
        # Capture test frame
        ret, _ = _capture_usb_image(camera)
        if not ret:
            camera.release()
            verbose_log('Couldn\'t get frame from {}'.format(camera_path))
            camera_port += 1
            continue
        break
    if not ret:
        _log_no_image()
        return
    verbose_log('First test frame captured.')
    # Let camera adjust
    failed_attempts = 0
    for _ in range(discard_frames):
        if not camera.grab():
            verbose_log('Could not get frame.')
            failed_attempts += 1
        if failed_attempts >= max_attempts:
            break
        sleep(0.1)

    # Take a photo
    verbose_log('Taking photo...')
    ret, image = _capture_usb_image(camera)

    # Close the camera
    camera.release()

    # Output
    if ret:  # an image has been returned by the camera
        verbose_log('Photo captured.')
        save_image(image)
    else:  # no image has been returned by the camera
        _log_no_image()

def rpi_camera_photo():
    'Take a photo using the Raspberry Pi Camera.'
    tempfile = upload_path('temporary')
    verbose_log('Taking photo with Raspberry Pi camera...')
    retcode = rpi_photo_call(tempfile)
    if retcode == 0:
        verbose_log('Image captured.')
        image = cv2.imread(tempfile)
        os.remove(tempfile)
        save_image(image)
    else:
        log('Raspberry Pi Camera not detected.', 'error')

def take_photo():
    'Take a photo.'
    CAMERA = get_camera_selection()

    if 'NONE' in CAMERA:
        log(CAMERA_DISABLED_MSG, 'error')
    elif 'RPI' in CAMERA:
        rpi_camera_photo()
    else:
        usb_camera_photo()

if __name__ == '__main__':
    take_photo()

ext: py

reason: import numpy


text:
"""https://adventofcode.com/2021/day/6"""

import numpy as np

data = open("day-06/input", encoding='utf-8').read()
df = np.array([int(x) for x in data.split(",")])

for day in range(80):
    df -= 1
    new = np.sum(df < 0)
    df = np.where(df < 0, 6, df)
    if new:
        df = np.hstack([df, np.full(new, 8)])
    print(f"After {day+1} days:", df)

print(df.shape)
# 362346


text:
#
# Gomoku图形界面：可与人工进行对弈。不是太友好，运行不流畅，AI考虑时要等待。
#
# 采用pygame写的五子棋对弈程序：通过判断每个点的价值。
#
#

from mytoolkit import print_time, load_config
import logging.config
# logging设置只能执行一次，要确保最先执行，在其它包含文件之前，否则包含文件会WARNING及以上才会记录。
logging.config.dictConfig(load_config('./conf/train_config.yaml')['train_logging'])
# 目的得到当前程序名，便于定位。
_logger = logging.getLogger(__name__)

from time import sleep
import pygame
from pygame.locals import *
from random import randint
import numpy as np
from game import Board
from mcts_alphaZero import MCTSPlayer
from policy_value_net_keras import PolicyValueNet               # Keras

level = 15
grade = 10
MAX = 1008611

def Scan(chesspad, color):
    shape = [[[0 for high in range(5)] for col in range(15)] for row in range(15)]
    # 扫描每一个点，然后在空白的点每一个方向上做出价值评估！！
    for i in range(15):
        for j in range(15):
            # 如果此处为空 那么就可以开始扫描周边
            if chesspad[i][j] == 0:
                m = i
                n = j
                # 如果上方跟当前传入的颜色参数一致，那么加分到0位！
                while n - 1 >= 0 and chesspad[m][n - 1] == color:
                    n -= 1
                    shape[i][j][0] += grade
                if n-1>=0 and chesspad[m][n - 1] == 0:
                    shape[i][j][0] += 1
                if n-1 >= 0 and chesspad[m][n - 1] == -color:
                    shape[i][j][0] -= 2
                m = i
                n = j

                # 如果下方跟当前传入的颜色参数一致，那么加分到0位！
                while (n + 1 < level  and chesspad[m][n + 1] == color):
                    n += 1
                    shape[i][j][0] += grade
                if n + 1 < level  and chesspad[m][n + 1] == 0:
                    shape[i][j][0] += 1
                if n + 1 < level  and chesspad[m][n + 1] == -color:
                    shape[i][j][0] -= 2
                m = i
                n = j

                # 如果左边跟当前传入的颜色参数一致，那么加分到1位！
                while (m - 1 >= 0 and chesspad[m - 1][n] == color):
                    m -= 1
                    shape[i][j][1] += grade
                if m - 1 >= 0 and chesspad[m - 1][n] == 0:
                    shape[i][j][1] += 1
                if m - 1 >= 0 and chesspad[m - 1][n] == -color:
                    shape[i][j][1] -= 2
                m = i
                n = j

                # 如果右边跟当前传入的颜色参数一致，那么加分到1位！
                while (m + 1 < level  and chesspad[m + 1][n] == color):
                    m += 1
                    shape[i][j][1] += grade
                if m + 1 < level  and chesspad[m + 1][n] == 0:
                    shape[i][j][1] += 1
                if m + 1 < level  and chesspad[m + 1][n] == -color:
                    shape[i][j][1] -= 2
                m = i
                n = j

                # 如果左下方跟当前传入的颜色参数一致，那么加分到2位！
                while (m - 1 >= 0 and n + 1 < level  and chesspad[m - 1][n + 1] == color):
                    m -= 1
                    n += 1
                    shape[i][j][2] += grade
                if m - 1 >= 0 and n + 1 < level  and chesspad[m - 1][n + 1] == 0:
                    shape[i][j][2] += 1
                if m - 1 >= 0 and n + 1 < level  and chesspad[m - 1][n + 1] == -color:
                    shape[i][j][2] -= 2
                m = i
                n = j

                # 如果右上方跟当前传入的颜色参数一致，那么加分到2位！
                while (m + 1 < level  and n - 1 >= 0 and chesspad[m + 1][n - 1] == color):
                    m += 1
                    n -= 1
                    shape[i][j][2] += grade
                if m + 1 < level  and n - 1 >= 0 and chesspad[m + 1][n - 1] == 0:
                    shape[i][j][2] += 1
                if m + 1 < level  and n - 1 >= 0 and chesspad[m + 1][n - 1] == -color:
                    shape[i][j][2] -= 2
                m = i
                n = j

                # 如果左上方跟当前传入的颜色参数一致，那么加分到3位！
                while (m - 1 >= 0 and n - 1 >= 0 and chesspad[m - 1][n - 1] == color):
                    m -= 1
                    n -= 1
                    shape[i][j][3] += grade
                if m - 1 >= 0 and n - 1 >= 0 and chesspad[m - 1][n - 1] == 0:
                    shape[i][j][3] += 1
                if m - 1 >= 0 and n - 1 >= 0 and chesspad[m - 1][n - 1] == -color:
                    shape[i][j][3] -= 2
                m = i
                n = j

                # 如果右下方跟当前传入的颜色参数一致，那么加分到3位！
                while m + 1 < level  and n + 1 < level  and chesspad[m + 1][n + 1] == color:
                    m += 1
                    n += 1
                    shape[i][j][3] += grade
                if m + 1 < level  and n + 1 < level  and chesspad[m + 1][n + 1] == 0:
                    shape[i][j][3] += 1
                if m + 1 < level  and n + 1 < level  and chesspad[m + 1][n + 1] == -color:
                    shape[i][j][3] -= 2
    return shape

def Sort(shape):
    for i in shape:
        for j in i:
            for x in range(5):
                for w in range(3, x - 1, -1):
                    if j[w - 1] < j[w]:
                        temp = j[w]
                        j[w - 1] = j[w]
                        j[w] = temp
    print("This Time Sort Done !")
    return shape

def Evaluate(shape):
    for i in range(level):
        for j in range(level):

            if shape[i][j][0] == 4:
                return i, j, MAX
            shape[i][j][4] = shape[i][j][0]*1000 + shape[i][j][1]*100 + shape[i][j][2]*10 + shape[i][j][3]
    max_x = 0
    max_y = 0
    max = 0
    for i in range(15):
        for j in range(15):
            if max < shape[i][j][4]:
                max = shape[i][j][4]
                max_x = i
                max_y = j
    print("the max is "+ str(max) + " at ( "+ str(max_x)+" , "+str(max_y)+" )")
    return max_x, max_y, max

class Chess(object):
    def __init__(self):
        self.a = [[0 for high in range(15)] for col in range(15)]
        # print(self.a)

    def fall(self, x, y, color):
        if (x < 0 or x > level - 1 or y < 0 or y > level - 1):
            return
        self.a[x][y] = color
        if Judge(x, y, color, self.a, 4):
            if color < 0:
                print("The Winner is White!!")
            else:
                print("The Winner is Black!!")

    def isEmpty(self, m, n):
        if self.a[m][n] != 0:
            return False
        else:
            return True

def Judge(x, y, color, CHESSLOCATION, length):
    count1, count2, count3, count4 = 0, 0, 0, 0
    # 横向判断
    i = x - 1
    while (i >= 0):
        if color == CHESSLOCATION[i][y]:
            count1 += 1
            i -= 1
        else:
            break
    i = x + 1
    while i < level:
        if CHESSLOCATION[i][y] == color:
            count1 += 1
            i += 1
        else:
            break

    # 纵向判断
    j = y - 1
    while (j >= 0):
        if CHESSLOCATION[x][j] == color:
            count2 += 1
            j -= 1
        else:
            break
    j = y + 1
    while j < level:
        if CHESSLOCATION[x][j] == color:
            count2 += 1
            j += 1
        else:
            break

    # 正对角线判断
    i, j = x - 1, y - 1
    while (i >= 0 and j >= 0):
        if CHESSLOCATION[i][j] == color:
            count3 += 1
            i -= 1
            j -= 1
        else:
            break
    i, j = x + 1, y + 1
    while (i < level and j < level):
        if CHESSLOCATION[i][j] == color:
            count3 += 1
            i += 1
            j += 1
        else:
            break
    # 反对角线判断
    i, j = x + 1, y - 1
    while (i < level and j >= 0):
        if CHESSLOCATION[i][j] == color:
            count4 += 1
            i += 1
            j -= 1
        else:
            break
    i, j = x - 1, y + 1
    while (i > 0 and j < level):
        if CHESSLOCATION[i][j] == color:
            count4 += 1
            i -= 1
            j += 1
        else:
            break

    if count1 >= length or count2 >= length or count3 >= length or count4 >= length:
        return True
    else:
        return False

def Autoplay(ch, m, n):
    a1 = [1,-1,1,-1,1,-1,0,0]
    b1 = [1,-1,-1,1,0,0,1,-1]
    rand = randint(0,7)
    while m+a1[rand]>=0 and m+a1[rand]<level and n+b1[rand]>=0 and n+b1[rand]<level and ch[m+a1[rand]][n+b1[rand]]!=0 :
        rand = randint(0,7)
    return m + a1[rand], n+b1[rand]

def BetaGo(ch, m, n, color, times):
    if times < 2:
        return Autoplay(ch, m, n)
    else:
        shape_P = Scan(ch, -color)
        shape_C = Scan(ch,color)
        shape_P = Sort(shape_P)
        shape_C = Sort(shape_C)
        max_x_P, max_y_P, max_P = Evaluate(shape_P)
        max_x_C, max_y_C, max_C = Evaluate(shape_C)

        if max_P > max_C and max_C<MAX:
            return max_x_P, max_y_P
        else:
            return max_x_C, max_y_C

def satrtGUI(board, mcts_player):
    pygame.init()
    bg = 'bg.png'
    white_image = 'white.png'
    black_image = 'black.png'

    screen = pygame.display.set_mode((750, 750), 0, 32)
    background = pygame.image.load(bg).convert()
    white = pygame.image.load(white_image).convert_alpha()
    black = pygame.image.load(black_image).convert_alpha()
    white = pygame.transform.smoothscale(white, (int(white.get_width() * 1.5), int(white.get_height() * 1.5)))
    black = pygame.transform.smoothscale(black, (int(black.get_width() * 1.5), int(black.get_height() * 1.5)))

    screen.blit(background, (0, 0))
    font = pygame.font.SysFont("黑体", 40)

    pygame.event.set_blocked([1, 4, KEYUP, JOYAXISMOTION, JOYBALLMOTION, JOYBUTTONDOWN, JOYBUTTONUP, JOYHATMOTION])
    pygame.event.set_allowed([MOUSEBUTTONDOWN, MOUSEBUTTONUP, 12, KEYDOWN])

    # 点坐标：
    dot_list = [(int(25 + i * 50 - white.get_width() / 2), int(25 + j * 50 - white.get_height() / 2)) for i in range(level) for
                j in range(level)]
    # print(dot_list)

    start_player = 0            # 0 or 1
    board.init_board(start_player)

    color = -1
    times = 0
    flag = False
    play_steps = []

    while not flag:
        for event in pygame.event.get():
            if event.type == QUIT:
                exit()
            elif event.type == MOUSEBUTTONDOWN:
                x, y = pygame.mouse.get_pos()
                if 25 <= x <= 725 and 25 <= y <= 725 and ((x - 25) % 50 <= level or (x - 25) % 50 >= 0) and (
                        (y - 25) % 50 <= level or (y - 25) % 50 >= 0):

                    color = -1 * color
                    m = int(round((x - 25) / 50))
                    n = int(round((y - 25) / 50))

                    move = board.location_to_move([14-n, m])
                    print("you press: [%d %d]-%d, original board(%d %d)" % (m, n, move, 14-n, m))

                    if move not in board.availables:
                    # if not ch.isEmpty(m, n):
                        print("Black OverWrite~~: [%d, %d]" % (m, n))
                        continue
                    # ch.fall(m, n, color)

                    board.do_move(move)
                    play_steps.append(move)
                    # 画黑棋：
                    # screen.blit(black, dot_list[level * m + n])
                    pygame.draw.circle(screen, (10, 0, 0), np.array(dot_list[level * m + n])+20, 20)

                    end, winner = board.game_end()
                    if end:
                    # if Judge(m, n, color, ch.a, 4):
                        _logger.warning("{'steps':[%s],'winner':%d,'start_player':1}" % (
                            ",".join(map(str, play_steps)), 1))
                        print("player%d win!" % winner)
                        screen.blit(font.render('GAME OVER,Black is win!', True, (110, 210, 30)), (80, 650))
                        break

                    color = -1 * color
                    sleep(0.1)

                    # 采用模型与AI对战：双方轮流走棋，对弈步骤都记录在board中。
                    move = mcts_player.get_action(board)
                    [x1, y1] = board.move_to_location(move)
                    x = y1
                    y = 14-x1
                    print("AI response: [%d %d]-%d, original board(%d %d)" % (x, y, move, x1, y1))
                    # x, y = BetaGo(ch.a, m, n, color, times)
                    times += 1
                    print("Predict:" + str(x) + " and " + str(y))
                    board.do_move(move)
                    play_steps.append(move)
                    # ch.fall(x, y, color)
                    # screen.blit(white, dot_list[level * x + y])
                    pygame.draw.circle(screen, (210, 200, 220), np.array(dot_list[level * x + y])+20, 20)

                    # if Judge(x, y, color, ch.a, 4):
                    end, winner = board.game_end()
                    if end:
                        win = 'GAME OVER, player%d is win!' % winner

                        _logger.warning("{'steps':[%s],'winner':%d,'start_player':1}" % (
                            ",".join(map(str, play_steps)), winner))

                        screen.blit(font.render(win, True, (217, 20, 30)), (80, 650))
                        break

        pygame.display.update()
        if flag:
            sleep(5)

if __name__ == '__main__':
    # now = Chess()
    model_file = 'current_policy.model'
    conf = load_config('./conf/train_config.yaml')

    _logger.info("AI model: %s" % model_file)
    best_policy = PolicyValueNet(conf["board_width"], conf["board_height"], model_file)

    # 采用训练的AI模型作为对手：n_playout越大，水平越高，速度明显快很多！
    mcts_player = MCTSPlayer(best_policy.policy_value_fn,
                         c_puct=5,
                         n_playout=600)

    board = Board(width=conf["board_width"], height=conf["board_height"], n_in_row=conf["n_in_row"])

    satrtGUI(board, mcts_player)


text:
#!/usr/bin/env python

from keras.layers import Input, Dense
from keras.models import Model
import keras.backend as K

import tensorflow as tf
sess = tf.Session()
K.set_session(sess)

#input = tf.placeholder(tf.float32, shape=(None, 2)) #Input(shape=[2])
input = Input(shape=[2])
probs = Dense(1, activation='linear', init='normal')(input)

model = Model(input, probs)
model.summary()

sess.run(tf.initialize_all_variables())

with sess.as_default():
    print sess.run([probs], feed_dict={input: [[1,1]], K.learning_phase(): 0})
    print sess.run([probs], feed_dict={input: [[2,1]], K.learning_phase(): 0})
    print sess.run([probs], feed_dict={input: [[1,0]], K.learning_phase(): 0})

    grads = tf.gradients(probs, input)
    print grads

batch_size = 1
weights = model.trainable_weights # weight tensors

print "====================="
print "WEIGHTS:"

#weights = [weight for weight in weights if model.get_layer(weight.name[:-2]).trainable] # filter down weights tensors to only ones which are trainable
#gradients = model.optimizer.get_gradients(model.total_loss, weights) # gradient tensors
if K._BACKEND == 'tensorflow':
  gradients = K.gradients(model.output, model.trainable_weights)
  gradients = [g / float(batch_size) for g in gradients]  # since TF sums over the batch
elif K._BACKEND == 'theano':
  import theano.tensor as T
  gradients = T.jacobian(model.output.flatten(), model.trainable_weights)
  gradients = [K.mean(g, axis=0) for g in gradients]
else:
  raise RuntimeError('Unknown Keras backend "{}".'.format(K._BACKEND))

for (wt, g) in zip(weights, gradients):
    print "--- wts/grads:"
    print sess.run(wt)

print "-------------------------"
print "Layers", model.layers
print "-------------------------"
print "Weights", weights
print "-------------------------"
print "Gradients", gradients
# ==> [dense_1_W, dense_1_b]

import keras.backend as K

'''
input_tensors = [model.inputs[0], # input data
                 model.sample_weights[0], # how much to weight each sample by
                 model.targets[0], # labels
                 K.learning_phase(), # train or test mode
]

get_gradients = K.function(inputs=input_tensors, outputs=gradients)

from keras.utils.np_utils import to_categorical

inputs0 = [[[1, 2]], # X
          [1], # sample weights
          [[1]], # y
          0 # learning phase in TEST mode
]
inputs1 = [[[1, 1]], # X
          [1], # sample weights
          [[1]], # y
          0 # learning phase in TEST mode
]
inputs2 = [[[1, 2]], # X
          [1], # sample weights
          [[0]], # y
          0 # learning phase in TEST mode
]

show_weights = [weight.value() for weight in weights]
#show_weights = [weight.initialized_value() for weight in weights]
print zip(show_weights, get_gradients(inputs0))
print zip(show_weights, get_gradients(inputs1))
print zip(show_weights, get_gradients(inputs2))
# ==> [(dense_1_W, array([[-0.42342907],
#                          [-0.84685814]], dtype=float32)),
#       (dense_1_b, array([-0.42342907], dtype=float32))]

print "================================"
print "================================"
print "================================"

'''

def collect_trainable_weights(layer):
    '''Collects all `trainable_weights` attributes,
    excluding any sublayers where `trainable` is set the `False`.
    '''
    trainable = getattr(layer, 'trainable', True)
    if not trainable:
        return []
    weights = []
    if layer.__class__.__name__ == 'Sequential':
        for sublayer in layer.flattened_layers:
            weights += collect_trainable_weights(sublayer)
    elif layer.__class__.__name__ == 'Model':
        for sublayer in layer.layers:
            weights += collect_trainable_weights(sublayer)
    elif layer.__class__.__name__ == 'Graph':
        for sublayer in layer._graph_nodes.values():
            weights += collect_trainable_weights(sublayer)
    else:
        weights += layer.trainable_weights
    # dedupe weights
    weights = list(set(weights))
    weights.sort(key=lambda x: x.name)
    return weights

# from https://github.com/fchollet/keras/issues/3062
import keras
def get_trainable_params(model):
    params = []
    for layer in model.layers:
        params += collect_trainable_weights(layer)
    return params

network_params = get_trainable_params(model)
#param_grad = tf.gradients(cost, network_params)

#import tensorflow as tf
# first define the model that takes state as input and outputs actions
#action = model(Input)

#x = tf.placeholder(tf.float32, shape=(None, 2))
#param_grad = tf.gradients(model.output, network_params)#, grad_ys=dQ_da)
#param_grad = tf.gradients(action, network_params, grad_ys=dQ_da)
# where dQ_da is grad of critic wrt action
#sess.run(param_grad, {x: [[1,0]]})

#print model.output
#print param_grad
#print grads([[1,2]])

reason: import theano


text:
import pandas as pd
import torch
import numpy as np
from util import Tool
import pdb
import re

class BaseRobot:
    def __init__(self, name):
        self.name = name

    def make_decision(self, data):
        pass

class TableRobot(BaseRobot):
    def __init__(self, name, base, n_after):
        super().__init__(name)
        self.base = base
        self.n = n_after
        # hour = time.gmtime()[3]
        # self.load(hour)
        # self._get_decision_table()

    def load(self, hour):
        base = self.base
        n = self.n
        fn = 'workspace/data/model/table_tr_'+str(base)+'_'+str(n)+'_'+str(hour)+'.csv'
        self.table = pd.read_csv(fn)

    def get_strict_decision_table(self):
        tmp = self.table.copy()
        base = self.base
        val_col = 'c' + str(base)
        tmp = tmp[tmp[val_col] != 0]

        tmp_pos = tmp.copy()
        tmp_pos.loc[tmp_pos[val_col]==1, ] = -1
        tmp_pos[val_col] = tmp_pos[val_col]/tmp_pos[val_col].abs()
        tmp_pos = tmp_pos[tmp_pos[val_col] > 0]

        tmp_neg = tmp.copy()
        tmp_neg.loc[tmp_neg[val_col]==-1, ] = 1
        tmp_neg[val_col] = tmp_neg[val_col]/tmp_neg[val_col].abs()
        tmp_neg = tmp_neg[tmp_neg[val_col] < 0]
        self.s_table = pd.concat([tmp_pos, tmp_neg]).reset_index(drop=True)
        del tmp
        return self.d_table

    def _get_decision_table(self):
        tmp = self.table.copy()
        base = self.base
        col = 'c' + str(base)
        tmp = tmp[tmp[col] != 0]

        tmp[col] = tmp[col]/tmp[col].abs()
        by = list(tmp.columns)[:base+1]
        tmp = tmp.groupby(by).sum()
        self.d_table = tmp.reset_index()
        del tmp
        return self.d_table

    def _trend_table(self, data, pre_decision):
        data = data[:3]
        if data[0] == 2:
            return 'wait'
        elif data[0] == 0:
            return 'buy'
        elif data[0] == 4:
            return 'sell'
        elif pre_decision == 'buy' and data[0] == 3:
            return 'close'
        elif pre_decision == 'sell' and data[0] == 1:
            return 'close'
        elif data[:2] == [1, 0]:
            return 'buy'
        elif data[:2] == [3, 4]:
            return 'sell'
        else:
            return 'wait'

    def make_decision(self, index, pre_decision):

        index = index[-self.base:]
        index = np.array(index).argsort()
        # index = pd.Series(index)
        # index = index.rank()-1
        # print(index)
        prob = self.search_table(index)
        # print(prob)
        data = prob.iloc[:, self.base+1:].values
        data = torch.Tensor(data)
        # print(data)
        max_, index = torch.max(data, dim=0)
        index = index.tolist()
        return self._trend_table(index, pre_decision)

    def brake(self, price, pre_decision, high, low):
        if pre_decision == 'sell' and (price-low)>10**-4:
            return 'close'
        elif pre_decision == 'buy' and (high-price)>10**-4:
            return 'close'
        else:
            return pre_decision

    def search_table(self, index):
        ret = self.table.copy()
        cnt = 0
        for i in index.tolist():
            col = 'c'+str(cnt)
            ret = ret[ret[col]==i]
            cnt += 1

        return ret

class ProbRobot(BaseRobot):

    def __init__(self, name, prob_table=None):
        super().__init__(name)
        if prob_table is None:
            fn = 'workspace/data/model/prob_in_trade_adjustment.csv'
            self.prob_table = pd.read_csv(fn)
        else:
            self.prob_table = prob_table
        self.di = None

    def set_opening_price(self, price):
        # a daily opening price, for forex we set it as end  price
        # in the last day
        self.open_price = price

    def calc_avg_price(self, history, period=100):
        if len(history) < period:
            self.avg = None
        else:
            data = history[-period:]
            self.avg = np.mean(data)

    def set_direction(self, di=None):
        self.di = di

    def make_decision(self, price):
        # search the best choice for buy and for sell
        table = self.prob_table
        last_data_num = table['data_num'].unique()[-1]
        data = table[table['data_num'] == last_data_num]
        data = data.fillna('nan')
        data = data.drop('data_num', axis=1)
        data = data.groupby(['prof', 'loss', 'open_di',\
                'fl_interval', 'mean_di']).mean()
        data = data.reset_index()
        buy_keywords, sell_keywords = self.gen_keyword(price)
        if buy_keywords is None:
            buy_best = None
        else:
            buy_table = self._search_table_by_keywords(buy_keywords)
            # pdb.set_trace()
            buy_best = buy_table.sort_values(['long_ex_prof'], ascending=False).head(1)

        if sell_keywords is None:
            sell_best = None
        else:
            sell_table = self._search_table_by_keywords(sell_keywords)
            sell_best = sell_table.sort_values(['short_ex_prof'], ascending=False).head(1)

        # compare buy and sell
        if buy_best is None and sell_best is None:
            return None
        elif buy_best is None:
            selected = 'sell'
        elif sell_best is None:
            selected = 'buy'
        else:
            try:
                buy_expected = buy_table['long_ex_prof'].tolist()[0]
                sell_expected = sell_table['short_ex_prof'].tolist()[0]
            except Exception:
                return None
            if buy_expected > sell_expected and buy_expected > 0:
                selected = 'buy'
            elif buy_expected <= sell_expected and sell_expected > 0:
                selected = 'sell'
            else:
                return None

        if selected == 'buy':
            prof = buy_best['prof'].tolist()[0]
            loss = buy_best['loss'].tolist()[0]
            return ['buy', prof, loss]
        if selected == 'sell':
            prof = sell_best['prof'].tolist()[0]
            loss = sell_best['loss'].tolist()[0]
            return ['sell', prof, loss]

    def _search_table_by_keywords(self, keywords):
        ret = pd.DataFrame()
        for keyword in keywords:
            table = self.prob_table.copy()
            for key in keyword:
                table = table[table[key] == keyword[key]]

            ret = pd.concat([ret, table])
        return ret

    def gen_keyword(self, curr_price):
        # generate some set of keywords and check if
        # the relative operation has a positive expected
        # profit
        open_diff = curr_price - self.open_price
        if self.avg is None:
            return None, None
        mean_diff = curr_price - self.avg

        # for buy
        if open_diff > 0:
            open_di = True
        elif open_diff < 0:
            open_di = False
        else:
            open_di = 'nan'
        if mean_diff > 0:
            mean_di = True
        elif mean_diff < 0:
            mean_di = False
        else:
            mean_di = 'nan'
        o_params = ['nan', open_di]
        m_params = ['nan', mean_di]
        o_params = list(dict.fromkeys(o_params))
        m_params = list(dict.fromkeys(m_params))
        params = [o_params, m_params]
        buy_keywords = Tool.get_all_combined(params)
        names = ['open_di', 'mean_di']
        buy_keywords = Tool.list2dict(buy_keywords, names)

        # for sell
        if open_diff > 0:
            open_di = False
        elif open_diff < 0:
            open_di = True
        else:
            open_di = 'nan'
        if mean_diff > 0:
            mean_di = False
        elif mean_diff < 0:
            mean_di = True
        else:
            mean_di = 'nan'

        o_params = ['nan', open_di]
        m_params = ['nan', mean_di]
        o_params = list(dict.fromkeys(o_params))
        m_params = list(dict.fromkeys(m_params))
        params = [o_params, m_params]
        sell_keywords = Tool.get_all_combined(params)
        names = ['open_di', 'mean_di']
        sell_keywords = Tool.list2dict(sell_keywords, names)

        if self.di == 'buy':
            sell_keywords = None
        elif self.di == 'sell':
            buy_keywords = None

        return buy_keywords, sell_keywords


text:
import re
import pandas as pd
import numpy as np
from sklearn import metrics

train = pd.read_csv('/Users/ccc/src/dev/titanic/data/train.csv')
test = pd.read_csv('/Users/ccc/src/dev/titanic/data/test.csv')

def breakup_names(df):
    """Breakup full name into surname, title, first name"""
    df[['surname','given_name']] = df['Name'].str.split(",", expand= True)
    df[['title', 'first_name']] = df['given_name'].str.split(n=1, expand= True)
    df = df.drop(columns=['Name', 'given_name'])
    return df

def predict_age(df1, df2):
    """Predict age via regression since we have so many NaN"""

train = breakup_names(train)
test = breakup_names(test)

max_stars_repo_name: camille-chanel/titanic

max_issues_repo_name: camille-chanel/titanic

max_forks_repo_name: camille-chanel/titanic


text:
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
fMRI Design Matrix creation functions.
"""

import numpy as np
from nipy.modalities.fmri import formula, utils, hrf

##########################################################
# Paradigm handling
##########################################################

class Paradigm(object):
    """
    Simple class to hanle the experimental paradigm in one session
    """

    def __init__(self, index=None, onset=None, amplitude=None):
        """
        Parameters
        ----------
        index: array of shape (n_events), type = int, optional
               index of the events (the index of their exprimental condition)
        onset: array of shape (n_events), type = float, optional,
               onset time (in s.) of the events
        amplitude: array of shape (n_events), type = float, optional,
                   amplitude of the events (if applicable)
        """
        self.index = index
        self.onset = onset
        self.amplitude = amplitude
        if index is not None:
            self.n_events = len(onset)
            self.index = np.ravel(np.array(index)).astype(np.int)
        if onset is not None:
            if len(onset) != self.n_events:
                raise ValueError, 'inconsistant definition of index and onsets'
            self.onset = np.ravel(np.array(onset))
        if amplitude is not None:
            if len(amplitude) != self.n_events:
                raise ValueError, 'inconsistant definition of amplitude'
            self.amplitude = np.ravel(np.array(amplitude))
        self.type= 'event'

class EventRelatedParadigm(Paradigm):
    """ Class to handle event-related paradigms
    """

    def __init__(self, index=None, onset=None, amplitude=None):
        """
        Parameters
        ----------
        index: array of shape (n_events), type = int, optional
               index of the events (the index of their exprimental condition)
        onset: array of shape (n_events), type = float, optional
               onset time (in s.) of the events
        amplitude: array of shape (n_events), type = float, optional,
                   amplitude of the events (if applicable)
        """
        Paradigm.__init__(self, index, onset, amplitude)

class BlockParadigm(Paradigm):
    """ Class to handle block paradigms
    """

    def __init__(self, index=None, onset=None, duration=None, amplitude=None):
        """
        Parameters
        ----------
        index: array of shape (n_events), type = int, optional
               index of the events (the index of their exprimental condition)
        onset: array of shape (n_events), type = float, optional
               onset time (in s.) of the events
        amplitude: array of shape (n_events), type = float, optional,
                   amplitude of the events (if applicable)

        """
        Paradigm.__init__(self, index, onset, amplitude)
        self.duration = duration
        self.type = 'block'
        if duration is not None:
            if len(duration) != self.n_events:
                raise ValueError, 'inconsistant definition of duration'
            self.duration = np.ravel(np.array(duration))

def load_protocol_from_csv_file(path, session=None):
    """
    Read a (.csv) paradigm file consisting of values yielding
    (occurence time, (duration), event ID, modulation)
    and returns a paradigm instance or a dictionary of paradigm instances

    Parameters
    ----------
    path: string,
          path to a .csv file that describes the paradigm
    session: int, optional
             session number.
             by default the output is a dictionary
             of session-level dictionaries indexed by session

    Returns
    -------
    paradigm, paradigm instance (if session is provided), or
              dictionary of paradigm instances otherwise,
              the resulting session-by-session paradigm

    Note
    ----
    It is assumed that the csv file contains the following columns:
    (session id, condition id, onset),
    plus possibly (duration) and (amplitude)
    If all the durations are 0, the paradigm will be treated as event-related
    If only some of the durations are zero, there will probably we trouble

    fixme
    -----
    would be much clearer if amplitude was put before duration in the .csv
    """
    import csv
    csvfile = open(path)
    dialect = csv.Sniffer().sniff(csvfile.read())
    csvfile.seek(0)
    reader = csv.reader(open(path, "rb"),dialect)

    # load the csv as a protocol array
    protocol = []
    for row in reader:
        protocol.append([float(row[j]) for j in range(len(row))])
    protocol = np.array(protocol)

    def read_session(protocol, session):
        """ return a paradigm instance corresponding to session
        """
        ps = (protocol[:,0] == session)
        if np.sum(ps)==0:
            return None
        if protocol.shape[1] > 4:
            paradigm = BlockParadigm(protocol[ps, 1], protocol[ps, 2],
                                     protocol[ps, 3],  protocol[ps, 4])
        elif protocol.shape[1] == 4:
            amplitude = np.ones(np.sum(ps))
            paradigm = BlockParadigm(protocol[ps, 1], protocol[ps, 2],
                                     protocol[ps, 3], amplitude)
        else:
            amplitude = np.ones(np.sum(ps))
            paradigm = EventRelatedParadigm(protocol[ps, 1], protocol[ps, 2],
                                            amplitude)
        if (protocol.shape[1] > 3) and (protocol[ps,3]==0).all():
            paradigm = EventRelatedParadigm(protocol[ps, 1], protocol[ps, 2],
                                            protocol[ps, 4])

        return paradigm

    sessions = np.unique(protocol[:,0])
    if session is None:
        paradigm = {}
        for s in sessions:
            paradigm[s] = read_session(protocol, s)
    else:
        paradigm = read_session(protocol, session)

    return paradigm

##########################################################
# Design Matrices
##########################################################

class DesignMatrix(object):
    """
    Class to handle design matrices
    """

    def __init__(self, frametimes=None, paradigm=None, hrf_model='Canonical',
               drift_model='Cosine', hfcut=128, drift_order=1, fir_delays=[0],
               fir_duration=1., cond_ids=None, add_regs=None, add_reg_names=None):
        """
        Parameters
        ----------
        frametimes: array of shape(nbframes), optional
                    the timing of the scans
        paradigm: Paradigm instance, optional
                  descritpion of the experimental paradigm
        hrf_model, string, optional,
                   that specifies the hemodynamic reponse function
                   it can be 'Canonical', 'Canonical With Derivative' or 'FIR'
        drift_model, string
                     specifies the desired drift model,
                     to be chosen among 'Polynomial', 'Cosine', 'Blank'
        hfcut=128:  float ,
                    cut frequency of the low-pass filter
        drift_order=1, int,
                       order of the dirft model (in case it is polynomial)
        fir_delays=[0], optional,
                        array of shape(nb_onsets) or list
                        in case of FIR design, yields the array of delays
                        used in the FIR model
        fir_duration=1., float, duration of the FIR block;
                         in general it should be equal to the tr
        cond_ids=None, list of strin of length (ncond),
                       ids of the experimental conditions.
                       If None this will be called 'c0',..,'cn'
        add_regs=None, array of shape(naddreg, nbframes)
                       additional user-supplied regressors
        add_reg_names=None: list of (naddreg) regressor names
                            if None, while naddreg>0, these will be termed
                            'reg_%i',i=0..naddreg-1
        """
        self.frametimes = frametimes
        self.paradigm = paradigm
        self.hrf_model = hrf_model
        self.drift_model = drift_model
        self.hfcut = hfcut
        self.drift_order = drift_order
        self.fir_delays = fir_delays
        self.fir_duration = fir_duration
        self.cond_ids = cond_ids
        self.estimated = False

        # todo : check arguments
        if add_regs==None:
            self.n_add_regs = 0
        else:
            # check that regressor specification is correct
            if add_regs.shape[0] == np.size(add_regs):
                add_regs = np.reshape(add_regs, (np.size(1, add_regs)))
            if add_regs.shape[0] != np.size(frametimes):
                raise ValueError, 'incorrect specification of additional regressors'
            self.n_add_regs = add_regs.shape[1]
        self.add_regs = add_regs

        if  add_reg_names == None:
            self.add_reg_names = ['reg%d'%k for k in range(self.n_add_regs)]
        elif len(add_reg_names)!= self.n_add_regs:
             raise ValueError, 'Incorrect number of additional regressors names \
                               was provided'
        else:
            self.add_reg_names = add_reg_names

    def estimate(self):
        """
        Numerical estimation of self
        """
        self.drift = set_drift(self.drift_model, self.frametimes,
                               self.drift_order, self.hfcut)
        if self.paradigm==None:
            self.n_conditions = 0
            self.n_main_regressors = 0
            self.formula = self.drift
            self.names = []
        else:
            self.n_conditions = int(self.paradigm.index.max()+1)
            if self.cond_ids==None:
                self.cond_ids = ['c%d'%k for k in range(self.n_conditions)]
            self.conditions, self.names = convolve_regressors(
                self.paradigm, self.hrf_model, self.cond_ids, self.fir_delays,
                self.fir_duration)
            self.n_main_regressors = len(self.names)
            self.formula = self.conditions + self.drift

        # sample the matrix
        self.matrix = build_dmtx(self.formula, self.frametimes).T
        # FIXME: ugly  workaround the fact that NaN can occur when trials
        # are earlier than scans (!)
        self.matrix[np.isnan(self.matrix)] = 0

        # add user-supplied regressors
        if self.n_add_regs>0:
            self.matrix = np.hstack((self.matrix[:,:self.n_main_regressors],
                                     self.add_regs,
                                     self.matrix[:,self.n_main_regressors:]))

            # add the corresponding names
            self.names = self.names + self.add_reg_names

        # Force the design matrix to be full rank at working precision
        self.matrix, self.design_cond = full_rank(self.matrix)

        # complete the names with the drift terms
        for k in range(len(self.drift.terms)-1):
            self.names.append('drift_%d'%(k+1))
        self.names.append('constant')
        self.estimated = True

    def write_csv(self, path):
        """
        write oneselfs as a csv

        Parameters
        ----------
        path: string, path of the resulting csv file
        """
        import csv
        if self.estimated==False:
           self.estimate()

        fid = open(path, "wb")
        writer = csv.writer(fid)
        writer.writerow(self.names)
        writer.writerows(self.matrix)
        fid.close()

    def read_from_csv(self, path):
        """
        load self.matrix and self.names  from a csv file
        Parameter
        ---------
        path: string,
            path of the .csv file that includes the matriox and related information

        fixme
        -----
        needs to check that this is coherent with the information of self ?
        """
        import csv
        csvfile = open(path)
        dialect = csv.Sniffer().sniff(csvfile.read())
        csvfile.seek(0)
        reader = csv.reader(open(path, "rb"),dialect)

        boolfirst = True
        design = []
        for row in reader:
            if boolfirst:
                names = [row[j] for j in range(len(row))]
                boolfirst = False
            else:
                design.append([row[j] for j in range(len(row))])

            x = np.array([[float(t) for t in xr] for xr in design])
        self.matrix = x
        self.names = names

    def show(self, rescale=True, ax=None):
        """
        Vizualization of a design matrix

        Parameters
        ----------
        rescale: bool, optional
                 rescale columns for visualization or not
        ax: figure handle, optional

        Returns
        -------
        ax, figure handle

        """
        if self.estimated==False:
            self.estimate()

        x = self.matrix.copy()

        if rescale:
            x = x/np.sqrt(np.sum(x**2,0))

        import matplotlib.pylab as mp
        if ax is None:
            mp.figure()
            ax = mp.subplot(1, 1, 1)

        ax.imshow(x, interpolation='Nearest', aspect='auto')
        ax.set_label('conditions')
        ax.set_ylabel('scan number')

        if self.names is not None:
            ax.set_xticks(range(len(self.names)))
            ax.set_xticklabels( self.names, rotation=60, ha='right')

        #mp.subplots_adjust(top=0.99, bottom=0.25)
        return ax

def dmtx_light(frametimes, paradigm=None, hrf_model='Canonical',
               drift_model='Cosine', hfcut=128, drift_order=1, fir_delays=[0],
               fir_duration=1., cond_ids=None, add_regs=None,
               add_reg_names=None, path=None):
    """
    Make a design matrix while avoiding framework

    Parameters
    ----------
    frametimes: array of shape(nbframes),
        the timing of the scans
    paradigm: Paradigm instance, optional
              descritpion of the experimental paradigm
              if paradigm==None, then no condition is included
    hrf_model, string, optional,
             that specifies the hemodynamic reponse function
             it can be 'Canonical', 'Canonical With Derivative' or 'FIR'
    drift_model, string that specifies the desired drift model,
                to be chosen among 'Polynomial', 'Cosine', 'Blank'
    hfcut=128  float , cut frequency of the low-pass filter
    drift_order=1, int, order of the dirft model (in case it is polynomial)
    fir_delays=[0], optional, array of shape(nb_onsets) or list
                    in case of FIR design, yields the array of delays
                    used in the FIR model
    fir_duration=1., float, duration of the FIR block;
                     in general it should be equal to the tr
    cond_ids=None, list of strin of length (ncond),
                     ids of the experimental conditions.
                If None this will be called 'c0',..,'cn'
    add_regs=None, array of shape(naddreg, nbframes)
                   additional user-supplied regressors
    add_reg_names=None, list of (naddreg) regressor names
                        if None, while naddreg>0, these will be termed
                        'reg_%i',i=0..naddreg-1
    path: string, optional
         if not None, the matrix is written as a .csv file at the given path

    Returns
    -------
    dmtx array of shape(nreg, nbframes):
        the sampled design matrix
    names list of strings of len (nreg)
        the names of the columns of the design matrix
    """
    DM = DesignMatrix(frametimes, paradigm, hrf_model, drift_model, hfcut,
                      drift_order, fir_delays, fir_duration, cond_ids,
                      add_regs, add_reg_names)
    DM.estimate()
    if path is not None:
        DM.write_csv(path)
    return DM.matrix, DM.names

def _polydrift(order, tmax):
    """
    Create a polynomial drift formula

    Parameters
    ----------
    order, int, number of polynomials in the drift model
    tmax, float maximal time value used in the sequence
          this is used to normalize porperly the columns

    Returns
    -------
    pol a formula that contains all the polynomial drift
    plus a constant regressor
    """
    t = formula.Term('t')
    pt = []
    # fixme : ideally  this should be orthonormalized
    for k in range(order):
        pt.append(formula.define('poly_drift_%d'%(k+1),t**(k+1)/tmax**(k+1)))
    pt.append(formula.define('constant',1.0+0*t))
    pol =  formula.Formula(pt)
    return pol

def _cosinedrift(hfcut, tmax, tsteps):
    """
    Create a cosine drift formula

    Parameters
    ----------
    hfcut, float , cut frequency of the low-pass filter
    tmax, float  maximal time value used in the sequence
    tsteps, int,  number of TRs in the sequence

    Returns
    -------
    cos  a formula that contains all the polynomial drift
    plus a constant regressor
    """
    t = formula.Term('t')
    pt = []
    order = int(np.floor(2 * float(tmax) / float(hfcut)) + 1)
    for k in range(1,order):
        u = np.sqrt(2.0/tmax) * utils.sympy_cos(np.pi*(t/tmax+ 0.5/tsteps)*k )
        pt.append(formula.define('cosine_drift_%d'%(k+1),u))
    pt.append(formula.define('constant',1.0+0*t))
    cos =  formula.Formula(pt)
    return cos

def _blankdrift():
    """
    Create the blank drift formula

    Returns
    -------
    df  a formula that contains a constant regressor
    """
    t = formula.Term('t')
    pt = [formula.define('constant',1.0+0*t)]
    df =  formula.Formula(pt)
    return df

def set_drift(DriftModel, frametimes, order=1, hfcut=128.):
    """
    Create the drift formula

    Parameters
    ----------
    DriftModel: string,
                to be chosen among 'Polynomial', 'Cosine', 'Blank'
                that specifies the desired drift model
    frametimes: array of shape(ntimes),
                list of values representing the desired TRs
    order: int, optional,
        order of the dirft model (in case it is polynomial)
    hfcut: float, optional,
        frequency cut in case of a cosine model

    Returns
    -------
    df, the resulting drift formula
    """
    if DriftModel=='Polynomial':
        d = _polydrift(order, frametimes.max())
    elif DriftModel=='Cosine':
        d = _cosinedrift(hfcut, frametimes.max(), frametimes.size)
    elif DriftModel=='Blank':
        d  = _blankdrift()
    else:
        raise NotImplementedError,"unknown drift model"
    return d

def convolve_regressors(paradigm, hrf_model, names=None, fir_delays=[0],
    fir_duration = 1.):
    """
    Creation of  a formula that represents
    the convolution of the conditions onset witha  certain hrf model

    Parameters
    ----------
    paradigm: paradigm instance
    hrf_model, string that can be 'Canonical',
               'Canonical With Derivative' or 'FIR'
               that specifies the hemodynamic reponse function
    names=None, list of strings corresponding to the condition names
                if names==None, these are create as 'c1',..,'cn'
                meaning 'condition 1'.. 'condition n'
    fir_delays=[0], optional, array of shape(nb_onsets) or list
                    in case of FIR design, yields the array of delays
                    used in the FIR model
    fir_duration=1., float, duration of the FIR block;
                     in general it should eb equal to the tr

    Returns
    -------
    f a formula object that contains the convolved regressors
      as functions of time
    names list of strings corresponding to the condition names
          the output names depend on teh hrf model used
          if 'Canonical' then this is identical to the input names
          if 'Canonical With Derivative', then two names are produced for
             input name 'name': 'name' and 'name_derivative'

    fixme:
    normalization of the columns of the design matrix ?
    """
    ncond = int(paradigm.index.max()+1)
    if names==None:
        names=["c%d" % k for k in range(ncond)]
    else:
        if len(names)<ncond:
            raise ValueError, 'the number of names is less than the \
                  number of conditions'
        else:
            ncond = len(names)
    listc = []
    hnames = []
    typep = paradigm.type

    for nc in range(ncond):
        onsets =  paradigm.onset[paradigm.index==nc]
        nos = np.size(onsets)
        if paradigm.amplitude is not None:
            values = paradigm.amplitude[paradigm.index==nc]
        else:
            values = np.ones(nos)
        if nos>0:
            if typep=='event':
                if hrf_model=="Canonical":
                    c = formula.define(names[nc],
                                       utils.events(onsets, values, f=hrf.glover))
                    listc.append(c)
                    hnames.append(names[nc])
                elif hrf_model=="Canonical With Derivative":
                    c1 = formula.define(names[nc],
                                        utils.events(onsets, values, f=hrf.glover))
                    c2 = formula.define(names[nc]+"_derivative",
                                        utils.events(onsets, values, f=hrf.dglover))
                    listc.append(c1)
                    listc.append(c2)
                    hnames.append(names[nc])
                    hnames.append(names[nc]+"_derivative")
                elif hrf_model=="FIR":
                    for i,ft in enumerate(fir_delays):
                        lnames = names[nc] + "_delay_%d"%i
                        changes = np.hstack((onsets+ft, onsets+ft+fir_duration))
                        ochanges = np.argsort(changes)
                        lvalues = np.hstack((values, np.zeros(nos)))
                        changes = changes[ochanges]
                        lvalues = lvalues[ochanges]

                        c = formula.define(lnames, utils.step_function(changes, lvalues))

                        listc.append(c)
                        hnames.append(lnames)
                else:
                    raise NotImplementedError,'unknown hrf model'
            elif typep=='block':
                offsets =  onsets + paradigm.duration[paradigm.type==nc]
                changes = np.hstack((onsets, offsets))
                values = np.hstack((values, -values))

                if hrf_model=="Canonical":
                    c = utils.events(changes, values, f=hrf.iglover)
                    listc.append(c)
                    hnames.append(names[nc])
                elif hrf_model=="Canonical With Derivative":
                    c1 = utils.events(changes,values, f=hrf.iglover)
                    c2 = utils.events(changes,values, f=hrf.glover)
                    listc.append(c1)
                    listc.append(c2)
                    hnames.append(names[nc])
                    hnames.append(names[nc]+"_derivative")
                elif hrf_model=="FIR":
                    raise NotImplementedError,\
                          'block design are not compatible with FIR at the moment'
                else:
                    raise NotImplementedError,'unknown hrf model'

    # create the formula
    p = formula.Formula(listc)

    return p, hnames

def build_dmtx(form, frametimes):
    """
    This is a work arount to control the order of the regressor
    in the design matrix construction

    Parameters
    ----------
    form, the formula that describes the design matrix
    frametimes, array of shape (nb_time_samples), the time sampling grid

    Returns
    -------
    X array of shape (nrows,nb_time_samples) the sdesign matrix
    """
    # fixme : workaround to control matrix columns order
    t = formula.make_recarray(frametimes, 't')
    X = []
    for ft in form.terms:
        lf = formula.Formula([ft])
        X.append(lf.design(t, return_float=True))
    X = np.array(X)
    return X

def full_rank(X, cmax=1e15):
    """
    This function possibly adds a scalar matrix to X
    to guarantee that the condition number is smaller than a given threshold.

    Parameters
    ----------
    X array of shape(nrows,ncols)
    cmax=1.e-15, float tolerance for condition number

    Returns
    -------
    X array of shape(nrows,ncols) after regularization
    cmax=1.e-15, float tolerance for condition number
    """
    U, s, V = np.linalg.svd(X,0)
    sM = s.max()
    sm = s.min()
    c = sM/sm
    if c < cmax:
        return X, c
    print 'Warning: matrix is singular at working precision, regularizing...'
    lda = (sM-cmax*sm)/(cmax-1)
    s = s + lda
    X = np.dot(U, np.dot(np.diag(s), V))
    return X, cmax

max_stars_repo_name: fperez/nipy

max_issues_repo_name: fperez/nipy

max_forks_repo_name: fperez/nipy


text:
#!/usr/bin/python
# -*- coding: utf-8 -*-
'''
 * Copyright (c) 2017, Matthias Julius Kannwischer
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
'''

import csv
import math
import numpy as np
from os import listdir
from os.path import isfile, join
import matplotlib.pyplot as plt

def parse_csv(file):
    ifile = open(file, "rb")
    reader = csv.reader(ifile, delimiter=';')
    values = []
    for row in reader:
        values.append(int(row[6]))
    ifile.close()
    return values

def parse_all(path):
    values = [];
    for f in listdir(path):
        file = join(path, f)
        if isfile(file) and file.endswith(".log"):
            values.extend(parse_csv(file))
    return values

def hist_txt(x):
    return "mean= %0.3f\nmax= %d\nmin= %d\nmedian= %d" % (np.mean(x), max(x), min(x), np.median(x))

def plot_hist(x, fname):
    f = plt.figure(figsize=(11,5))
    plt.hist(x, range(min(x), max(x)+1), normed=True)
    plt.ylabel("relative frequency")
    plt.xlabel("number of required faulty signatures")
    #plt.figtext(.7, .7, hist_txt(x))

    table_data = [["mean", "%.1f" %np.mean(x)],["max", max(x)],["min", min(x)],["median","%d"%np.median(x)]]
    table = plt.table(cellText=table_data, colWidths = [0.05,0.05],
          loc='center right')

    table.scale(1.8,1.8)
    plt.tight_layout()
    f.savefig(fname)

def scatter_txt(x,y):
    txt = ""
    for p, d in zip(x, y):
        txt += "\np=%02d mean=%.3f" % (p, d)
    return txt

def plot_scatter(x, y, fname):
    f = plt.figure(figsize=(11,5))
    plt.plot(x,y, "x", ms=18.0)
    plt.ylabel("number of required faulty signatures")
    plt.xlabel("number of forgery trials p")
    plt.ylim(ymin=0)
    plt.ylim(ymax=math.ceil(max(y)+2))
    plt.tight_layout()
    plt.figtext(.1, .2, scatter_txt(x, y))
    #plt.show()
    f.savefig(fname)

def plot_boxplot(x, y, fname):
    f = plt.figure(figsize=(11,5))
    plt.ylabel("number of required faulty signatures")
    plt.xlabel("number of forgery trials p")
    plt.boxplot(y, labels=x, positions=x,showfliers=False)
    #plt.figtext(.8, .65,scatter_txt(x,[np.mean(a) for a in y]));
    table_data = zip(x,map(lambda y: "%.1f"%np.mean(y), y))
    table_cols = ["p", "mean"]
    table = plt.table(cellText=table_data,
                  colLabels=table_cols, loc='center right',bbox=[0.78, 0.55, 0.20, 0.4])

    plt.tight_layout()
    f.savefig(fname)

plt.style.use('classic')

data_P_1_N_32 = parse_all("../data/P_1/N_32")
plot_hist(data_P_1_N_32,"FA_P_1_N_32.pdf")

data_P_1_N_64 = parse_all("../data/P_1/N_64")
plot_hist(data_P_1_N_64,"FA_P_1_N_64.pdf")

data_P_1 = data_P_1_N_32[:1000]
data_P_2 = parse_all("../data/P_NOT1/P_2")[:1000]
data_P_4 = parse_all("../data/P_NOT1/P_4")[:1000]
data_P_8 = parse_all("../data/P_NOT1/P_8")[:1000]
data_P_16 = parse_all("../data/P_NOT1/P_16")[:1000]
data_P_24 = parse_all("../data/P_NOT1/P_24")[:1000]
data_P_32 = parse_all("../data/P_NOT1/P_32")[:1000]

p = [1,2,4,8,16,24,32]
d = [np.mean(data_P_1), np.mean(data_P_2), np.mean(data_P_4), np.mean(data_P_8), np.mean(data_P_16), np.mean(data_P_24), np.mean(data_P_32)]
plot_scatter(p,d, "FA_P_1_2_4_8_16_32.pdf")

d = [data_P_1, data_P_2, data_P_4, data_P_8, data_P_16, data_P_24, data_P_32]
plot_boxplot(p, d, "FA_P_boxplot.pdf")

max_stars_repo_name: mkannwischer/xmss-fault-attack

max_issues_repo_name: mkannwischer/xmss-fault-attack

max_forks_repo_name: mkannwischer/xmss-fault-attack


text:
import numpy as np
import tensorflow as tf

from fast_bfmatcher.matchers import Matcher, MatchResult

if tf.__version__ < "2.0.0":
    if tf.__version__ >= "1.14.0":
        if not tf.compat.v1.executing_eagerly():
            msg = "TFL2BFMatcher requires eager mode to be enabled! Call tf.enable_eager_execution() first."
            print(msg)
            raise ImportError(msg)
    else:
        msg = f"TF Matcher Requires minimum version of Tensorflow is 2.0.0, got: {tf.__version__}"
        print(msg)
        raise ImportError(msg)

@tf.function(input_signature=(tf.TensorSpec(shape=[None, None], dtype=tf.float32),))
def find_row_col_min_values_tf(X: tf.Tensor):
    row_indices = tf.argmin(X, 1)
    row_values = tf.reduce_min(X, 1)
    col_values = tf.reduce_min(X, 0)
    return row_indices, row_values, col_values

@tf.function(
    input_signature=(
        tf.TensorSpec(shape=[None, None], dtype=tf.float32),
        tf.TensorSpec(shape=[None, None], dtype=tf.float32),
    )
)
def l2_distance_matrix_tf(X: tf.Tensor, Y: tf.Tensor) -> tf.Tensor:
    sqnorm1 = tf.reduce_sum(tf.square(X), 1, keepdims=True)
    sqnorm2 = tf.reduce_sum(tf.square(Y), 1, keepdims=True)
    innerprod = tf.matmul(X, Y, transpose_a=False, transpose_b=True)
    return sqnorm1 + tf.transpose(sqnorm2) - 2 * innerprod

class TFL2CCBFMatcher(tf.Module, Matcher):
    def __init__(self, name: str = None):
        super(TFL2CCBFMatcher, self).__init__(name=name)
        self.dtype = tf.float32

    @classmethod
    def cast_input(cls, X: np.ndarray, dtype: np.dtype) -> np.ndarray:
        if X.dtype != dtype:
            X = tf.cast(X, dtype)
        return X

    @classmethod
    def l2_distance_matrix(cls, X: np.ndarray, Y: np.ndarray):
        D = l2_distance_matrix_tf(X, Y)
        return D.numpy()

    @classmethod
    def find_row_col_min_values(cls, X: np.ndarray):
        row_indices, row_values, col_values = find_row_col_min_values_tf(X)
        return row_indices.numpy(), row_values.numpy(), col_values.numpy()

    @tf.function(
        input_signature=(
            tf.TensorSpec(shape=[None, None], dtype=tf.float32),
            tf.TensorSpec(shape=[None, None], dtype=tf.float32),
        )
    )
    def match_fn(self, X: tf.Tensor, Y: tf.Tensor):

        dist_mat = l2_distance_matrix_tf(X, Y)

        row_matches = tf.argmin(dist_mat, 1)
        col_matches = tf.argmin(dist_mat, 0)

        num_rows = tf.cast(tf.shape(row_matches)[0], dtype=row_matches.dtype)

        inverse_row_indices = tf.gather(col_matches, row_matches)
        row_indices = tf.range(0, num_rows, dtype=row_matches.dtype)

        cross_checked = tf.equal(row_indices, inverse_row_indices)
        rows = row_indices[cross_checked]
        cols = row_matches[cross_checked]

        indices = tf.transpose(tf.stack([rows, cols]))
        distances = tf.gather_nd(dist_mat, indices)
        distances = tf.sqrt(tf.maximum(0.0, distances))
        return indices, distances

    def match(self, X: np.ndarray, Y: np.ndarray) -> MatchResult:
        X, Y = self.cast_inputs(X, Y)
        indices, distances = self.match_fn(X, Y)
        return MatchResult(indices.numpy(), distances.numpy())

class TFL2RTBFMatcher(TFL2CCBFMatcher):
    def __init__(self, ratio: float = 0.7, name: str = None):
        super(TFL2CCBFMatcher, self).__init__(name)
        self.ratio = ratio

    @tf.function(
        input_signature=(
            tf.TensorSpec(shape=[None, None], dtype=tf.float32),
            tf.TensorSpec(shape=[None, None], dtype=tf.float32),
        )
    )
    def match_fn(self, X: tf.Tensor, Y: tf.Tensor):

        dist_mat = l2_distance_matrix_tf(X, Y)

        neg_distances, indices = tf.nn.top_k(-dist_mat, 2)
        distances = -neg_distances
        ratio2 = self.ratio ** 2
        is_good = distances[:, 0] < ratio2 * distances[:, 1]

        num_rows = tf.shape(dist_mat)[0]
        row_indices = tf.range(0, num_rows, dtype=indices.dtype)

        rows = row_indices[is_good]
        cols = indices[is_good][:, 0]
        distances = distances[is_good][:, 0]

        indices = tf.transpose(tf.stack([rows, cols]))
        distances = tf.sqrt(tf.maximum(0.0, distances))
        return indices, distances

    def match(self, X: np.ndarray, Y: np.ndarray) -> MatchResult:
        X, Y = self.cast_inputs(X, Y)
        indices, distances = self.match_fn(X, Y)
        return MatchResult(indices.numpy(), distances.numpy())

max_stars_repo_name: kmkolasinski/fast-bfmatcher

max_issues_repo_name: kmkolasinski/fast-bfmatcher

max_forks_repo_name: kmkolasinski/fast-bfmatcher


text:
import numpy as np
import math

def hartmann(params):
    X = np.array([params['x'+str(i)] for i in range(6)]).reshape(1, -1)
    alpha = np.array([1.0, 1.2, 3.0, 3.2]).reshape(1, -1)
    A = np.array([
        [10, 3, 17, 3.5, 1.7, 8],
        [0.05, 10, 17, 0.1, 8, 14],
        [3, 3.5, 1.7, 10, 17, 8],
        [17, 8, 0.05, 10, 0.1, 14],
    ])
    P = np.array([
        [1312, 1696, 5569, 124, 8283, 5886],
        [2329, 4135, 8307, 3736, 1004, 9991],
        [2348, 1451, 3522, 2883, 3047, 6650],
        [4047, 8828, 8732, 5743, 1091, 381],
    ]) * 1e-4
    exponents = np.zeros(4)
    for i in range(4):
        for j in range(6):
            exponents[i] -= A[i, j] * np.power(X[0, j] - P[i, j], 2)
    ans = 0
    for i in range(4):
        ans -= alpha[0, i] * np.exp(exponents[i])

    return ans

# Write a function like this called 'main'
def main(job_id, params):
    print 'Anything printed here will end up in the output directory for job #%d' % job_id
    print params
    return hartmann(params)

max_stars_repo_name: piboauthors/PiBO-Spearmint

max_issues_repo_name: piboauthors/PiBO-Spearmint

max_forks_repo_name: piboauthors/PiBO-Spearmint