text:
export BrickIMUAcceleration
struct BrickIMUAcceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUMagneticField
struct BrickIMUMagneticField
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUAngularVelocity
struct BrickIMUAngularVelocity
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUAllData
struct BrickIMUAllData
    acc_x::Integer
    acc_y::Integer
    acc_z::Integer
    mag_x::Integer
    mag_y::Integer
    mag_z::Integer
    ang_x::Integer
    ang_y::Integer
    ang_z::Integer
    temperature::Integer
end

export BrickIMUOrientation
struct BrickIMUOrientation
    roll::Integer
    pitch::Integer
    yaw::Integer
end

export BrickIMUQuaternion
struct BrickIMUQuaternion
    x::Real
    y::Real
    z::Real
    w::Real
end

export BrickIMUSPITFPBaudrateConfig
struct BrickIMUSPITFPBaudrateConfig
    enable_dynamic_baudrate::Bool
    minimum_dynamic_baudrate::Integer
end

export BrickIMUSPITFPErrorCount
struct BrickIMUSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickIMUProtocol1BrickletName
struct BrickIMUProtocol1BrickletName
    protocol_version::Integer
    firmware_version::Vector{Integer}
    name::String
end

export BrickIMUIdentity
struct BrickIMUIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickIMU
"""
Full fledged AHRS with 9 degrees of freedom
"""
mutable struct BrickIMU <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickIMU(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.brick_imu")
        deviceInternal = package.BrickIMU(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the calibrated acceleration from the accelerometer for the
x, y and z axis.

If you want to get the acceleration periodically, it is recommended
to use the :cb:`Acceleration` callback and set the period with
:func:`Set Acceleration Period`.
"""
function get_acceleration(device::BrickIMU)
    return device.deviceInternal.get_acceleration()
end

export get_magnetic_field
"""
    $(SIGNATURES)

Returns the calibrated magnetic field from the magnetometer for the
x, y and z axis.

If you want to get the magnetic field periodically, it is recommended
to use the :cb:`Magnetic Field` callback and set the period with
:func:`Set Magnetic Field Period`.
"""
function get_magnetic_field(device::BrickIMU)
    return device.deviceInternal.get_magnetic_field()
end

export get_angular_velocity
"""
    $(SIGNATURES)

Returns the calibrated angular velocity from the gyroscope for the
x, y and z axis in °/14.375s (you have to divide by 14.375 to
get the value in °/s).

If you want to get the angular velocity periodically, it is recommended
to use the :cb:`Angular Velocity` callback and set the period with
:func:`Set Angular Velocity Period`.
"""
function get_angular_velocity(device::BrickIMU)
    return device.deviceInternal.get_angular_velocity()
end

export get_all_data
"""
    $(SIGNATURES)

Returns the data from :func:`Get Acceleration`, :func:`Get Magnetic Field`
and :func:`Get Angular Velocity` as well as the temperature of the IMU Brick.

If you want to get the data periodically, it is recommended
to use the :cb:`All Data` callback and set the period with
:func:`Set All Data Period`.
"""
function get_all_data(device::BrickIMU)
    return device.deviceInternal.get_all_data()
end

export get_orientation
"""
    $(SIGNATURES)

Returns the current orientation (roll, pitch, yaw) of the IMU Brick as Euler
angles. Note that Euler angles always experience a
`gimbal lock <https://en.wikipedia.org/wiki/Gimbal_lock>`__.

We recommend that you use quaternions instead.

The order to sequence in which the orientation values should be applied is
roll, yaw, pitch.

If you want to get the orientation periodically, it is recommended
to use the :cb:`Orientation` callback and set the period with
:func:`Set Orientation Period`.
"""
function get_orientation(device::BrickIMU)
    return device.deviceInternal.get_orientation()
end

export get_quaternion
"""
    $(SIGNATURES)

Returns the current orientation (x, y, z, w) of the IMU as
`quaternions <https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation>`__.

You can go from quaternions to Euler angles with the following formula::

 xAngle = atan2(2*y*w - 2*x*z, 1 - 2*y*y - 2*z*z)
 yAngle = atan2(2*x*w - 2*y*z, 1 - 2*x*x - 2*z*z)
 zAngle =  asin(2*x*y + 2*z*w)

This process is not reversible, because of the
`gimbal lock <https://en.wikipedia.org/wiki/Gimbal_lock>`__.

It is also possible to calculate independent angles. You can calculate
yaw, pitch and roll in a right-handed vehicle coordinate system according to
DIN70000 with::

 yaw   =  atan2(2*x*y + 2*w*z, w*w + x*x - y*y - z*z)
 pitch = -asin(2*w*y - 2*x*z)
 roll  = -atan2(2*y*z + 2*w*x, -w*w + x*x + y*y - z*z))

Converting the quaternions to an OpenGL transformation matrix is
possible with the following formula::

 matrix = [[1 - 2*(y*y + z*z),     2*(x*y - w*z),     2*(x*z + w*y), 0],
           [    2*(x*y + w*z), 1 - 2*(x*x + z*z),     2*(y*z - w*x), 0],
           [    2*(x*z - w*y),     2*(y*z + w*x), 1 - 2*(x*x + y*y), 0],
           [                0,                 0,                 0, 1]]

If you want to get the quaternions periodically, it is recommended
to use the :cb:`Quaternion` callback and set the period with
:func:`Set Quaternion Period`.
"""
function get_quaternion(device::BrickIMU)
    return device.deviceInternal.get_quaternion()
end

export get_imu_temperature
"""
    $(SIGNATURES)

Returns the temperature of the IMU Brick.
"""
function get_imu_temperature(device::BrickIMU)
    return device.deviceInternal.get_imu_temperature()
end

export leds_on
"""
    $(SIGNATURES)

Turns the orientation and direction LEDs of the IMU Brick on.
"""
function leds_on(device::BrickIMU)
    device.deviceInternal.leds_on()
end

export leds_off
"""
    $(SIGNATURES)

Turns the orientation and direction LEDs of the IMU Brick off.
"""
function leds_off(device::BrickIMU)
    device.deviceInternal.leds_off()
end

export are_leds_on
"""
    $(SIGNATURES)

Returns *true* if the orientation and direction LEDs of the IMU Brick
are on, *false* otherwise.
"""
function are_leds_on(device::BrickIMU)
    return device.deviceInternal.are_leds_on()
end

export set_acceleration_range
"""
    $(SIGNATURES)

Not implemented yet.
"""
function set_acceleration_range(device::BrickIMU, range)
    device.deviceInternal.set_acceleration_range(range)
end

export get_acceleration_range
"""
    $(SIGNATURES)

Not implemented yet.
"""
function get_acceleration_range(device::BrickIMU)
    return device.deviceInternal.get_acceleration_range()
end

export set_magnetometer_range
"""
    $(SIGNATURES)

Not implemented yet.
"""
function set_magnetometer_range(device::BrickIMU, range)
    device.deviceInternal.set_magnetometer_range(range)
end

export get_magnetometer_range
"""
    $(SIGNATURES)

Not implemented yet.
"""
function get_magnetometer_range(device::BrickIMU)
    return device.deviceInternal.get_magnetometer_range()
end

export set_convergence_speed
"""
    $(SIGNATURES)

Sets the convergence speed of the IMU Brick. The convergence speed
determines how the different sensor measurements are fused.

If the orientation of the IMU Brick is off by 10° and the convergence speed is
set to 20°/s, it will take 0.5s until the orientation is corrected. However,
if the correct orientation is reached and the convergence speed is too high,
the orientation will fluctuate with the fluctuations of the accelerometer and
the magnetometer.

If you set the convergence speed to 0, practically only the gyroscope is used
to calculate the orientation. This gives very smooth movements, but errors of the
gyroscope will not be corrected. If you set the convergence speed to something
above 500, practically only the magnetometer and the accelerometer are used to
calculate the orientation. In this case the movements are abrupt and the values
will fluctuate, but there won't be any errors that accumulate over time.

In an application with high angular velocities, we recommend a high convergence
speed, so the errors of the gyroscope can be corrected fast. In applications with
only slow movements we recommend a low convergence speed. You can change the
convergence speed on the fly. So it is possible (and recommended) to increase
the convergence speed before an abrupt movement and decrease it afterwards
again.

You might want to play around with the convergence speed in the Brick Viewer to
get a feeling for a good value for your application.
"""
function set_convergence_speed(device::BrickIMU, speed)
    device.deviceInternal.set_convergence_speed(speed)
end

export get_convergence_speed
"""
    $(SIGNATURES)

Returns the convergence speed as set by :func:`Set Convergence Speed`.
"""
function get_convergence_speed(device::BrickIMU)
    return device.deviceInternal.get_convergence_speed()
end

export set_calibration
"""
    $(SIGNATURES)

There are several different types that can be calibrated:

.. csv-table::
 :header: "Type", "Description", "Values"
 :widths: 10, 30, 110

 "0",    "Accelerometer Gain", "``[mul x, mul y, mul z, div x, div y, div z, 0, 0, 0, 0]``"
 "1",    "Accelerometer Bias", "``[bias x, bias y, bias z, 0, 0, 0, 0, 0, 0, 0]``"
 "2",    "Magnetometer Gain",  "``[mul x, mul y, mul z, div x, div y, div z, 0, 0, 0, 0]``"
 "3",    "Magnetometer Bias",  "``[bias x, bias y, bias z, 0, 0, 0, 0, 0, 0, 0]``"
 "4",    "Gyroscope Gain",     "``[mul x, mul y, mul z, div x, div y, div z, 0, 0, 0, 0]``"
 "5",    "Gyroscope Bias",     "``[bias xl, bias yl, bias zl, temp l, bias xh, bias yh, bias zh, temp h, 0, 0]``"

The calibration via gain and bias is done with the following formula::

 new_value = (bias + orig_value) * gain_mul / gain_div

If you really want to write your own calibration software, please keep
in mind that you first have to undo the old calibration (set bias to 0 and
gain to 1/1) and that you have to average over several thousand values
to obtain a usable result in the end.

The gyroscope bias is highly dependent on the temperature, so you have to
calibrate the bias two times with different temperatures. The values ``xl``,
``yl``, ``zl`` and ``temp l`` are the bias for ``x``, ``y``, ``z`` and the
corresponding temperature for a low temperature. The values ``xh``, ``yh``,
``zh`` and ``temp h`` are the same for a high temperatures. The temperature
difference should be at least 5°C. If you have a temperature where the
IMU Brick is mostly used, you should use this temperature for one of the
sampling points.

.. note::
 We highly recommend that you use the Brick Viewer to calibrate your
 IMU Brick.
"""
function set_calibration(device::BrickIMU, typ, data)
    device.deviceInternal.set_calibration(typ, data)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration for a given type as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickIMU, typ)
    return device.deviceInternal.get_calibration(typ)
end

export set_acceleration_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Acceleration` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_acceleration_period(device::BrickIMU, period)
    device.deviceInternal.set_acceleration_period(period)
end

export get_acceleration_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Acceleration Period`.
"""
function get_acceleration_period(device::BrickIMU)
    return device.deviceInternal.get_acceleration_period()
end

export set_magnetic_field_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Magnetic Field` callback is
triggered periodically. A value of 0 turns the callback off.
"""
function set_magnetic_field_period(device::BrickIMU, period)
    device.deviceInternal.set_magnetic_field_period(period)
end

export get_magnetic_field_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Magnetic Field Period`.
"""
function get_magnetic_field_period(device::BrickIMU)
    return device.deviceInternal.get_magnetic_field_period()
end

export set_angular_velocity_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Angular Velocity` callback is
triggered periodically. A value of 0 turns the callback off.
"""
function set_angular_velocity_period(device::BrickIMU, period)
    device.deviceInternal.set_angular_velocity_period(period)
end

export get_angular_velocity_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Angular Velocity Period`.
"""
function get_angular_velocity_period(device::BrickIMU)
    return device.deviceInternal.get_angular_velocity_period()
end

export set_all_data_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`All Data` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_all_data_period(device::BrickIMU, period)
    device.deviceInternal.set_all_data_period(period)
end

export get_all_data_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set All Data Period`.
"""
function get_all_data_period(device::BrickIMU)
    return device.deviceInternal.get_all_data_period()
end

export set_orientation_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Orientation` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_orientation_period(device::BrickIMU, period)
    device.deviceInternal.set_orientation_period(period)
end

export get_orientation_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Orientation Period`.
"""
function get_orientation_period(device::BrickIMU)
    return device.deviceInternal.get_orientation_period()
end

export set_quaternion_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Quaternion` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_quaternion_period(device::BrickIMU, period)
    device.deviceInternal.set_quaternion_period(period)
end

export get_quaternion_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Quaternion Period`.
"""
function get_quaternion_period(device::BrickIMU)
    return device.deviceInternal.get_quaternion_period()
end

export orientation_calculation_on
"""
    $(SIGNATURES)

Turns the orientation calculation of the IMU Brick on.

As default the calculation is on.

.. versionadded:: 2.0.2\$nbsp;(Firmware)
"""
function orientation_calculation_on(device::BrickIMU)
    device.deviceInternal.orientation_calculation_on()
end

export orientation_calculation_off
"""
    $(SIGNATURES)

Turns the orientation calculation of the IMU Brick off.

If the calculation is off, :func:`Get Orientation` will return
the last calculated value until the calculation is turned on again.

The trigonometric functions that are needed to calculate the orientation
are very expensive. We recommend to turn the orientation calculation
off if the orientation is not needed, to free calculation time for the
sensor fusion algorithm.

As default the calculation is on.

.. versionadded:: 2.0.2\$nbsp;(Firmware)
"""
function orientation_calculation_off(device::BrickIMU)
    device.deviceInternal.orientation_calculation_off()
end

export is_orientation_calculation_on
"""
    $(SIGNATURES)

Returns *true* if the orientation calculation of the IMU Brick
is on, *false* otherwise.

.. versionadded:: 2.0.2\$nbsp;(Firmware)
"""
function is_orientation_calculation_on(device::BrickIMU)
    return device.deviceInternal.is_orientation_calculation_on()
end

export set_spitfp_baudrate_config
"""
    $(SIGNATURES)

The SPITF protocol can be used with a dynamic baudrate. If the dynamic baudrate is
enabled, the Brick will try to adapt the baudrate for the communication
between Bricks and Bricklets according to the amount of data that is transferred.

The baudrate will be increased exponentially if lots of data is sent/received and
decreased linearly if little data is sent/received.

This lowers the baudrate in applications where little data is transferred (e.g.
a weather station) and increases the robustness. If there is lots of data to transfer
(e.g. Thermal Imaging Bricklet) it automatically increases the baudrate as needed.

In cases where some data has to transferred as fast as possible every few seconds
(e.g. RS485 Bricklet with a high baudrate but small payload) you may want to turn
the dynamic baudrate off to get the highest possible performance.

The maximum value of the baudrate can be set per port with the function
:func:`Set SPITFP Baudrate`. If the dynamic baudrate is disabled, the baudrate
as set by :func:`Set SPITFP Baudrate` will be used statically.

.. versionadded:: 2.3.5\$nbsp;(Firmware)
"""
function set_spitfp_baudrate_config(device::BrickIMU, enable_dynamic_baudrate, minimum_dynamic_baudrate)
    device.deviceInternal.set_spitfp_baudrate_config(enable_dynamic_baudrate, minimum_dynamic_baudrate)
end

export get_spitfp_baudrate_config
"""
    $(SIGNATURES)

Returns the baudrate config, see :func:`Set SPITFP Baudrate Config`.

.. versionadded:: 2.3.5\$nbsp;(Firmware)
"""
function get_spitfp_baudrate_config(device::BrickIMU)
    return device.deviceInternal.get_spitfp_baudrate_config()
end

export get_send_timeout_count
"""
    $(SIGNATURES)

Returns the timeout count for the different communication methods.

The methods 0-2 are available for all Bricks, 3-7 only for Master Bricks.

This function is mostly used for debugging during development, in normal operation
the counters should nearly always stay at 0.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_send_timeout_count(device::BrickIMU, communication_method)
    return device.deviceInternal.get_send_timeout_count(communication_method)
end

export set_spitfp_baudrate
"""
    $(SIGNATURES)

Sets the baudrate for a specific Bricklet port.

If you want to increase the throughput of Bricklets you can increase
the baudrate. If you get a high error count because of high
interference (see :func:`Get SPITFP Error Count`) you can decrease the
baudrate.

If the dynamic baudrate feature is enabled, the baudrate set by this
function corresponds to the maximum baudrate (see :func:`Set SPITFP Baudrate Config`).

Regulatory testing is done with the default baudrate. If CE compatibility
or similar is necessary in your applications we recommend to not change
the baudrate.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function set_spitfp_baudrate(device::BrickIMU, bricklet_port, baudrate)
    device.deviceInternal.set_spitfp_baudrate(bricklet_port, baudrate)
end

export get_spitfp_baudrate
"""
    $(SIGNATURES)

Returns the baudrate for a given Bricklet port, see :func:`Set SPITFP Baudrate`.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_spitfp_baudrate(device::BrickIMU, bricklet_port)
    return device.deviceInternal.get_spitfp_baudrate(bricklet_port)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Brick side. All
Bricklets have a similar function that returns the errors on the Bricklet side.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_spitfp_error_count(device::BrickIMU, bricklet_port)
    return device.deviceInternal.get_spitfp_error_count(bricklet_port)
end

export enable_status_led
"""
    $(SIGNATURES)

Enables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function enable_status_led(device::BrickIMU)
    device.deviceInternal.enable_status_led()
end

export disable_status_led
"""
    $(SIGNATURES)

Disables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function disable_status_led(device::BrickIMU)
    device.deviceInternal.disable_status_led()
end

export is_status_led_enabled
"""
    $(SIGNATURES)

Returns *true* if the status LED is enabled, *false* otherwise.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function is_status_led_enabled(device::BrickIMU)
    return device.deviceInternal.is_status_led_enabled()
end

export get_protocol1_bricklet_name
"""
    $(SIGNATURES)

Returns the firmware and protocol version and the name of the Bricklet for a
given port.

This functions sole purpose is to allow automatic flashing of v1.x.y Bricklet
plugins.
"""
function get_protocol1_bricklet_name(device::BrickIMU, port)
    return device.deviceInternal.get_protocol1_bricklet_name(port)
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has an
accuracy of ±15%. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickIMU)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Brick. Calling this function
on a Brick inside of a stack will reset the whole stack.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickIMU)
    device.deviceInternal.reset()
end

export write_bricklet_plugin
"""
    $(SIGNATURES)

Writes 32 bytes of firmware to the bricklet attached at the given port.
The bytes are written to the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_bricklet_plugin(device::BrickIMU, port, offset, chunk)
    device.deviceInternal.write_bricklet_plugin(port, offset, chunk)
end

export read_bricklet_plugin
"""
    $(SIGNATURES)

Reads 32 bytes of firmware from the bricklet attached at the given port.
The bytes are read starting at the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function read_bricklet_plugin(device::BrickIMU, port, offset)
    return device.deviceInternal.read_bricklet_plugin(port, offset)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Brick is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position is the position in the stack from '0' (bottom) to '8' (top).

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickIMU)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickIMU, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using DBI, PostgreSQL

connect(Postgres) do conn
  stmt = prepare(conn,"SELECT * FROM \"MediaType\"");
  res = execute(stmt);
  for row in res
	@printf "ID: %2d, Media: %s\n" row[1] row[2]
  end
  finish(stmt);
end

ext: jl

max_stars_repo_path: Module 3/Chapter08/postgres-dbi.jl

max_issues_repo_path: Module 3/Chapter08/postgres-dbi.jl

max_forks_repo_path: Module 3/Chapter08/postgres-dbi.jl


text:
module Auth

# use both cookie and header to validate user incase some server stripped out
# the cookie e.g. Google server then use the header to stored the verified user
# with the validate token

using StructTypes, Dates, JWTs, HTTP
using ..Model

const JWT_AUTH_KEYS = Ref{JWKSet}()

function init(authkeysfile)
    JWT_AUTH_KEYS[] = JWKSet(authkeysfile)
    refresh!(JWT_AUTH_KEYS[])
    return
end

# generating a new token for a user
const DATE_FORMAT = DateFormat("e, dd u yyyy HH:MM:SS G\\MT") # Wed, 21 Oct 2015 07:28:00 GMT
const JWT_TOKEN_COOKIE_NAME = "X-MusicAlbums-Jwt-Token"

function addtoken!(resp::HTTP.Response, user::User)
    exp = Dates.now(Dates.UTC) + Dates.Hour(12)
    payload = Dict("iss"=>"MusicAlbums.jl", "exp"=>Dates.datetime2unix(exp), "sub"=>"managing albums", "aud"=>user.username, "uid"=>user.id)
    jwt = JWT(; payload=payload)
    keyid = first(first(JWT_AUTH_KEYS[].keys))
    sign!(jwt, JWT_AUTH_KEYS[], keyid)
    # handle the case where Google cloud stripped out all cookies to their server automatically
    # check not just the cookie but header to see if it valid
    HTTP.setheader(resp, "Set-Cookie" => "$JWT_TOKEN_COOKIE_NAME=$(join([jwt.header, jwt.payload, jwt.signature], '.')); Expires=$(Dates.format(exp, DATE_FORMAT))")
    HTTP.setheader(resp, JWT_TOKEN_COOKIE_NAME => join([jwt.header, jwt.payload, jwt.signature], '.'))
    return resp
end

struct Unauthenticated <: Exception end

# parsing token from HTTP request
# and validate the user
# then add the validate token to the cookie and header
function User(req::HTTP.Request)
    if HTTP.hasheader(req, "Cookie")
        cookies = filter(x->x.name == JWT_TOKEN_COOKIE_NAME, HTTP.cookies(req))
        if !isempty(cookies) && !isempty(cookies[1].value)
            jwt = JWT(; jwt=cookies[1].value)
            verified = false                  # start to validate the user cookie
            for kid in JWT_AUTH_KEYS[].keys
                validate!(jwt, JWT_AUTH_KEYS[], kid[1])
                verified |= isverified(jwt)
            end
            if verified
                parts = claims(jwt)  # user is good add token
                return User(parts["uid"], parts["aud"])
            end
        end
    elseif HTTP.hasheader(req, JWT_TOKEN_COOKIE_NAME)
        jwt = JWT(; jwt=String(HTTP.header(req, JWT_TOKEN_COOKIE_NAME)))
        verified = false                # in case cookie is stripped out use the header
        for kid in JWT_AUTH_KEYS[].keys
            validate!(jwt, JWT_AUTH_KEYS[], kid[1])
            verified |= isverified(jwt)
        end
        if verified
            parts = claims(jwt)   # ssame thing with the cookie
            return User(parts["uid"], parts["aud"])
        end
    end
    throw(Unauthenticated())
end

end # module

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
using Slurm
using Base.Test

# write your own tests here

include("hostlist.jl")


text:
abstract type AbstractNBodySystem end
abstract type AbstractInteraction end
abstract type AbstractSummationAlgorithm{T} end

struct Pairwise{isparallel} <: AbstractSummationAlgorithm{isparallel} end
struct PairwiseSymmetric{isparallel} <: AbstractSummationAlgorithm{isparallel} end
struct PairList{isparallel} <: AbstractSummationAlgorithm{isparallel} end

function number_of_particles(::AbstractNBodySystem, ::AbstractInteraction) end
function interactions(::AbstractNBodySystem) end
function pairlist(::AbstractNBodySystem, ::AbstractInteraction) end
function initial_conditions(::AbstractNBodySystem) end
function atype(::AbstractNBodySystem) end
function force!(dv, ::AbstractInteraction, ::AbstractNBodySystem) end
function force!(dv, ::AbstractInteraction, ::AbstractNBodySystem, i, j) end
function particle_mass(::AbstractNBodySystem, i) end
function summation_alg(::AbstractInteraction) end
function issparse(::AbstractInteraction) end

boundary_conditions!(s) = CallbackSet()
thermostats!(s) = CallbackSet()
neighbours(s) = CallbackSet()


text:
# this is used to call defns and thms quickly. just an idea but might build out more

function unitary()
    print("""A complex square matrix U is unitary if U†U == UU† == Id

    Properties of unitary matricies:

    * diagonalizable (SVD)
    * norm-preserving
    * invertible
    """)
end

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
generic_sphere(ref_dirs) = ref_dirs ./ ([norm(v) for v in ref_dirs])

function DTLZ_hypersphere!(fx, x; α = 1)
    m = length(fx)
    for i = 1:m
        if i < m
            fx[i] *= prod(cos.((π*0.5) * x[1:m-i] .^ α ))
        end
        if i > 1
            fx[i] *= sin((π*0.5) * x[1+m - i] .^ α)
        end
    end

    return fx
end

function DTLZ_g1(x, m)
    y = view(x, m:length(x)) .- 0.5
    return 100*( length(y) + sum( y.^2 - cos.(20π*y) ))
end

function DTLZ_g2(x, m)
    return sum( (view(x, m:length(x)) .- 0.5).^2 )
end

"""
    DTLZ1(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ1 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- convex
- multifrontal
"""
function DTLZ1(m=3, ref_dirs = gen_ref_dirs(m, 12))
    f(x,m=m) = begin
        g = DTLZ_g1(x, m)
        D = length(x)

        fx = fill(0.5*(1 + g), m)
        for i in 1:m
            fx[i] *= prod(x[1:m-i])
            if i > 1
                fx[i] *= 1 - x[1+m - i]
            end
        end

        return fx, [0.0], [0.0]
    end

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = 0.5ref_dirs
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return f, bounds, pareto_set
end

"""
    DTLZ2(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ2 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
"""
function DTLZ2(m=3, ref_dirs = gen_ref_dirs(m, 12))
    f(x,m=m) = begin
        g = DTLZ_g2(x, m)
        fx = fill(1.0 + g, m)
        DTLZ_hypersphere!(fx, x)

        return fx, [0.0], [0.0]
    end

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return f, bounds, pareto_set
end

"""
    DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ3 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- multifrontal
"""
function DTLZ3(m = 3, ref_dirs = gen_ref_dirs(m, 12))
    f(x,m=m) = begin
        g = DTLZ_g1(x, m)
        fx = fill(1 + g, m)
        DTLZ_hypersphere!(fx, x)

        return fx, [0.0], [0.0]
    end

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return f, bounds, pareto_set
end

"""
    DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))

DTLZ4 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions
- `ref_dirs` number of Pareto solutions (default: Das and Dennis' method).

Main properties:
- nonconvex
- unifrontal
"""
function DTLZ4(m = 3, ref_dirs = gen_ref_dirs(m, 12))
    f(x,m=m) = begin
        g = DTLZ_g2(x, m)
        fx = fill(1.0 + g, m)
        DTLZ_hypersphere!(fx, x; α = 100)

        return fx, [0.0], [0.0]
    end

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    pf = generic_sphere(ref_dirs)
    pareto_set = [ generateChild(zeros(0), (fx, [0.0], [0.0])) for fx in pf ]

    return f, bounds, pareto_set
end

"""
    DTLZ5(m = 3)

DTLZ5 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions

"""
function DTLZ5(m = 3)
    f(x,m=m) = begin
        g = DTLZ_g2(x, m)
        fθ = fill(1.0 + g, m)
        θ = @. 1 / (2*(1 + g)) * (1 + 2g * x[1:m-1] )
        θ[1] = x[1]

        DTLZ_hypersphere!(fθ, θ)

        return fθ, [0.0], [0.0]
    end

    D = 10 + m - 1

    bounds = Array([zeros(D) ones(D)]')

    if m == 3
        ref_dirs = gen_ref_dirs(2, 100)
        X = fill(0.5, length(ref_dirs), D)

        for i in eachindex(ref_dirs)
            X[i,1:2] = ref_dirs[i]
        end
        pareto_set = [ generateChild(X[i,:], f(X[i,:])) for i in eachindex(ref_dirs) ]

    else
        pareto_set = []
    end

    return f, bounds, pareto_set
end

"""
    DTLZ6(m = 3)

DTLZ6 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `m` number of objective functions

"""
function DTLZ6(m = 3)
    f(x,m=m) = begin
        g = sum(x[m:end] .^ 0.1)
        fθ = fill(1.0 + g, m)
        θ = @. 1 / (2*(1 + g)) * (1 + 2g * x[1:m-1] )
        θ[1] = x[1]

        DTLZ_hypersphere!(fθ, θ)

        return fθ, [0.0], [0.0]
    end

    D = 10

    bounds = Array([zeros(D) ones(D)]')

    if m == 3
        ref_dirs = gen_ref_dirs(2, 100)
        X = fill(0.0, length(ref_dirs), D)

        for i in eachindex(ref_dirs)
            X[i,1:2] = ref_dirs[i]
        end
        pareto_set = [ generateChild(X[i,:], f(X[i,:])) for i in eachindex(ref_dirs) ]

    else
        pareto_set = []
    end

    return f, bounds, pareto_set
end


text:
# different parameterizations of the spin system FID parameters.

struct SpinSysParamsType2{T}
    κs_λ::Vector{Vector{T}} # a multiplier for each (spin group, partition element).
    κs_β::Vector{Vector{T}} # a vector coefficient for each (spin group). vector length: number of spins in the spin group.
    d::Vector{Vector{T}} # a multiplier for each (spin group, partition element).
end

function SpinSysParamsType2(x::T) where T
    return SpinSysParamsType2(Vector{Vector{T}}(undef,0), Vector{Vector{T}}(undef, 0), Vector{Vector{T}}(undef, 0))
end

struct SpinSysParamsType1{T}
    κs_λ::Vector{T} # a common multiplier for each spin group. length: number of spin groups.
    κs_β::Vector{Vector{T}} # a vector coefficient for each (spin group). vector length: number of spins in the spin group.
    d::Vector{T} # a common multiplier for each spin group. length: number of spin groups.
end

function SpinSysParamsType1(x::T) where T
    return SpinSysParamsType1(Vector{T}(undef,0), Vector{Vector{T}}(undef, 0), Vector{T}(undef, 0))
end

function constructorSSFID(x::SpinSysParamsType1{T}, y...)::SpinSysParamsType1{T} where T
    return SpinSysParamsType1(y...)
end

# function constructorSSFID(x::SpinSysFIDType2{T}, y...)::SpinSysFIDType2{T} where T
#     return SpinSysFIDType2(y...)
# end

struct FIDModelType{T,SST} # parameters for surrogate model.

    # non-singlet spin systems.
    qs::Vector{Vector{Function}} # spin group, partition element index.

    # κs_λ::Vector{T} # multipliers for each spin group.
    # κs_β::Vector{Vector{T}} # coefficients for each (spin group, partition element).
    # d::Vector{T}
    ss_params::SST

    # singlets.
    κs_λ_singlets::Vector{T}
    β_singlets::Vector{T}
    d_singlets::Vector{T}

    # misc.
    Δcs_max::Vector{T}
    λ0::T
end

mutable struct καFIDModelType{T,SST}
    κs_α::Vector{Vector{T}} # spin group, partition element index.
    κs_α_singlets::Vector{T}
    core::FIDModelType{T,SST}
end

function καFIDModelType(core::FIDModelType{T,SST}) where {T,SST}

    N_spins = length(core.qs)
    κs_α = Vector{Vector{T}}(undef, N_spins)
    for i = 1:N_spins
        κs_α[i] = ones(T, length(core.qs[i]))
    end

    κs_α_singlets = ones(T, length(core.d_singlets))

    return καFIDModelType(κs_α, κs_α_singlets, core)
end

######

struct SHType{T} # output of the SH simulation.

    # resonance components in non-singlet spin systems.
    αs::Vector{Vector{T}}
    Ωs::Vector{Vector{T}}

    Δc_m_compound::Vector{Vector{Vector{T}}}
    part_inds_compound::Vector{Vector{Vector{Int}}}
    Δc_bar::Vector{Vector{Vector{T}}}

    N_spins_sys::Vector{Int}

    # resonance components in singlets spin systems.
    αs_singlets::Vector{T}
    Ωs_singlets::Vector{T}

    # misc input for the SH simulation.
    fs::T
    SW::T
    ν_0ppm::T
end

"""
Remove `Δc_m_compound` and `part_inds_compound` to reduce the memory footprint of Vector{CompoundFIDType} types.
"""
function removeauxinfo!(As::Vector{SHType{T}}) where T

    for n = 1:length(As)
        A = As[n]
        for i = 1:length(A.Δc_m_compound)
            resize!(A.Δc_m_compound[i], 0)
            resize!(A.part_inds_compound[i], 0)
        end
    end

    return nothing
end

function varinfocompound(A)

    property_names = propertynames(A)

    for j = 1:length(property_names)
        display(property_names[j])
        asdf = getfield(A, property_names[j]);

        info_string = varinfo(r"asdf")
        display(info_string)
    end
end


text:
# Kronecker power utils
export kron_id,
       kron_sandwich,
       kron_sum,
       kron_pow

# error bounds
export error_bound_apriori,
       convergence_radius_apriori,
       error_bound_pseries,
       convergence_radius_pseries,
       error_bound_specabs,
       convergence_radius_specabs

# linearization
export build_matrix


text:
#=
tuple:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module Tuple
    # tuple
    t = (1, 2, 3)

    # named tuple
    t_n = (a = 1, b = 2, c = 3)

    # array
    array = rand(4, 3)

    # methods
    function get_type(tuple)
        return typeof(tuple)
    end

    function get_val(tuple, idx)
        return tuple[idx]
    end

    function get_tuple(x...)
        return x
    end

    function main()
        println("Type of tuple, t = ", get_type(t))
        println("Type of tuple, t_n = ", get_type(t_n))
        println("Value of t[1] = ", get_val(t, 1))
        println("Value of t_n[1] = ", t_n[1])
        println("Value of t_n.a = ", t_n.a)
        println("Value of t_n[:a] = ", t_n[:a])
        println("Keys of t_n = ", keys(t_n))
        println("Values of t_n = ", values(t_n))
        println("rand(4, 3) = ", array)
        println("Size of rand(4, 3) = ", size(array))
        println(get_tuple(1, 2, 3))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Tuple
    Tuple.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
####Newsvendor Factories

#simple discrete quantile
function nv_quantile(pk, supp_k, s)
	if s < 0
		return supp_k[1]
	end

	cdf_ = 0.
    for i = 1:length(pk)
        cdf_ += pk[i]
        if cdf_ >= s
        	return supp_k[i]
        end
    end
    return supp_k[end]
end

#returns a d x K matrix of univariate functions
function getNewsVendorCosts(supps, s, K)
	return [x-> supps[i, k] > x ? s/(1-s) * (supps[i, k] - x) : (x - supps[i, k])
				for i =1:size(supps, 1), k =1:K]
end

#returns an K array of functions which train S-SAA
#hyperparams = (p0, alpha)
function genSSAAtrainers(supps, s)
	#Generic computation of the sth quantile
	function x_k(mhat_k, k, s, p0, alpha, supps, palpha)
	    Nhat_k = sum(mhat_k)
	    JS._shrink!(mhat_k, p0, alpha, Nhat_k, palpha)
		nv_quantile(palpha, supps[:, k], s)
	end
	in_place_pk = zeros(size(supps, 1)) #allocate the shared workspace once here
	return [(mhat_k, (p0, alpha))-> x_k(mhat_k, k, s, p0, alpha, supps, in_place_pk) for k = 1:size(supps, 2)]
end

function genKSTrainers(supps, s, K, method)
	#Computes KS Order quantity assuming Gamma small enough
	function x_k(mhat_k, k, s, Gamma::Number)
	    Nhat_k = sum(mhat_k)
	    if Nhat_k == 0  #when there's no data, return maximum
	    	return supps[end, k]
	    end
	    @assert Gamma <= 1-s "KS Newsvendor not support Gamma > 1-s: $Gamma > $(1-s)"
	    phat = mhat_k./Nhat_k
	    qlo = nv_quantile(phat, supps[:, k], ceil(Nhat_k * (s - Gamma))/Nhat_k)
	    qhigh = nv_quantile(phat, supps[:, k], floor(Nhat_k * (s + Gamma) + 1)/Nhat_k)
	    return s * qhigh + (1 - s) * qlo
	end

	#does K-fold cross-validation to determine Gamma
	cs = getNewsVendorCosts(supps, s, K)
	function x_k(mhat_k, k, s, Gamma_grid, numFolds)
	    Nhat_k = round(Int, sum(mhat_k))
	    #divy the data into cross-val sets
    	splits = mod.(1:Nhat_k, numFolds) .+ 1
    	shuffle!(splits)
    	cv_data = [zero(mhat_k) for fold = 1:numFolds]

    	ix_split = 1
        for i = 1:length(mhat_k)
            for j = 1:round(Int, mhat_k[i])
                cv_data[splits[ix_split]][i] += 1
                ix_split += 1
            end
        end

	    #for each set train and evaluate
	    out = zero(Gamma_grid)
	    for (ix, Gamma) in enumerate(Gamma_grid)
		    for fold = 1:numFolds
		    	train_data = sum(cv_data[setdiff(1:numFolds, fold)])
		    	Nhat_train = sum(train_data)
		    	x = x_k(train_data, k, s, Gamma)
	    		costs = [c(x) for c in cs[:, k]]
	    		out[ix] += dot(cv_data[fold], costs) / (Nhat_k - Nhat_train)
	    	end
	    	out[ix] /= numFolds
	    end

	    #identify who is best and retrain
	    Gammastar = Gamma_grid[argmin(out)]
	    #println(k, "\t", Gammastar)
	    x_k(mhat_k, k, s, Gammastar)
	end

	if method == :aPriori
		xs  = [(mhat_k, Gamma)-> x_k(mhat_k, k, s, Gamma) for k = 1:K]
	elseif method == :crossVal
		xs = [(mhat_k, (Gamma_grid, numFolds)) -> x_k(mhat_k, k, s, Gamma_grid, numFolds) for k = 1:K]
	else
		throw("Method must be one of :aPriori or :crossVal")
	end

	xs
end


text:
using DependencyParser.DependencyParsing

@testset "Model" begin
  @testset "initialize untrained with embeddings" begin
    connlu_file = "test/fixtures/connlu_test.txt"
    embeddings_file = "test/fixtures/embeddings_test.vec"
    connlu_sentences = load_connlu_file(connlu_file)
    system = ArcEager.ArcEagerSystem()
    settings = Settings(embeddings_size = 4, hidden_size = 2, batch_size = 3)

  model = Model(settings, system, embeddings_file, connlu_sentences)
    @test length(model.embeddings[1, :]) == settings.embeddings_size
    @test length(model.hidden_layer.weight[:, 1]) == settings.hidden_size
    @test length(model.hidden_layer.weight[1, :]) == settings.batch_size * settings.embeddings_size
    @test length(model.output_layer.weight[:, 1]) == ArcEager.transitions_number(system)
  end

  @testset "Write to file" begin
    filename = "tmp/modelfile.txt"
    settings = Settings(embeddings_size = 4, hidden_size = 2, batch_size = 3)
    labels_num = 10
    model = build_model(settings)

    words_count = length(model.word_ids)
    tags_count = length(model.tag_ids)
    labels_count = length(model.label_ids)
    total_ids_count = words_count + tags_count + labels_count

    total_lines_count = total_ids_count + settings.hidden_size + labels_num + 3

    write_to_file!(model, filename)

    word = first(model.word_ids)
    tag = first(model.tag_ids)
    label = first(model.label_ids)

    # manualy count from model and settings
    correct_embeddings_info_line = "6 6 3 4"
    correct_layers_info_line = "3 2 10"

    word_correct_line = "$(word[begin]) $(join(model.embeddings[word[end], :], " "))"
    tag_correct_line = "$(tag[begin]) $(join(model.embeddings[tag[end], :], " "))"
    label_correct_line = "$(label[begin]) $(join(model.embeddings[label[end], :], " "))"

    hidden_weight_correct_line = join(model.hidden_layer.weight[begin, :], " ")
    hidden_bias_correct_line = join(model.hidden_layer.bias, " ")
    softmax_weight_correct_line = join(model.output_layer.weight[begin, :], " ")

    lines = readlines(filename)

    @test length(lines) == total_lines_count
    @test lines[begin] == correct_embeddings_info_line
    @test lines[begin + 1] == word_correct_line
    @test lines[begin + length(model.word_ids) + 1] == tag_correct_line
    @test lines[begin + length(model.word_ids) + length(model.tag_ids) + 1] == label_correct_line
    @test lines[begin + length(model.word_ids) + length(model.tag_ids) + length(model.label_ids) + 1] == correct_layers_info_line
    @test lines[begin + length(model.word_ids) + length(model.tag_ids) + length(model.label_ids) + 2] == hidden_weight_correct_line
    @test lines[begin + length(model.word_ids) + length(model.tag_ids) + length(model.label_ids) + settings.hidden_size + 2] == hidden_bias_correct_line
    @test lines[begin + length(model.word_ids) + length(model.tag_ids) + length(model.label_ids) + settings.hidden_size + 3] == softmax_weight_correct_line

    rm(filename)
  end

  @testset "Read from file" begin
    filename = "test/fixtures/modelfile_test.txt"

    model = Model(filename)

    # manualy calculated from fixture model
    word_pair = Pair{String, Integer}("here", 1)
    tag_pair = Pair{String, Integer}("TAG2", 7)
    label_pair = Pair{String, Integer}("LABEL2", 13)
    word_embedding = Float64[1, 2, 3, 4]
    tag_embedding = Float64[3, 2, 4, 1]
    label_embedding = Float64[4, 2, 1, 3]
    hidden_size = 2
    batch_size = 3
    labels_num = 10
    hidden_layer_weight_vector = Float64[1,2,3,4,5,6,7,8,9,10,9,8]
    hidden_layer_bias = Float64[1, 2]
    softmax_layer_weight_vector = [3, 4]

    @test length(model.word_ids) == 6
    @test length(model.tag_ids) == 6
    @test length(model.label_ids) == 3
    @test haskey(model.word_ids, word_pair[begin]) && model.word_ids[word_pair[begin]] == word_pair[end]
    @test haskey(model.tag_ids, tag_pair[begin]) && model.tag_ids[tag_pair[begin]] == tag_pair[end]
    @test haskey(model.label_ids, label_pair[begin]) && model.label_ids[label_pair[begin]] == label_pair[end]
    @test model.embeddings[word_pair[end], :] == word_embedding
    @test model.embeddings[tag_pair[end], :] == tag_embedding
    @test model.embeddings[label_pair[end], :] == label_embedding
    @test length(model.hidden_layer.bias) == hidden_size
    @test length(model.output_layer.weight[begin, :]) == hidden_size
    @test length(model.output_layer.weight[:, begin]) == labels_num
    @test model.hidden_layer.bias == hidden_layer_bias
    @test model.hidden_layer.weight[begin, :] == hidden_layer_weight_vector
    @test model.output_layer.weight[begin, :] == softmax_layer_weight_vector
  end

  @testset "Model prediction" begin
    parser = build_dep_parser()
    config = build_configuration()
    model = parser.model
    batch = build_correct_batch(parser.model, config)

    prediction = predict(model, batch)

    @test length(prediction) == ArcEager.transitions_number(parser.parsing_system)
    @test round(sum(prediction)) == 1
  end

  @testset "Model update" begin
    parser = build_dep_parser()
    config = build_configuration()
    gold_state = build_updated_gold_state(config, system=parser.parsing_system)
    model = parser.model

    batch = build_correct_batch(parser.model, config)
    gold = ArcEager.transition_costs(gold_state)

    loss_before_update = loss_function(model, batch, gold)
    update_model!(model, batch, gold)
    loss_after_update = loss_function(model, batch, gold)

    @test loss_after_update < loss_before_update
  end

  @testset "Test form batch" begin
    parser = build_dep_parser()
    config = build_configuration()
    model = parser.model

    correct_batch = build_correct_batch(model, config)
    @test form_batch(model, parser.settings, config) == correct_batch
  end
end

# Temp used for test

# parser = build_dep_parser()
# config = build_configuration()
# gold_state = build_updated_gold_state(config, system=parser.parsing_system)
# model = parser.model

# batch = build_correct_batch(parser.model, config)
# costs = ArcEager.transition_costs(gold_state)
# gold = ArcEager.transition_costs(gold_state) |> gold -> map(c -> -c, gold) |> softmax

# model.hidden_layer.weight
# model.hidden_layer.bias

# model.chain

# params = Flux.params(model.chain)

# params

# loss_before_update = loss_function(model, batch, gold)
# update_model!(model, batch, costs)
# update_model!_2(model, batch, costs)

# predict(model, batch)
# gold

model_test = DependencyParser.DependencyParsing.Model("tmp/test12.txt")

words = ["He", "wrote", "her", "a", "letter", "."]
tags = ["PRP", "VBN", "PRP", "DT", "NN", "JJ"]

words_with_tags = zip(words, tags) |> collect

parser = DependencyParser.DependencyParsing.DepParser(DependencyParser.DependencyParsing.Settings(), model_test, DependencyParser.DependencyParsing.ArcEager.ArcEagerSystem())

sentence = DependencyParser.DependencyParsing.Sentence(words_with_tags)
config = DependencyParser.DependencyParsing.Configuration(sentence)

parser(words_with_tags)

batch = form_batch(parser.model, Settings(), config)
pr = DependencyParser.DependencyParsing.predict(parser.model, batch)
transition = predict_transition(parser, config)
execute_transition(parser, config, transition)
config

using TextAnalysis
using TextModels

text = StringDocument("Any use of the work other than as authorized under this license or copyright law is prohibited")
words = tokens(text)
pos = TextModels.PerceptronTagger(true)
words_with_tags = pos(words) |> words -> map(word -> (String(word[begin]), String(word[end])), words)
parser(words_with_tags)

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
#Sparse Ensemble Kalman Inversion: specific structures and function definitions
using Convex, SCS
using SparseArrays

"""
    SparseInversion <: Process

A sparse ensemble Kalman Inversion process
"""
struct SparseInversion{FT <: AbstractFloat, IT <: Int} <: Process
    "upper limit of l1-norm"
    γ::FT
    "the flag to enable thresholding"
    threshold_eki::Bool
    "a small value to threshold results that are close to zero"
    threshold_value::FT
    "a small value for regularization to enhance robustness of convex optimization"
    reg::FT
    "a list of index to indicate the parameters included in the evaluation of l1-norm"
    uc_idx::Vector{IT}
end

"""
    sparse_qp(
        ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT,IT}},
        v_j::Vector{FT},
        cov_vv_inv::Array{FT, 2},
        H_u::SparseArrays.SparseMatrixCSC{FT},
        H_g::SparseArrays.SparseMatrixCSC{FT},
        y_j::Vector{FT};
        H_uc::SparseArrays.SparseMatrixCSC{FT} = H_u,
    ) where {FT, IT}

Solving quadratic programming problem with sparsity constraint.
"""
function sparse_qp(
    ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT, IT}},
    v_j::Vector{FT},
    cov_vv_inv::Array{FT, 2},
    H_u::SparseArrays.SparseMatrixCSC{FT},
    H_g::SparseArrays.SparseMatrixCSC{FT},
    y_j::Vector{FT};
    H_uc::SparseArrays.SparseMatrixCSC{FT} = H_u,
) where {FT, IT}

    P = H_g' * (ekp.obs_noise_cov \ H_g) + cov_vv_inv
    P = 0.5 * (P + P')
    q = -(cov_vv_inv * v_j + H_g' * (ekp.obs_noise_cov \ y_j))
    N_params = size(H_uc)[1]
    P1 = vcat(
        hcat(P, fill(FT(0), size(P)[1], N_params)),
        hcat(fill(FT(0), N_params, size(P)[1]), fill(FT(0), N_params, N_params)),
    )
    q1 = vcat(q, fill(FT(0), N_params, 1))
    H_uc_abs = 1.0 * I(N_params)
    G = hcat(vcat(H_uc, -1.0 * H_uc), vcat(-1.0 * H_uc_abs, -1.0 * H_uc_abs))
    h = fill(FT(0), 2 * N_params, 1)
    G1 = vcat(G, hcat(fill(FT(0), 1, size(P)[1]), fill(FT(1), 1, N_params)))
    h1 = vcat(h, ekp.process.γ)
    x = Variable(size(P1)[1])
    problem = minimize(0.5 * quadform(x, P1; assume_psd = true) + q1' * x, [G1 * x <= h1])
    solve!(problem, () -> SCS.Optimizer(verbose = false))

    return hcat(H_u, fill(FT(0), size(H_u)[1], N_params)) * evaluate(x)
end

"""
    update_ensemble!(
        ekp::EnsembleKalmanProcess{FT, IT, <:SparseInversion{FT,IT}},
        g::Array{FT,2};
        cov_threshold::FT=0.01,
        Δt_new=nothing
    ) where {FT, IT}

Updates the ensemble according to which type of Process we have. Model outputs `g` need to be a `N_obs × N_ens` array (i.e data are columms).
"""
function update_ensemble!(
    ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT, IT}},
    g::Array{FT, 2};
    cov_threshold::FT = 0.01,
    Δt_new = nothing,
    deterministic_forward_map = true,
) where {FT, IT}

    # Update follows eqns. (4) and (5) of Schillings and Stuart (2017)

    #catch works when g non-square
    if !(size(g)[2] == ekp.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g do not match, try transposing g or check ensemble size"))
    end

    # u: N_par × N_ens
    # g: N_obs × N_ens
    u = get_u_final(ekp)
    N_obs = size(g, 1)

    cov_init = cov(u, dims = 2)
    cov_ug = cov(u, g, dims = 2, corrected = false) # [N_par × N_obs]
    cov_gg = cov(g, g, dims = 2, corrected = false) # [N_obs × N_obs]
    cov_uu = cov(u, u, dims = 2, corrected = false) # [N_par × N_par]

    if !isnothing(Δt_new)
        push!(ekp.Δt, Δt_new)
    elseif isnothing(Δt_new) && isempty(ekp.Δt)
        push!(ekp.Δt, FT(1))
    else
        push!(ekp.Δt, ekp.Δt[end])
    end

    v = hcat(u', g')

    # Scale noise using Δt
    scaled_obs_noise_cov = ekp.obs_noise_cov / ekp.Δt[end]
    noise = rand(MvNormal(zeros(N_obs), scaled_obs_noise_cov), ekp.N_ens)

    # Add obs_mean (N_obs) to each column of noise (N_obs × N_ens) if
    # G is deterministic
    y = deterministic_forward_map ? (ekp.obs_mean .+ noise) : (ekp.obs_mean .+ zero(noise))

    # N_obs × N_obs \ [N_obs × N_ens]
    # --> tmp is [N_obs × N_ens]
    tmp = (cov_gg + scaled_obs_noise_cov) \ (y - g)
    u += (cov_ug * tmp) # [N_par × N_ens]

    # store new parameters (and model outputs)
    push!(ekp.u, DataContainer(u, data_are_columns = true))
    push!(ekp.g, DataContainer(g, data_are_columns = true))

    # Sparse EKI
    cov_vv = vcat(hcat(cov_uu, cov_ug), hcat(cov_ug', cov_gg))
    H_u = hcat(1.0 * I(size(u)[1]), fill(FT(0), size(u)[1], size(g)[1]))
    H_g = hcat(fill(FT(0), size(g)[1], size(u)[1]), 1.0 * I(size(g)[1]))

    H_uc = H_u[ekp.process.uc_idx, :]

    cov_vv_inv = cov_vv \ (1.0 * I(size(cov_vv)[1]))

    # Loop over ensemble members to impose sparsity
    Threads.@threads for j in 1:(ekp.N_ens)
        # Solve a quadratic programming problem
        u[:, j] = sparse_qp(ekp, v[j, :], cov_vv_inv, H_u, H_g, y[:, j], H_uc = H_uc)

        # Threshold the results if needed
        if ekp.process.threshold_eki
            u[ekp.process.uc_idx, j] =
                u[ekp.process.uc_idx, j] .* (abs.(u[ekp.process.uc_idx, j]) .> ekp.process.threshold_value)
        end

        # Add small noise to constrained elements of u
        u[ekp.process.uc_idx, j] +=
            rand(MvNormal(zeros(size(ekp.process.uc_idx)[1]), ekp.process.reg * I(size(ekp.process.uc_idx)[1])))
    end

    # Store error
    compute_error!(ekp)

    # Check convergence
    cov_new = cov(get_u_final(ekp), dims = 2)
    cov_ratio = det(cov_new) / det(cov_init)
    if cov_ratio < cov_threshold
        @warn string(
            "New ensemble covariance determinant is less than ",
            cov_threshold,
            " times its previous value.",
            "\nConsider reducing the EK time step.",
        )
    end
end


text:
using ExoplanetsSysSim
using StatsFuns
using JLD
using DataFrames
#import ExoplanetsSysSim.StellarTable.df
#import ExoplanetsSysSim.StellarTable.usable
import Compat: UTF8String, ASCIIString

## simulation_parameters
function setup_sim_param_model(args::Vector{ASCIIString} = Array{ASCIIString}(0) )   # allow this to take a list of parameter (e.g., from command line)
  sim_param = SimParam()
  #add_param_fixed(sim_param,"max_planets_in_sys",20)
  #add_param_fixed(sim_param,"max_tranets_in_sys",8)
  add_param_fixed(sim_param,"max_planets_in_sys",2)
  add_param_fixed(sim_param,"max_tranets_in_sys",2)
  add_param_fixed(sim_param,"num_targets_sim_pass_one",150969)                      # Note this is used for the number of stars in the simulations, not necessarily related to number of Kepler targets
  add_param_fixed(sim_param,"num_kepler_targets",150969)                            # Note this is used for the number of Kepler targets for the observational catalog
  #add_param_fixed(sim_param,"generate_star",generate_star_dumb)
  add_param_fixed(sim_param,"generate_planetary_system", ExoplanetsSysSim.generate_planetary_system_simple)
  #add_param_fixed(sim_param,"generate_planetary_system", ExoplanetsSysSim.generate_planetary_system_uncorrelated_incl)

  # add_param_fixed(sim_param,"generate_kepler_target",ExoplanetsSysSim.generate_kepler_target_simple)
  add_param_fixed(sim_param,"generate_kepler_target",generate_kepler_target_from_table)
  #=
  add_param_fixed(sim_param,"generate_num_planets",ExoplanetsSysSim.generate_num_planets_poisson)
  add_param_active(sim_param,"log_eta_pl",log(2.0))
  =#
  #=
  add_param_fixed(sim_param,"generate_num_planets",generate_num_planets_poisson_mixture)
  add_param_active(sim_param,"frac_zero_planet",0.3)
  add_param_active(sim_param,"frac_one_planet",0.4)
  add_param_active(sim_param,"frac_multi_planet",0.3)
  add_param_active(sim_param,"eta_pl",4.0)
  =#
  add_param_fixed(sim_param,"generate_num_planets",generate_num_planets_categorical)
  add_param_active(sim_param,"fracs_num_planets",ones(get_int(sim_param,"max_planets_in_sys")+1)/(1+get_int(sim_param,"max_planets_in_sys")))

  add_param_fixed(sim_param,"generate_planet_mass_from_radius",ExoplanetsSysSim.generate_planet_mass_from_radius_powerlaw)
  add_param_fixed(sim_param,"mr_power_index",2.0)
  add_param_fixed(sim_param,"mr_const",1.0)
  #add_param_fixed(sim_param,"generate_period_and_sizes",ExoplanetsSysSim.generate_period_and_sizes_log_normal)
  #add_param_active(sim_param,"mean_log_planet_radius",log(2.0*ExoplanetsSysSim.earth_radius))
  #add_param_active(sim_param,"sigma_log_planet_radius",log(2.0))
  ##add_param_fixed(sim_param,"sigma_log_planet_radius",log(2.0))
  #add_param_active(sim_param,"mean_log_planet_period",log(10.0))
  #add_param_active(sim_param,"sigma_log_planet_period",log(2.0))
  ##add_param_fixed(sim_param,"sigma_log_planet_period",log(2.0))
  #add_param_fixed(sim_param,"generate_period_and_sizes", ExoplanetsSysSim.generate_period_and_sizes_power_law)
  #add_param_active(sim_param,"power_law_P",0.3)
  #add_param_active(sim_param,"power_law_r",-2.44)
  add_param_fixed(sim_param,"min_period",1.0)
  add_param_fixed(sim_param,"max_period",400.0)
  add_param_fixed(sim_param,"min_radius",0.5*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"max_radius",10.0*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"generate_e_omega",ExoplanetsSysSim.generate_e_omega_rayleigh)
  add_param_fixed(sim_param,"sigma_hk",0.03)
  #add_param_fixed(sim_param,"sigma_hk_one",0.3)
  #add_param_fixed(sim_param,"sigma_hk_multi",0.03)
  add_param_fixed(sim_param,"sigma_incl",0.0)   # degrees; 0 = coplanar w/ generate_kepler_target_simple; ignored by generate_planetary_system_uncorrelated_incl
  add_param_fixed(sim_param,"calc_target_obs_sky_ave",ExoplanetsSysSim.calc_target_obs_sky_ave)
  add_param_fixed(sim_param,"calc_target_obs_single_obs",ExoplanetsSysSim.calc_target_obs_single_obs)
  add_param_fixed(sim_param,"read_target_obs",ExoplanetsSysSim.simulated_read_kepler_observations)
  add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_fixed_noise)
  #add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_diagonal)
  # add_param_fixed(sim_param,"rng_seed",1234)   # If you want to be able to reproduce simulations

  #set_inactive(sim_param,["eta_pl","power_law_P","power_law_r"])
  add_param_fixed(sim_param,"star_table_setup",setup_star_table_christiansen)
  add_param_fixed(sim_param,"stellar_catalog","q1_q17_dr25_stellar.jld")
  p_lim_arr_num = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
  r_lim_arr_num = [1., 1.25, 1.5, 1.75, 2.]
  add_param_fixed(sim_param, "p_lim_arr", p_lim_arr_num)
  add_param_fixed(sim_param, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
  #add_param_fixed(sim_param,"generate_num_planets",generate_num_planets_from_rate_table)
  add_param_fixed(sim_param,"generate_period_and_sizes", generate_period_and_sizes_from_rate_table)
  #=
  rate_tab_init = [0.027 0.149 0.354 1.37 2.56 1.96 2.85 4.77 0.0001 ;
                   0.02 0.063 0.432 0.9 1.72 1.46 2.4 5.42 3.63 ;
                   0.0074 0.032 0.389 0.545 1.32 1.65 1.96 1.69 1.17 ;
                   0.0001 0.017 0.156 0.542 0.777 1.35 0.965 1.22 2.49]*0.01

  add_param_active(sim_param,"obs_par", rate_tab_init)
  =#
  return sim_param
end

## planetary_system
using Distributions
function generate_num_planets_categorical(s::Star, sim_param::SimParam)
  const n::Int64 = get_int(sim_param,"max_tranets_in_sys")
  fracs::Array{Float64,1} = get(sim_param,"fracs_num_planets", ones(n+1)/(n+1))
  fracs = abs.(fracs)
  fracs ./= sum(fracs)
  n = rand(Distributions.Categorical(fracs))-1
end

function generate_num_planets_poisson_mixture(s::Star, sim_param::SimParam)
  frac_zero::Float64 = min(max(get_real(sim_param,"frac_zero_planet"), 0.0),1.0)
  frac_one::Float64  = min(max(get_real(sim_param,"frac_one_planet"), 0.0),1.0)
  frac_multi::Float64  = min(max(get_real(sim_param,"frac_multi_planet"),0.0), 1.0)
  frac_sum = frac_zero+frac_one+frac_multi
  frac_zero /= frac_sum
  frac_one /= frac_sum
  frac_multi /= frac_sum
  u = rand()
  if u <= frac_zero
     return 0
  elseif u<= frac_zero+frac_one
     return 1
  else
     const lambda::Float64 = get_real(sim_param,"eta_pl")
     const max_planets_in_sys::Int64 = get_int(sim_param,"max_planets_in_sys")
     return ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_planets_in_sys,min_planets_in_sys=2)
  end
end

function generate_num_planets_from_rate_table(s::Star, sim_param::SimParam)
  const max_planets_in_sys::Int64 = get_int(sim_param,"max_planets_in_sys")
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})
  lambda = sum_kbn(rate_tab)
  #println("# lambda= ", lambda)
  ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_planets_in_sys)
end

function generate_period_and_sizes_from_rate_table(s::Star, sim_param::SimParam; num_pl::Integer = 1)
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})

  limitP::Array{Float64,1} = get_any(sim_param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(sim_param, "r_lim_arr", Array{Float64,1})

  @assert ((length(limitP)-1) == size(rate_tab, 2))
  @assert ((length(limitRp)-1) == size(rate_tab, 1))
  Plist = Array{Float64}(num_pl)
  Rplist = Array{Float64}(num_pl)
  rate_tab_1d = reshape(rate_tab,length(rate_tab))
  #logmaxcuml = logsumexp(rate_tab_1d)
  #cuml = cumsum_kbn(exp(rate_tab_1d-logmaxcuml))
  maxcuml = sum(rate_tab_1d)
  cuml = cumsum_kbn(rate_tab_1d/maxcuml)

  for n in 1:num_pl
    rollp = rand()
    idx = findfirst(x -> x > rollp, cuml)
    i_idx = (idx-1)%size(rate_tab,1)+1
    j_idx = floor(Int64,(idx-1)//size(rate_tab,1))+1
    ### TODO: Keep uniform log sampling here?
    Rp = exp(rand()*(log(limitRp[i_idx+1]/limitRp[i_idx]))+log(limitRp[i_idx]))
    P = exp(rand()*(log(limitP[j_idx+1]/limitP[j_idx]))+log(limitP[j_idx]))
    #push!(Plist, P)
    #push!(Rplist, Rp)
    if !(0<P<Inf)
       println("# BAD P: ",P," should be in ", limitP[j_idx], " - ", limitP[j_idx+1])
    end
    Plist[n] = P
    Rplist[n] = Rp
  end
  return Plist, Rplist
end

## stellar_table
function setup_star_table_christiansen(sim_param::SimParam; force_reread::Bool = false)
  #global df
  df = ExoplanetsSysSim.StellarTable.df
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return df
     #return data
  end
  stellar_catalog_filename = convert(ASCIIString,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(ASCIIString,get(sim_param,"stellar_catalog","q1_q17_dr25_stellar.csv")) ) )
  df = setup_star_table_christiansen(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  set_star_table(df)
  return df
end

function setup_star_table_christiansen(filename::ASCIIString; force_reread::Bool = false)
  #global df, usable
  df = ExoplanetsSysSim.StellarTable.df
  usable = ExoplanetsSysSim.StellarTable.usable
  if ismatch(r".jld$",filename)
  try
    data = load(filename)
    df::DataFrame = data["stellar_catalog"]
    usable::Array{Int64,1} = data["stellar_catalog_usable"]
    set_star_table(df, usable)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld format."))
  end
  else
  try
    df = readtable(filename)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end

  has_mass = ! (isna(df[:mass]) | isna(df[:mass_err1]) | isna(df[:mass_err2]))
  has_radius = ! (isna(df[:radius]) | isna(df[:radius_err1]) | isna(df[:radius_err2]))
  has_dens = ! (isna(df[:dens]) | isna(df[:dens_err1]) | isna(df[:dens_err2]))
  has_rest = ! (isna(df[:rrmscdpp04p5]) | isna(df[:dataspan]) | isna(df[:dutycycle]))
  in_Q1Q12 = []
  for x in df[:st_quarters]
    subx = string(x)
    subx = ("0"^(17-length(subx)))*subx
    indQ = search(subx, '1')
    if ((indQ < 1) | (indQ > 12))
      push!(in_Q1Q12, false)
    else
      push!(in_Q1Q12, true)
    end
  end
  is_FGK = []
  for x in 1:length(df[:teff])
    if ((df[x,:teff] > 4000.0) & (df[x,:teff] < 7000.0) & (df[x,:logg] > 4.0))
      push!(is_FGK, true)
    else
      push!(is_FGK, false)
    end
  end
  is_usable = has_radius & is_FGK & has_mass & has_rest #& in_Q1Q12 # & has_dens
  #if contains(filename,"q1_q12_christiansen.jld")
  if contains(filename,"q1_q12_christiansen")   # TODO: Ask Danely what he's trying to do here.
    is_usable = is_usable & in_Q1Q12
  end
  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  # TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp04p5, :dataspan, :dutycycle ]
  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  set_star_table(df, usable)
  end
  return df
end

## summary_statistics

# Compile indices for N-tranet systems
function calc_summary_stats_idx_n_tranets!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  if haskey(css.cache,"idx_n_tranets")
     return css.cache["idx_n_tranets"]
  end
     max_tranets_in_sys = get_int(param,"max_tranets_in_sys")
     idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
     for n in 1:max_tranets_in_sys-1
       idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target )
     end
     idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target )
     css.cache["idx_n_tranets"] = idx_n_tranets
  return idx_n_tranets
end

# Count total number of tranets using lists of indices for N-tranet systems
function calc_summary_stats_num_tranets!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  if haskey(css.stat,"num_tranets")
     return css.stat["num_tranets"]
  elseif haskey(css.cache,"num_tranets")
     css.stat["num_tranets"] = css.cache["num_tranets"]
     return css.stat["num_tranets"]
  end
     idx_n_tranets = calc_summary_stats_idx_n_tranets!(css,cat_obs,param)
     num_tranets = 0
     for n in 1:length(idx_n_tranets)
         num_tranets += n*length(idx_n_tranets[n])
     end
     css.stat["num_tranets"] = num_tranets
  return num_tranets
end

function calc_summary_stats_num_targets!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam ; trueobs_cat::Bool = false)
  if !trueobs_cat
    css.stat["num targets"] = get_int(param,"num_targets_sim_pass_one")
  else
    css.stat["num targets"] = get_int(param,"num_kepler_targets")
  end
end

function calc_summary_stats_num_n_tranet_systems!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css,cat_obs,param)
  #max_tranets_in_sys = get_int(param,"max_tranets_in_sys")
  num_n_tranet_systems = map(n->length(idx_n_tranets[n]), 1:length(idx_n_tranets) )
  #for n in 1:length(idx_n_tranets)
  #  num_n_tranet_systems[n] = length(idx_n_tranets[n])
  #end
  css.stat["num n-tranet systems"] = num_n_tranet_systems
end

function calc_summary_stats_duration_ratios_neighbors!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  if haskey(css.stat,"duration_ratio_list")
     return css.stat["duration_ratio_list"]
  elseif haskey(css.cache,"duration_ratio_list")
     return css.cache["duration_ratio_list"]
  end
  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css,cat_obs,param)
  @assert length(idx_n_tranets) >= 1

  # Calculate how many duration ratios there will be & allocate storage
  num_ratios = 0
  for i in 2:length(idx_n_tranets)
     num_ratios += length(idx_n_tranets[i])*(i-1)
  end
  duration_ratio_list = Array{Float64}(num_ratios)

  k = 0
  for n in 2:length(idx_n_tranets)         # Loop over number of tranets in system
    for i in idx_n_tranets[n]              # Loop over systems with n tranets
       period_in_sys = Array{Float64}(n)
       duration_in_sys = Array{Float64}(n)
       for j in 1:n                        # Loop over periods within a system
         period_in_sys[j] = cat_obs.target[i].obs[j].period
         duration_in_sys[j] = cat_obs.target[i].obs[j].duration
       end
       perm = sortperm(period_in_sys)
       for j in 1:(n-1)                       # Loop over period ratios within a system
          inv_period_ratio = period_in_sys[perm[j+1]]/period_in_sys[perm[j]]
          if 1<inv_period_ratio<Inf
             k += 1
             duration_ratio_list[k] = duration_in_sys[perm[j]]/duration_in_sys[perm[j+1]] * inv_period_ratio^(1//3)
          end
       end
    end
  end
  resize!(duration_ratio_list,k)
  css.stat["duration_ratio_list"] = duration_ratio_list

  return duration_ratio_list
end

function calc_summary_stats_period_radius_ratios_neighbors_internal!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  if haskey(css.stat,"period_ratio_list") && haskey(css.stat,"radius_ratio_list")
     return (css.stat["period_ratio_list"], css.stat["radius_ratio_list"])
  elseif haskey(css.cache,"period_ratio_list") && haskey(css.cache,"radius_ratio_list")
     return (css.cache["period_ratio_list"], css.cache["radius_ratio_list"])
  end
  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css,cat_obs,param)
  @assert length(idx_n_tranets) >= 1

  # Calculate how many period ratios there will be & allocate storage
  num_ratios = 0
  for i in 2:length(idx_n_tranets)
     num_ratios += length(idx_n_tranets[i])*(i-1)
  end
  period_ratio_list = Array{Float64}(num_ratios)
  radius_ratio_list = Array{Float64}(num_ratios)

  k = 0
  for n in 2:length(idx_n_tranets)         # Loop over number of tranets in system
    period_in_sys = Array{Float64}(n)
    #radius_in_sys = Array{Float64}(n)
    depth_in_sys = Array{Float64}(n)
    for i in idx_n_tranets[n]              # Loop over systems with n tranets
       for j in 1:n                        # Loop over periods within a system
         period_in_sys[j] = cat_obs.target[i].obs[j].period
         #radius_in_sys[j] = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius
         depth_in_sys[j] = cat_obs.target[i].obs[j].depth
       end
       perm = sortperm(period_in_sys)
       for j in 1:(n-1)                       # Loop over period ratios within a system
          k = k+1
          period_ratio_list[k] = period_in_sys[perm[j]]/period_in_sys[perm[j+1]]
          #radius_ratio_list[k] = radius_in_sys[perm[j]]/radius_in_sys[perm[j+1]]
          radius_ratio_list[k] = sqrt(depth_in_sys[perm[j]]/depth_in_sys[perm[j+1]])
       end
    end
  end
  css.cache["period_ratio_list"] = period_ratio_list
  css.cache["radius_ratio_list"] = radius_ratio_list

  return (period_ratio_list,radius_ratio_list)
end

function calc_summary_stats_period_radius_ratios_neighbors!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  (period_ratio_list,radius_ratio_list) = calc_summary_stats_period_radius_ratios_neighbors!(css,cat_obs,param)
  css.stat["period_ratio_list"] = period_ratio_list
  css.stat["radius_ratio_list"] = radius_ratio_list
  return (period_ratio_list, radius_ratio_list)
end
function calc_summary_stats_period_ratios_neighbors!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  period_ratio_list = calc_summary_stats_period_radius_ratios_neighbors_internal!(css,cat_obs,param)[1]
  css.stat["period_ratio_list"] = period_ratio_list
  return period_ratio_list
end
function calc_summary_stats_radius_ratios_neighbors!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  radius_ratio_list = calc_summary_stats_period_radius_ratios_neighbors_internal!(css,cat_obs,param)[2]
  css.stat["radius_ratio_list"] = radius_ratio_list
  return radius_ratio_list
end

function calc_summary_stats_mean_std_log_period_depth!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  # Allocate arrays to store values for each tranet
  num_tranets  = calc_summary_stats_num_tranets!(css, cat_obs, param)
  period_list = zeros(num_tranets)
  depth_list = zeros(num_tranets)
  #weight_list = ones(num_tranets)

  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css, cat_obs, param)
  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")
  @assert max_tranets_in_sys >= 1
   i = 1   # tranet id
   for targ in cat_obs.target                        # For each target
     for j in 1:min(length(targ.obs),max_tranets_in_sys)          # For each tranet around that target (but truncated if too many tranets in one system)
         #println("# i= ",i," j= ",j)
         period_list[i] = targ.obs[j].period
         depth_list[i] = targ.obs[j].depth
         #weight_list[i] = 1.0
         i = i+1
      end
   end

  css.cache["P list"] = period_list                                     # We can store whole lists, e.g., if we want to compute K-S distances
  css.cache["depth list"] = depth_list
  #css.cache["weight list"] = weight_list

  idx_good = Bool[ period_list[i]>0.0 && depth_list[i]>0.0 for i in 1:length(period_list) ]
  log_period_list = log10(period_list[idx_good])
  log_depth_list = log10(depth_list[idx_good])
  css.stat["mean log10 P"]  =  mean_log_P = mean(log_period_list)
  css.stat["mean log10 depth"]  =  mean_log_depth = mean(log_depth_list)
  css.stat["std log10 P"]  = std_log_P = stdm(log_period_list,mean_log_P)
  css.stat["std log10 depth"]  = std_log_depth = stdm(log_depth_list,mean_log_depth)

  return (mean_log_P, std_log_P, mean_log_depth, std_log_depth)
end

function calc_summary_stats_cuml_period_depth!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam)
  # Allocate arrays to store values for each tranet
  num_tranets  = calc_summary_stats_num_tranets!(css, cat_obs, param)
  period_list = zeros(num_tranets)
  depth_list = zeros(num_tranets)
  #weight_list = ones(num_tranets)

  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css, cat_obs, param)
  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")
  @assert max_tranets_in_sys >= 1
  i = 0   # tranet id
   for targ in cat_obs.target                        # For each target
     for j in 1:min(length(targ.obs),max_tranets_in_sys)          # For each tranet around that target (but truncated if too many tranets in one system)
         i = i+1
         #println("# i= ",i," j= ",j)
         period_list[i] = targ.obs[j].period
         depth_list[i] = targ.obs[j].depth
         #weight_list[i] = 1.0
      end
   end
  resize!(period_list,i)
  resize!(depth_list,i)
  css.stat["P list"] = period_list                                     # We can store whole lists, e.g., if we want to compute K-S distances
  css.stat["depth list"] = depth_list
  #css.cache["weight list"] = weight_list
  #println("# P list = ",period_list)
  return (period_list,depth_list)
end

function calc_summary_stats_obs_binned_rates!(css::CatalogSummaryStatistics, cat_obs::KeplerObsCatalog, param::SimParam; trueobs_cat::Bool = false)
  num_tranets  = calc_summary_stats_num_tranets!(css, cat_obs, param)
  idx_n_tranets = calc_summary_stats_idx_n_tranets!(css, cat_obs, param)

  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet
  css.cache["idx_tranets"] = idx_tranets                                   # We can save lists of indices to summary stats for pass 2, even though we won't use these for computing a distance or probability

  #=
  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target )
  css.cache["idx_n_tranets"] = idx_n_tranets

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  if ( length( find(x::KeplerTargetObs-> length(x.obs) > max_tranets_in_sys, cat_obs.target ) ) > 0)   # Make sure max_tranets_in_sys is at least big enough for observed systems
    warn("Observational data has more transiting planets in one systems than max_tranets_in_sys allows.")
  end
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations
  css.cache["num_tranets"] = num_tranets

  num_sys_tranets = zeros(max_tranets_in_sys)                           # Since observed data, don't need to calculate probabilities.
  for n in 1:max_tranets_in_sys                                         # Make histogram of N-tranet systems
    num_sys_tranets[n] = length(idx_n_tranets[n])
  end
  css.stat["num_sys_tranets"] = num_sys_tranets
  css.stat["planets detected"] = num_tranets
  =#

  period_list = zeros(num_tranets)
  radius_list = zeros(num_tranets)
  weight_list = ones(num_tranets)

  n = 1    # tranet id
  for i in 1:length(cat_obs.target)
    for j in 1:num_planets(cat_obs.target[i])
      period_list[n] = cat_obs.target[i].obs[j].period
      radius_list[n] = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius
      n = n+1
    end
  end

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})
  @assert length(limitP)>=2 && length(limitRp)>=2

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #css.stat["planets detected"] = sum(np_bin)
  css.stat["planets table"] = np_bin

  return css
end

function calc_summary_stats_model(cat_obs::KeplerObsCatalog, param::SimParam; trueobs_cat::Bool = false)
  css = CatalogSummaryStatistics()
  calc_summary_stats_num_targets!(css,cat_obs,param,trueobs_cat=trueobs_cat)
  calc_summary_stats_num_tranets!(css,cat_obs,param)
  calc_summary_stats_num_n_tranet_systems!(css,cat_obs,param)
  calc_summary_stats_cuml_period_depth!(css,cat_obs,param)
  #calc_summary_stats_obs_binned_rates!(css,cat_obs,param)
  #calc_summary_stats_mean_std_log_period_depth!(css,cat_obs,param)
  calc_summary_stats_period_ratios_neighbors!(css,cat_obs,param)
  calc_summary_stats_duration_ratios_neighbors!(css,cat_obs,param)
  return css
end

## abc_distance
function calc_distance_model_ks(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  d1 = calc_distance_num_planets(summary1,summary2,sim_param,verbose=verbose)
  d2 = calc_distance_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d3 = calc_distance_num_planets_binned(summary1,summary2,sim_param,verbose=verbose)
  #d4 = calc_distance_mean_std_log_period_depth(summary1,summary2,sim_param,verbose=verbose)
  d5 = calc_distance_ks_period(summary1,summary2,sim_param,verbose=verbose)
  d6 = calc_distance_ks_depth(summary1,summary2,sim_param,verbose=verbose)
  d7 = calc_distance_ks_period_ratios(summary1,summary2,sim_param,verbose=verbose)
  d8 = calc_distance_ks_duration_ratios(summary1,summary2,sim_param,verbose=verbose)
  return vcat(d1, d2, d5, d6, d7, d8)
end

function calc_distance_model_kl(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  d1 = calc_distance_kl_num_planets(summary1,summary2,sim_param,verbose=verbose)
  d2 = calc_distance_kl_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d2 = calc_distance_hellinger_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d3 = calc_distance_num_planets_binned(summary1,summary2,sim_param,verbose=verbose)
  #d4 = calc_distance_mean_std_log_period_depth(summary1,summary2,sim_param,verbose=verbose)
  d5 = calc_distance_kl_period(summary1,summary2,sim_param,verbose=verbose)
  d6 = calc_distance_kl_depth(summary1,summary2,sim_param,verbose=verbose)
  d7 = calc_distance_kl_period_ratios(summary1,summary2,sim_param,verbose=verbose)
  d8 = calc_distance_kl_duration_ratios(summary1,summary2,sim_param,verbose=verbose)
  return vcat(d1, d2, d5, d6, d7, d8)
end
calc_distance_model(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false) = calc_distance_model_kl(summary1,summary2,sim_param,verbose=verbose)

function calc_distance_num_planets(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = haskey(summary1.stat,"num_tranets") ? summary1.stat["num_tranets"] : summary1.stat["expected planets detected"]
    np2 = haskey(summary2.stat,"num_tranets") ? summary2.stat["num_tranets"] : summary2.stat["expected planets detected"]
    #println("np1 = ",np1,", np2 = ",np2)

    dist_np = dist_L1_abs(np1/summary1.stat["num targets"], np2/summary2.stat["num targets"])
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", d = ",dist_np)
  return dist_np
end

function calc_distance_num_planets_binned(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = haskey(summary1.stat,"planets table") ? summary1.stat["planets table"] : summary1.stat["expected planets table"]
    np2 = haskey(summary2.stat,"planets table") ? summary2.stat["planets table"] : summary2.stat["expected planets table"]
    #println("np1 = ",np1,", np2 = ",np2)

    dist_np_bin = zeros(length(np1))
    for n in 1:length(np1)
      dist_np_bin[n] = dist_L1_abs(np1[n]/summary1.stat["num targets"], np2[n]/summary2.stat["num targets"])
      #println("True # [Bin ", n,"] = ",np1[n],", Expected # [Bin ", n,"] = ",np2[n])
    end
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", dist = ",dist_np_bin)
  return dist_np_bin
end

function calc_distance_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    d[n] = dist_L1_abs(summary1.stat["num n-tranet systems"][n]/summary1.stat["num targets"], summary2.stat["num n-tranet systems"][n]/summary2.stat["num targets"])
  end
  return d
end

function calc_distance_mean_std_log_period_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    d1 = dist_L1_abs(summary1.stat["mean log10 P"], summary2.stat["mean log10 P"])
    d2 = dist_L1_abs(summary1.stat["mean log10 depth"], summary2.stat["mean log10 depth"])
    d3 = dist_L1_abs(summary1.stat["std log10 P"], summary2.stat["std log10 P"])
    d4 = dist_L1_abs(summary1.stat["std log10 depth"], summary2.stat["std log10 depth"])
    return [d1, d2, d3, d4]
end

# compute supremum of differences between empirical cdfs.
# Borrowed from JuliaStats/HypothesisTests.jl
function ksstats{T<:Real, S<:Real}(x::AbstractVector{T}, y::AbstractVector{S})
    n_x, n_y = length(x), length(y)
    sort_idx = sortperm([x; y])
    pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)
end

function calc_distance_ks_period(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["P list"]
  samp2 = summary2.stat["P list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["depth list"]
  samp2 = summary2.stat["depth list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_period_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["period_ratio_list"]
  samp2 = summary2.stat["period_ratio_list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_duration_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["duration_ratio_list"]
  samp2 = summary2.stat["duration_ratio_list"]
  return ksstats(samp1,samp2)[5]
end

# Function for Relative Entropy / K-L divergence
include("kde.jl")

function calc_distance_kl_num_planets(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = (haskey(summary1.stat,"num_tranets") ? summary1.stat["num_tranets"] : summary1.stat["expected planets detected"])
    np2 = (haskey(summary2.stat,"num_tranets") ? summary2.stat["num_tranets"] : summary2.stat["expected planets detected"])
    ntarg1 = summary1.stat["num targets"]
    ntarg2 = summary2.stat["num targets"]
    mintarg = min(ntarg1,ntarg2)
    alpha1 = 1+np1
    beta1 = 1+ntarg1
    #mu1 = alpha1/beta1
    #sigma1 = mu1 /sqrt(alpha1)
    #dist = alpha+lgamma(alpha)+(1-alpha)*digamma(alpha) -log(beta)
    alpha2 = 1+np2
    beta2 = 1+ntarg2
    #dist -= alpha+lgamma(alpha)+(1-alpha)*digamma(alpha) -log(beta)
    #mu2 = alpha2/beta2
    #sigma2 = mu2 /sqrt(alpha1)
    #dist = 0.5*( (sigma2/sigma1)^2 - 1 + 2*log(sigma1/sigma2) + (mu2-mu1)^2/sigma1^2 )  # Normal approximation, substitute denom for difference in rates
    #= Gamma PDF
    dist = abs( (alpha2-alpha1)*digamma(alpha1) - lgamma(alpha1) + lgamma(alpha2)  ) # Gamma approximation, same beta, why need abs?
    dist += abs( (alpha1-alpha2)*digamma(alpha2) - lgamma(alpha2) + lgamma(alpha1)  )  # Gamma approximation, same beta, why need abs?
    print("# a1=",alpha1, " a2=",alpha2, " digamma=",digamma(alpha1), " ", digamma(alpha2)," gamma term=",lgamma(alpha1)-lgamma(alpha2))
    if ntarg1!=ntarg2
       print("# ntarg=",ntarg1," ",ntarg2," dist(before betas)=", dist)
       dist += alpha2*(log(beta1)-log(beta2)) + alpha1*(beta2-beta1)/beta1  # Gammma approximation, beta terms
       dist += alpha1*(log(beta2)-log(beta1)) + alpha2*(beta1-beta2)/beta2  # Gammma approximation, beta terms
       println(" dist (after betas)=",dist)
    end
    =#
    ## #= Poisson
    #dist = abs( alpha2-alpha1+alpha1*log(alpha1/alpha2) )
    #dist += abs( alpha2-alpha1+alpha1*log(alpha1/alpha2) )
    rate1 = alpha1/beta1*mintarg
    rate2 = alpha2/beta2*mintarg
    dist = (rate1-rate2)*log(rate1/rate2)
    ## =#
    dist /= 2
    return dist
end

function calc_distance_kl_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  #= categorical distribution
  f1sum = sum(summary1.stat["num n-tranet systems"]) # summary1.stat["num targets"]
  f2sum = sum(summary2.stat["num n-tranet systems"]) # summary2.stat["num targets"]
  if !(f1sum>0 && f2sum>0) return 0.0  end
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    f1 = summary1.stat["num n-tranet systems"][n]/f1sum
    f2 = summary2.stat["num n-tranet systems"][n]/f2sum
    m = (f1+f2)/2
    if m>zero(m)
    if f1>zero(f1)
       d[n] += 0.5*f1*log(f1/m)
    end
    if f2>zero(f2)
       d[n] += 0.5*f2*log(f2/m)
    end
#   else
#       d += Inf
    end
  end
  =#
  # Poisson distributions for each
    ntarg1 = summary1.stat["num targets"]
    ntarg2 = summary2.stat["num targets"]
    #mintarg = min(ntarg1,ntarg2)
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    f1 = (1+summary1.stat["num n-tranet systems"][n])/(1+ntarg1) # *(1+mintarg)
    f2 = (1+summary2.stat["num n-tranet systems"][n])/(1+ntarg2) # *(1+mintarg)
    d[n] = (f1-f2)*log(f1/f2)
  end
  return d
end
function calc_distance_hellinger_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  f1sum = sum(summary1.stat["num n-tranet systems"]) # summary1.stat["num targets"]
  f2sum = sum(summary2.stat["num n-tranet systems"]) # summary2.stat["num targets"]
  d = 1
  for n in 1:max_tranets_in_sys
    f1 = summary1.stat["num n-tranet systems"][n]/f1sum
    f2 = summary2.stat["num n-tranet systems"][n]/f2sum
    d -= sqrt(f1*f2)
  end
  #d = sqrt(d)
  return d
end

function calc_distance_kl_period(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  #println("# P list 1: n=",length(summary1.stat["P list"])," min=",minimum(summary1.stat["P list"]), " max=",maximum(summary1.stat["P list"]))
  #println("# P list 2: n=",length(summary2.stat["P list"])," min=",minimum(summary2.stat["P list"]), " max=",maximum(summary2.stat["P list"]))
  samp1 = log(summary1.stat["P list"] )
  samp2 = log(summary2.stat["P list"] )
  calc_kl_distance_ab(samp1,samp2,log(0.5),log(320.) )
end

function calc_distance_kl_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = log(summary1.stat["depth list"] )
  samp2 = log(summary2.stat["depth list"] )
  calc_kl_distance_ab(samp1,samp2,log(0.000025),log(0.025) )
end

function calc_distance_kl_period_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  min_ratios_to_compute_distances = 3
  distance_when_not_enough_ratios = 10.0
  samp1 = summary1.stat["period_ratio_list"]
  samp2 = summary2.stat["period_ratio_list"]
  if length(samp1)<min_ratios_to_compute_distances || length(samp2)<min_ratios_to_compute_distances
     return distance_when_not_enough_ratios
  end
  calc_kl_distance_ab(samp1,samp2,0.0,1.0)
end

function calc_distance_kl_duration_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  min_ratios_to_compute_distances = 3
  distance_when_not_enough_ratios = 10.0
  samp1 = log( summary1.stat["duration_ratio_list"] )
  samp2 = log( summary2.stat["duration_ratio_list"] )
  if length(samp1)<min_ratios_to_compute_distances || length(samp2)<min_ratios_to_compute_distances
     return distance_when_not_enough_ratios
  end
  calc_kl_distance_ab(samp1,samp2,-log(10.0),log(10.0) )
end

## eval_model
function test_model()
  global sim_param_closure = setup_sim_param_model()
  cat_phys = generate_kepler_physical_catalog(sim_param_closure)
  cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param_closure)
  global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs,sim_param_closure)
  global cat_phys_try_closure = generate_kepler_physical_catalog(sim_param_closure)
  global cat_obs_try_closure  = observe_kepler_targets_sky_avg(cat_phys_try_closure,sim_param_closure)
  global summary_stat_try_closure  = calc_summary_stats_sim_pass_one_demo(cat_obs_try_closure,cat_phys_try_closure,sim_param_closure)
  summary_stat_try_closure   = calc_summary_stats_sim_pass_two_demo(cat_obs_try_closure,cat_phys_try_closure,summary_stat_try_closure,sim_param_closure)
  param_guess = make_vector_of_sim_param(sim_xparam_closure)
  evaluate_model_scalar_ret( param_guess)
end


text:
using LaserTypes, Test
using StaticArrays
using Unitful

@testset "Poynting vector" begin
    laser = setup_laser(GaussLaser, :SI_unitful, profile=ConstantProfile())
    r = SVector{3}(0,0,0)u"m"
    t = 0u"s"
    s = S(r, t, laser)

    @testset "Units" begin
        @test dimension(eltype(s)) == dimension(u"W/m^2")
    end
    @testset "Generic" begin
        @test isreal(s)
        @test ustrip(s[1]) isa Real
    end
end


text:
using PointsOnASphere
using Test

@testset "Constructor" begin
	r,θ,ϕ = 2,pi/2,pi/3

	@testset "Point2D" begin
		n = Point2D(θ,ϕ)
	    @test n.θ == θ
	    @test n.ϕ == ϕ

	    @test Point2D(n) === n

	    n2 = Point2D((θ,ϕ))
	    @test n2 == n

	    n3 = Point2D(Point3D(1,θ,ϕ))
	    @test n3 == n

	    @test_throws ArgumentError Point2D(θ,ϕ + 2pi)
	    @test_throws ArgumentError Point2D(θ,ϕ - 2pi)
	    @test_throws ArgumentError Point2D(θ+π,ϕ)
	    @test_throws ArgumentError Point2D(-θ,ϕ)
	end

	@testset "Point3D" begin
		x = Point3D(r,θ,ϕ)
		@test x.r == r
		@test x.θ == θ
	    @test x.ϕ == ϕ

	    @test Point3D(x) === x

	    x2 = Point3D(r,(θ,ϕ))
	    @test x2 == x

	    x3 = Point3D(r,Point2D(θ,ϕ))
	    @test x3 == x

	    x4 = Point3D((θ,ϕ))
	    @test x4.r == 1
	    @test x4.θ == θ
	    @test x4.ϕ == ϕ

	    x5 = Point3D(Point2D(θ,ϕ))
	    @test x4 == x5

	    x6 = Point3D((r,θ,ϕ))
	    @test x6 == x

	    @test_throws ArgumentError Point3D(r,θ,ϕ + 2pi)
	    @test_throws ArgumentError Point3D(r,θ,ϕ - 2pi)
	    @test_throws ArgumentError Point3D(r,θ+π,ϕ)
	    @test_throws ArgumentError Point3D(r,-θ,ϕ)
	    @test_throws ArgumentError Point3D(-r,θ,ϕ)
	end
end

@testset "show" begin
	r,θ,ϕ = 2,pi/2,pi/3

	io = IOBuffer()
	pt = Point2D(θ,ϕ)
	show(io,pt)
	strexp = "(θ=$(round(pt.θ,sigdigits=3)), ϕ=$(round(pt.ϕ,sigdigits=3)))"
	@test String(take!(io)) == strexp

	show(io,MIME"text/plain"(),pt)
	strexp = "(θ=$(pt.θ), ϕ=$(pt.ϕ))"
	@test String(take!(io)) == strexp

	pt = Point3D(1,θ,ϕ)
	show(io,pt)
	strexp = "(r=$(round(pt.r,sigdigits=3)), θ=$(round(pt.θ,sigdigits=3)), ϕ=$(round(pt.ϕ,sigdigits=3)))"
	@test String(take!(io)) == strexp

	show(io,MIME"text/plain"(),pt)
	strexp = "(r=$(pt.r), θ=$(pt.θ), ϕ=$(pt.ϕ))"
	@test String(take!(io)) == strexp
end


text:
using KhepriPOVRay
using Test

@testset "KhepriPOVRay.jl" begin
    # Write your tests here.
end


text:
__precompile__()

module Bio

include("IO.jl")
include("Ragel.jl")
include("StringFields.jl")
include("util/Util.jl")
include("seq/Seq.jl")
include("intervals/Intervals.jl")
include("align/Align.jl")
include("util/tokenize.jl")
include("util/indexing.jl")
include("var/Var.jl")
include("phylo/Phylo.jl")
include("structure/Structure.jl")
include("tools/Tools.jl")
include("util/windows.jl")
include("precompile.jl")

end  # module Bio


text:
"""
    SomeVertexInCircumsphere(simplex1, r2, c2)

Decides whether some vertices of simplex1 lies inside the circumsphere of simplex2
(defined by the radius and centroid of simplex1, r2 and c2). Each column of simplex1 is a vertex.

Arguments
---------
c2::Vector{Float64} Column vector. Centroid of simplex2.
r2::Float64 Radius of simplex 2

Returns either 0 or 1; 1 if some of the vertices of simplex1 are contained in the circumsphere
of simplex 2, 0 otherwise.
"""
function SomeVertexInCircumsphere(simplex1, r2, c2)
    # The dimension
    n = size(simplex1, 1)

    i = 1
    some_vertex_in_circumsphere = false

    while i <= (n + 1) && !some_vertex_in_circumsphere
        # Difference between the i-th vertex of simplex1 and centroid of simplex2
        ith_vertex = simplex1[:, i] - c2
        tmp = heaviside0(r2^2 - transpose(ith_vertex) * ith_vertex) # Radius times norm

        if tmp == 1
            some_vertex_in_circumsphere = true
        end

        i = i + 1
    end

    return some_vertex_in_circumsphere
end