text:
module AccuracyBenchmark

using DataFrames
using Distributions
import FITSIO
using StaticArrays
import WCS
import CSV

import ..Config
import ..DeterministicVI
import ..Model
import ..ParallelRun
using ..SDSSIO
using ..Coordinates: angular_separation, match_coordinates

const ARCSEC_PER_DEGREE = 3600
const SDSS_ARCSEC_PER_PIXEL = 0.396
const SDSS_DATA_DIR = joinpath(Pkg.dir("Celeste"), "test", "data")
const STRIPE82_RCF = SDSSIO.RunCamcolField(4263, 5, 119)
const COADD_CATALOG_FITS = joinpath(Pkg.dir("Celeste"), "test", "data", "coadd_for_4263_5_119.fit")
const GALAXY_ONLY_COLUMNS = [:gal_frac_dev, :gal_axis_ratio, :gal_radius_px, :gal_angle_deg]

struct BenchmarkFitsFileNotFound <: Exception
    filename::String
end

struct MatchException <: Exception
    msg::String
end

struct MissingColumnsError <: Exception
    missing_columns::Vector{Symbol}
end

################################################################################
# Load/store a catalog DF from/to a CSV.
################################################################################

CATALOG_COLUMNS = Set([
    :ra,
    :dec,
    :is_star,
    :flux_r_nmgy,
    :color_ug,
    :color_gr,
    :color_ri,
    :color_iz,
    :gal_frac_dev,
    :gal_axis_ratio,
    :gal_radius_px,
    :gal_angle_deg,
])

STDERR_COLUMNS = Set([
    :log_flux_r_stderr,
    :color_ug_stderr,
    :color_gr_stderr,
    :color_ri_stderr,
    :color_iz_stderr,
])

BAD_COADD_OBJID = Set([
    # this object is actually multiple objects,
    # see http://legacysurvey.org/viewer/jpeg-cutout/?ra=0.5636&dec=0.4445&zoom=16&layer=decals-dr3
    8647474692482203853,

    # bright neighbor not accounted for,
    # see http://skyserver.sdss.org/dr10/en/tools/chart/navi.aspx?ra=0.551730166101009&dec=0.48085411481715&scale=0.2
    8647474692482203816,

    #  possibly a quasar; update: decals labels this light source a galaxy!
    # see http://legacysurvey.org/viewer?ra=0.5531&dec=0.4530&zoom=16&layer=decals-dr3&sources
    8647474692482204612,

    # Qusar: http://skyserver.sdss.org/dr10/en/get/SpecById.ashx?id=435863481728657408
    8647474692482204147,
])

function assert_columns_are_present(catalog_df::DataFrame, required_columns::Set{Symbol})
    missing_columns = setdiff(required_columns, Set(names(catalog_df)))
    if !isempty(missing_columns)
        throw(MissingColumnsError([missing_columns...]))
    end
end

function read_catalog(csv_file::String)
    catalog_df = CSV.read(csv_file, rows_for_type_detect=100)
    assert_columns_are_present(catalog_df, CATALOG_COLUMNS)
    catalog_df
end

function write_catalog(filename::String, catalog_df::DataFrame;
                       append_hash=false)
    assert_columns_are_present(catalog_df, CATALOG_COLUMNS)

    if append_hash
        # Serialize the data frame into an array of bytes.
        # (CSV.write(::IO, ...) currently broken, so we use a temp file.
        tmp = tempname()
        CSV.write(tmp, catalog_df)
        data = open(tmp) do f
            read(f)
        end
        rm(tmp)

        # Hash the bytes and add the hash string to the filename.
        hash_string = hex(hash(data))[1:10]
        base, extension = splitext(filename)
        filename = @sprintf("%s_%s%s", base, hash_string, extension)

        # Write out the file
        open(filename, "w") do f
            write(f, data)
        end
    else
        CSV.write(filename, catalog_df)
    end

    return filename
end

################################################################################
# Read various catalogs to a common catalog DF format
################################################################################

"""
Convert between SDSS asinh mags and SDSS flux (nMgy).
See http://www.sdss.org/dr12/algorithms/magnitudes/#asinh.
"""
const ASINH_SOFTENING_PARAMETERS = [
    1.4e-10, # for band 1 = u
    0.9e-10,
    1.2e-10,
    1.8e-10,
    7.4e-10, # for band 5 = z
]

function mag_to_flux(mags::AbstractFloat, band_index::Int)
    b = ASINH_SOFTENING_PARAMETERS[band_index]
     1e9 * 2 * b * sinh(-log(10) / 2.5 * mags - log(b))
end

function flux_to_mag(flux_nmgy::AbstractFloat, band_index::Int)
    b = ASINH_SOFTENING_PARAMETERS[band_index]
    -2.5 / log(10) * (asinh(flux_nmgy * 1e-9 / (2*b)) + log(b))
end

function color_from_fluxes(flux1::AbstractFloat, flux2::AbstractFloat)
    if flux1 <= 0 || flux2 <= 0
        missing
    else
        log(flux2 / flux1)
    end
end

function color_from_mags(mags1::AbstractFloat, band1::Int, mags2::AbstractFloat, band2::Int)
    color_from_fluxes(mag_to_flux(mags1, band1), mag_to_flux(mags2, band2))
end

canonical_angle(gal_angle_deg) = gal_angle_deg - floor(gal_angle_deg / 180) * 180

STRIPE_82_CATALOG_KEYS = [
    :objid, :rerun, :run, :camcol, :field, :flags,
    :ra, :dec, :probpsf,
    :psfmag_u, :psfmag_g, :psfmag_r, :psfmag_i, :psfmag_z,
    :devmag_u, :devmag_g, :devmag_r, :devmag_i, :devmag_z,
    :expmag_u, :expmag_g, :expmag_r, :expmag_i, :expmag_z,
    :fracdev_r,
    :devab_r, :expab_r,
    :devphi_r, :expphi_r,
    :devrad_r, :exprad_r,
    :flags, :is_saturated,
]

function load_fits_table_as_data_frame(filename, extension_index, field_names)
    result_df = DataFrame()
    fits = FITSIO.FITS(filename)
    try
        for key in field_names
            result_df[key] = read(fits[extension_index], string(key))
        end
        return result_df
    finally
        close(fits)
    end
end

function load_stripe82_fits_catalog_as_data_frame(filename, extension_index)
    load_fits_table_as_data_frame(filename, extension_index, STRIPE_82_CATALOG_KEYS)
end

"""
Load Stripe 82 objects into a DataFrame. `fits_filename` should be a FITS file created by running a
CasJobs (skyserver.sdss.org/casjobs/) query on the Stripe82 database. See
https://github.com/jeff-regier/Celeste.jl/wiki/About-SDSS-and-Stripe-82.
"""
function load_coadd_catalog(fits_filename)
    raw_df = load_stripe82_fits_catalog_as_data_frame(fits_filename, 2)

    usedev = raw_df[:fracdev_r] .> 0.5  # true=> use dev, false=> use exp
    dev_or_exp(dev_column, exp_column) = ifelse.(usedev, raw_df[dev_column], raw_df[exp_column])
    is_star = [x != 0 for x in raw_df[:probpsf]]
    function star_or_galaxy(star_column, galaxy_dev_column, galaxy_exp_column)
        ifelse.(is_star, raw_df[star_column], dev_or_exp(galaxy_dev_column, galaxy_exp_column))
    end

    mag_u = star_or_galaxy(:psfmag_u, :devmag_u, :expmag_u)
    mag_g = star_or_galaxy(:psfmag_g, :devmag_g, :expmag_g)
    mag_r = star_or_galaxy(:psfmag_r, :devmag_r, :expmag_r)
    mag_i = star_or_galaxy(:psfmag_i, :devmag_i, :expmag_i)
    mag_z = star_or_galaxy(:psfmag_z, :devmag_z, :expmag_z)

    result = DataFrame()
    result[:objid] = [string(objid) for objid in raw_df[:objid]]
    result[:ra] = raw_df[:ra]
    result[:dec] = raw_df[:dec]
    result[:is_star] = is_star

    flux_r = mag_to_flux.(mag_r, 3)
    result[:flux_r_nmgy] = ifelse.(flux_r .> 0, flux_r, missing)

    result[:color_ug] = color_from_mags.(mag_u, 1, mag_g, 2)
    result[:color_gr] = color_from_mags.(mag_g, 2, mag_r, 3)
    result[:color_ri] = color_from_mags.(mag_r, 3, mag_i, 4)
    result[:color_iz] = color_from_mags.(mag_i, 4, mag_z, 5)

    # gal shape -- fracdev
    result[:gal_frac_dev] = raw_df[:fracdev_r]

    # Note that the SDSS photo pipeline doesn't constrain the de Vaucouleur
    # profile parameters and exponential disk parameters (A/B, angle, scale)
    # to be the same, whereas Celeste does. Here, we pick one or the other
    # from SDSS, based on fracdev - we'll get the parameters corresponding
    # to the dominant component. Later, we limit comparison to objects with
    # fracdev close to 0 or 1 to ensure that we're comparing apples to apples.

    result[:gal_axis_ratio] = dev_or_exp(:devab_r, :expab_r)

    # gal effective radius (re)
    re_arcsec = dev_or_exp(:devrad_r, :exprad_r)
    re_pixel = re_arcsec ./ SDSS_ARCSEC_PER_PIXEL
    result[:gal_radius_px] = convert(Vector{Float64}, re_pixel)

    # gal angle (degrees)
    raw_phi = dev_or_exp(:devphi_r, :expphi_r)
    result[:gal_angle_deg] = canonical_angle.(raw_phi)

    is_saturated = raw_df[:is_saturated] .!= 0
    result = result[.!is_saturated, :]
    bad_rows = [x in BAD_COADD_OBJID for x in result[:, :objid]]
    result = result[.!bad_rows, :]

    # for stars, ensure galaxy-only fields are "missing"
    for col in GALAXY_ONLY_COLUMNS
        result[col] = convert(Vector{Union{Missing, Float64}}, result[col])
        result[result[:is_star], col] = missing
    end

    return result
end

function object_dict_to_data_frame(objects::Dict)
    result_df = DataFrame()
    for (key, values) in objects
        result_df[Symbol(key)] = values
    end
    result_df
end

"""
Load the SDSS photoObj catalog used to initialize celeste, and reformat column
names to match what the rest of the scoring code expects.
"""
function load_primary(rcf::SDSSIO.RunCamcolField, stagedir::String)
    dataset = SDSSDataSet(stagedir)
    raw_df = object_dict_to_data_frame(SDSSIO.read_photoobj(dataset, rcf))

    usedev = raw_df[:frac_dev] .> 0.5  # true=> use dev, false=> use exp
    dev_or_exp(dev_column, exp_column) = ifelse.(usedev, raw_df[dev_column], raw_df[exp_column])
    function star_or_galaxy(star_column, galaxy_dev_column, galaxy_exp_column)
        ifelse.(raw_df[:is_star], raw_df[star_column], dev_or_exp(galaxy_dev_column, galaxy_exp_column))
    end

    flux_u = star_or_galaxy(:psfflux_u, :devflux_u, :expflux_u)
    flux_g = star_or_galaxy(:psfflux_g, :devflux_g, :expflux_g)
    flux_r = star_or_galaxy(:psfflux_r, :devflux_r, :expflux_r)
    flux_i = star_or_galaxy(:psfflux_i, :devflux_i, :expflux_i)
    flux_z = star_or_galaxy(:psfflux_z, :devflux_z, :expflux_z)

    result = DataFrame()
    result[:objid] = raw_df[:objid]
    result[:ra] = raw_df[:ra]
    result[:dec] = raw_df[:dec]
    result[:is_star] = raw_df[:is_star]

    result[:flux_r_nmgy] = flux_r

    result[:color_ug] = color_from_fluxes.(flux_u, flux_g)
    result[:color_gr] = color_from_fluxes.(flux_g, flux_r)
    result[:color_ri] = color_from_fluxes.(flux_r, flux_i)
    result[:color_iz] = color_from_fluxes.(flux_i, flux_z)

    result[:gal_frac_dev] = raw_df[:frac_dev]
    result[:gal_axis_ratio] = dev_or_exp(:ab_dev, :ab_exp)

    # gal effective radius (re)
    re_arcsec = dev_or_exp(:theta_dev, :theta_exp)
    re_pixel = re_arcsec ./ SDSS_ARCSEC_PER_PIXEL
    result[:gal_radius_px] = convert(Vector{Float64}, re_pixel)

    # gal angle (degrees)
    raw_phi = dev_or_exp(:phi_dev, :phi_exp)
    result[:gal_angle_deg] = canonical_angle.(raw_phi)

    # primary is better at indicating oversaturation than coadd
    is_saturated = flux_to_mag.(raw_df[:psfflux_r], 3) .< 16
    result = result[.!is_saturated, :]

    return result
end

function fluxes_from_colors(flux_r_nmgy::Float64, colors::AbstractVector)
    @assert length(colors) == 4
    color_ratios = exp.(colors)
    fluxes = similar(colors, 5)
    fluxes[3] = flux_r_nmgy
    fluxes[4] = fluxes[3] * color_ratios[3]
    fluxes[5] = fluxes[4] * color_ratios[4]
    fluxes[2] = fluxes[3] / color_ratios[2]
    fluxes[1] = fluxes[2] / color_ratios[1]
    fluxes
end

function get_median_fluxes(variational_params::Vector{Float64}, source_type::Int64)
    fluxes_from_colors(
        exp(variational_params[Model.ids.flux_loc[source_type]]),
        variational_params[Model.ids.color_mean[:, source_type]],
    )
end

function variational_parameters_to_data_frame_row(variational_params::Vector{Float64})
    ids = Model.ids
    result = DataFrame()
    result[:ra] = variational_params[ids.pos[1]]
    result[:dec] = variational_params[ids.pos[2]]
    result[:is_star] = variational_params[ids.is_star[1, 1]]
    result[:gal_frac_dev] = variational_params[ids.gal_frac_dev]
    result[:gal_axis_ratio] = variational_params[ids.gal_axis_ratio]
    result[:gal_radius_px] = (
        variational_params[ids.gal_radius_px] * sqrt(variational_params[ids.gal_axis_ratio])
    )
    result[:gal_angle_deg] = canonical_angle(180 / pi * variational_params[ids.gal_angle])

    star_galaxy_index = (result[1, :is_star] > 0.5 ? 1 : 2)
    fluxes = get_median_fluxes(variational_params, star_galaxy_index)
    result[:flux_r_nmgy] = fluxes[3]
    result[:color_ug] = color_from_fluxes(fluxes[1], fluxes[2])
    result[:color_gr] = color_from_fluxes(fluxes[2], fluxes[3])
    result[:color_ri] = color_from_fluxes(fluxes[3], fluxes[4])
    result[:color_iz] = color_from_fluxes(fluxes[4], fluxes[5])

    result[:log_flux_r_stderr] = sqrt(variational_params[ids.flux_scale[star_galaxy_index]])
    result[:color_ug_stderr] = sqrt(variational_params[ids.color_var[1, star_galaxy_index]])
    result[:color_gr_stderr] = sqrt(variational_params[ids.color_var[2, star_galaxy_index]])
    result[:color_ri_stderr] = sqrt(variational_params[ids.color_var[3, star_galaxy_index]])
    result[:color_iz_stderr] = sqrt(variational_params[ids.color_var[4, star_galaxy_index]])

    result
end

"""
Convert Celeste results to a dataframe.
"""
function celeste_to_df(results::Vector{ParallelRun.OptimizedSource})
    rows = []
    for result in results
        if !result.is_sky_bad
            row = variational_parameters_to_data_frame_row(result.vs)
            push!(rows, row)
        end
    end
    vcat(rows...)
end

################################################################################
# Generate a random catalog from the Celeste prior
################################################################################

# This is the ratio of stars derived from the catalog used the generate the prior; the value 0.99 is
# currently used in inference to account for the extra flexibility of the galaxy model
const PRIOR_PROBABILITY_OF_STAR = 0.28

"""
Draw a random source from the Celeste prior, returning a 1-row data frame.
"""
function draw_source_params(prior)
    is_star = (rand(Bernoulli(PRIOR_PROBABILITY_OF_STAR)) == 1)
    source_type_index = is_star ? 1 : 2
    flux_r_nmgy = exp(rand(
        Normal(prior.flux_mean[source_type_index], sqrt(prior.flux_var[source_type_index]))
    ))

    color_mixture_weights = prior.k[:, source_type_index]
    NUM_COLOR_COMPONENTS = length(color_mixture_weights)
    color_components = MvNormal[
        MvNormal(prior.color_mean[:, k, source_type_index], prior.color_cov[:, :, k, source_type_index])
        for k in 1:NUM_COLOR_COMPONENTS
    ]
    colors = rand(MixtureModel(color_components, color_mixture_weights))

    if !is_star
        gal_radius_px = exp(rand(
            Normal(prior.gal_radius_px_mean, sqrt(prior.gal_radius_px_var))
        ))
        gal_angle_deg = rand(Uniform(0, 180))
        gal_axis_ratio = rand(Beta(2, 2))
        gal_frac_dev = rand(Beta(0.5, 0.5))
    else
        gal_radius_px = -1
        gal_angle_deg = -1
        gal_axis_ratio = -1
        gal_frac_dev = -1
    end

    # We should change ra/dec ranges below if we start using another field
    @assert STRIPE82_RCF.run == 4263 && STRIPE82_RCF.camcol == 5 && STRIPE82_RCF.field == 119

    # Use approximate size of SDSS field in degrees
    ra = rand(Uniform(0.443, 0.606))
    dec = rand(Uniform(0.411, 0.635))

    DataFrame(
        ra=ra,
        dec=dec,
        is_star=is_star,
        flux_r_nmgy=flux_r_nmgy,
        color_ug=colors[1],
        color_gr=colors[2],
        color_ri=colors[3],
        color_iz=colors[4],
        gal_frac_dev=gal_frac_dev,
        gal_axis_ratio=gal_axis_ratio,
        gal_radius_px=gal_radius_px,
        gal_angle_deg=gal_angle_deg,
    )
end

"""
Draw sources at random from Celeste prior, returning a catalog DF.
"""
function generate_catalog_from_celeste_prior(num_sources::Int64, seed::Int64)
    srand(seed)
    prior = Model.load_prior()
    result = vcat(
        [draw_source_params(prior) for index in 1:num_sources]...
    )

    # for stars, ensure galaxy-only fields are "missing"
    for col in GALAXY_ONLY_COLUMNS
        result[col] = convert(Vector{Union{Missing, Float64}}, result[col])
        result[result[:is_star], col] = missing
    end

    return result
end

################################################################################
# Support for running Celeste on test imagery
################################################################################

## Load a multi-extension FITS imagery file

struct FitsImage
    pixels::Matrix{Float32}
    header::FITSIO.FITSHeader
    wcs::WCS.WCSTransform
end

function read_fits(filename::String)
    if !isfile(filename)
        throw(BenchmarkFitsFileNotFound(filename))
    end

    fits = FITSIO.FITS(filename)
    try
        map(fits) do extension
            pixels = read(extension)
            header = FITSIO.read_header(extension)
            wcs = WCS.from_header(FITSIO.read_header(fits[1], String))[1]
            FitsImage(pixels, header, wcs)
        end
    finally
        close(fits)
    end
end

## Create Celeste images from FITS imagery

function make_simple_psf(psf_sigma_px)
    alphaBar = [1.; 0.]
    xiBar = [0.; 0.]
    tauBar = [psf_sigma_px^2 0.; 0. psf_sigma_px^2]
    [
        Model.PsfComponent(
            alphaBar[k],
            StaticArrays.SVector{2, Float64}(xiBar),
            StaticArrays.SMatrix{2, 2, Float64, 4}(tauBar)
        )
        for k in 1:2
    ]
end

function make_psf_from_header(header::FITSIO.FITSHeader)
    if haskey(header, "CLSIGMA")
        return make_simple_psf(header["CLSIGMA"])
    end

    num_components = header["PSFNCOMP"]
    map(1:num_components) do index
        get(field) = header["PSF$(index)_$(field)"]
        alpha_bar = get("A")
        xi_bar = Float64[get("X1"), get("X2")]
        tau_bar = Float64[get("T11") get("T12"); get("T12") get("T22")]
        Model.PsfComponent(
            alpha_bar,
            StaticArrays.SVector{2, Float64}(xi_bar),
            StaticArrays.SMatrix{2, 2, Float64, 4}(tau_bar),
        )
    end
end

function serialize_psf_to_header(psf::Vector{Model.PsfComponent}, header::FITSIO.FITSHeader)
    header["PSFNCOMP"] = length(psf)
    for (index, component) in enumerate(psf)
        @assert component.tauBar[1, 2] == component.tauBar[2, 1]
        function set(field, value)
            header["PSF$(index)_$(field)"] = value
        end
        set("A", component.alphaBar)
        set("X1", component.xiBar[1])
        set("X2", component.xiBar[2])
        set("T11", component.tauBar[1, 1])
        set("T12", component.tauBar[1, 2])
        set("T22", component.tauBar[2, 2])
    end
end

function make_image(pixels::Matrix{Float32}, band_index::Int,
                    wcs::WCS.WCSTransform, psf::Vector{Model.PsfComponent},
                    sky_level_nmgy::Real, nelec_per_nmgy::Real)
    nx, ny = size(pixels)

    # Render the PSF on a grid, to be used as a (spatially constant) PSF map.
    psfstamp = Model.render_psf(psf, (51, 51))

    Model.Image(
        pixels,
        band_index,
        wcs,
        psf,
        fill(Float32(sky_level_nmgy), nx, ny),
        fill(Float32(nelec_per_nmgy), nx),
        Model.ConstantPSFMap(psfstamp)
    )
end

function make_images(band_extensions::Vector{FitsImage})
    map(enumerate(band_extensions)) do pair
        band_index, extension = pair
        make_image(
            extension.pixels,
            band_index,
            extension.wcs,
            make_psf_from_header(extension.header),
            extension.header["CLSKY"]::Float64,
            extension.header["CLIOTA"]::Float64
        )
    end
end

## Create an initialization catalog for Celeste

function typical_band_fluxes(is_star::Bool)
    source_type_index = is_star ? 1 : 2
    prior_parameters::Model.PriorParams = Model.load_prior()
    # this is the mode. brightness is log normal.
    flux_r = exp(
        prior_parameters.flux_mean[source_type_index] - prior_parameters.flux_var[source_type_index]
    )
    # Band relative intensities are a mixture of lognormals. Which mixture component has the most
    # weight?
    dominant_component = indmax(prior_parameters.k[:, source_type_index])
    # What are the most typical log relative intensities for that component?
    colors = (
        prior_parameters.color_mean[:, dominant_component, source_type_index]
        - diag(prior_parameters.color_cov[:, :, dominant_component, source_type_index])
    )
    fluxes_from_colors(flux_r, colors)
end

function make_catalog_entry(x_position_world_coords::Float64,
                            y_position_world_coords::Float64)
    Model.CatalogEntry(
        [x_position_world_coords, y_position_world_coords],
        false, # is_star
        typical_band_fluxes(true),
        typical_band_fluxes(false),
        0.1, # gal_frac_dev
        0.7, # gal_axis_ratio
        pi / 4, # gal_angle
        4., # gal_radius_px
    )
end

ensure_small_flux(value) = (ismissing(value) || value <= 0) ? 1e-6 : value

na_to_default(value, default) = ismissing(value) ? default : value

function make_catalog_entry(row::DataFrameRow)
    colors = Union{Float64, Missing}[row[:color_ug], row[:color_gr],
                                     row[:color_ri], row[:color_iz]]
    fluxes = fluxes_from_colors(row[:flux_r_nmgy], colors)
    fluxes = ensure_small_flux.(fluxes)
    gal_axis_ratio = na_to_default(row[:gal_axis_ratio], 0.8)
    Model.CatalogEntry(
        [row[:ra], row[:dec]],
        row[:is_star] > 0.5,
        fluxes,
        fluxes,
        na_to_default(row[:gal_frac_dev], 0.5),
        gal_axis_ratio,
        na_to_default(row[:gal_angle_deg], 0.) / 180.0 * pi,
        na_to_default(row[:gal_radius_px] / sqrt(gal_axis_ratio), 2.)
        )
end

function make_initialization_catalog(catalog::DataFrame, use_full_initialzation::Bool)
    position_offset_width = SDSS_ARCSEC_PER_PIXEL / ARCSEC_PER_DEGREE # 1 pixel, in degrees
    map(eachrow(catalog)) do row
        if use_full_initialzation
            make_catalog_entry(row)
        else
            make_catalog_entry(
                row[:ra] + position_offset_width,
                row[:dec] - 0.5 * position_offset_width
            )
        end
    end
end

################################################################################
# Support for generating imagery using Synthetic.jl
################################################################################

struct ImageGeometry
    height_px::Int64
    width_px::Int64
    world_coordinate_origin::Tuple{Float64, Float64} # (ra, dec)
    ra_degrees_per_pixel::Float64
    dec_degrees_per_pixel::Float64
end

function get_image_geometry(catalog_data::DataFrame; field_expand_arcsec=20.0)
    min_ra_deg = minimum(catalog_data[:ra])
    max_ra_deg = maximum(catalog_data[:ra])
    min_dec_deg = minimum(catalog_data[:dec])
    max_dec_deg = maximum(catalog_data[:dec])

    width_arcsec = (max_ra_deg - min_ra_deg) * ARCSEC_PER_DEGREE + 2 * field_expand_arcsec
    height_arcsec = (max_dec_deg - min_dec_deg) * ARCSEC_PER_DEGREE + 2 * field_expand_arcsec
    width_px = convert(Int64, round(width_arcsec / SDSS_ARCSEC_PER_PIXEL))
    height_px = convert(Int64, round(height_arcsec / SDSS_ARCSEC_PER_PIXEL))

    dec_degrees_per_pixel = SDSS_ARCSEC_PER_PIXEL / ARCSEC_PER_DEGREE
    ra_degrees_per_pixel = dec_degrees_per_pixel / cosd(min_dec_deg)

    ImageGeometry(
        height_px,
        width_px,
        (
            min_ra_deg - field_expand_arcsec / ARCSEC_PER_DEGREE,
            min_dec_deg - field_expand_arcsec / ARCSEC_PER_DEGREE,
        ),
        ra_degrees_per_pixel,
        dec_degrees_per_pixel,
    )
end

function make_template_images(
    catalog_data::DataFrame, psf_sigma_px::Float64, band_sky_level_nmgy::Vector{Float64},
    band_nelec_per_nmgy::Vector{Float64}
)
    @assert length(band_sky_level_nmgy) == 5
    @assert length(band_nelec_per_nmgy) == 5
    geometry = get_image_geometry(catalog_data)
    println("  Image dimensions $(geometry.height_px) H x $(geometry.width_px) W px")
    wcs = WCS.WCSTransform(
        2, # dimensions
        # reference pixel coordinates...
        crpix=[1., 1.],
        # ...and corresponding reference world coordinates
        crval=[geometry.world_coordinate_origin[1], geometry.world_coordinate_origin[2]],
        # this WCS is a simple linear transformation
        ctype=["RA---TAN", "DEC--TAN"],
        cunit=["deg", "deg"],
        # these are [du/dx du/dy; dv/dx dv/dy]. (u, v) = world coords, (x, y) = pixel coords.
        pc=[0. geometry.ra_degrees_per_pixel; geometry.dec_degrees_per_pixel 0.],
    )
    map(1:5) do band
        make_image(
            zeros(Float32, (geometry.height_px, geometry.width_px)),
            band,
            wcs,
            make_simple_psf(psf_sigma_px),
            band_sky_level_nmgy[band],
            band_nelec_per_nmgy[band],
        )
    end
end

## Parse a FITS header from a raw string representation

# It's silly that we have to do this, but the Julia FITSIO library doesn't support this (it can only
# read a header from a file), while the Julia WCS library will only generate a header in raw string
# form. (The Julia FITSIO library also won't write a raw string header to a file, so scrap that idea
# :).) Fortunately, the FITS header format is a simple fixed-width textual format with only a few
# data types. See https://fits.gsfc.nasa.gov/fits_primer.html.

function parse_fits_header_value(raw_value::AbstractString)
    if raw_value == "T"
        true
    elseif raw_value == "F"
        false
    elseif raw_value[1] == '\''
        convert(String, raw_value[2:(end - 1)])
    elseif ismatch(r"^[0-9]+$", raw_value)
        parse(Int64, raw_value)
    else
        parse(Float64, raw_value)
    end
end

function parse_fits_header_from_string(header_string::AbstractString)
    header = FITSIO.FITSHeader(String[], [], String[])
    start_position = 1
    while start_position < length(header_string)
        field_string = header_string[start_position:(start_position + 79)]
        field_name = strip(field_string[1:8])
        if length(field_name) != 0 && field_string[9:10] == "= "
            field_parts = split(field_string[11:80], "/", limit=2)
            raw_value = field_parts[1]
            header[field_name] = parse_fits_header_value(strip(raw_value))
            if length(field_parts) == 2
                comment = strip(field_parts[2])
                if length(comment) > 0
                    FITSIO.set_comment!(header, field_name, convert(String, comment))
                end
            end
        end
        start_position += 80
    end
    header
end

function save_images_to_fits(filename::String, images::Vector{<:Model.Image})
    println("Writing images to $filename...")
    fits_file = FITSIO.FITS(filename, "w")
    for band_image in images
        header = parse_fits_header_from_string(WCS.to_header(band_image.wcs))
        serialize_psf_to_header(band_image.psf, header)
        header["CLSKY"] = mean(band_image.sky.sky_small) * mean(band_image.sky.calibration)
        FITSIO.set_comment!(header, "CLSKY", "Mean sky background per pixel, nMgy")
        header["CLIOTA"] = mean(band_image.nelec_per_nmgy)
        FITSIO.set_comment!(header, "CLIOTA", "Gain, nelec per nMgy")
        write(
            fits_file,
            band_image.pixels,
            header=header,
            name="CELESTE_FIELD_$(band_image.b)",
        )
    end
    close(fits_file)
end

################################################################################
# Score a set of predictions against a ground truth catalog
################################################################################

COLOR_COLUMNS = [:color_ug, :color_gr, :color_ri, :color_iz]

ABSOLUTE_ERROR_COLUMNS = vcat(
    [:gal_frac_dev, :gal_axis_ratio, :gal_radius_px],
    COLOR_COLUMNS,
)

function degrees_to_diff(a, b)
    angle_between = abs.(a - b) .% 180
    min.(angle_between, 180 .- angle_between)
end

"""
Given two results data frame, one containing ground truth (i.e Coadd)
and one containing predictions (i.e., either Primary of Celeste),
compute an a data frame containing each prediction's error.
(It's not an average of the errors, it's each error.)
Let's call the return type of this function an \"error data frame\".
"""
function get_error_df(truth::DataFrame, predicted::DataFrame)
    errors = DataFrame()

    predicted_galaxy = predicted[:is_star] .< .5
    true_galaxy = truth[:is_star] .< .5
    errors[:missed_stars] = ifelse.(.!true_galaxy, predicted_galaxy, missing)
    errors[:missed_galaxies] = ifelse.(true_galaxy, .!predicted_galaxy, missing)

    errors[:position] =
        (ARCSEC_PER_DEGREE / SDSS_ARCSEC_PER_PIXEL) .*
        angular_separation.(truth[:ra],
                            truth[:dec],
                            predicted[:ra],
                            predicted[:dec])

    # compare flux in both mags and nMgy for now
    errors[:flux_r_mag] = abs.(
        flux_to_mag.(truth[:flux_r_nmgy], 3)
        .- flux_to_mag.(predicted[:flux_r_nmgy], 3)
    )
    errors[:flux_r_nmgy] = abs.(
        truth[:flux_r_nmgy] .- predicted[:flux_r_nmgy]
    )
    errors[:gal_angle_deg] = degrees_to_diff(truth[:gal_angle_deg], predicted[:gal_angle_deg])

    for column_symbol in ABSOLUTE_ERROR_COLUMNS
        errors[column_symbol] = abs.(truth[column_symbol] .- predicted[column_symbol])
    end
    for color_column in COLOR_COLUMNS
        # to match up with Stripe82Score, which used differences of mags
        errors[color_column] .*= 2.5 / log(10)
    end

    errors
end

function is_good_row(truth_row::DataFrameRow, error_row::DataFrameRow, column_name::Symbol)
    if ismissing(error_row[column_name]) || isnan(error_row[column_name])
        return false
    elseif !ismissing(truth_row[:gal_radius_px]) && truth_row[:gal_radius_px] > 20
        return false
    end

    if column_name in [:gal_axis_ratio, :gal_radius_px, :gal_angle_deg,
                       :gal_frac_dev]
        has_mixture_weight = !ismissing(truth_row[:gal_frac_dev])
        if has_mixture_weight && (0.05 < truth_row[:gal_frac_dev] < 0.95)
            return false
        end
    end
    if column_name == :gal_angle_deg
        if !ismissing(truth_row[:gal_axis_ratio]) && truth_row[:gal_axis_ratio] > .6
            return false
        end
    end
    return true
end

function filter_rows(truth::DataFrame, errors::DataFrame, column_name::Symbol)
    map(zip(eachrow(truth), eachrow(errors))) do rows
        (truth_row, error_row) = rows
        is_good_row(truth_row, error_row, column_name)
    end
end

function score_column(errors::AbstractVector)
    DataFrame(
        N=length(errors),
        first=mean(errors),
    )
end

function score_column(first_errors::AbstractVector,
                      second_errors::AbstractVector)
    @assert length(first_errors) == length(second_errors)
    scores = score_column(first_errors)
    scores[:second] = mean(second_errors)
    diffs = first_errors .- second_errors
    scores[:diff] = mean(diffs)
    scores[:diff_sd] = std(abs.(diffs)) / sqrt(length(diffs))
    scores
end

function get_scores_df(
    truth::DataFrame, first_errors::DataFrame, second_errors::Nullable{DataFrame}
)
    score_rows = DataFrame[]
    for column_name in names(first_errors)
        good_row = filter_rows(truth, first_errors, column_name)
        if !isnull(second_errors)
            good_row .&= filter_rows(truth, get(second_errors), column_name)
        end
        if sum(good_row) <= 1
            continue
        else
            if isnull(second_errors)
                row = score_column(first_errors[good_row, column_name])
            else
                row = score_column(
                    first_errors[good_row, column_name],
                    get(second_errors)[good_row, column_name],
                )
            end

            row[:field] = column_name
            push!(score_rows, row)
        end
    end
    vcat(score_rows...)
end

"""
    match_catalogs(truth, predictions)

Return subset of rows in `truth` and `predictions`, such that each entry in
`truth` has a corresponding match in all `predictions`. A "match" is the
closest object within `tol` degrees. The output data frames all have the same
number of rows, and rows correspond to matching objects.
"""
function match_catalogs(truth::DataFrame, prediction_dfs::Vector{DataFrame};
                        tol::Float64=(0.396 / 3600.))

    matched = trues(nrow(truth))
    idxs = Vector{Int}[]
    for prediction in prediction_dfs
        # remove missings before feeding to match_coordinates
        @assert sum(ismissing.(truth[:ra])) == 0
        @assert sum(ismissing.(truth[:dec])) == 0
        @assert sum(ismissing.(prediction[:ra])) == 0
        @assert sum(ismissing.(prediction[:dec])) == 0
        idx, dists = match_coordinates(Vector{Float64}(truth[:ra]),
                                       Vector{Float64}(truth[:dec]),
                                       Vector{Float64}(prediction[:ra]),
                                       Vector{Float64}(prediction[:dec]))
        matched .&= dists .< tol
        push!(idxs, idx)
    end

    # For each row in `truth`, `matched` is now true if and only if all
    # predictions have a close match.
    #
    # `idxs[i]` is same length as `truth` and contains the index in
    # `prediction_dfs[i]` that is the closest match to each row in
    # `truth`.
    matched_truth = truth[matched, :]
    matched_predictions = [prediction[idx[matched], :]
                           for (prediction, idx) in zip(prediction_dfs, idxs)]

    return matched_truth, matched_predictions
end

function score_predictions(truth::DataFrame, prediction_dfs::Vector{DataFrame})
    matched_truth, matched_prediction_dfs = match_catalogs(truth, prediction_dfs)
    error_dfs = [get_error_df(matched_truth, predictions) for predictions in matched_prediction_dfs]
    @assert length(prediction_dfs) <= 2
    get_scores_df(
        matched_truth,
        error_dfs[1],
        length(error_dfs) > 1 ? Nullable(error_dfs[2]) : Nullable{DataFrame}(),
    )
end

function get_uncertainty_df(truth::DataFrame, predictions::DataFrame)
    assert_columns_are_present(predictions, STDERR_COLUMNS)
    matched_truth, matched_prediction_dfs = match_catalogs(truth, [predictions])
    matched_predictions = matched_prediction_dfs[1]

    valid_rows = (matched_truth[:flux_r_nmgy] .> 0)
    matched_truth = matched_truth[valid_rows, :]
    matched_predictions = matched_predictions[valid_rows, :]

    get_errors(column, map_fn) =
        map_fn.(matched_predictions[column]) .- map_fn.(matched_truth[column])
    get_errors(column) = get_errors(column, x -> x)
    errors = [
        get_errors(:flux_r_nmgy, log),
        get_errors(:color_ug),
        get_errors(:color_gr),
        get_errors(:color_ri),
        get_errors(:color_iz),
    ]
    std_errs = [
        matched_predictions[:log_flux_r_stderr],
        matched_predictions[:color_ug_stderr],
        matched_predictions[:color_gr_stderr],
        matched_predictions[:color_ri_stderr],
        matched_predictions[:color_iz_stderr],
    ]
    names = [:log_flux_r_nmgy, :color_ug, :color_gr,
             :color_ri, :color_iz]

    matched_truth[:log_flux_err] = errors[1]
    matched_truth[:log_flux_stderr] = std_errs[1]
    matched_truth[:flux_r_celeste] = matched_predictions[:flux_r_nmgy]
    CSV.write("stderr.csv", matched_truth)

    mapreduce(vcat, zip(names, errors, std_errs)) do values
        name, error, std_err = values
        DataFrame(name=fill(name, length(error)),
                  error=error,
                  posterior_std_err=std_err)
    end
end

function score_uncertainty(uncertainty_df::DataFrame)
    mapreduce(vcat, groupby(uncertainty_df, :name)) do group_df
        abs_error_sds = abs.(group_df[:error] ./ group_df[:posterior_std_err])
        abs_error_sds = abs_error_sds[.!ismissing.(abs_error_sds)]
        DataFrame(
            field=group_df[1, :name],
            within_half_sd=mean(abs_error_sds .<= 1/2),
            within_1_sd=mean(abs_error_sds .<= 1),
            within_2_sd=mean(abs_error_sds .<= 2),
            within_3_sd=mean(abs_error_sds .<= 3),
        )
    end
end

function plot_image(img)
    xs = img.pixels'

    xs -= 550
    black = minimum(filter(.!isnan, xs))
    xs = map((x)->isnan(x) ? black : x, xs)
    xs = min.(xs, 10_000)
    xs = log.(xs + 100)
    #cutoffs = quantile(xs2[:], 0:0.01:1)
    #xs3 = map(x->findfirst(x .< cutoffs), xs2)
    xs -= log(black + 100)
    xs /= log(10_100)
    return xs
end

end # module AccuracyBenchmark


text:
module GRAMPC

using BitFlags
using GRAMPC_jll
using LinearAlgebra

const CallbackArg = Union{Symbol, Function}
const Callback    = Union{Function, Nothing}
const MassMatrix  = Union{Array{Float64}, UniformScaling}

# Define the library error handling
include( "cinterface/errorHandling.jl" )

# Constants used in this package
include( "constants.jl" )

# Must come before the ctypes because this defines the wrapper struct
include( "cinterface/wrapperFunctions.jl" )

# Most types used in this package
include( "cinterface/ctypes.jl" )
include( "types.jl" )

# Low-level problem translation for the C wrapper
include( "cinterface/options.jl" )
include( "cinterface/interface.jl" )

# High-level Julia part of the interface
include( "callbacks/dynamics.jl" )
include( "callbacks/cost.jl" )
include( "callbacks/constraints.jl" )
include( "settings.jl" )

end

ext: jl


text:
"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.

# Examples
```jldoctest
julia> m = Chain(x -> x^2, x -> x+1);

julia> m(5) == 26
true

julia> m = Chain(Dense(10, 5), Dense(5, 2));

julia> x = rand(10);

julia> m(x) == m[2](m[1](x))
true
```
"""
struct Chain{T<:Tuple}
  layers::T
  Chain(xs...) = new{typeof(xs)}(xs)
end

@forward Chain.layers Base.getindex, Base.length, Base.first, Base.last,
  Base.iterate, Base.lastindex

functor(::Type{<:Chain}, c) = c.layers, ls -> Chain(ls...)

applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(tail(fs), first(fs)(x))

(c::Chain)(x) = applychain(c.layers, x)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

# This is a temporary and naive implementation
# it might be replaced in the future for better performance
# see issue https://github.com/FluxML/Flux.jl/issues/702
# Johnny Chen -- @johnnychen94
# only slightly changed to better handle interaction with Zygote @dsweber2
"""
    activations(c::Chain, input)

Calculate the forward results of each layers in Chain `c` with `input` as model input.
"""
function activations(c::Chain, input)
    extraChain(c.layers, input)
end

function extraChain(fs::Tuple, x)
    res = first(fs)(x)
    return (res, extraChain(Base.tail(fs), res)...)
end

extraChain(::Tuple{}, x) = ()

"""
    Dense(in, out, σ=identity; bias=true, init=glorot_uniform)
    Dense(W::AbstractMatrix, [bias, σ])

Create a traditional `Dense` layer, whose forward pass is given by:

    y = σ.(W * x .+ bias)

The input `x` should be a vector of length `in`, or batch of vectors represented
as an `in × N` matrix, or any array with `size(x,1) == in`.
The out `y` will be a vector  of length `out`, or a batch with
`size(y) == (out, size(x)[2:end]...)`

Keyword `bias=false` will switch off trainable bias for the layer.
The initialisation of the weight matrix is `W = init(out, in)`, calling the function
given to keyword `init`, with default [`glorot_uniform`](@doc Flux.glorot_uniform).
The weight matrix and/or the bias vector (of length `out`) may also be provided explicitly.

# Examples
```jldoctest
julia> d = Dense(5, 2)
Dense(5, 2)

julia> d(rand(Float32, 5, 64)) |> size
(2, 64)

julia> d(rand(Float32, 5, 1, 1, 64)) |> size  # treated as three batch dimensions
(2, 1, 1, 64)

julia> d1 = Dense(ones(2, 5), false, tanh)  # using provided weight matrix
Dense(5, 2, tanh; bias=false)

julia> d1(ones(5))
2-element Vector{Float64}:
 0.9999092042625951
 0.9999092042625951

julia> Flux.params(d1)  # no trainable bias
Params([[1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]])
```
"""
struct Dense{F, M<:AbstractMatrix, B}
  weight::M
  bias::B
  σ::F
  function Dense(W::M, bias = true, σ::F = identity) where {M<:AbstractMatrix, F}
    b = create_bias(W, bias, size(W,1))
    new{F,M,typeof(b)}(W, b, σ)
  end
end

function Dense(in::Integer, out::Integer, σ = identity;
               initW = nothing, initb = nothing,
               init = glorot_uniform, bias=true)

  W = if initW !== nothing
    Base.depwarn("keyword initW is deprecated, please use init (which similarly accepts a funtion like randn)", :Dense)
    initW(out, in)
  else
    init(out, in)
  end

  b = if bias === true && initb !== nothing
    Base.depwarn("keyword initb is deprecated, please simply supply the bias vector, bias=initb(out)", :Dense)
    initb(out)
  else
    bias
  end

  return Dense(W, b, σ)
end

@functor Dense

function (a::Dense)(x::AbstractVecOrMat)
  W, b, σ = a.weight, a.bias, a.σ
  return σ.(W*x .+ b)
end

(a::Dense)(x::AbstractArray) =
  reshape(a(reshape(x, size(x,1), :)), :, size(x)[2:end]...)

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.weight, 2), ", ", size(l.weight, 1))
  l.σ == identity || print(io, ", ", l.σ)
  l.bias == Zeros() && print(io, "; bias=false")
  print(io, ")")
end

"""
    Diagonal(α, β)
    Diagonal(size::Integer...)

Create an element-wise linear layer, which performs

    y = α .* x .+ β

The learnable arrays are initialised `α = ones(Float32, size)` and
`β = zeros(Float32, size)`.

Used by [`LayerNorm`](@ref).
"""
struct Diagonal{T}
  α::T
  β::T
end

function Diagonal(sz::Integer...; initα = nothing, initβ = nothing)
  α = if initα !== nothing
    Base.depwarn("keyword initα is deprecated, please simply supply the desired vectors", :Diagonal)
    initα(sz...)
  else
    ones(sz...)
  end
  β = if initβ !== nothing
    Base.depwarn("keyword initβ is deprecated, please simply supply the desired vectors", :Diagonal)
    initβ(sz...)
  else
    zeros(sz...)
  end
  Diagonal(α, β)
end

@functor Diagonal

(a::Diagonal)(x) = a.α .* x .+ a.β

function Base.show(io::IO, l::Diagonal)
  print(io, "Diagonal(", join(size(l.α), ", "), ")")
end

"""
    Maxout(over)

The [Maxout](https://arxiv.org/abs/1302.4389) layer has a number of
internal layers which all receive the same input. It returns the elementwise
maximum of the internal layers' outputs.

Maxout over linear dense layers satisfies the univeral approximation theorem.
"""
struct Maxout{FS<:Tuple}
    over::FS
end

"""
    Maxout(f, n_alts)

Construct a Maxout layer over `n_alts` instances of the layer given by `f`.
The function takes no arguments and should return some callable layer.
Conventionally, this is a linear dense layer.

# Examples

This constructs a `Maxout` layer over 4 internal dense linear layers, each
identical in structure (784 inputs, 128 outputs):
```jldoctest
julia> insize = 784;

julia> outsize = 128;

julia> Maxout(()->Dense(insize, outsize), 4);
```
"""
function Maxout(f, n_alts)
  over = Tuple(f() for _ in 1:n_alts)
  return Maxout(over)
end

@functor Maxout

function (mo::Maxout)(input::AbstractArray)
    mapreduce(f -> f(input), (acc, out) -> max.(acc, out), mo.over)
end

"""
    SkipConnection(layer, connection)

Create a skip connection which consists of a layer or `Chain` of consecutive
layers and a shortcut connection linking the block's input to the output
through a user-supplied 2-argument callable. The first argument to the callable
will be propagated through the given `layer` while the second is the unchanged,
"skipped" input.

The simplest "ResNet"-type connection is just `SkipConnection(layer, +)`.
Here is a more complicated example:
```jldoctest
julia> m = Conv((3,3), 4 => 7, pad=(1,1));

julia> x = ones(Float32, 5, 5, 4, 10);

julia> size(m(x)) == (5, 5, 7, 10)
true

julia> sm = SkipConnection(m, (mx, x) -> cat(mx, x, dims=3));

julia> size(sm(x)) == (5, 5, 11, 10)
true
```
"""
struct SkipConnection{T,F}
  layers::T
  connection::F  #user can pass arbitrary connections here, such as (a,b) -> a + b
end

@functor SkipConnection

function (skip::SkipConnection)(input)
  skip.connection(skip.layers(input), input)
end

function Base.show(io::IO, b::SkipConnection)
  print(io, "SkipConnection(", b.layers, ", ", b.connection, ")")
end

"""
    Bilinear(in1, in2, out, σ=identity; bias=true, init=glorot_uniform)
    Bilinear(W::AbstractArray, [bias, σ])

Creates a Bilinear layer, which operates on two inputs at the same time.
Its output, given vectors `x` & `y`, is another vector `z` with,
for all `i ∈ 1:out`:

    z[i] = σ(x' * W[i,:,:] * y + bias[i])

If `x` and `y` are matrices, then each column of the output `z = B(x, y)` is of this form,
with `B` a Bilinear layer.

If `y` is not given, it is taken to be equal to `x`, i.e. `B(x) == B(x, x)`

The two inputs may also be provided as a tuple, `B((x, y)) == B(x, y)`,
which is accepted as the input to a `Chain`.

The initialisation works as for [`Dense`](@ref) layer, with `W = init(out, in1, in2)`.
By default the bias vector is `zeros(Float32, out)`, option `bias=false` will switch off
trainable bias. Either of these may be provided explicitly.

# Examples
```jldoctest
julia> x, y = randn(Float32, 5, 32), randn(Float32, 5, 32);

julia> B = Flux.Bilinear(5, 5, 7);

julia> B(x) |> size  # interactions based on one input
(7, 32)

julia> B(x,y) == B((x,y))  # two inputs, may be given as a tuple
true

julia> sc = SkipConnection(
                Chain(Dense(5, 20, tanh), Dense(20, 9, tanh)),
                Flux.Bilinear(9, 5, 3, bias=false),
            );  # used as the recombinator, with skip as the second input

julia> sc(x) |> size
(3, 32)

julia> Flux.Bilinear(rand(4,8,16), false, tanh)  # first dim of weight is the output
Bilinear(8, 16, 4, tanh, bias=false)
```
"""
struct Bilinear{F,A,B}
  weight::A
  bias::B
  σ::F
  function Bilinear(W::A, bias = true, σ::F = identity) where {A<:AbstractArray, F}
    ndims(A) == 3 || throw(ArgumentError("expected a 3-array of weights"))
    b = create_bias(W, bias, size(W,1))
    new{F,A,typeof(b)}(W, b, σ)
  end
end

@functor Bilinear

function Bilinear(in1::Integer, in2::Integer, out::Integer, σ = identity;
                  init = glorot_uniform, bias = true)
  Bilinear(init(out, in1, in2), bias, σ)
end

function (a::Bilinear)(x::AbstractMatrix, y::AbstractMatrix)
  W, b, σ = a.weight, a.bias, a.σ

  d_z, d_x, d_y = size(W)
  d_x == size(x,1) && d_y == size(y,1) || throw(DimensionMismatch("number of rows in data must match W"))
  size(x,2) == size(y,2) || throw(DimensionMismatch("Data inputs must agree on number of columns, got $(size(x,2)) and $(size(y,2))"))

  # @einsum Wy[o,i,s] := W[o,i,j] * y[j,s]
  Wy = reshape(reshape(W, (:, d_y)) * y, (d_z, d_x, :))

  # @einsum Z[o,s] := Wy[o,i,s] * x[i,s]
  Wyx = batched_mul(Wy, reshape(x, (d_x, 1, :)))
  Z = reshape(Wyx, (d_z, :))

  # @einsum out[o,s] := σ(Z[o,i] + b[o])
  σ.(Z .+ b)
end

(a::Bilinear)(x::AbstractVecOrMat) = a(x, x)
(a::Bilinear)(x::AbstractVector, y::AbstractVector) = vec(a(reshape(x, :,1), reshape(y, :,1)))
(a::Bilinear)(x::NTuple{2, AbstractArray}) = a(x[1], x[2])

function Base.show(io::IO, l::Bilinear)
  print(io, "Bilinear(", size(l.weight, 2), ", ", size(l.weight, 3), ", ", size(l.weight, 1))
  l.σ == identity || print(io, ", ", l.σ)
  l.bias == Flux.Zeros() && print(io, ", bias=false")
  print(io, ")")
end

"""
    Parallel(connection, layers...)

Create a 'Parallel' layer that passes an input array to each path in
`layers`, reducing the output with `connection`.

Called with one input `x`, this is equivalent to `reduce(connection, [l(x) for l in layers])`.
If called with multiple inputs, they are `zip`ped with the layers, thus `Parallel(+, f, g)(x, y) = f(x) + g(y)`.

# Examples

```jldoctest
julia> model = Chain(Dense(3, 5),
                     Parallel(vcat, Dense(5, 4), Chain(Dense(5, 7), Dense(7, 4))),
                     Dense(8, 17));

julia> size(model(rand(3)))
(17,)

julia> model = Parallel(+, Dense(10, 2), Dense(5, 2))
Parallel(+, Dense(10, 2), Dense(5, 2))

julia> size(model(rand(10), rand(5)))
(2,)
```
"""
struct Parallel{F, T}
  connection::F
  layers::T
end

Parallel(connection, layers...) = Parallel(connection, layers)

@functor Parallel

(m::Parallel)(x::AbstractArray) = mapreduce(f -> f(x), m.connection, m.layers)
(m::Parallel)(xs::Vararg{<:AbstractArray}) = mapreduce((f, x) -> f(x), m.connection, m.layers, xs)
(m::Parallel)(xs::Tuple) = m(xs...)

Base.getindex(m::Parallel, i::Integer) = m.layers[i]
Base.getindex(m::Parallel, i::AbstractVector) = Parallel(m.connection, m.layers[i]...)

function Base.show(io::IO, m::Parallel)
  print(io, "Parallel(", m.connection, ", ")
  join(io, m.layers, ", ")
  print(io, ")")
end


text:
"""
FunctionBlocks
"""

export QuantizerBlock

mutable struct QuantizerBlock <: AbstractFunctionBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    quantizationinterval::AbstractSignalLine

    function QuantizerBlock(::Type{Tv} = Float64; quantizationinterval = 0) where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}(),
            signal(quantizationinterval, Tv),
            ) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::QuantizerBlock)
    println(io, "Quantizer()")
end

"""
to expr
"""

function _toexpr(blk::QuantizerBlock)
    quantizationinterval = _toexpr(blk.quantizationinterval)
    inport = _toexpr(blk.inport[1])
    :($quantizationinterval * round($inport / $quantizationinterval))
end


text:
module PendulumWithSpring

using ModiaLang
using Unitful

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

simulationModel = nothing
get_simulationModel() = simulationModel

Pendulum = Model(
    m = 1.0,
    g = 9.81,

    world = Object3D(feature=Scene()),
    body  = Object3D(feature = Solid(massProperties=MassProperties(mass=:m, centerOfMass=[0.5, 0.0, 0.0]))),
    rev   = RevoluteWithFlange(obj1=:world, obj2=:body)
)

PendulumWithSpr = Model(
    pendulum = buildModia3D(Pendulum | Map(Lx=0.5, m=2, rev=Map(phi=Var(init=1.0)))),
    spring   = Spring | Map(c=100.0u"N*m/rad"),
	support  = Fixed,

    connect = :[ (pendulum.rev.flange, spring.flange_b)
                 (spring.flange_a, support.flange)
               ]
)

pendulumWithSpring = @instantiateModel(PendulumWithSpr, aliasReduction=false, unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 1.0
requiredFinalStates = [-0.1166492890518372, -15.325157691194002]
simulate!(pendulumWithSpring, stopTime=stopTime, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(pendulumWithSpring, ["pendulum.rev.flange.phi", "pendulum.rev.variables[1]"], figure=1)

end


text:
# 2017.07.21: Copied from the IJCNN codebase
#

module BayesChpOrig

# using Metadata

export bayes_chp_demo, bayes_detector, pccf_bayes_detector, extract_changes

function constant_hazard(r, lambda)
    p = 1.0 / lambda * ones(size(r))
end

"Non-standardized Student's t-distribution"
function studentpdf(x, mu, var, nu)
    ## As in the Matlab code by author:
    # c = exp(gammaln(nu/2 + 0.5) - gammaln(nu/2)) .* (nu.*pi.*var).^(-0.5);
    # p = c .* (1 + (1./(nu.*var)).*(x-mu).^2).^(-(nu+1)/2);
    c = exp.(lgamma.(nu/2 + 0.5) - lgamma.(nu/2)) .* (nu .* pi .* var).^(-0.5)
    p = c .* (1 + (1.0 ./ (nu .* var)) .* (x-mu).^2).^(-(nu+1)/2)
end

function extract_changes_from_output_naive(seq)
    chps = []
    for i in 2:length(seq)
        if seq[i] < seq[i-1]
            push!(chps, i)
        end
    end
    chps
end

"Function to extract changepoints locations from the bayes_detector() output"
function extract_changes(seq)
    #
    # seq[j] contains current estimate
    # of the location of the last change at the time moment j
    #
    n = length(seq)
    chps = zeros(n)
    for i in 1:n
        chps[i] = i - seq[i]
    end
    out = unique(chps)
    out = map(x->convert(Int,x), out)
end
function test_extract_changes()
    inp1 = [0,1,2,3,4,3]
    rt = extract_changes(inp1)
    print("... test changes extractor ", rt)
end

"""
Dynamic PCCF-based Bayesian detector.

h ~ 1/lambda => larger lambda correspondes to the smaller h

threshold : threshold for PCCF
h_min     : lambda value when PCCF <= threshold;
h_max     : lambda value when PCCF > threshold
h_common  : is either equal to h_min or h_max - the value to be used after max_d moment
max_d     : Max delay, after which we don't use PCCF's predictions

- V.0.1 (04 Nov 2016)
    Added max_d parameter to specify till which
    moment PCCF predictions will be used.

"""
function pccf_bayes_detector(in_sig, lambda_t, h_min, h_max, threshold; max_d = Nullable{Int64}(), h_common = Nullable{Float64}())
    @assert length(in_sig) == length(lambda_t)
    if ! isnull(h_common)
        @assert h_common == h_min || h_common == h_max
    end
    T       = length(in_sig)
    mu0     = 0.0
    kappa0  = 1
    alpha0  = 1
    beta0   = 1
    CP      = [0]
    R       = zeros(T+1, T+1)
    R[1, 1] = 1
    muT     = mu0
    kappaT  = kappa0
    alphaT  = alpha0
    betaT   = beta0
    maxes   = zeros(T+1)

    if isnull(max_d)
        max_d = T
    end

    for t=1:T
        predprobs = studentpdf(in_sig[t], muT, betaT .* (kappaT+1)./ (alphaT.*kappaT), 2*alphaT)

        if t <= max_d
            if lambda_t[t] <= threshold
                lambda = h_min
            else
                lambda = h_max
            end

        else
            lambda = h_common
        end

        H             = constant_hazard([1:t;], lambda)
        R[2:t+1, t+1] = R[1:t, t] .* predprobs .* (1-H)
        R[1, t+1]     = sum(R[1:t, t] .* predprobs .* H)
        R[:, t+1]     = R[:, t+1] / sum(R[:, t+1])
        muT0          = vcat(mu0   , (kappaT .* muT + in_sig[t]) ./ (kappaT + 1))
        kappaT0       = vcat(kappa0, kappaT + 1)
        alphaT0       = vcat(alpha0, alphaT + 0.5)
        betaT0        = vcat(beta0 , kappaT + (kappaT .* (in_sig[t] - muT).^2) ./ (2* (kappaT + 1)) )
        muT           = muT0
        kappaT        = kappaT0
        alphaT        = alphaT0
        betaT         = betaT0
        maxes[t]      = findmax(R[:,t])[2]
    end
    R = -log.(R)
    maxes = map(x -> convert(Int, x), maxes)
    maxes, R
end

#--------------------------------------------------

"Bayesian detector"
function bayes_detector(in_sig; lambda=200)
    T      = length(in_sig)
    # lambda = 200
    mu0    = 0.0
    kappa0 = 1
    alpha0 = 1
    beta0  = 1
    #X = zeros(T)
    CP = [0]
    R = zeros(T+1, T+1)
    R[1, 1] = 1
    #
    muT    = mu0
    kappaT = kappa0
    alphaT = alpha0
    betaT  = beta0
    #
    maxes = zeros(T+1)
    #
    for t=1:T
        # muT - is a vector of all possible hypos
        # We calc. predictive distribution for each hypothesis
        predprobs = studentpdf(in_sig[t], muT, betaT .* (kappaT+1)./ (alphaT.*kappaT), 2*alphaT)

        H = constant_hazard([1:t;], lambda)

        R[2:t+1, t+1] = R[1:t, t] .* predprobs .* (1-H)
        R[1, t+1] = sum(R[1:t, t] .* predprobs .* H)
        R[:, t+1] = R[:, t+1] / sum(R[:, t+1])

        ## Update statistics
        muT0    = vcat(mu0   , (kappaT .* muT + in_sig[t]) ./ (kappaT + 1))
        kappaT0 = vcat(kappa0, kappaT + 1)
        alphaT0 = vcat(alpha0, alphaT + 0.5)
        betaT0  = vcat(beta0 , kappaT + (kappaT .* (in_sig[t] - muT).^2) ./ (2* (kappaT + 1)) )

        muT    = muT0
        kappaT = kappaT0
        alphaT = alphaT0
        betaT  = betaT0

        maxes[t] = findmax(R[:,t])[2]
    end
    R = -log.(R)
    maxes = map(x->convert(Int,x), maxes)
    maxes, R
end

function bayes_chp_demo()
    sig    = readdlm("./data/sig.dat")
    T      = length(sig) # 1000
    lambda = 200
    mu0    = 0.0
    kappa0 = 1
    alpha0 = 1
    beta0  = 1
    #X = zeros(T)
    CP = [0]

    R = zeros(T+1, T+1)
    R[1, 1] = 1

    muT    = mu0
    kappaT = kappa0
    alphaT = alpha0
    betaT  = beta0

    maxes = zeros(T+1)

    for t=1:T
        # muT - is a vector of all possible hypos
        # We calc. predictive distribution for each hypothesis
        predprobs = studentpdf(sig[t], muT, betaT .* (kappaT+1)./ (alphaT.*kappaT), 2*alphaT)

        H = constant_hazard([1:t;], lambda)

        R[2:t+1, t+1] = R[1:t, t] .* predprobs .* (1-H)
        R[1, t+1] = sum(R[1:t, t] .* predprobs .* H)
        R[:, t+1] = R[:, t+1] / sum(R[:, t+1])

        ## Update statistics
        muT0    = vcat(mu0   , (kappaT .* muT + sig[t]) ./ (kappaT + 1))
        kappaT0 = vcat(kappa0, kappaT + 1)
        alphaT0 = vcat(alpha0, alphaT + 0.5)
        betaT0  = vcat(beta0 , kappaT + (kappaT .* (sig[t] - muT).^2) ./ (2* (kappaT + 1)) )

        muT    = muT0
        kappaT = kappaT0
        alphaT = alphaT0
        betaT  = betaT0

        maxes[t] = findmax(R[:,t])[2]
    end
    R = -log.(R)
    writedlm(outdir * "Rmatrix.txt", R, ",")
    writedlm(outdir * "Rmaxes.txt", maxes, ",")
    print(" write into ", outdir)
end

end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
struct GalileoE1BBase{C <: AbstractMatrix} <: AbstractGNSS{C}
    codes::C
end

"""
GalileoE1B is in theory CBOC. Here it is approximated as BOCcos.
"""
struct GalileoE1B{C <: AbstractMatrix} <: AbstractGNSSBOCcos{C, 1, 1}
    system::GalileoE1BBase{C}
end

get_system_string(s::GalileoE1B) = "GalileoE1B"

function read_from_documentation(raw_code)
    raw_code_without_spaces = replace(replace(raw_code, " " => ""), "\n" => "")
    code_hex_array = map(x -> parse(UInt16, x, base = 16), collect(raw_code_without_spaces))
    code_bit_string = string(string.(code_hex_array, base = 2, pad = 4)...)
    map(x -> parse(Int8, x, base = 2), collect(code_bit_string)) .* Int8(2) .- Int8(1)
end

function read_galileo_e1b_codes()
    read_in_codes(
        Int8,
        joinpath(dirname(pathof(GNSSSignals)), "..", "data", "codes_galileo_e1b.bin"),
        50,
        4092
    )
end

function GalileoE1BBase(;use_gpu = Val(false))
    _GalileoE1BBase(use_gpu)
end

function _GalileoE1BBase(::Val{false})
    GalileoE1BBase(Int16.(read_galileo_e1b_codes()))
end
function _GalileoE1BBase(::Val{true})
    GalileoE1BBase(CuMatrix{Float32}(read_galileo_e1b_codes()))
end

function GalileoE1B(;use_gpu = Val(false))
    GalileoE1B(GalileoE1BBase(use_gpu = use_gpu))
end

"""
$(SIGNATURES)

Get code length of GNSS system GalileoE1B.
```julia-repl
julia> get_code_length(GalileoE1B)
```
"""
@inline function get_code_length(galileo_e1b_base::GalileoE1BBase)
    4092
end

"""
$(SIGNATURES)

Get shortest code length of GNSS system GalileoE1B.
```julia-repl
julia> get_shortest_code_length(GalileoE1B)
```
"""
@inline function get_secondary_code_length(galileo_e1b_base::GalileoE1BBase)
    1
end

"""
$(SIGNATURES)

Get center frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_center_frequency(GalileoE1B)
```
"""
@inline function get_center_frequency(galileo_e1b_base::GalileoE1BBase)
    1_575_420_000Hz
end

"""
$(SIGNATURES)

Get code frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_code_frequency(GalileoE1B)
```
"""
@inline function get_code_frequency(galileo_e1b_base::GalileoE1BBase)
    1023_000Hz
end

"""
$(SIGNATURES)

Get data frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_data_frequency(GalileoE1B)
```
"""
@inline function get_data_frequency(galileo_e1b_base::GalileoE1BBase)
    250Hz
end


text:
using MKL
using Base.LinAlg.BLAS
using Gadfly
using DataFrames
using Colors

include("../utility.jl")

# BENCH PARAMETERS
path = Pkg.dir("MKL") * "/benchmark/b2/"

# Define

tup = (
    ("sqrt(X)",   "VJulia" ,    t,  sqrt        ,   (a,)         ),
    ("sqrt(X)",   "LJulia" ,    t,  lj_sqrt     ,   (n,a,1,b,1)  ),
    ("sqrt(X)",   "iLJulia",    t,  ilj_sqrt    ,   (a,b)        ),
    ("sqrt(X)",   "MJulia" ,    t,  map!        ,   (a,b)        ),
    ("sqrt(X)",   "BJulia" ,    t,  broadcast!  ,   (a,b)        ),
    ("sqrt(X)",   "BLAS"   ,    t,  nanf        ,   ()           ),
    ("sqrt(X)",   "MKL"    ,    t,  vj_dot      ,   (a,b)        ),
    )

# BENCH

bench_b2

function bench_sqrt(N::Vector{Int}, T::DataType, timing=0.01, gcdisable=true; α=0.05, print_header=true, print_eval=false, print_stats=false, print_detail=false)

    df = DataFrame( Generic = ASCIIString[],
                    N = Int[],
                    VJulia = Float64[],     # vectorized Julia
                    LJulia = Float64[],     # loop Julia
                    iLJulia = Float64[],    # @inbounds loop Julia
                    MJulia = Float64[],     # map! Julia
                    BJulia = Float64[],     # broadcast! Julia
                    MKL = Float64[])        # MKL

    for n in N

        a = rand(T,n)
        dest = rand(T,n)
        factor = 1e9/n

        t_vj = mean(ibench(sqrt,(a,),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])
        t_lj = mean(ibench(lj_sqrt,(dest,a),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])
        t_ilj = mean(ibench(ilj_sqrt,(dest,a),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])
        t_mj = mean(ibench(map!, (sqrt,dest,a),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])
        t_bj = mean(ibench(broadcast!, (sqrt,dest,a),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])

        t_mkl = mean(ibench(mkl_sqrt!, (n,a,dest),timing, gcdisable ; α=α, print_header= print_header, print_stats=print_stats, print_detail=print_detail)[1])

        push!(df,["sqrt(X)",n, t_vj*factor, t_lj*factor, t_ilj*factor, t_mj*factor, t_bj*factor, t_mkl]*factor)
        # push!(df,["sqrt(X)",n, t_vj/t_mkl, t_lj/t_mkl, t_ilj/t_mkl, t_mj/t_mkl, t_bj/t_mkl, t_mkl/t_mkl])
        # push!(res,[t_vj/t_mkl t_lj/t_mkl t_ilj/t_mkl t_mj/t_mkl t_bj/t_mkl t_mkl])
    end
    df
end

T = Float32
N = [
        10,20,30,40,50,100,250,500,750,1000,
        2500,5000,7500,10000,
        25000,50000,75000,100000,
        # 250_000,500_000,750_000,1_000_000,
        # 2500000,5000000,7500000,10000000,
        ]

using RDatasets
println(dataset("Zelig", "approval"))

df = bench_sin(N, T, 0.01, true)
println(df)
colors = colormap("Grays",10)[3:7]
colors[4]=colorant"red"
push!(colors,colorant"deepskyblue")

p = Gadfly.plot(
        layer(df, x="N", y="VJulia", Geom.line, Theme(default_color=colors[1])),
        layer(df, x="N", y="LJulia", Geom.line, Theme(default_color=colors[2])),
        layer(df, x="N", y="iLJulia", Geom.line, Theme(default_color=colors[3])),
        layer(df, x="N", y="MJulia", Geom.line, Theme(default_color=colors[4])),
        layer(df, x="N", y="BJulia", Geom.line, Theme(default_color=colors[5])),
        layer(df, x="N", y="MKL", Geom.line,Theme(default_color=colors[6])),
        Scale.x_log10(minvalue=10,maxvalue  =100_000),
        Scale.y_log10,
        Guide.xlabel("n-dim"),
        Guide.ylabel("relative CPU time"),
        Guide.title("relative CPU time against MKL ($T)"),
        # Guide.manual_color_key("Legend", ["Julia", "Julia", "MKL"], ["red", "gray", "deepskyblue"]),
        Guide.manual_color_key("Legend",
                ["Julia vectorized", "Julia loop", "Julia inbounds loop", "Julia map!", "Julia broadcast!", "MKL"],
                colors),
        )

filepath = string(path, "sqrt_$T.png")
println(filepath)
draw(PNG(filepath, 20cm, 10cm), p)
p


text:
import LLVM: refcheck

LLVM.@checked struct Logic
    ref::API.EnzymeLogicRef
    function Logic()
        ref = API.CreateLogic()
        new(ref)
    end
end
Base.unsafe_convert(::Type{API.EnzymeLogicRef}, logic::Logic) = logic.ref
LLVM.dispose(logic::Logic) = API.FreeLogic(logic)

# typedef bool (*CustomRuleType)(int /*direction*/, CTypeTree * /*return*/,
#                                CTypeTree * /*args*/, size_t /*numArgs*/,
#                                LLVMValueRef)=T


text:
using Transits
using Documenter

setup = quote
    using Transits
end

DocMeta.setdocmeta!(Transits, :DocTestSetup, setup; recursive = true)

makedocs(;
    modules=[Transits],
    authors="Miles Lucas <mdlucas@hawaii.edu> and contributors",
    repo="https://github.com/juliaastro/Transits.jl/blob/{commit}{path}#L{line}",
    sitename="Transits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaastro.github.io/Transits.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Introduction" => "introduction.md",
        "Getting Started" => "gettingstarted.md",
        "Benchmarks" => "bench.md",
        "API/Reference" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaAstro/Transits.jl",
    push_preview=true,
    devbranch="main"
)


text:
## import packages
include("load_packages.jl")

experiment = "3_calculation_time_company"
capacity_steps      = 5
min_cap             = 10000
max_cap             = 50000
warehouse_versions  = 1
buffer_addition     = 1
buffer_size         = 0.2
capacity_diffsteps  = 3
capacity_diffmax    = 0.3

total_rows = capacity_steps * warehouse_versions * (capacity_diffsteps) * (buffer_addition + 1)
total_cols = warehouse_versions + 1
base_cap   = min_cap:(max_cap-min_cap)/(capacity_steps-1):max_cap
base_ware  = 2:warehouse_versions+1

benchmark = Matrix{Int64}(undef,total_rows,total_cols) .= 0
sku_list  = Vector{Int64}(undef,total_rows) .= 0

current_row = 1
for buffstep = 1:buffer_addition + 1
    for cstep in base_cap
        if buffstep == 1
            capacity = cstep
        else
            capacity = cstep * (1+buffer_size)
        end
        for warestep in base_ware
            for capvarstep = 1:capacity_diffsteps
                vararray = Vector{Float64}(undef,warestep) .= 0
                if capvarstep == 1
                    vararray .= 1/length(vararray)
                else
                    capacity_diff = capacity_diffmax/(capacity_diffsteps-capvarstep+1)
                    step = 1
                    for i = (1+capacity_diff):-(capacity_diff*2)/(length(vararray)-1):(1-capacity_diff)
                        vararray[step] = i/length(vararray)
                        step += 1
                    end
                end
                for warehouse = 1:warestep
                    if warehouse == 1
                        benchmark[current_row,warehouse] = ceil(vararray[warehouse] * capacity)
                    elseif  warehouse != warestep
                        benchmark[current_row,warehouse] = round(vararray[warehouse] * capacity)
                    else
                        benchmark[current_row,warehouse] = capacity - sum(benchmark[current_row,:])
                    end
                end
                sku_list[current_row] = sum(benchmark[current_row,:])
                global current_row += 1
            end
        end
    end
end
writedlm("capacity/capacity_$experiment.csv", benchmark, ';')
writedlm("capacity/skus_$experiment.csv", sku_list)

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
immutable RateParameter
end

@withkw type Rate
    param::RateParameter = RateParameter()
    N::Int = 100
    x::Vector{Float} = 0.5randn(N)
    r::Vector{Float} = tanh(x)
    g::Vector{Float} = zeros(N)
    I::Vector{Float} = zeros(N)
    records::Dict = Dict()
end

@replace function integrate!(p::Rate, param::RateParameter, dt::Float)
    @inbounds for i = 1:N
        x[i] += dt * (-x[i] + g[i] + I[i])
        r[i] = tanh(x[i]) #max(0, x[i])
    end
end


text:
############################################################
########### Incremental Pruning Solver #####################
############################################################

"""
    PruneSolver <: Solver

POMDP solver type using the incremental pruning algorithm.
"""
mutable struct PruneSolver <: Solver
    max_iterations::Int64
    tolerance::Float64
    optimizer_factory::OptimizerFactory
end

"""
    PruneSolver(; max_iterations, tolerance)

Initialize an incremental pruning solver with the `max_iterations` limit and desired `tolerance`.
"""
function PruneSolver(;max_iterations::Int64=10, tolerance::Float64=1e-3, optimizer_factory::OptimizerFactory=JuMP.with_optimizer(GLPK.Optimizer))
    return PruneSolver(max_iterations, tolerance, optimizer_factory)
end

"""
    AlphaVec

Alpha vector type of paired vector and action.
"""
struct AlphaVec
    alpha::Vector{Float64} # alpha vector
    action::Any # action associated wtih alpha vector
end

"""
    AlphaVec(vector, action_index)

Create alpha vector from `vector` and `action_index`.
"""
AlphaVec() = AlphaVec([0.0], 0)

# define alpha vector equality
==(a::AlphaVec, b::AlphaVec) = (a.alpha,a.action) == (b.alpha, b.action)
Base.hash(a::AlphaVec, h::UInt) = hash(a.alpha, hash(a.action, h))

"""
    create_policy(prune_solver, pomdp)

Create AlphaVectorPolicy for `prune_solver` using immediate rewards from `pomdp`.
"""
function create_policy(solver::PruneSolver, pomdp::POMDP)
    ns = length(states(pomdp))
    na = length(actions(pomdp))
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    sa_reward = LazyCachedSAR(pomdp)
    alphas = [[sa_reward(S[i],A[j]) for i in 1:ns] for j in 1:na]
    AlphaVectorPolicy(pomdp, alphas, A)
end

"""
    xsum(A,B)

Compute the cross-sum of alpha vectors `A` and `B`.

# Examples
```julia-repl
julia> a, b, c = AlphaVec([1.0, 2.0], 1), AlphaVec([3.0, 4.0], 1), AlphaVec([5.0, 6.0], 1)
julia> xsum(Set([a]), Set([b, c]))
Set(AlphaVec[AlphaVec([6.0, 8.0], 1), AlphaVec([4.0, 6.0], 1)])
```
"""
function xsum(A::Set{AlphaVec}, B::Set{AlphaVec})
    X = Set{AlphaVec}()
    for a in A, b in B
        @assert a.action == b.action "action mismatch"
        newVec = AlphaVec(a.alpha + b.alpha, a.action)
        push!(X, newVec)
    end
    X
end

"""
    xsum(A,B)

Compute the cross-sum of vectors `A` and `B`.

# Examples
```julia-repl
julia> xsum(Set([[1.0, 2.0]]), Set([[3.0, 4.0], [5.0, 6.0]]))
Set(Array{Float64,1}[[4.0, 6.0], [6.0, 8.0]])
```
"""
function xsum(A::Set{Array{Float64,1}}, B::Set{Array{Float64,1}})
    X = Set{Array{Float64,1}}()
    for a in A, b in B
        push!(X, a + b)
    end
    X
end

"""
    dominate(α, A)

The set of vectors in `A` dominated by `α`.
"""
function dominate(α::Array{Float64,1}, A::Set{Array{Float64,1}}, optimizer_factory::OptimizerFactory)
    ns = length(α)
    αset = Set{Array{Float64,1}}()
    push!(αset, α)
    Adiff = setdiff(A,αset)
    L = Model(optimizer_factory)
    @variable(L, x[1:ns])
    @variable(L, δ)
    @objective(L, Max, δ)
    @constraint(L, sum(x) == 1)
    @constraint(L, x[1:ns] .<= 1)
    @constraint(L, x[1:ns] .>= 0)
    # dx = size(x)
    # da = size(α)
    # println("Size of x: $dx")
    # println("Size of α: $da")
    for ap in Adiff
        # dap = size(ap)
        # println("Size of ap: $dap")
        @constraint(L, dot(x, α) >= δ + dot(x, ap))
    end
    JuMP.optimize!(L)
    sol_status = JuMP.termination_status(L)
    if sol_status == :Infeasible
        return :Perp
        # println("Perp")
    else
        xval = JuMP.value.(x)
        dval = JuMP.value.(δ)
        if dval > 0
            return xval
            # println("xval")
        else
            return :Perp
            # println("Perp")
        end
    end
end # dominate
@deprecate dominate(α::Array{Float64,1}, A::Set{Array{Float64,1}}) dominate(α, A, JuMP.with_optimizer(GLPK.Optimizer))

"""
    filtervec(F)

The set of vectors in `F` that contribute to the value function.
"""
function filtervec(F::Set{Array{Float64,1}}, optimizer_factory::OptimizerFactory)
    ns = length(sum(F))
    W = Set{Array{Float64,1}}()
    for i = 1: ns
        if !isempty(F)
            # println("i: $i  ")
            w = Array{Float64,1}()
            fsmax = -Inf
            for f in F
                # println("f: $f")
                if f[i] > fsmax
                    fsmax = f[i]
                    w = f
                    end
            end
            wset = Set{Array{Float64,1}}()
            push!(wset, w)
            # println("w: $w")
            push!(W,w)
            setdiff!(F,wset)
        end
    end
    while !isempty(F)
        # println("F Before: $F")
        ϕ = pop!(F)
        # println("F After: $F")
        # println("ϕ: $ϕ")
        # println("W: $W")
        x = dominate(ϕ, W, optimizer_factory)
        if x != :Perp
            push!(F, ϕ)
            w = Array{Float64,1}()
            fsmax = -Inf
            for f in F
                if dot(x, f) > fsmax
                    fsmax = dot(x, f)
                    w = f
                    end
            end
            wset = Set{Array{Float64,1}}()
            push!(wset, w)
            push!(W,w)
            setdiff!(F,wset)
        end
    end
    temp = [Float64[]]
    setdiff!(W,temp)
    W
end # filtervec
@deprecate filtervec(F::Set{Array{Float64,1}}) filtervec(F, JuMP.with_optimizer(GLPK.Optimizer))

"""
    filtervec(F)

The set of alpha vectors in `F` that contribute to the value function.
"""
function filtervec(A::Set{AlphaVec}, optimizer_factory::OptimizerFactory)
    ns = [length(av.alpha) for av in A][1]
    W = Set{AlphaVec}()
    for i = 1:ns
        if !isempty(A)
            # println("i: $i  ")
            w = AlphaVec()
            fsmax = -Inf
            for f in A
                # println("f: $f")
                if f.alpha[i] > fsmax
                    fsmax = f.alpha[i]
                    w = f
                    end
            end
            wset = Set{AlphaVec}()
            push!(wset, w)
            # println("w: $w")
            push!(W,w)
            setdiff!(A,wset)
        end
    end
    while !isempty(A)
        # println("F Before: $F")
        ϕ = pop!(A)
        ϕa = ϕ.alpha
        Wa = Set([w.alpha for w in W])
        # println("F After: $F")
        # println("ϕ: $ϕ")
        # println("W: $W")
        x = dominate(ϕa, Wa, optimizer_factory)
        if x != :Perp
            push!(A, ϕ)
            w = AlphaVec()
            fsmax = -Inf
            for f in A
                if dot(x, f.alpha) > fsmax
                    fsmax = dot(x, f.alpha)
                    w = f
                    end
            end
            wset = Set{AlphaVec}()
            push!(wset, w)
            push!(W,w)
            setdiff!(A,wset)
        end
    end
    W
end # filtervec
@deprecate filtervec(A::Set{AlphaVec}) filtervec(A, JuMP.with_optimizer(GLPK.Optimizer))

"""
    incprune(Sz)

Standard incremental pruning of the set of vectors `Sz`.
"""
function incprune(SZ::Vector{Set{Vector{Float64}}}, optimizer_factory::OptimizerFactory)
    W = filtervec(xsum(SZ[1], SZ[2]), optimizer_factory)
    for i = 3:length(SZ)
        W = filtervec(xsum(W, SZ[i]), optimizer_factory)
    end
    W
end # incprune
@deprecate incprune(SZ::Vector{Set{Vector{Float64}}}) incprune(SZ, JuMP.with_optimizer(GLPK.Optimizer))

"""
    dpval(α, a, z, pomdp, sa_reward=LazyCachedSAR(pomdp))

Dynamic programming backup value of `α` for action `a` and observation `z` in `pomdp`.
"""
function dpval(α::Array{Float64,1}, a, z, prob::POMDP, sa_reward::StateActionReward=LazyCachedSAR(prob))
    S = ordered_states(prob)
    A = ordered_actions(prob)
    ns = length(states(prob))
    nz = length(observations(prob))
    γ = discount(prob)
    τ = Array{Float64,1}(undef, ns)
    for (sind,s) in enumerate(S)
        dist_t = transition(prob,s,a)
        exp_sum = 0.0
        for (spind, sp) in enumerate(S)
            dist_o = observation(prob, s, a, sp)
            pt = pdf(dist_t,sp)
            po = pdf(dist_o,z)
            exp_sum += α[spind] * po * pt
        end
        τ[sind] = (1 / nz) * sa_reward(s,a) + γ * exp_sum
    end
    τ
end

"""
    dpupdate(F, pomdp, sa_reward, optimizer_factory)

Dynamic programming update of `pomdp` for the set of alpha vectors `F`.
"""
function dpupdate(F::Set{AlphaVec}, prob::POMDP, sa_reward::StateActionReward, optimizer_factory::OptimizerFactory)
    alphas = [avec.alpha for avec in F]
    A = ordered_actions(prob)
    Z = ordered_observations(prob)
    na = length(actions(prob))
    nz = length(observations(prob))
    Sp = Set{AlphaVec}()
    # tcount = 0
    Sa = Set{AlphaVec}()
    for (aind, a) in enumerate(A)
        Sz = Vector{Set{Vector{Float64}}}(undef, nz)
        for (zind, z) in enumerate(Z)
            # tcount += 1
            # println("DP Update Inner Loop: $tcount")
            V = Set(dpval(α, a, z, prob, sa_reward) for α in alphas)
            # println("V: $V")
            Sz[zind] = filtervec(V, optimizer_factory)
        end
        Sa = Set([AlphaVec(α,a) for α in incprune(Sz, optimizer_factory)])
        union!(Sp,Sa)
    end
    filtervec(Sp, optimizer_factory)
end
@deprecate dpupdate(F::Set{AlphaVec}, prob::POMDP) dpupdate(F, prob, LazyCachedSAR(prob), JuMP.with_optimizer(GLPK.Optimizer))

"""
    diffvalue(Vnew, Vold, pomdp, sa_reward, optimizer_factory)

Maximum difference between new alpha vectors `Vnew` and old alpha vectors `Vold` in `pomdp`.
"""
function diffvalue(Vnew::Vector{AlphaVec}, Vold::Vector{AlphaVec}, pomdp::POMDP, sa_reward::StateActionReward, optimizer_factory::OptimizerFactory)
    ns = length(states(pomdp)) # number of states in alpha vector
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    Anew = [avec.alpha for avec in Vnew]
    Aold = [avec.alpha for avec in Vold]
    dmax = -Inf # max difference
    for avecnew in Anew
        L = Model(optimizer_factory)
        @variable(L, x[1:ns])
        @variable(L, t)
        @objective(L, Max, t)
        @constraint(L, x .>= 0)
        @constraint(L, x .<= 1)
        @constraint(L, sum(x) == 1)
        for avecold in Aold
            @constraint(L, (avecnew - avecold)' * x >= t)
        end
        JuMP.optimize!(L)
        dmax = max(dmax, JuMP.objective_value(L))
    end
    rmin = minimum(sa_reward(s, a) for s in S, a in A) # minimum reward
    if rmin < 0 # if negative rewards, find max difference from old to new
        for avecold in Aold
            L = Model(optimizer_factory)
            @variable(L, x[1:ns])
            @variable(L, t)
            @objective(L, Max, t)
            @constraint(L, x .>= 0)
            @constraint(L, x .<= 1)
            @constraint(L, sum(x) == 1)
            for avecnew in Anew
                @constraint(L, (avecold - avecnew)' * x >= t)
            end
            JuMP.optimize!(L)
            dmax = max(dmax, JuMP.objective_value(L))
        end
    end
    dmax
end
@deprecate diffvalue(Vnew::Vector{AlphaVec},Vold::Vector{AlphaVec},pomdp::POMDP) diffvalue(Vnew, Vold, pomdp, LazyCachedSAR(pomdp), JuMP.with_optimizer(GLPK.Optimizer))

"""
    solve(solver::PruneSolver, pomdp)

AlphaVectorPolicy for `pomdp` caluclated by the incremental pruning algorithm.
"""
function solve(solver::PruneSolver, prob::POMDP)
    # println("Solver started...")
    ϵ = solver.tolerance
    replimit = solver.max_iterations
    policy = create_policy(solver, prob)
    avecs = [AlphaVec(policy.alphas[i], policy.action_map[i]) for i in 1:length(policy.action_map)]
    Vold = Set(avecs)
    Vnew = Set{AlphaVec}()
    del = Inf
    reps = 0
    sa_reward = LazyCachedSAR(prob)
    while del > ϵ && reps < replimit
        reps += 1
        Vnew = dpupdate(Vold, prob, sa_reward, solver.optimizer_factory)
        del = diffvalue(collect(Vnew), collect(Vold), prob, sa_reward, solver.optimizer_factory)
        Vold = Vnew
    end
    alphas_new = [v.alpha for v in Vnew]
    actions_new = [v.action for v in Vnew]
    policy = AlphaVectorPolicy(prob, alphas_new, actions_new)
    return policy
end

@POMDP_require solve(solver::PruneSolver, pomdp::POMDP) begin
    P = typeof(pomdp)
    S = state_type(P)
    A = action_type(P)
    @req discount(::P)
    @subreq ordered_states(pomdp)
    @subreq ordered_actions(pomdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req state_index(::P,::S)
    as = actions(pomdp)
    @req length(::typeof(as))
    ss = states(pomdp)
    @req length(::typeof(ss))
    @req iterator(::typeof(as))
    @req iterator(::typeof(ss))
    s = first(iterator(ss))
    a = first(iterator(as))
    dist = transition(pomdp, s, a)
    D = typeof(dist)
    @req iterator(::D)
    @req pdf(::D,::S)
end


text:
using ITensors
using ChainRulesCore: ZeroTangent

# Write a custom `rrule` for this.
function inner_circuit(ϕ::ITensor, U::Vector{ITensor}, ψ::ITensor)
  Uψ = ψ
  for u in U
    s = commoninds(u, Uψ)
    s′ = s'
    Uψ = replaceinds(u * Uψ, s′ => s)
  end
  return (ϕ * Uψ)[]
end

name(g::Tuple{String,Vararg}) = g[1]
gate_sites(g::Tuple{<:Any,Tuple{Vararg{Int}},Vararg}) = g[2]
params(g::Tuple{<:Any,<:Any,<:NamedTuple}) = g[3]
params(g::Tuple{<:Any,<:Any}) = NamedTuple()

function gate(g::Tuple, s::Vector{<:Index})
  return gate(name(g), params(g), s[collect(gate_sites(g))])
end

function gate(gn::String, params::NamedTuple, s::Vector{<:Index})
  return gate(gate(gn, params), s)
end

function gate(gn, params::NamedTuple)
  return gate(gn; params...)
end

function gate(gn::String, params::NamedTuple)
  return gate(Val{Symbol(gn)}(), params)
end

function gate(g::Matrix, s::Vector{<:Index})
  s = reverse(s)
  return itensor(g, s'..., dag(s)...)
end

function buildcircuit(gates::Vector, s::Vector{<:Index})
  return [gate(g, s) for g in gates]
end

# Gate definitions
function gate(gn::Val; params...)
  return error("Gate $gn not defined")
end

function gate(::Val{:Ry}; θ)
  return [
    cos(θ / 2) -sin(θ / 2)
    sin(θ / 2) cos(θ / 2)
  ]
end

function gate(::Val{:CX})
  return [
    1 0 0 0
    0 1 0 0
    0 0 0 1
    0 0 1 0
  ]
end

# XXX: For some reason Zygote needs these definitions?
Base.reverse(z::ZeroTangent) = z
Base.adjoint(::Tuple{Nothing}) = nothing
Base.adjoint(::Tuple{Nothing,Nothing}) = nothing


text:
export ToeplitzOperator, HankelOperator

mutable struct ToeplitzOperator{T<:Number} <: Operator{T}
    negative::Vector{T}
    nonnegative::Vector{T}
end

ToeplitzOperator(V::Vector{T},W::Vector{Q}) where {T<:Number,Q<:Number} =
    ToeplitzOperator{promote_type(T,Q)}(V,W)
ToeplitzOperator(V::AbstractVector,W::AbstractVector) =
    ToeplitzOperator(collect(V),collect(W))

convert(::Type{Operator{TT}},T::ToeplitzOperator) where {TT} =
    ToeplitzOperator(convert(Vector{TT},T.negative),convert(Vector{TT},T.nonnegative))

for op in (:(Base.real), :(Base.imag))
    @eval $op(T::ToeplitzOperator) = ToeplitzOperator($op(T.negative), $op(T.nonnegative))
end

function SymToeplitzOperator(V::Vector)
    W=V[2:end]
    V=copy(V)
    V[1]*=2
    ToeplitzOperator(W,V)
end

for OP in (:domainspace,:rangespace)
    @eval $OP(T::ToeplitzOperator) = ℓ⁰
end

getindex(T::ToeplitzOperator,k::Integer,j::Integer) =
    toeplitz_getindex(T.negative,T.nonnegative,k,j)

function toeplitz_getindex(negative::AbstractVector{T},nonnegative::AbstractVector{T},k::Integer,j::Integer) where T
    if 0<k-j≤length(negative)
        negative[k-j]
    elseif 0≤j-k≤length(nonnegative)-1
        nonnegative[j-k+1]
    else
        zero(T)
    end
end

function toeplitz_getindex(cfs::AbstractVector{T},k::Integer,j::Integer) where T
    if k==j && !isempty(cfs)
        2cfs[1]
    elseif 0<k-j≤length(cfs)-1
        cfs[k-j+1]
    elseif 0<j-k≤length(cfs)-1
        cfs[j-k+1]
    else
        zero(T)
    end
end

function BandedMatrix(S::SubOperator{T,ToeplitzOperator{T},Tuple{UnitRange{Int},UnitRange{Int}}}) where T
    ret = BandedMatrix(Zeros, S)

    kr,jr=parentindices(S)

    neg=parent(S).negative
    pos=parent(S).nonnegative

    toeplitz_axpy!(1.0,neg,pos,kr,jr,ret)
end

bandwidths(T::ToeplitzOperator)=(length(T.negative),length(T.nonnegative)-1)

# slice of a ToeplitzOPerator is a ToeplitzOperator

function Base.getindex(T::ToeplitzOperator,kr::InfRanges,jr::InfRanges)
    sh=first(jr)-first(kr)
    st=step(jr)
    @assert st==step(kr)
    if sh ≥0
        ToeplitzOperator([reverse!(T.nonnegative[1:sh]);T.negative],T.nonnegative[sh+1:st:end])
    else
        ToeplitzOperator(T.negative[-sh+1:st:end],[reverse!(T.negative[1:-sh]);T.nonnegative])
    end
end

## Hankel Operator

mutable struct HankelOperator{T<:Number} <: Operator{T}
    coefficients::Vector{T}
end

for OP in (:domainspace,:rangespace)
    @eval $OP(T::HankelOperator) = ℓ⁰
end

HankelOperator(V::AbstractVector)=HankelOperator(collect(V))

HankelOperator(f::Fun)=HankelOperator(f.coefficients)

@eval convert(::Type{Operator{TT}},T::HankelOperator) where {TT}=HankelOperator(convert(Vector{TT},T.coefficients))

function hankel_getindex(v::AbstractVector,k::Integer,j::Integer)
   if k+j-1 ≤ length(v)
        v[k+j-1]
    else
        zero(eltype(v))
    end
end

getindex(T::HankelOperator,k::Integer,j::Integer) =
    hankel_getindex(T.coefficients,k,j)

function BandedMatrix(S::SubOperator{T,HankelOperator{T},Tuple{UnitRange{Int},UnitRange{Int}}}) where T
    ret=BandedMatrix(Zeros, S)

    kr,jr=parentindices(S)
    cfs=parent(S).coefficients

    hankel_axpy!(1.0,cfs,kr,jr,ret)
end

bandwidths(T::HankelOperator) = (max(0,length(T.coefficients)-1),max(0,length(T.coefficients)-1))

## algebra

function Base.maximum(T::ToeplitzOperator)
    if isempty(T.negative)
        maximum(T.nonnegative)
    elseif isempty(T.nonnegative)
        maximum(T.negative)
    else
        max(maximum(T.negative),maximum(T.nonnegative))
    end
end

-(T::ToeplitzOperator)=ToeplitzOperator(-T.negative,-T.nonnegative)
*(c::Number,T::ToeplitzOperator)=ToeplitzOperator(c*T.negative,c*T.nonnegative)

-(H::HankelOperator)=HankelOperator(-H.coefficients)
*(c::Number,H::HankelOperator)=HankelOperator(c*H.coefficients)

## inv

function Base.inv(T::ToeplitzOperator)
    @assert length(T.nonnegative)==1
    ai=\(T,[1.0];maxlength=100000)
    ToeplitzOperator(ai[2:end],ai[1:1])
end

function Fun(T::ToeplitzOperator)
   if length(T.nonnegative)==1
      Fun(Taylor(),[T.nonnegative;T.negative])
    elseif length(T.negative)==0
        Fun(Hardy{false}(),T.nonnegative)
    else
        Fun(Laurent(Circle()),interlace(T.nonnegative,T.negative))
    end
end


text: run(`../build.sh`)


text:
# Julia wrapper for header: /usr/include/scip/pub_relax.h
# Automatically generated using Clang.jl wrap_c

function SCIPrelaxComp(elem1, elem2)
    ccall((:SCIPrelaxComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPrelaxCompName(elem1, elem2)
    ccall((:SCIPrelaxCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPrelaxGetData(relax)
    ccall((:SCIPrelaxGetData, libscip), Ptr{SCIP_RELAXDATA}, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxSetData(relax, relaxdata)
    ccall((:SCIPrelaxSetData, libscip), Cvoid, (Ptr{SCIP_RELAX}, Ptr{SCIP_RELAXDATA}), relax, relaxdata)
end

function SCIPrelaxGetName(relax)
    ccall((:SCIPrelaxGetName, libscip), Cstring, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetDesc(relax)
    ccall((:SCIPrelaxGetDesc, libscip), Cstring, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetPriority(relax)
    ccall((:SCIPrelaxGetPriority, libscip), Cint, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetFreq(relax)
    ccall((:SCIPrelaxGetFreq, libscip), Cint, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetSetupTime(relax)
    ccall((:SCIPrelaxGetSetupTime, libscip), Cdouble, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetTime(relax)
    ccall((:SCIPrelaxGetTime, libscip), Cdouble, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxGetNCalls(relax)
    ccall((:SCIPrelaxGetNCalls, libscip), Clonglong, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxIsInitialized(relax)
    ccall((:SCIPrelaxIsInitialized, libscip), UInt32, (Ptr{SCIP_RELAX},), relax)
end

function SCIPrelaxMarkUnsolved(relax)
    ccall((:SCIPrelaxMarkUnsolved, libscip), Cvoid, (Ptr{SCIP_RELAX},), relax)
end


text:
"""
    XDR (machine-independent binary) RData format IO stream wrapper.
"""
type XDRIO{T<:IO} <: RDAIO
    sub::T             # underlying IO stream
    buf::Vector{UInt8} # buffer for strings

    XDRIO(io::T) = new(io, Array(UInt8, 1024))
    @compat (::Type{XDRIO}){T <: IO}(io::T) = new{T}(io, Array(UInt8, 1024))
end

readint32(io::XDRIO) = ntoh(read(io.sub, Int32))
readuint32(io::XDRIO) = ntoh(read(io.sub, UInt32))
readfloat64(io::XDRIO) = ntoh(read(io.sub, Float64))

readintorNA(io::XDRIO) = readint32(io)
readintorNA(io::XDRIO, n::RVecLength) = map!(ntoh, read(io.sub, Int32, n))

# this method have Win32 ABI issues, see JuliaStats/RData.jl#5
# R's NA is silently converted to NaN when the value is loaded in the register(?)
#readfloatorNA(io::XDRIO) = readfloat64(io)

readfloatorNA(io::XDRIO, n::RVecLength) = reinterpret(Float64, map!(ntoh, read(io.sub, UInt64, n)))

readuint8(io::XDRIO, n::RVecLength) = readbytes(io.sub, n)

function readnchars(io::XDRIO, n::Int32)  # a single character string
    readbytes!(io.sub, io.buf, n)
    convert(RString, unsafe_string(pointer(io.buf), n))
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function format(obj::VectorFormat, arg0::JVector)
    return jcall(obj, "format", JString, (JVector,), arg0)
end

function format(obj::VectorFormat, arg0::JVector, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (JVector, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_available_locales(::Type{VectorFormat})
    return jcall(VectorFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_format(obj::VectorFormat)
    return jcall(obj, "getFormat", NumberFormat, ())
end

function get_prefix(obj::VectorFormat)
    return jcall(obj, "getPrefix", JString, ())
end

function get_separator(obj::VectorFormat)
    return jcall(obj, "getSeparator", JString, ())
end

function get_suffix(obj::VectorFormat)
    return jcall(obj, "getSuffix", JString, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function parse(obj::VectorFormat, arg0::JString)
    return jcall(obj, "parse", JVector, (JString,), arg0)
end

function parse(obj::VectorFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", JVector, (JString, ParsePosition), arg0, arg1)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
# Test that TimeZones.jl can be safely used in a multithreaded environment.
# Note that the number of threads being used cannot be changed dynamically, so
# this test file spawns a new julia process running with multiple threads.

using Test

const program = """
using TimeZones
using Test

@assert Threads.nthreads() > 1 "This system does not support multiple threads, so the thread-safety tests cannot be run."

@testset "Multithreaded TimeZone brute force test" begin
    function create_zdt(year, month, day, tz_name)
        ZonedDateTime(DateTime(year, month, day), TimeZone(tz_name))
    end
    function cycle_zdts()
        return [
            try
                create_zdt(year, month, day, tz_name)
            catch e
                # Ignore ZonedDateTimes that aren't valid
                e isa Union{ArgumentError,AmbiguousTimeError,NonExistentTimeError} || rethrow()
                nothing
            end
            for year in 2000:2020
            for month in 1:5
            for day in 10:15
            for tz_name in timezone_names()
        ]
    end

    outputs = Channel(Inf)
    @sync begin
        for _ in 1:15
            Threads.@spawn begin
                put!(outputs, cycle_zdts())
            end
        end
    end
    close(outputs)

    tzs = collect(outputs)

    # Test that every Task produced the same result
    allsame(x) = all(y -> y == first(x), x)
    @test allsame(tzs)
end

#----------------------------------------------------

@testset "Interleaved compile() and TimeZone construction" begin
    @sync for i in 1:20
        if (i % 5 == 0)
            TimeZones.TZData.compile()
        end
        Threads.@spawn begin
            TimeZone("US/Eastern", TimeZones.Class(:LEGACY))
        end
    end
end
"""

@info "Running Thread Safety tests"
@testset "Multithreaded TimeZone construction" begin
    withenv("JULIA_NUM_THREADS" => "8") do
        run(`$(Base.julia_cmd()) --proj -E $(program)`)
    end
end


text:
using DataStructures
using CSV

"Parase genus assignment from a dada results file"
function parse_dada_results(fd::String)
    centroid_genus_map = Dict{String,String}()
    ct = 0
    nact = 0
    totalnr = 0
    genuses = Set{String}()
    genuses_counts = Dict{String,Int64}()
    for l in readlines(fd)
        ct += 1
        if ct > 1
            parts = split(replace(l, "\""=>""), "\t")
            rn = split(parts[1], ";")[1]
            size = parse(Int64, split(split(parts[1], ";")[2], "=")[2])
            totalnr += size
            gn = replace(parts[7], " "=>".")
            gn = replace(gn, "/"=>"-")
            if gn == "NA"
                nact += 1
                gn = "$gn"*"_"*"$nact"
            else
                if gn in genuses
                    genuses_counts[gn] += 1
                else
                    push!(genuses,gn)
                    genuses_counts[gn] = 1
                end
                gnnr = genuses_counts[gn]
                gn = "$gn"*"_"*"$gnnr"
            end
            centroid_genus_map[rn] = gn
        end
    end
    return(centroid_genus_map)
end

function fasta_to_dict(sfasta::String;remove_abundance=false)::OrderedDict{String,LongSequence{BioSequences.DNAAlphabet{4}}}
    out = OrderedDict{String,LongSequence{BioSequences.DNAAlphabet{4}}}()
    record = FASTA.Record()
    reader = FASTA.Reader(open(sfasta, "r"))
    while !eof(reader)
        read!(reader, record)
        id = identifier(record)
        if remove_abundance
            id = split(id, ";")[1]
        end
        seq = sequence(record)
        out[id]=seq
    end
    close(reader)
    return(out)
end

"Read salmon output"
function read_salmon(inf::String;remove_abundance=true) #by defaulyt remove ;size= fragment
    din = DataFrame(CSV.File(inf))
    if remove_abundance
        din.Name = map(x->split(x, ";")[1], din.Name)
    end
    out = Dict(zip(din.Name, din.TPM))
    return(out)
end

function dict_to_fasta(inp::OrderedDict{String,BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}, sfasta::String)
    writefasta(sfasta, inp)
end

#test aligned part percentage identity and length fraction compared to the s2 sequence
function check_pair(s1::BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}, s2::BioSequences.LongSequence{BioSequences.DNAAlphabet{4}})
    scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=-1)
    res = pairalign(LocalAlignment(), s1, s2, scoremodel)
    aln = alignment(res)
    pairs = collect(aln)
    l = count_aligned(aln)
    match = count_matches(aln)
    mismatch = BioAlignments.count_mismatches(aln)
    identityf = round(match/l, digits=3)
    coveragef = Base.maximum([(match + mismatch)/length(s1), ((match + mismatch)/length(s2))])
    coveragef = round(coveragef, digits=3)
    return(identityf, coveragef)
end

#read mapping of gjc file - "tracking of clusterings file"
function parse_mapping(mapfile::String)
    pairs = Array{Tuple{String,String},1}()
    r2c = Dict{String,String}() #read to cluster representative
    c2r = Dict{String,Array{String,1}}() #cluster rep to reads map
    for l in readlines(mapfile)
        parts = split(l)
        push!(pairs, (parts[1], parts[2]))
        r2c[String(parts[1])] = String(parts[2])
    end
    sort!(pairs, lt=(x,y)->isless(x[2], y[2]))
    cur_r2 = ""
    len = length(pairs)
    ct = 0
    cl_ct = 0
    r1s = Array{String,1}()
    cluster_4_work = Array{Tuple{Array{String,1},String},1}()
    @showprogress 1 "Parsing mapping..."  for pair in pairs
        ct += 1
        if ct == 1
            cur_r2 = pair[2]
        end
        if pair[2] != cur_r2 && ct !=1
            cl_ct += 1
            push!(cluster_4_work, (r1s, cur_r2))
            r1s = Array{String,1}()
            cur_r2 = pair[2]
        end
        c2r[cur_r2] = r1s
        push!(r1s, pair[1])
    end
    push!(cluster_4_work, (r1s,cur_r2))
    c2r[cur_r2] = r1s
    return(cluster_4_work, r2c, c2r)
end


text:
using CSV, DataFrames, Dates, LinearAlgebra, Distributions, DelimitedFiles, SharedArrays

import Base.Threads.@spawn
# Custom package
using Jevo

# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

gap = 10
l_0 = 20
n = 4
N = 100
fl = 0.6l_0/2N
f0 = 200l_0/2N
κ_arr = [0, 1, 2, 3, 4, 5, 10, 20, 40]
F = Jevo.num_fermi(n, l_0, gap, f0, fl)
emat = gap/l_0 * (ones(n, n) - Matrix{Float64}(I, n, n))

function make_histogram(data)
    x = sort(unique(data))
    y = [sum(data .== t) for t in x]
    return x, y/sum(y)
end

# Run one rep
function Run(N, l, emat, F, κ, l_0, gap, steps)
    ad_steps = Int(floor((steps*(1+κ)/2)))
    E_arr = zeros(Float64, ad_steps)
    L_arr = zeros(Float64, ad_steps)
    pop = Jevo.mono_pop(N=N, l=l)
    Jevo.initiate!(pop, emat)
    for r in 1:ad_steps
        Jevo.bp_substitution!(pop, emat, F)
        if rand() < κ/N
            Jevo.driver_mutation!(pop)
        end
        if rand() < 1/(N * (1+κ))
            Jevo.l_substitution!(pop, emat, F)
        end
        E_arr[r] = Jevo.get_energy(pop, emat)
        L_arr[r] = length(pop.seqs)
        if E_arr[r]/(L_arr[r] * gap/l_0) > (Jevo.γ_1(l, n, l_0) + 0.15)
            Jevo.initiate!(pop, emat)
        end
    end
    return E_arr, L_arr
end

x_list = Array{Array{Float64, 1}, 1}()
y_list = Array{Array{Float64, 1}, 1}()
κ_list = Array{Array{Float64, 1}, 1}()
# Predefine output arrays
for i in 1:length(κ_arr)
    push!(x_list, Float64[])
    push!(y_list, Float64[])
    push!(κ_list, Float64[])
end

Threads.@threads for i in 1:length(κ_arr)
    E, L = Run(N, 130, emat, F, κ_arr[i], l_0, gap, 10^8)
    x, y = make_histogram(L)
    x_list[i] = x
    y_list[i] = y
    κ_list[i] = ones(length(x)) .* κ_arr[i]
    println("Run $i done.")
end

df = DataFrame(p_l=[(y_list...)...], l=[(x_list...)...], kappa=[(κ_list...)...])

CSV.write(date * "_results.csv", df)

open(date*"_METADATA.txt", "a") do io
    write(io, "gap=$gap\n")
    write(io, "l_0=$l_0\n")
    write(io, "f0=$f0\n")
    write(io, "fl=$fl\n")
    write(io, "kappa=$κ_arr\n")
    write(io, "n=$n\n")
    write(io, "N=$N\n")
end