text:
#=
Copyright 2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

# These are more complex case folding functions, and maybe belong in a separate UTF8Str.jl package

# Note: these only check for cases in Unicode where 2 byte sequences
# could expand to 3 byte sequences.  In the standard Unicode tables,
# that is the only expansion that occurs for upper or lower case

function _lower_utf8(beg, off, len)
    # Note, the final length may be larger or smaller
    buf, out = _allocate(UInt8, len)
    unsafe_copyto!(out, beg, off)
    fin = beg + len
    pnt = beg + off
    outend = out + len
    out += off
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0x80
            set_codeunit!(out, ch + (_isupper_a(ch) << 5))
            out += 1
        elseif ch < 0xc4
            ch = (ch << 6) | (get_codeunit(pnt += 1) & 0x3f)
            out = output_utf8_2byte!(out, ch + (_isupper_l(ch) << 5))
        elseif ch < 0xe0
            # 2 byte
            c16 = get_utf8_2byte(pnt += 1, ch)
            if _isupper_u(c16)
                c16 = _lowercase_u(c16)
                # Check if still 2 byte, could increase to 3 byte, decrease to 1 byte
                if c16 < 0x80
                    set_codeunit!(out, c16%UInt8)
                    out += 1
                elseif c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    # Check to see if we need to resize
                    diff = (outend - out - 3) - (fin - pnt - 1)
                    if diff < 0
                        outend -= diff
                        resize!(buf, outend - out)
                        out = pointer(buf)
                        outend = out + sizeof(buf)
                    end
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_2byte!(out, c16)
            end
        elseif ch < 0xf0
            # 3 byte
            c16 = get_utf8_3byte(pnt += 2, ch)
            if _isupper_u(c16)
                c16 = _lowercase_u(c16)
                # Check if still 3 byte, could drop to 2 byte
                if c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_3byte!(out, c16)
            end
        else
            # 4 byte
            c32 = get_utf8_4byte(pnt += 3, ch)
            _isupper_u(c32) && (c32 = _lowercase_u(c32))
            out = output_utf8_4byte!(out, c32)
        end
        pnt += 1
    end
    out < outend && (buf = resize!(buf, out - pointer(buf)))
    Str(UTF8CSE, buf)
end

function _upper_utf8(beg, off, len)
    # Note, the final length may be larger or smaller
    buf, out = _allocate(UInt8, len)
    unsafe_copyto!(out, beg, off)
    fin = beg + len
    pnt = beg + off
    outend = out + len
    out += off
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0x80
            set_codeunit!(out, ch - (_islower_a(ch)<<5))
            out += 1
        elseif ch < 0xc4
            ch = (ch << 6) | (get_codeunit(pnt += 1) & 0x3f)
            c16 = ch == 0xb5 ? 0x39c : (ch == 0xff ? 0x178 : (ch - _can_upper_l(ch)<<5)%UInt16)
            out = output_utf8_2byte!(out, c16)
        elseif ch < 0xe0
            # 2 byte
            c16 = get_utf8_2byte(pnt += 1, ch)
            if _islower_u(c16)
                c16 = _uppercase_u(c16)
                # Check if still 2 byte, could increase to 3 byte, or decrease to 1 byte
                if c16 < 0x80
                    set_codeunit!(out, c16%UInt8)
                    out += 1
                elseif c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    # Increasing from 2 to 3 bytes, check to see if we need to resize
                    diff = (outend - out - 3) - (fin - pnt - 1)
                    if diff < 0
                        outend -= diff
                        resize!(buf, outend - out)
                        out = pointer(buf)
                        outend = out + sizeof(buf)
                    end
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_2byte!(out, c16)
            end
        elseif ch < 0xf0
            # 3 byte
            c16 = get_utf8_3byte(pnt += 2, ch)
            if _islower_u(c16)
                c16 = _uppercase_u(c16)
                # Check if still 3 byte, uppercase form could drop to 2 byte
                if c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_3byte!(out, c16)
            end
        else
            # 4 byte
            c32 = get_utf8_4byte(pnt += 3, ch)
            _islower_u(c32) && (c32 = _uppercase_u(c32))
            out = output_utf8_4byte!(out, c32)
        end
        pnt += 1
    end
    out < outend && (buf = resize!(buf, out - pointer(buf)))
    Str(UTF8CSE, buf)
end

function lowercase(str::Str{UTF8CSE})
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            (ch < 0x80
             ? _isupper_a(ch)
             : (ch < 0xc4
                ? _isupper_l((ch << 6) | (get_codeunit(pnt += 1) & 0x3f))
                : _isupper_u(ch >= 0xf0
                             ? get_utf8_4byte(pnt += 3, ch)
                             : (ch < 0xe0
                                ? get_utf8_2byte(pnt += 1, ch)
                                : get_utf8_3byte(pnt += 2, ch))%UInt32))) &&
                                    return _lower_utf8(beg, prv-beg, ncodeunits(str))
            pnt += 1
        end
        str
    end
end

function uppercase(str::Str{UTF8CSE})
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            (ch < 0x80
             ? _islower_a(ch)
             : (ch < 0xc4
                ? _can_upper_only_latin((ch << 6) | (get_codeunit(pnt += 1) & 0x3f))
                : _islower_u(ch >= 0xf0
                             ? get_utf8_4byte(pnt += 3, ch)
                             : (ch < 0xe0
                                ? get_utf8_2byte(pnt += 1, ch)
                                : get_utf8_3byte(pnt += 2, ch))%UInt32))) &&
                                    return _upper_utf8(beg, prv-beg, ncodeunits(str))
            pnt += 1
        end
        str
    end
end


text:
# This file was generated, do not modify it. # hide
#hideall
for f in readdir(joinpath("src", "pages"))
    println(f)
end

ext: jl

max_stars_repo_name: abhishalya/julialang-judoc

max_issues_repo_name: abhishalya/julialang-judoc

max_forks_repo_name: abhishalya/julialang-judoc


text:
# This file is a part of RadiationSpectra.jl, licensed under the MIT License (MIT).
using Test

using Distributions
using StatsBase
using IntervalSets
using ValueShapes

# using BAT
using RadiationSpectra

@testset "Package RadiationSpectra" begin
    T = Float64
    N = 10^6
    true_pars = T[N, 0.5, 1]
    h = fit(Histogram, rand(Normal(true_pars[2:3]...), N), nbins = 400)

    @testset "Fit - Optimizer" begin
        fitted_dist, opt_result_nt = RadiationSpectra.fit(RadiationSpectra.NormalPeakUvD, h)
        @test isapprox(true_pars[1], fitted_dist.A, rtol = 1e-2)
        @test isapprox(true_pars[2], fitted_dist.UvNormal.μ, rtol = 1e-2)
        @test isapprox(true_pars[3], fitted_dist.UvNormal.σ, rtol = 1e-2)
        @test isapprox(0, fitted_dist.bgll, rtol = true_pars[1] * 1e-2)
        @test isapprox(0, fitted_dist.bglr, rtol = true_pars[1] * 1e-2)
    end

    @testset "Custom Spectrum Density" begin
        struct CustomSpectrumDensity{T}  <: RadiationSpectra.UvSpectrumDensity{T}
            A::T
            σ::T
            μ::T
            offset::T
        end
        RadiationSpectra.evaluate(d::CustomSpectrumDensity, x) = d.A * pdf(Normal(d.μ, d.σ), x) + d.offset
        # Constructor
        function CustomSpectrumDensity(nt::NamedTuple{(:A, :μ, :σ, :offset)})
            T = promote_type(typeof.(values(nt))...)
            CustomSpectrumDensity(T(nt.A), T(nt.σ), T(nt.μ), T(nt.offset))
        end
        p0 = (A = 0.9*N, μ = 0.0, σ = 2.0, offset = 1.0)
        lower_bounds = (A = 0.0, μ = -5, σ = 0.01, offset = 0.0)
        upper_bounds = (A = 10.0*N, μ = 5, σ = 200.0, offset = 500.0)

        fitted_dist, backend_result = fit(CustomSpectrumDensity, h, p0, lower_bounds, upper_bounds)
        @test isapprox(true_pars[1], fitted_dist.A, rtol = 1e-2)
        @test isapprox(true_pars[2], fitted_dist.μ, rtol = 1e-2)
        @test isapprox(true_pars[3], fitted_dist.σ, rtol = 1e-2)
        @test isapprox(0, fitted_dist.offset, atol = N / 1000)

        # Fix one parameter in the fit:
        p0 = (A = 0.9*N, μ = 0.0, σ = 2.0, offset = 0)
        shape = NamedTupleShape(
            A = ScalarShape{Real}(),
            μ = ScalarShape{Real}(),
            σ = ScalarShape{Real}(),
            offset = ConstValueShape{Real}(p0.offset),
        )
        fitted_dist, backend_result = fit(CustomSpectrumDensity, h, p0, lower_bounds, upper_bounds, shape)
        @test isapprox(true_pars[1], fitted_dist.A, rtol = 1e-2)
        @test isapprox(true_pars[2], fitted_dist.μ, rtol = 1e-2)
        @test isapprox(true_pars[3], fitted_dist.σ, rtol = 1e-2)
        @test fitted_dist.offset == p0.offset
    end

    @testset "Auto Calibration" begin
        h_uncal = RadiationSpectra.get_example_spectrum()

        gamma_lines = T[609.312, 911.204, 1120.287, 1460.830, 1764.494, 2614.533]
        h_cal, h_deconv, peakPositions, threshold, c, c_precal = RadiationSpectra.calibrate_spectrum(h_uncal, gamma_lines, min_n_peaks = 30, σ = 2.0)

        # @info "Calibration constant: c = $(c)"
        c_true = 0.011259105696794384
        @testset "Auto calibration" begin
            @test abs(c - c_true) / c_true < 0.001
        end
    end

    # @testset "BAT Backend" begin
    #     # ENV["JULIA_DEBUG"] = "BAT"
    #     fitted_dist, bat_samples_result = RadiationSpectra.fit(RadiationSpectra.NormalPeakUvD, float(h), backend = :BAT, strict = false)
    #     @test isapprox(true_pars[1], fitted_dist.A, rtol = 1e-2)
    #     @test isapprox(true_pars[2], fitted_dist.UvNormal.μ, rtol = 1e-2)
    #     @test isapprox(true_pars[3], fitted_dist.UvNormal.σ, rtol = 1e-2)
    # end
end # testset


text:
# Conversions
# -----------

# no-op conversions

for CV in CVparametric
    @eval begin
        convert{T}(::Type{$CV{T}}, c::$CV{T}) = c
        convert(::Type{$CV}, c::$CV) = c
    end
end
convert(::RGB24, c::RGB24) = c
convert(::ARGB32, c::ARGB32) = c

# changing the datatype without changing the underlying colorspace
# Note: this assumes that the underlying type constructor doesn't swap the layout!
for CV in CVparametric
    @eval begin
        convert{T}(::Type{$CV{T}}, c::$CV) = $CV{T}(convert(T, getfield(c, 1)), convert(T, getfield(c, 2)), convert(T, getfield(c, 3)))
    end
end

# Conversions where the datatype is not specified
for (ElementClass, Colorspace) in [(Fractional, CVfractional), (FloatingPoint, CVfloatingpoint)]
    for CV in Colorspace
        @eval begin
            # preserves the datatype of the original space
            convert{T<:$ElementClass}(::Type{$CV}, c::ColorValue{T}) = convert($CV{T}, c)
            # fallback is Float64 (needed for RGB24)
            convert{T}(::Type{$CV}, c::ColorValue{T}) = convert($CV{Float64}, c)
        end
    end
end
convert{T}(::Type{RGB}, c::ColorValue{T}) = _convert(RGB{T}, c)
convert{T}(::Type{RGB{T}}, c::ColorValue) = _convert(RGB{T}, c)
convert(::Type{RGB}, c::RGB24) = _convert(RGB{Float64}, c)

# Everything to RGB
# -----------------

function _convert{CV<:AbstractRGB}(::Type{CV}, c::HSV)
    h = c.h / 60
    i = floor(h)
    f = h - i
    if round(Int, i) & 1 == 0
        f = 1 - f
    end
    m = c.v * (1 - c.s)
    n = c.v * (1 - c.s * f)
    i = round(Int, i)
    if i == 6 || i == 0; CV(c.v, n, m)
    elseif i == 1;       CV(n, c.v, m)
    elseif i == 2;       CV(m, c.v, n)
    elseif i == 3;       CV(m, n, c.v)
    elseif i == 4;       CV(n, m, c.v)
    else;                CV(c.v, m, n)
    end
end

function _convert{CV<:AbstractRGB}(::Type{CV}, c::HSL)
    function qtrans(u, v, hue)
        if     hue > 360; hue -= 360
        elseif hue < 0;   hue += 360
        end

        if     hue < 60;  u + (v - u) * hue / 60
        elseif hue < 180; v
        elseif hue < 240; u + (v - u) * (240 - hue) / 60
        else;             u
        end
    end

    v = c.l <= 0.5 ? c.l * (1 + c.s) : c.l + c.s - (c.l * c.s)
    u = 2 * c.l - v

    if c.s == 0; CV(c.l, c.l, c.l)
    else;        CV(qtrans(u, v, c.h + 120),
                        qtrans(u, v, c.h),
                        qtrans(u, v, c.h - 120))
    end
end

const M_XYZ_RGB = [ 3.2404542 -1.5371385 -0.4985314
                   -0.9692660  1.8760108  0.0415560
                    0.0556434 -0.2040259  1.0572252 ]

correct_gamut{CV<:AbstractRGB}(c::CV) = CV(clamp01(c.r), clamp01(c.g), clamp01(c.b))
clamp01{T<:Fractional}(v::T) = ifelse(v < zero(T), zero(T), ifelse(v > one(T), one(T), v))

function srgb_compand(v)
    v <= 0.0031308 ? 12.92v : 1.055v^(1/2.4) - 0.055
end

function _convert{CV<:AbstractRGB}(::Type{CV}, c::XYZ)
    ans = M_XYZ_RGB * [c.x, c.y, c.z]
    correct_gamut(CV(srgb_compand(ans[1]),
                     srgb_compand(ans[2]),
                     srgb_compand(ans[3])))
end

_convert{CV<:AbstractRGB}(::Type{CV}, c::xyY)    = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::Lab)    = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::LCHab)  = _convert(CV, convert(Lab{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::Luv)    = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::LCHuv)  = _convert(CV, convert(Luv{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::DIN99)  = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::DIN99o) = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::DIN99d) = _convert(CV, convert(XYZ{eltype(c)}, c))
_convert{CV<:AbstractRGB}(::Type{CV}, c::LMS)    = _convert(CV, convert(XYZ{eltype(c)}, c))

_convert{CV<:AbstractRGB{Ufixed8}}(::Type{CV}, c::RGB24) = CV(Ufixed8(c.color&0x00ff0000>>>16,0), Ufixed8(c.color&0x0000ff00>>>8,0), Ufixed8(c.color&0x000000ff,0))
_convert{CV<:AbstractRGB}(::Type{CV}, c::RGB24) = CV((c.color&0x00ff0000>>>16)/255, ((c.color&0x0000ff00)>>>8)/255, (c.color&0x000000ff)/255)

# Everything to HSV
# -----------------

function convert{T}(::Type{HSV{T}}, c::AbstractRGB)
    c_min = min(c.r, c.g, c.b)
    c_max = max(c.r, c.g, c.b)
    if c_min == c_max
        return HSV(zero(T), zero(T), c_max)
    end

    if c_min == c.r
        f = c.g - c.b
        i = convert(T, 3)
    elseif c_min == c.g
        f = c.b - c.r
        i = convert(T, 5)
    else
        f = c.r - c.g
        i = convert(T, 1)
    end

    HSV{T}(60 * (i - f / (c_max - c_min)),
        (c_max - c_min) / c_max,
        c_max)
end

convert{T}(::Type{HSV{T}}, c::ColorValue) = convert(HSV{T}, convert(RGB{T}, c))

# Everything to HSL
# -----------------

function convert{T}(::Type{HSL{T}}, c::AbstractRGB)
    c_min = min(c.r, c.g, c.b)
    c_max = max(c.r, c.g, c.b)
    l = (c_max + c_min) / 2

    if c_max == c_min
        return HSL(zero(T), zero(T), l)
    end

    if l < 0.5; s = (c_max - c_min) / (c_max + c_min)
    else;       s = (c_max - c_min) / (convert(T, 2) - c_max - c_min)
    end

    if c_max == c.r
        h = (c.g - c.b) / (c_max - c_min)
    elseif c_max == c.g
        h = convert(T, 2) + (c.b - c.r) / (c_max - c_min)
    else
        h = convert(T, 4) + (c.r - c.g) / (c_max - c_min)
    end

    h *= 60
    if h < 0
        h += 360
    elseif h > 360
        h -= 360
    end

    HSL{T}(h,s,l)
end

convert{T}(::Type{HSL{T}}, c::ColorValue) = convert(HSL{T}, convert(RGB{T}, c))

# Everything to XYZ
# -----------------

function invert_rgb_compand(v)
    v <= 0.04045 ? v/12.92 : ((v+0.055) /1.055)^2.4
end

const M_RGB_XYZ =
    [ 0.4124564  0.3575761  0.1804375
      0.2126729  0.7151522  0.0721750
      0.0193339  0.1191920  0.9503041 ]

function convert{T}(::Type{XYZ{T}}, c::AbstractRGB)
    v = [invert_rgb_compand(c.r),
         invert_rgb_compand(c.g),
         invert_rgb_compand(c.b)]
    ans = M_RGB_XYZ * v
    XYZ{T}(ans[1], ans[2], ans[3])
end

convert{T}(::Type{XYZ{T}}, c::HSV) = convert(XYZ{T}, convert(RGB{T}, c))
convert{T}(::Type{XYZ{T}}, c::HSL) = convert(XYZ{T}, convert(RGB{T}, c))

function convert{T}(::Type{XYZ{T}}, c::xyY)
    X = c.Y*c.x/c.y
    Z = c.Y*(1-c.x-c.y)/c.y
    XYZ{T}(X, c.Y, Z)
end

const xyz_epsilon = 216 // 24389
const xyz_kappa   = 24389 // 27

function convert{T}(::Type{XYZ{T}}, c::Lab, wp::XYZ)
    fy = (c.l + 16) / 116
    fx = c.a / 500 + fy
    fz = fy - c.b / 200

    fx3 = fx^3
    fz3 = fz^3

    x = fx3 > xyz_epsilon ? fx3 : (116fx - 16) / xyz_kappa
    y = c.l > xyz_kappa * xyz_epsilon ? ((c. l+ 16) / 116)^3 : c.l / xyz_kappa
    z = fz3 > xyz_epsilon ? fz3 : (116fz - 16) / xyz_kappa

    XYZ{T}(x*wp.x, y*wp.y, z*wp.z)
end

convert{T}(::Type{XYZ{T}}, c::Lab)   = convert(XYZ{T}, c, WP_DEFAULT)
convert{T}(::Type{XYZ{T}}, c::LCHab) = convert(XYZ{T}, convert(Lab{T}, c))
convert{T}(::Type{XYZ{T}}, c::DIN99) = convert(XYZ{T}, convert(Lab{T}, c))
convert{T}(::Type{XYZ{T}}, c::DIN99o) = convert(XYZ{T}, convert(Lab{T}, c))

convert{T<:Ufixed}(::Type{XYZ{T}}, c::LCHab) = convert(XYZ{T}, convert(Lab{eltype(c)}, c))
convert{T<:Ufixed}(::Type{XYZ{T}}, c::DIN99) = convert(XYZ{T}, convert(Lab{eltype(c)}, c))
convert{T<:Ufixed}(::Type{XYZ{T}}, c::DIN99o) = convert(XYZ{T}, convert(Lab{eltype(c)}, c))

function xyz_to_uv(c::XYZ)
    d = c.x + 15c.y + 3c.z
    d==0 && return (d, d)
    u = 4c.x / d
    v = 9c.y / d
    return (u, v)
end

function convert{T}(::Type{XYZ{T}}, c::Luv, wp::XYZ = WP_DEFAULT)
    (u_wp, v_wp) = xyz_to_uv(wp)

    a = (52 * (c.l==0 ? zero(T) : c.l / (c.u + 13 * c.l * u_wp)) - 1) / 3
    y = c.l > xyz_kappa * xyz_epsilon ? wp.y * ((c.l + 16) / 116)^3 : wp.y * c.l / xyz_kappa
    b = -5y
    d = y * (39 * (c.l==0 ? zero(T) : c.l / (c.v + 13 * c.l * v_wp)) - 5)
    x = d==b ? zero(T) : (d - b) / (a + 1//3)
    z = a * x + b + zero(T)

    XYZ{T}(x, y, z)
end

convert{T}(::Type{XYZ{T}}, c::LCHuv) = convert(XYZ{T}, convert(Luv{T}, c))
convert{T<:Ufixed}(::Type{XYZ{T}}, c::LCHuv) = convert(XYZ{T}, convert(Luv{eltype(c)}, c))

function convert{T}(::Type{XYZ{T}}, c::DIN99d)

    # Go back to C-h space
    # FIXME: Clean this up (why is there no atan2d?)
    h = rad2deg(atan2(c.b,c.a)) + 50
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end

    C = sqrt(c.a^2 + c.b^2)

    # Intermediate terms
    G = (exp(C/22.5)-1)/0.06
    f = G*sind(h - 50)
    ee = G*cosd(h - 50)

    l = (exp(c.l/325.22)-1)/0.0036
    # a = ee*cosd(50) - f/1.14*sind(50)
    a = ee*0.6427876096865393 - f/1.14*0.766044443118978
    # b = ee*sind(50) - f/1.14*cosd(50)
    b = ee*0.766044443118978 - f/1.14*0.6427876096865393

    adj = convert(XYZ, Lab(l, a, b))

    XYZ{T}((adj.x + 0.12*adj.z)/1.12, adj.y, adj.z)

end

function convert{T}(::Type{XYZ{T}}, c::LMS)
    ans = CAT02_INV * [c.l, c.m, c.s]
    XYZ{T}(ans[1], ans[2], ans[3])
end

convert{T}(::Type{XYZ{T}}, c::RGB24) = convert(XYZ{T}, convert(RGB{T}, c))

# Everything to xyY
# -----------------

function convert{T}(::Type{xyY{T}}, c::XYZ)

    x = c.x/(c.x + c.y + c.z)
    y = c.y/(c.x + c.y + c.z)

    xyY{T}(x, y, convert(typeof(x), c.y))

end

convert{T}(::Type{xyY{T}}, c::ColorValue) = convert(xyY{T}, convert(XYZ{T}, c))

# Everything to Lab
# -----------------

convert{T}(::Type{Lab{T}}, c::AbstractRGB) = convert(Lab{T}, convert(XYZ{T}, c))
convert{T}(::Type{Lab{T}}, c::HSV) = convert(Lab{T}, convert(RGB{T}, c))
convert{T}(::Type{Lab{T}}, c::HSL) = convert(Lab{T}, convert(RGB{T}, c))

function convert{T}(::Type{Lab{T}}, c::XYZ, wp::XYZ)
    function f(v)
        v > xyz_epsilon ? cbrt(v) : (xyz_kappa * v + 16) / 116
    end

    fx, fy, fz = f(c.x / wp.x), f(c.y / wp.y), f(c.z / wp.z)
    Lab{T}(116fy - 16, 500(fx - fy), 200(fy - fz))
end

convert{T}(::Type{Lab{T}}, c::XYZ{T}) = convert(Lab{T}, c, WP_DEFAULT)

function convert{T}(::Type{Lab{T}}, c::LCHab)
    hr = deg2rad(c.h)
    Lab{T}(c.l, c.c * cos(hr), c.c * sin(hr))
end

function convert{T}(::Type{Lab{T}}, c::DIN99)

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99 space
    cc = sqrt(c.a^2 + c.b^2)

    # NOTE: This is calculated in degrees, against the standard, to save
    # computation steps later.
    if (c.a > 0 && c.b >= 0)
        h = atand(c.b/c.a)
    elseif (c.a == 0 && c.b > 0)
        h = 90
    elseif (c.a < 0)
        h = 180+atand(c.b/c.a)
    elseif (c.a == 0 && c.b < 0)
        h = 270
    elseif (c.a > 0 && c.b <= 0)
        h = 360 + atand(c.b/c.a)
    else
        h = 0
    end

    # Temporary variable for chroma
    g = (e^(0.045*cc*kch*ke)-1)/0.045

    # Temporary redness
    ee = g*cosd(h)

    # Temporary yellowness
    f = g*sind(h)

    # CIELAB a*b*
    # ciea = ee*cosd(16) - (f/0.7)*sind(16)
    ciea = ee*0.9612616959383189 - (f/0.7)*0.27563735581699916
    # cieb = ee*sind(16) + (f/0.7)*cosd(16)
    cieb = ee*0.27563735581699916 + (f/0.7)*0.9612616959383189

    # CIELAB L*
    ciel = (e^(c.l*ke/105.51)-1)/0.0158

    Lab{T}(ciel, ciea, cieb)
end

function convert{T}(::Type{Lab{T}}, c::DIN99o)

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99o space
    co = sqrt(c.a^2 + c.b^2)

    # hue angle h99o
    h = atan2(c.b, c.a)

    # revert rotation by 26°
    ho= rad2deg(h)-26

    # revert logarithmic chroma compression
    g = (e^(co*kch*ke/23.0)-1)/0.075

    # Temporary redness
    eo = g*cosd(ho)

    # Temporary yellowness
    fo = g*sind(ho)

    # CIELAB a*b* (revert b* axis compression)
    # ciea = eo*cosd(26) - (fo/0.83)*sind(26)
    ciea = eo*0.898794046299167 - (fo/0.83)*0.4383711467890774
    # cieb = eo*sind(26) + (fo/0.83)*cosd(26)
    cieb = eo*0.4383711467890774 + (fo/0.83)*0.898794046299167

    # CIELAB L* (revert logarithmic lightness compression)
    ciel = (e^(c.l*ke/303.67)-1)/0.0039

    Lab{T}(ciel, ciea, cieb)
end

convert{T}(::Type{Lab{T}}, c::ColorValue) = convert(Lab{T}, convert(XYZ{T}, c))

# Everything to Luv
# -----------------

convert{T}(::Type{Luv{T}}, c::AbstractRGB) = convert(Luv{T}, convert(XYZ{T}, c))
convert{T}(::Type{Luv{T}}, c::HSV) = convert(Luv{T}, convert(RGB{T}, c))
convert{T}(::Type{Luv{T}}, c::HSL) = convert(Luv{T}, convert(RGB{T}, c))

function convert{T}(::Type{Luv{T}}, c::XYZ, wp::XYZ = WP_DEFAULT)
    (u_wp, v_wp) = xyz_to_uv(wp)
    (u_, v_) = xyz_to_uv(c)

    y = c.y / wp.y

    l = y > xyz_epsilon ? 116 * cbrt(y) - 16 : xyz_kappa * y
    u = 13 * l * (u_ - u_wp) + zero(T)
    v = 13 * l * (v_ - v_wp) + zero(T)

    Luv{T}(l, u, v)
end

function convert{T}(::Type{Luv{T}}, c::LCHuv)
    hr = deg2rad(c.h)
    Luv{T}(c.l, c.c * cos(hr), c.c * sin(hr))
end

convert{T}(::Type{Luv{T}}, c::ColorValue) = convert(Luv{T}, convert(XYZ{T}, c))

# Everything to LCHuv
# -------------------

function convert{T}(::Type{LCHuv{T}}, c::Luv)
    h = rad2deg(atan2(c.v, c.u))
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end
    LCHuv{T}(c.l, sqrt(c.u^2 + c.v^2), h)
end

convert{T}(::Type{LCHuv{T}}, c::ColorValue) = convert(LCHuv{T}, convert(Luv{T}, c))

# Everything to LCHab
# -------------------

function convert{T}(::Type{LCHab{T}}, c::Lab)
    h = rad2deg(atan2(c.b, c.a))
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end
    LCHab{T}(c.l, sqrt(c.a^2 + c.b^2), h)
end

convert{T}(::Type{LCHab{T}}, c::ColorValue) = convert(LCHab{T}, convert(Lab{T}, c))

# Everything to DIN99
# -------------------

function convert{T}(::Type{DIN99{T}}, c::Lab)

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99 L
    l99 = (1/ke)*105.51*log(1+0.0158*c.l)

    # Temporary value for redness and yellowness
    # ee = c.a*cosd(16) + c.b*sind(16)
    ee = c.a*0.9612616959383189 + c.b*0.27563735581699916
    # f = -0.7*c.a*sind(16) + 0.7*c.b*cosd(16)
    f = -0.7*c.a*0.27563735581699916 + 0.7*c.b*0.9612616959383189

    # Temporary value for chroma
    g = sqrt(ee^2 + f^2)

    # Hue angle
    # Calculated in degrees, against the specification.
    if (ee > 0 && f >= 0)
        h = atand(f/ee)
    elseif (ee == 0 && f > 0)
        h = 90
    elseif (ee < 0)
        h = 180+atand(f/ee)
    elseif (ee == 0 && f < 0)
        h = 270
    elseif (ee > 0 && f <= 0)
        h = 360 + atand(f/ee)
    else
        h = 0
    end

    # DIN99 chroma
    cc = log(1+0.045*g)/(0.045*kch*ke)

    # DIN99 chromaticities
    a99, b99 = cc*cosd(h), cc*sind(h)

    DIN99{T}(l99, a99, b99)

end

convert{T}(::Type{DIN99{T}}, c::ColorValue) = convert(DIN99{T}, convert(Lab{T}, c))

# Everything to DIN99d
# --------------------

function convert{T}(::Type{DIN99d{T}}, c::XYZ{T})

    # Apply tristimulus-space correction term
    adj_c = XYZ(1.12*c.x - 0.12*c.z, c.y, c.z)

    # Apply L*a*b*-space correction
    lab = convert(Lab, adj_c)
    adj_L = 325.22*log(1+0.0036*lab.l)

    # Calculate intermediate parameters
    # ee = lab.a*cosd(50) + lab.b*sind(50)
    ee = lab.a*0.6427876096865393 + lab.b*0.766044443118978
    # f = 1.14*(lab.b*cosd(50) - lab.a*sind(50))
    f = 1.14*(lab.b*0.6427876096865393 - lab.a*0.766044443118978)
    G = sqrt(ee^2+f^2)

    # Calculate hue/chroma
    C = 22.5*log(1+0.06*G)
    # FIXME: Clean this up (why is there no atan2d?)
    h = rad2deg(atan2(f,ee)) + 50
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end

    DIN99d{T}(adj_L, C*cosd(h), C*sind(h))

end

convert{T}(::Type{DIN99d{T}}, c::ColorValue) = convert(DIN99d{T}, convert(XYZ{T}, c))

# Everything to DIN99o
# -------------------

function convert{T}(::Type{DIN99o{T}}, c::Lab)

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99o L (logarithmic compression)
    l99 = 303.67/ke*log(1+0.0039*c.l)

    # Temporary value for redness and yellowness
    # including rotation by 26°
    # eo = c.a*cosd(26) + c.b*sind(26)
    eo = c.a*0.898794046299167 + c.b*0.4383711467890774
    # compression along the yellowness (blue-yellow) axis
    # fo = 0.83 * (c.b*cosd(26) - c.a*sind(26))
    fo = 0.83 * (c.b*0.898794046299167 - c.a*0.4383711467890774)

    # Temporary value for chroma
    go = sqrt(eo^2 + fo^2)
    ho = atan2(fo,eo)
    # rotation of the color space by 26°
    h  = rad2deg(ho) + 26

    # DIN99o chroma (logarithmic compression)
    cc = 23.0*log(1+0.075*go)/(kch*ke)

    # DIN99o chromaticities
    a99, b99 = cc*cosd(h), cc*sind(h)

    DIN99o{T}(l99, a99, b99)

end

convert{T}(::Type{DIN99o{T}}, c::ColorValue) = convert(DIN99o{T}, convert(Lab{T}, c))

# Everything to LMS
# -----------------

# Chromatic adaptation from CIECAM97s
const CAT97s = [ 0.8562  0.3372 -0.1934
                -0.8360  1.8327  0.0033
                 0.0357 -0.0469  1.0112 ]

const CAT97s_INV = inv(CAT97s)

# Chromatic adaptation from CIECAM02
const CAT02 = [ 0.7328 0.4296 -0.1624
               -0.7036 1.6975  0.0061
                0.0030 0.0136  0.9834 ]

const CAT02_INV = inv(CAT02)

function convert{T}(::Type{LMS{T}}, c::XYZ{T})
    ans = CAT02 * [c.x, c.y, c.z]
    LMS{T}(ans[1], ans[2], ans[3])
end

convert{T}(::Type{LMS{T}}, c::ColorValue) = convert(LMS{T}, convert(XYZ{T}, c))

# Everything to RGB24
# -------------------

convert(::Type{RGB24}, c::AbstractRGB{Ufixed8}) = RGB24(c.r, c.g, c.b)
convert(::Type{RGB24}, c::AbstractRGB) = RGB24(round(Uint32, 255*c.r)<<16 +
                                               round(Uint32, 255*c.g)<<8 +
                                               round(Uint32, 255*c.b))
to32(x::Ufixed8) = convert(Uint32, reinterpret(x))
convert(::Type{RGB24}, val::Uint32) = RGB24(val & 0x00ffffff)

convert(::Type{RGB24}, c::ColorValue) = convert(RGB24, convert(RGB{Ufixed8}, c))

# To Uint32
# ----------------

convert(::Type{Uint32}, c::RGB24) = c.color

convert(::Type{Uint32}, ac::ARGB32) = ac.color

#### AlphaColorValue

convert{C,T}(::Type{AlphaColorValue{C,T}}, c::AlphaColorValue{C,T}) = c
function convert{C,T,D,U}(::Type{AlphaColorValue{C,T}}, c::AbstractAlphaColorValue{D,U})
    AlphaColorValue{C,T}(convert(C, c.c), c.alpha)
end
for CV in CVAlpha
    @eval begin
        function convert{D,T}(AC::Type{$CV}, c::AbstractAlphaColorValue{D,T})
            AlphaColorValue(convert(colortype($CV), c.c), c.alpha)
        end
        convert(AC::Type{$CV}, c::ColorValue) = AlphaColorValue(convert(colortype($CV), c))
    end
end

convert{C,T}(::Type{AlphaColorValue{C,T}}, c::ColorValue) =
    AlphaColorValue{C,T}(convert(C, c), one(T))
convert{C<:ColorValue,D,T}(::Type{C}, c::AlphaColorValue{D, T}) = convert(C, c.c)

convert(::Type{ARGB32}, c::ARGB32) = c
convert{CV<:AbstractRGB{Ufixed8}}(::Type{ARGB32}, c::AbstractAlphaColorValue{CV,Ufixed8}) =
    ARGB32(c.c.r, c.c.g, c.c.b, c.alpha)
convert(::Type{ARGB32}, c::AbstractAlphaColorValue) =
    ARGB32(convert(RGB24, c.c).color | round(Uint32, 255*c.alpha)<<24)
convert(::Type{ARGB32}, c::ColorValue) = ARGB32(convert(RGB24, c).color | 0xff000000)
convert(::Type{ARGB32}, val::Uint32) = ARGB32(val)

function convert(::Type{RGBA{Ufixed8}}, ac::ARGB32)
    RGBA{Ufixed8}(Ufixed8(ac.color&0x00ff0000>>>16,0),
                  Ufixed8(ac.color&0x0000ff00>>>8,0),
                  Ufixed8(ac.color&0x000000ff,0),
                  Ufixed8(ac.color>>>24,0))
end
function convert{T}(::Type{RGBA{T}}, c::ARGB32)
    RGBA{T}((c.color&0x00ff0000>>>16)/255,
            (c.color&0x0000ff00>>>8)/255,
            (c.color&0x000000ff)/255,
            (c.color>>>24)/255)
end
function convert{C,T}(::Type{AlphaColorValue{C,T}}, c::ARGB32)
    convert(AlphaColorValue{C,T}, convert(RGBA{Ufixed8}, c))
end

### Equality
==(c1::AbstractRGB, c2::AbstractRGB) = c1.r == c2.r && c1.g == c2.g && c1.b == c2.b


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# Pyr5 elasticity with point load

nodes = Dict(
    1 => [-1.0,-1.0,-1.0],
    2 => [ 1.0,-1.0,-1.0],
    3 => [ 1.0, 1.0,-1.0],
    4 => [-1.0, 1.0,-1.0],
    5 => [ 0.0, 0.0, 1.0])

element1 = Element(Pyr5, (1, 2, 3, 4, 5))
baseQuad = Element(Quad4, (1, 2, 3, 4))
tipPoint = Element(Poi1, (5,))

update!((element1, baseQuad, tipPoint), "geometry", nodes)
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)

update!(tipPoint, "displacement traction force 1",  5.0)
update!(tipPoint, "displacement traction force 2", -7.0)
update!(tipPoint, "displacement traction force 3",  3.0)

elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
add_elements!(elasticity_problem, element1, baseQuad, tipPoint)

update!(baseQuad, "displacement 1", 0.0)
update!(baseQuad, "displacement 2", 0.0)
update!(baseQuad, "displacement 3", 0.0)
boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
add_elements!(boundary_problem, baseQuad)

analysis = Analysis(Linear)
add_problems!(analysis, elasticity_problem, boundary_problem)
run!(analysis)

xi, time = (0.0, 0.0, 1.0), 0.0
u = element1("displacement", xi, time)
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@debug("displacement at tip", xi, time, u, u_expected)
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
@test isapprox(u, u_expected)


text:
fuel(mass::Real) = max(0, fld(mass, 3) - 2)

const input = map(x->parse(Int, x), eachline("input"))

println("Part 1 answer: ", sum(fuel, input))

function totfuel(mass::Real)
    total = 0
    while mass > 0
        mass = fuel(mass)
        total += mass
    end
    total
end

println("Part 2 answer: ", sum(totfuel, input))

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
mutable struct F5Data{I, SΓ<:SigPolynomialΓ{I}}
    ctx::SΓ
    R # singular ring
    trace_sig_tail_tags::Vector{Symbol}
    remasks_left::Int
end

function f5data(I::Vector{P};
                mod_order=:POT,
                trace_sig_tail_tags=Symbol[],
                index_type=UInt32,
                mask_type=UInt32,
                pos_type=UInt32,
                max_remasks=3,
                use_macaulay_bound=false,
                kwargs...) where {P <: AA.MPolyElem}

    R = parent(first(I))
    char = AA.characteristic(R)
    coeffs = Nmod32Γ(char)
    if use_macaulay_bound
        deg_bound = mac_bound(I)
    else
        deg_bound = 0
    end
    ctx = idxsigpolynomialctx(coeffs, length(I), index_type=index_type,
                              mask_type=mask_type, pos_type=pos_type,
                              deg_bound=deg_bound; kwargs...)

    for (i, f) in enumerate(I)
        sig = ctx(i, R(1))
        ctx(sig, f)
    end
    F5Data{pos_type, typeof(ctx)}(ctx, R, trace_sig_tail_tags, max_remasks)
end


text:
using Optim

export CFADistributionTheta

"Confirmatory factor analysis distribution object."
type CFADistributionTheta
    Theta_X::SparseMatrixCSC
    A::SparseMatrixCSC
    Theta_L::AbstractMatrix
end
Base.length(d::CFADistributionTheta) = Int(size(d.A)[1] + nnz(d.A) + (size(d.A)[2]-1)*size(d.A)[2]/2)

"Likelihood of model given the data represented by the covariance S from N samples"
function Distributions.loglikelihood(d::CFADistributionTheta, S::AbstractMatrix, N::Int64)
    @assert minimum(d.Theta_X.nzval) > 0.0 "Invalid CFADistributionTheta! (Theta_X has elements <= 0)"
    @assert minimum(eig(d.Theta_L)[1]) > 0.0 "Invalid CFADistributionTheta! (Theta_L has eigen value <= 0)"

    # Woodbury transformation of original form allows us to only compute a KxK inversion
    Theta = d.Theta_X - d.Theta_X*d.A*inv(d.Theta_L + d.A'*d.Theta_X*d.A)*d.A'*d.Theta_X

    v = logdet(Theta) - trace(S*Theta)
    #println("v :", v, " ", minimum(d.Theta_X.nzval), " ", minimum(eig(Theta)[1]), " ", logdet(Theta))
    v
end

function Base.rand(d::CFADistributionTheta, N::Int)
    P,K = size(d.A)
    distL = MvNormal(zeros(K), inv(d.Theta_L))
    L = rand(distL, N)
    X = d.A*L
    for i in 1:P
        X[i,:] .+= transpose(rand(Normal(0, 1 / d.Theta_X.nzval[i]), N))
    end
    X
end
Base.rand(d::CFADistributionTheta) = rand(d::CFADistributionTheta, 1)

function normalize_Sigma_L!(d::CFADistributionTheta)
    Sigma_L = inv(d.Theta_L)
    scaling = sqrt(diag(Sigma_L))
    d.Theta_L[:,:] = inv(Base.cov2cor!(Sigma_L, scaling))
    vals = nonzeros(d.A)
    K = size(d.A)[2]
    for col = 1:K
        for j in nzrange(d.A, col)
            vals[j] *= scaling[col]
        end
    end
end

"Gradient of the CFA model."
type CFADistributionThetaGradient
    Theta_X::SparseMatrixCSC
    A::SparseMatrixCSC
    Theta_L::AbstractMatrix
end

"""
Pulls values from a flattened vector into the full matrix state.
"""
function vec2state!(d::Union{CFADistributionTheta,CFADistributionThetaGradient}, x::AbstractVector)
    P = length(d.Theta_X.nzval)
    d.Theta_X.nzval[:] = x[1:P]

    offset = length(d.A.nzval)
    d.A.nzval[:] = x[P+1:P+offset]

    pos = P+offset+1
    K = size(d.Theta_L)[1]
    for i in 1:K, j in i:K
        d.Theta_L[i,j] = d.Theta_L[j,i] = x[pos]
        pos += 1
    end
end
# function vec2state!(Sigma_X::SparseMatrixCSC, A::SparseMatrixCSC, Sigma_L::AbstractMatrix, x::AbstractVector)
#     P = length(Sigma_X.nzval)
#     Sigma_X.nzval[:] = x[1:P]
#
#     offset = length(A.nzval)
#     A.nzval[:] = x[P+1:P+offset]
#
#     pos = P+offset+1
#     K = size(Sigma_L)[1]
#     for i in 1:K, j in i+1:K
#         Sigma_L[i,j] = Sigma_L[j,i] = x[pos]
#         pos += 1
#     end
# end

"""
Pulls values from the full matrix state into a flattened vector.
"""
function state2vec!(x::AbstractVector, d::Union{CFADistributionTheta,CFADistributionThetaGradient})
    P = length(d.Theta_X.nzval)
    x[1:P] = d.Theta_X.nzval

    offest = length(d.A.nzval)
    x[P+1:P+offest] = d.A.nzval

    pos = P+offest+1
    K = size(d.Theta_L)[1]
    for i in 1:K, j in i:K
        x[pos] = d.Theta_L[i,j]
        pos += 1
    end
end

function dloglikelihood(d::CFADistributionTheta, S::AbstractMatrix, N::Int64)
    P,K = size(d.A)

    # Woodbury transformation of orignal form allows us to only compute a KxK inversion
    Theta = d.Theta_X - d.Theta_X*d.A*inv(d.Theta_L + d.A'*d.Theta_X*d.A)*d.A'*d.Theta_X

    commonPart = Theta - Theta*S*Theta

    # invert Sigma_X efficiently then compute gradient
    Theta_Xinv = deepcopy(d.Theta_X)
    Theta_Xinv.nzval[:] = 1 ./ Theta_Xinv.nzval
    dTheta_X = spdiagm(diag(Theta_Xinv*commonPart*Theta_Xinv))

    Theta_Linv = inv(d.Theta_L)
    dTheta_L = 2*Theta_Linv*d.A'*commonPart*d.A*Theta_Linv
    #dTheta_Linv = 2*d.A'*commonPart*d.A
    #for i in 1:K dTheta_Linv[i,i] = 0.0 end # we shouldn't change the diagonal of the covariance matrix
    #dTheta_L = Theta_Linv*dTheta_Linv*Theta_Linv
    for i in 1:size(dTheta_L)[1] dTheta_L[i,i] /= 2 end # diagonals should not be times 2

    dA_dense = -2*commonPart*d.A*Theta_Linv
    dA = deepcopy(d.A)
    rows = rowvals(dA)
    vals = nonzeros(dA)
    for col = 1:K
        for j in nzrange(dA, col)
            vals[j] = dA_dense[rows[j],col]
        end
    end

    CFADistributionThetaGradient(dTheta_X, dA, dTheta_L)
end

function Distributions.fit_mle(::Type{CFADistributionTheta}, A::SparseMatrixCSC, S::AbstractMatrix, N::Int64; iterations=1000, show_trace=true)
    P,K = size(A)

    # init our parameters
    d = CFADistributionTheta(speye(P), deepcopy(A), eye(K))
    x_init = zeros(Int(P + nnz(A) + (K+1)*K/2))
    state2vec!(x_init, d)

    function f(x)
        vec2state!(d, x)
        try
            return -loglikelihood(d, S, N)
        catch
            1e16 # invalid set of parameters
        end
    end

    function df!(x, g)
        vec2state!(d, x)
        gState = dloglikelihood(d, S, N)
        state2vec!(g, gState)
        g[:] = -g
    end

    res = optimize(f, df!,
        x_init,
        method = :l_bfgs,
        show_trace = show_trace,
        iterations = iterations
    )

    vec2state!(d, res.minimum)

    # rescale X and Theta_L so that Sigma_L is a correlation matrix
    normalize_Sigma_L!(d)
    d
end

function fit_map(pri::Normal, ::Type{CFADistributionTheta}, A::SparseMatrixCSC, S::AbstractMatrix, N::Int64; iterations=1000, show_trace=true)
    @assert mean(pri) == 0.0 "Prior must be zero mean!"
    rho = 1/(std(pri)*N)
    P,K = size(A)

    # init our parameters
    d = CFADistributionTheta(speye(P), deepcopy(A), eye(K))
    x_init = zeros(Int(P + nnz(A) + (K+1)*K/2))
    state2vec!(x_init, d)

    function f(x)
        vec2state!(d, x)

        # get likelihood
        ll = 0.0
        try
            ll = loglikelihood(d, S, N)
        catch
            return 1e16 # invalid set of parameters
        end

        # account for prior
        #ll -= rho*sum(d.Theta_X.^2)
        ll -= rho*sum(d.A.^2)
        ll -= rho*sum(d.Theta_L.^2)

        -ll
    end

    function df!(x, g)
        vec2state!(d, x)

        # get gradient of likelihood
        gState = dloglikelihood(d, S, N)

        # account for prior
        #gState.Theta_X .-= 2*rho*d.Theta_X
        gState.A .-= 2*rho*d.A
        gState.Theta_L .-= 2*rho*d.Theta_L

        state2vec!(g, gState)
        g[:] = -g
    end

    res = optimize(f, df!,
        x_init,
        method = :l_bfgs,
        show_trace = show_trace,
        iterations = iterations
    )

    vec2state!(d, res.minimum)
    d
end


text:
#===============
    STRUCTS
===============#

struct ray
    origin::rt_tuple
    direction::rt_tuple
end

abstract type RTObject end

struct intersection
    t::Float64
    object::RayTracer.RTObject
end

#===============
    FUNCTIONS
===============#

function position(r::ray,t::T)::rt_tuple where{T<:Real}
    return r.origin + t*r.direction
end

function hit(is::Array{intersection,1})
    index = -1
    closest_i = -1
    for (idx,i) in enumerate(is)
        if i.t>0 && (i.t<closest_i || index<0)
            index = idx
            closest_i = i.t
        end
    end
    if index<0
        return nothing
    else
        return is[index]
    end
end

function transform(r::ray,M::Array{T,2})::ray where {T<:Real}
    return ray(M*r.origin,M*r.direction)
end

function set_transform!(o::RayTracer.RTObject,M::Array{T,2}) where {T<:Real}
    o.transform = M*o.transform
end


text:
# Autogenerated wrapper script for Entwine_jll for powerpc64le-linux-gnu-cxx11
export entwine, libentwine

using PDAL_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Entwine")
JLLWrappers.@declare_executable_product(entwine)
JLLWrappers.@declare_library_product(libentwine, "libentwine.so.3")
function __init__()
    JLLWrappers.@generate_init_header(PDAL_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        entwine,
        "bin/entwine",
    )

    JLLWrappers.@init_library_product(
        libentwine,
        "lib/libentwine.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
for coeff_name in coeff_list
  include(lowercase(coeff_name) * ".jl")
end


text:
using MOBC, vOptGeneric, CPLEX, PyPlot, Suppressor
using Base.Test

function random_instance(range, n, rng)
	MT = MersenneTwister(rng)
	m = vModel(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
	p1 = rand(MT, range, n)
	p2 = rand(MT, range, n)
	w = [rand(MT, range, n) for i = 1:3]
	c = sum.(w).÷2

	@variable(m, x[1:length(p1)], Bin)
	@addobjective(m, Max, dot(x, p1))
	@addobjective(m, Max, dot(x, p2))
	@constraint(m, dot(x, w[1]) <= c[1])
	@constraint(m, dot(x, w[2]) <= c[2])
	@constraint(m, dot(x, w[3]) <= c[3])
	m
end

function hard_instance(solver=CplexSolver(CPX_PARAM_SCRIND = 0))
    n = 30
    p1 = [85 ,88 ,82 ,64 ,55 ,63 ,77 ,92 ,52 ,64 ,91 ,59 ,90 ,74 ,81 ,80 ,75 ,95 ,52 ,83 ,88 ,50 ,94 ,67 ,91 ,95 ,50 ,74 ,56 ,76]
    p2 = [56 ,58 ,78 ,55 ,69 ,86 ,92 ,84 ,88 ,60 ,62 ,99 ,91 ,62 ,56 ,96 ,75 ,71 ,77 ,87 ,80 ,87 ,72 ,99 ,91 ,74 ,94 ,50 ,98 ,100]
    w1 = [79 ,72 ,60 ,72 ,53 ,53 ,56 ,69 ,57 ,73 ,55 ,66 ,51 ,67 ,74 ,57 ,85 ,98 ,73 ,77 ,58 ,69 ,70 ,76 ,72 ,71 ,87 ,55 ,51 ,83]
    w2 = [67 ,67 ,96 ,88 ,83 ,52 ,60 ,77 ,60 ,75 ,79 ,66 ,71 ,73 ,65 ,58 ,69 ,100 ,67 ,99 ,59 ,69 ,81 ,95 ,59 ,79 ,86 ,56 ,65 ,73]
    w3 = [55 ,72 ,77 ,62 ,51 ,58 ,85 ,76 ,57 ,52 ,83 ,74 ,69 ,89 ,61 ,66 ,66 ,71 ,94 ,61 ,59 ,60 ,69 ,51 ,64 ,84 ,92 ,62 ,67 ,50]
    b = [1019, 1097, 1018]

    m = vModel(solver = solver)
    @variable(m, x[1:n], Bin)
    @addobjective(m, Max, dot(x,p1))
    @addobjective(m, Max, dot(x,p2))
    @constraint(m, dot(x, w1) <= b[1])
    @constraint(m, dot(x, w2) <= b[2])
    @constraint(m, dot(x, w3) <= b[3])

    m
end

function benchmark(itemrange, nrange, runpersize ; time_limit=300, parragh=true, epsilon=true, args...)
	rng = 0
	res_stidsen=Dict{Int, Vector{Float64}}() ; nb_nodes_stidsen=Dict{Int, Vector{Int}}()
	res_parragh=Dict{Int, Vector{Float64}}() ; nb_nodes_parragh=Dict{Int, Vector{Int}}()
	res_eps=Dict{Int, Vector{Float64}}()
	for n in nrange
		res_stidsen[n]=Float64[] ; res_parragh[n]=Float64[] ; res_eps[n]=Float64[]
		nb_nodes_stidsen[n]=Int[] ; nb_nodes_parragh[n]=Int[]
		for i = 1:runpersize
			m = random_instance(itemrange, n, rng)
			rng += 1
			global m_current = m
			m_copy = copy(m)
			status,tEPS,_ = @suppress @timed solve(m_copy, method=:epsilon, round_results=true);
			valBC, tBC, _ = @timed solve_stidsen(m ; time_limit=time_limit, args...);
			valPAR, tPAR, _ = 	if parragh
									@timed solve_parragh(m ; time_limit=time_limit, args...);
								else
									valBC, tBC, _
								end

			if (!(valBC.timeout) && length(getY_N(m_copy)) != length(valBC.YN)) || (!(valPAR.timeout) && length(valPAR.YN) != length(getY_N(m_copy)))
				figure(3)
				clf()
				plot(map(first, valPAR.YN), map(last,valPAR.YN), "b>", markersize="5", label="parragh")
				plot(map(first, valBC.YN), map(last, valBC.YN), "r<", markersize="5", label="stidsen")
				plot(map(first, getY_N(m)), map(last, getY_N(m)), "go", markersize="2", label="epsilon")
				legend()
				error()
			end

			push!(res_stidsen[n], valBC.timeout ? time_limit : tBC) ; push!(res_parragh[n], valPAR.timeout ? time_limit : tPAR) ; push!(res_eps[n], tEPS)
			!(valPAR.timeout) && push!(nb_nodes_parragh[n], valPAR.nodes)
			!(valBC.timeout) && push!(nb_nodes_stidsen[n], valBC.nodes)

			#figure(1) ; clf()
			#xlabel("n") ; ylabel("t(s)") ; title("time per instance")
			###parragh || epsilon && plot([],[],"gx",markersize=3,label="λ-scalarisation")
			###parragh && plot([],[],"bx",markersize=3,label="parragh")
			##epsilon && plot([],[],"k.",markersize=3,label="epsilon")
			#for k in keys(res_stidsen)
			#	plot(fill(k, length(res_stidsen[k])), res_stidsen[k], "gx", markersize=3)
				#parragh && plot(fill(k, length(res_stidsen[k])), res_parragh[k], "bx", markersize=3)
	#		#	epsilon && plot(fill(k, length(res_stidsen[k])), res_eps[k], "k.", markersize=3)
			#end
	#		legend()

			figure(2) ; clf()
			xlabel("n") ; ylabel("#noeuds") ; title("temps et nombre de noeuds en moyenne")
			k = sort(collect(keys(res_stidsen)))
			plot(k, [mean(nb_nodes_stidsen[a]) for a in k], "gs", markersize=3) ; plot(k, [mean(nb_nodes_stidsen[a]) for a in k], "g--", linewidth=1)
			if parragh
				plot(k, [mean(nb_nodes_parragh[a]) for a in k], "bs", markersize=3) ; plot(k, [mean(nb_nodes_parragh[a]) for a in k], "b--", linewidth=1)
			end
			twinx() ; ylabel("temps moyen (s)")
			plot([],[],"g--",markersize=3,label="#noeuds λ-scalarisation") #labels dont work after twinx(), had to put a dummy plot
			parragh && plot([],[],"b--",markersize=3,label="#noeuds parragh")
			plot(k, [mean(res_stidsen[a]) for a in k], "gx", markersize=3) ; plot(k, [mean(res_stidsen[a]) for a in k], "g-", linewidth=1, label="tmoy λ-scalarisation")
			if parragh
				plot(k, [mean(res_parragh[a]) for a in k], "bx", markersize=3) ; plot(k, [mean(res_parragh[a]) for a in k], "b-", linewidth=1, label="tmoy parragh")
			end
			if epsilon
				plot(k, [mean(res_eps[a])for a in k], "k-", linewidth=1, label="tmoy epsilon")
			end
			legend()
		end
		println("\n####################################\nn = $n")
		println("tmoy eps : $(mean(res_eps[n]))") ; println("tmoy stidsen : $(mean(res_stidsen[n]))") ; parragh && println("tmoy parragh : $(mean(res_parragh[n]))")
		println("moy : $(round(mean(nb_nodes_stidsen[n]),1)) noeuds par BC Stidsen") ; parragh && println("moy : $(round(mean(nb_nodes_parragh[n]),1)) noeuds par BC Parragh")
	end
	res_stidsen, res_parragh, res_eps, nb_nodes_stidsen, nb_nodes_parragh
end

function benchmark_stidsen(itemrange, nrange, runpersize ; time_limit=300,  epsilon=true)
	rng = 0
	config_1 = "avec cover cuts"
	config_2 = "sans cover cuts"
	res_1=Dict{Int, Vector{Float64}}() ; nb_nodes_1=Dict{Int, Vector{Int}}()
	res_2=Dict{Int, Vector{Float64}}() ; nb_nodes_2=Dict{Int, Vector{Int}}()
	res_eps=Dict{Int, Vector{Float64}}()
	for n in nrange
		res_1[n]=Float64[] ; res_2[n]=Float64[] ; res_eps[n]=Float64[]
		nb_nodes_1[n]=Int[] ; nb_nodes_2[n]=Int[]
		for i = 1:runpersize
			m = random_instance(itemrange, n, rng)
			rng += 1
			global m_current = m
			m_copy = copy(m)
			status,tEPS,_ = @suppress @timed solve(m_copy, method=:epsilon, round_results=true);
			val2, t2, _ = @timed solve_stidsen(m ; time_limit=time_limit, docovercuts=false, lift_covers=false, use_nsga=false, global_nsga=false, global_branch=false, preprocess=false, fill_triangles=false);
			val1, t1, _ = @timed solve_stidsen(m ; time_limit=time_limit, docovercuts=true, lift_covers=true, use_nsga=false, global_nsga=false, global_branch=false, preprocess=false, fill_triangles=false);

			if (!(val1.timeout) && !(val2.timeout) && length(val1.YN) != length(val2.YN))
				figure(3)
				clf()
				plot(map(first, val1.YN), map(last,val1.YN), "b>", markersize="5", label=config_1)
				plot(map(first, val1.YN), map(last, val1.YN), "r<", markersize="5", label=config_2)
				plot(map(first, getY_N(m)), map(last, getY_N(m)), "go", markersize="2", label="epsilon")
				legend()
				error()
			end

			push!(res_1[n], val1.timeout ? time_limit : t1) ; push!(res_2[n], val2.timeout ? time_limit : t2) ; push!(res_eps[n], tEPS)
			!(val1.timeout) && push!(nb_nodes_1[n], val1.nodes)
			!(val2.timeout) && push!(nb_nodes_2[n], val2.nodes)

			figure(1) ; clf()
			xlabel("n") ; ylabel("#noeuds") ; title("temps et nombre de noeuds en moyenne")
			k = sort(collect(keys(res_1)))
			plot(k, [mean(nb_nodes_1[a]) for a in k], "gs", markersize=3) ; plot(k, [mean(nb_nodes_1[a]) for a in k], "g--", linewidth=1)
			plot(k, [mean(nb_nodes_2[a]) for a in k], "bs", markersize=3) ; plot(k, [mean(nb_nodes_2[a]) for a in k], "b--", linewidth=1)

			twinx() ; ylabel("temps moyen (s)")
			plot([],[],"g--",markersize=3,label="#noeuds $config_1") #labels dont work after twinx(), had to put a dummy plot
			plot([],[],"b--",markersize=3,label="#noeuds $config_2")
			plot(k, [mean(res_1[a]) for a in k], "gx", markersize=3) ; plot(k, [mean(res_1[a]) for a in k], "g-", linewidth=1, label="tmoy $config_1")
			plot(k, [mean(res_2[a]) for a in k], "gx", markersize=3) ; plot(k, [mean(res_2[a]) for a in k], "b-", linewidth=1, label="tmoy $config_2")
			if epsilon
				plot(k, [mean(res_eps[a])for a in k], "k-", linewidth=1, label="tmoy epsilon")
			end
			legend()
		end
		println("\n####################################\nn = $n")
		println("tmoy eps : $(mean(res_eps[n]))") ; println("tmoy $config_1 : $(mean(res_1[n]))") ; println("tmoy $config_2 : $(mean(res_2[n]))")
		println("moy : $(round(mean(nb_nodes_1[n]),1)) noeuds $config_1") ; println("moy : $(round(mean(nb_nodes_2[n]),1)) noeuds $config_2")
	end
end

benchmark(1:20, 10:11, 1, use_nsga=true, global_branch=false, docovercuts=true, lift_covers=true)
# inst = hard_instance2(CplexSolver(CPX_PARAM_SCRIND = 0))
srand(0)
benchmark(1:20, 10:5:60, 10, use_nsga=true, global_branch=false, docovercuts=true, lift_covers=true, time_limit=30, parragh=false)


text:
abstract type Oracle end

function initialize_oracle end
function ask_oracle end

function initialize_oracle(oracle, data::Array{T, 2}, labels::Vector{Symbol}, params::Dict{Symbol, Any}=Dict{Symbol, Any}())::Oracle where T <: Real
    if isa(oracle, DataType) && oracle <: Oracle
        return oracle(data, labels, params)
    end
    throw(ErrorException("Unknown oracle $(oracle)."))
end


text:
"""
    estimatenoise(d::NMRData)

Estimate the rms noise level in the data and update `:noise` metadata.
"""
function estimatenoise!(d::NMRData)
    α = 0.25 # fraction of data to discard for noise estimation

    nsamples = 1000
    n0 = length(d)
    step = Int(ceil(n0/nsamples))

    y = sort(vec(data(d))[1:step:end])
    n = length(y)

    # select central subset of points
    i1 = ceil(Int,(α/2)*n)
    i2 = floor(Int,(1-α/2)*n)
    y = y[i1:i2]

    μ0 = mean(y)
    σ0 = std(y)
    a = y[1]
    b = y[end]
    #histogram(y)|>display

    # MLE of truncated normal distribution
    𝜙(x) = (1/sqrt(2π))*exp.(-0.5*x.^2)
    𝛷(x) = 0.5*erfc.(-x/sqrt(2))
    logP(x,μ,σ) = @. log(𝜙((x-μ)/σ) / (σ*(𝛷((b-μ)/σ) - 𝛷((a-μ)/σ))))
    ℒ(p) = -sum(logP(y, p...))

    p0 = [μ0, σ0]
    res = optimize(ℒ, p0)
    p = Optim.minimizer(res)
    d[:noise] = abs(p[2])
end

"""
    estimatenoise(spectra::Array{<:NMRData,1})

Estimate the rms noise level and update `:noise` metadata for a list of spectra.
"""
estimatenoise!(spectra::Array{<:NMRData,1}) = [estimatenoise!(spectrum) for spectrum in spectra]


text:
function sayhi()
    string = "Hello from STABLE."
    println(string)
    return string
end

"Copy and sort a dictionary by value"
function sortbykeys(d::Dict)
    return sort!(OrderedDict(d),byvalue=false)
end

"Copy and sort a dictionary by value"
function sortbyvals(d::Dict;rev=true)
    return sort!(OrderedDict(d),byvalue=true,rev=rev)
end

"Display the data from the input files"
function display_input(InputDict)
    for sc in keys(InputDict)
        println("\n$(sc):")
        for k in keys(InputDict[sc])
            println("$(k) => $(InputDict[sc][k])")
        end
    end
    return nothing
end

"Get a Run instance for a given year from a set of Instances"
function getRunInstance(run_instances_dc::OrderedDict,run_year::Integer)
    if !(run_year in keys(run_instances_dc))
        error("The run instance year $run_year is not in the provided Instances.")
    end

    RunInstance = run_instances_dc[run_year]
    return RunInstance
end

"Flattens all STABLE configuration settings into a Dieter settings `Dict`"
function create_Dieter_settings(ST_run::STABLErun)
    dtr_settings = Dict{Symbol,Any}()

    hoursInYear = Dieter.hoursInYear

    instance = ST_run.RunInstance
    scenario = instance.InstanceScenario
    global_settings = scenario.ScenarioSettings
    run_settings = ST_run.RunSettings

    dtr_settings[:hours] = global_settings.Hours
    dtr_settings[:timestep] = global_settings.Timestep
    dtr_settings[:interest] = global_settings.DiscountRate
    dtr_settings[:cost_scaling] = global_settings.CostScaling

    dtr_settings[:periods] = round(Int,hoursInYear*(2/global_settings.Timestep))

    dtr_settings[:scen] = scenario.ScenarioIdentifier
    dtr_settings[:scen_name] = scenario.ScenarioName
    dtr_settings[:scen_number] = scenario.ScenarioNumber
    dtr_settings[:tg_scen] = scenario.TGScenarioIdentifier
    dtr_settings[:tg_scen_number] = scenario.TGScenarioNumber
    dtr_settings[:weather_year] = scenario.WeatherYear
    dtr_settings[:trace_year] = scenario.TraceYear
    dtr_settings[:battery_improve_year] = scenario.BatteryImproveYear
    # dtr_settings[:grid_following_fraction] = scenario.GridFollowFraction
    dtr_settings[:min_res] = scenario.MinimumRenewShare
    dtr_settings[:min_res_system] = scenario.MinimumRenewShare
    dtr_settings[:tech_cost_scen] = scenario.TechCostScenario
    dtr_settings[:fuel_cost_scen] = scenario.FuelCostScenario
    dtr_settings[:vpp_scen] = scenario.VPPbatteryScenario
    dtr_settings[:h2_scen] = scenario.HydrogenScenario
    dtr_settings[:h2_g2p_scen] = scenario.HydrogenGasPeaker
    dtr_settings[:FixExistingCap_flag] = scenario.FixExistingCap
    dtr_settings[:NewGasAllowed_flag] = scenario.NewGasAllowed
    dtr_settings[:NewDistillateAllowed_flag] = scenario.NewDistillateAllowed
    dtr_settings[:NoNewFossilFuels_year] = scenario.NoNewFossilFuelsYear
    dtr_settings[:tx_scen] = scenario.TxScenario

    dtr_settings[:inst_id] = instance.InstanceIdentifier
    dtr_settings[:inst_year] = instance.InstanceYear

    # dtr_settings[:scen] changed to dtr_settings[:run_id]
    dtr_settings[:run_id] = ST_run.RunIdentifier
    dtr_settings[:run_time_start] = ST_run.RunTimestampStart
    dtr_settings[:run_time_finish] = ST_run.RunTimestampFinish
    dtr_settings[:run_status] = ST_run.RunStatus
    dtr_settings[:run_note] = ST_run.RunNote
    dtr_settings[:results_id] = ST_run.ResultsIdentifier
    dtr_settings[:results_dir] = ST_run.ResultsPath

    dtr_settings[:ev_flag] = run_settings.Flags["EV"]
    if dtr_settings[:ev_flag] == false
        dtr_settings[:ev] = missing
    else
        dtr_settings[:ev] = dtr_settings[:ev_flag]
    end

    dtr_settings[:heat_flag] = run_settings.Flags["HeatLoads"]
    if dtr_settings[:heat_flag] == false
        dtr_settings[:heat] = missing
    else
        dtr_settings[:heat] = dtr_settings[:heat_flag]
    end

    dtr_settings[:h2_flag] = run_settings.Flags["H2"]  # 0 (`false`) means H2 not included, positive number means H2 included
    if dtr_settings[:h2_scen] > 0   # Overwrite h2_flag settting by using the Scenario setting (may revise in future)
        dtr_settings[:h2_flag] = true
    else
        dtr_settings[:h2_flag] = false
    end
    if dtr_settings[:h2_flag]
        # If h2_scen is 1, we just consider a single scenario of P2G H2 demand. So a setting of 1 means G2P and GasStorages are not included.
        dtr_settings[:h2] = dtr_settings[:h2_scen]
    else
        dtr_settings[:h2] = missing
    end

    dtr_settings[:ae_capacity_factor_lb] = run_settings.Parameters["AE_Capacity_Factor_LB"]
    dtr_settings[:pem_capacity_factor_lb] = run_settings.Parameters["PEM_Capacity_Factor_LB"]
    dtr_settings[:h2_p2g_capacity_factor_ub] = run_settings.Parameters["H2_P2G_Capacity_Factor_UB"]
    dtr_settings[:capacity_factor_recipH2_ub] = run_settings.Parameters["RecipH2_Capacity_Factor_UB"]
    dtr_settings[:recipH2_timing] = run_settings.Parameters["RecipH2_Timing"]

    dtr_settings[:loss_factor_tx] = run_settings.Parameters["Tx_Loss_Factor"]
    dtr_settings[:lifetime_Tx] = run_settings.Parameters["Tx_Lifetime"]
    dtr_settings[:scaling_Tx] = run_settings.Parameters["Tx_Cost_Scaling"]
    dtr_settings[:tx_exp_upper_bound] = run_settings.Parameters["Tx_Exp_UB_Factor"]

    dtr_settings[:syncon_parameters] = run_settings.Parameters["SynCon_Parameters"]
    # Remove with update:
    dtr_settings[:cost_syncon] = run_settings.Parameters["SynCon_Parameters"]["SynConFlywheel"]["CapCost"]
    dtr_settings[:lifetime_syncon] = run_settings.Parameters["SynCon_Parameters"]["SynConFlywheel"]["Lifetime"]

    dtr_settings[:capacity_syncon_SA] = run_settings.Parameters["SynCon_SA_Capacity"]
    dtr_settings[:timing_syncon_SA] = run_settings.Parameters["SynCon_SA_Timing"]

    dtr_settings[:timing_central_west] = run_settings.Parameters["Central_West_Timing"]

    dtr_settings[:timing_snowy2] = run_settings.Parameters["Snowy2_Timing"]
    dtr_settings[:capacity_snowy2] = run_settings.Parameters["Snowy2_Capacity"]
    dtr_settings[:hours_storage_snowy2] = run_settings.Parameters["Snowy2_HoursStorage"]

    # # Legacy:
    # dtr_settings[:coal_adjust] = 1;
    # dtr_settings[:peak_factor] = 2.5;
    # dtr_settings[:carbon_budget] = 1000*CarbonBudgetDict[ScenarioName,ScenarioYear]

    return dtr_settings
end

# Functions to construct STABLE configuration data structures:

function InitGlobalSettings()
    return STABLEglobalSettings(8760,2,0.06,1)
end

function InitGlobalSettings(dc::OrderedDict)
    return STABLEglobalSettings(
            dc["Hours"],
            dc["Timestep"],
            dc["DiscountRate"],
            dc["CostScaling"]
            )
end
# InitGlobalSettings(dc::OrderedDict) = InitGlobalSettings(Dict(dc))

function InitScenarios(scen_input_dc::OrderedDict,ScenarioSettings::STABLEglobalSettings)
    Scenarios_Dict = OrderedDict{String,STABLEscenario}()
    for k in keys(scen_input_dc)
        scen_dc = scen_input_dc[k]
        Scenarios_Dict[k] =
            STABLEscenario(
                ScenarioSettings,
                scen_dc["ScenarioIdentifier"],
                scen_dc["ScenarioName"],
                scen_dc["ScenarioNumber"],
                scen_dc["TGScenarioIdentifier"],
                scen_dc["TGScenarioNumber"],
                scen_dc["WeatherYear"],
                scen_dc["TraceYear"],
                scen_dc["BatteryImproveYear"],
                scen_dc["MinimumRenewShare"],
                scen_dc["TechCostScenario"],
                scen_dc["FuelCostScenario"],
                scen_dc["VPPbatteryScenario"],
                scen_dc["HydrogenScenario"],
                scen_dc["HydrogenGasPeaker"],
                scen_dc["FixExistingCap"],
                scen_dc["NewGasAllowed"],
                scen_dc["NewDistillateAllowed"],
                scen_dc["NoNewFossilFuelsYear"],
                scen_dc["TxScenario"]
            )
    end
    return Scenarios_Dict
end

function InitInstances(inst_input_dc::OrderedDict,Scenarios_Dict::OrderedDict{String,STABLEscenario})
    Instances_Dict = OrderedDict{String,OrderedDict{Integer,STABLEinstance}}()
    for inst_id in keys(inst_input_dc)
        inst_dc = inst_input_dc[inst_id]
        if !(inst_dc["ScenarioIdentifier"] in keys(Scenarios_Dict))
            error("Scenario $(inst_id) not found in the given Scenario Dict")
        end
        Instances_Dict[inst_id] = OrderedDict{Integer,STABLEinstance}()

        if !isempty(inst_dc["InstanceYears"])
            CreateYears = inst_dc["InstanceYears"]
        else
            start_year = inst_dc["InstancesStart"]
            end_year = inst_dc["InstancesEnd"]
            year_step = inst_dc["InstancesStep"]
            CreateYears = collect((start_year):(year_step):(end_year))
        end

        for year in CreateYears
            Instances_Dict[inst_id][year] =
                STABLEinstance(
                    Scenarios_Dict[inst_dc["ScenarioIdentifier"]], # STABLEscenario
                    inst_id, # InstanceIdentifier
                    year     # InstanceYear
                )
        end

    end
    return Instances_Dict
end

function InitRunSettings(runsettings_input_dc::OrderedDict)
    RunSettingsDict = OrderedDict{String,STABLErunSettings}()
    for k in keys(runsettings_input_dc)
        runsettings_dc = runsettings_input_dc[k]
        RunSettingsDict[k] =
            STABLErunSettings(
                String(k),   # RunSettingsIdentifier
                runsettings_dc["Flags"],
                runsettings_dc["Parameters"]
            )
    end
    return RunSettingsDict
end


text:
module Params
using Parameters

default = @with_kw (
    bigm = 250, # Number of stocks
    bign = 1000, # Number of investors
    bigt = 1350, # Number of days
    bigk = 250, # Number of funds

    # Market parameters
    marketstartval = 100,
    drift = 0.000307888,
    marketvol = 0.01,

    # Stock parameters
    impactrange = 0.00075:0.0001:0.0075,
    betamean = 1,
    betastd = 0.3,
    stockstartval = 100,
    stockvolrange = 0.001:0.001:0.01,

    # Investor parameters
    reviewprobability = 1/63,
    perfwindow = 1:100,
    invcaprange = (50,150),
    thresholdmean = 0,
    thresholdstd = 0.05,

    # Fund parameters
    portfsizerange = 10:100,
    plotpath = "./plots/")

end #module


text:
# #############################################################################
# diagnostic.jl
# Authors: J. von Hardenberg (ISAC-CNR, Italy)
# #############################################################################
# Description
# Example of ESMValTool diagnostic written in Julia
#
# Modification history
#    20190807-vonhardenberg_jost written for v2.0
#    20191117-vonhardenberg_jost added more realistic writing of file and plot
# ############################################################################

import YAML
import JSON
using NetCDF
# Used to write output NetCDF file with original attributes
using RainFARM
using Statistics

using PyPlot
# Avoid plotting to screen
pygui(false)

# Provides the plotmap() function
include(joinpath(dirname(@__DIR__), "shared/external.jl"))

function provenance_record(infile)
    xprov = Dict("ancestors" => infile,
                 "authors" => ["vonhardenberg_jost", "arnone_enrico"],
                 "references" => ["zhang11wcc"],
                 "projects" => ["crescendo", "c3s-magic"],
                 "caption" => "Example diagnostic in Julia",
                 "statistics" => ["other"],
                 "realms" => ["atmos"],
                 "themes" => ["phys"],
                 "domains" => ["global"])
    return(xprov)
end

function compute_diagnostic(metadata, varname, diag_base, parameter,
                            work_dir, plot_dir)
    provenance = Dict()
    for (infile, value) in metadata
        dataset = value["dataset"]
        reference_dataset = value["reference_dataset"]
        start_year = value["start_year"]
        end_year = value["end_year"]
        exp = value["exp"]
        ensemble = value["ensemble"]
        println(diag_base, ": working on file ", infile)
        println(diag_base, ": calling diagnostic with following parameters")
        println(dataset, " ", reference_dataset, " ", start_year, " ",
                end_year, " ", exp, " ", ensemble, parameter)
        # Call here actual diagnostic
        println(diag_base, ": I am your Julia diagnostic")

        # Read the variable, lon and lat
        var = ncread(infile, varname)
        lon = ncread(infile, "lon")
        lat = ncread(infile, "lat")

        units = ncgetatt(infile, varname, "units")

        # Compute time average and add parameter
        varm = mean(var, dims = 3) .+ parameter

        # Output filename
        outfile = string(work_dir, "/", varname, "_", dataset, "_", exp, "_",
                         ensemble, "_", start_year, "-",
                         end_year, "_timmean.nc")

        # Use the RainFARM function write_netcdf2d to write variable to
        # output file copying original attributes from infile
        write_netcdf2d(outfile, varm, lon, lat, varname, infile)

        # Create provenance record for the output file
        xprov = provenance_record(infile)

        # Plot the field
        plotfile = string(plot_dir, "/", varname, "_", dataset, "_", exp, "_",
                         ensemble, "_", start_year, "-",
                         end_year, "_timmean.png")
        title = string("Mean ", varname, " ", dataset, " ", exp, " ", ensemble,
                       " ", start_year, "-", end_year)
        plotmap(lon, lat, var, title = title, proj = "robinson", clabel = units)
        savefig(plotfile)
        xprov["plot_file"] = plotfile
        provenance[outfile] = xprov
    end
    return provenance
end

function main(settings)

    metadata = YAML.load_file(settings["input_files"][1])
    climofiles = collect(keys(metadata))
    climolist = metadata[climofiles[1]]
    varname = climolist["short_name"]
    diag_base = climolist["diagnostic"]

    println(diag_base, ": starting routine")
    println(diag_base, ": creating work and plot directories")
    work_dir = settings["work_dir"]
    plot_dir = settings["plot_dir"]
    run_dir = settings["run_dir"]
    mkpath(work_dir)
    mkpath(run_dir)
    mkpath(plot_dir)
    cd(run_dir)

    # Reading an example parameter from the settings
    parameter = settings["parameter1"]

    # Compute the main diagnostic
    provenance = compute_diagnostic(metadata, varname, diag_base,
                                    parameter, work_dir, plot_dir)

    # setup provenance file
    provenance_file = string(run_dir, "/diagnostic_provenance.yml")

    # Write provenance file
    open(provenance_file, "w") do io
        JSON.print(io, provenance, 4)
    end
end

settings = YAML.load_file(ARGS[1])
main(settings)

max_stars_repo_name: JGCRI/ESMValTool

max_issues_repo_name: JGCRI/ESMValTool

max_forks_repo_name: JGCRI/ESMValTool


text:
####
#### Robust online (co-)variance estimators.
####

abstract type VarEstimator{T} end

# NOTE: this naive variance estimator is used only in testing
mutable struct NaiveVar{T} <: VarEstimator{T}
    n :: Int
    S :: Vector{Vector{T}}
end

NaiveVar(::Type{T}=Float64) where {T} = NaiveVar(0, Vector{Vector{T}}())
NaiveVar(::Type{T}, ::Int) where {T} = NaiveVar(T)
NaiveVar(::Int) = NaiveVar(Float64)

function add_sample!(nc::NaiveVar, s::AbstractVector)
    nc.n += 1
    push!(nc.S, s)
end

function reset!(nc::NaiveVar{T}) where {T<:AbstractFloat}
    nc.n = 0
    nc.S = Vector{Vector{T}}()
end

function get_var(nc::NaiveVar{T}) where {T<:AbstractFloat}
    @assert nc.n >= 2 "Cannot get variance with only one sample"
    return Statistics.var(nc.S)
end

# Ref： https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/welford_var_estimator.hpp
mutable struct WelfordVar{T<:AbstractFloat, AT<:AbstractVector{T}} <: VarEstimator{T}
    n :: Int
    μ :: AT
    M :: AT
end

WelfordVar(::Type{T}, d::Int) where {T} = WelfordVar(0, zeros(T, d), zeros(T, d))
WelfordVar(d::Int) = WelfordVar(Float64, d)

function reset!(wv::WelfordVar{T, AT}) where {T<:AbstractFloat, AT<:AbstractVector{T}}
    wv.n = 0
    wv.μ .= zero(T)
    wv.M .= zero(T)
end

function add_sample!(wv::WelfordVar, s::AbstractVector)
    wv.n += 1
    @unpack μ, M, n = wv
    for i in eachindex(s)
        δ = s[i] - μ[i]
        μ[i] += δ / n
        M[i] += δ * (s[i] - μ[i])
    end
end

# https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/var_adaptation.hpp
function get_var(wv::WelfordVar{T, AT}) where {T<:AbstractFloat, AT<:AbstractVector{T}}
    n, M = T(wv.n), wv.M
    @assert n >= 2 "Cannot get covariance with only one sample"
    return (n / ((n + 5) * (n - 1))) .* M .+ T(1e-3) * (5 / (n + 5))
end

abstract type CovEstimator{T} end

# NOTE: This naive covariance estimator is used only in testing.
mutable struct NaiveCov{T} <: CovEstimator{T}
    n :: Int
    S :: Vector{Vector{T}}
end
NaiveCov(::Type{T}=Float64) where {T} = NaiveCov(0, Vector{Vector{T}}())

function add_sample!(nc::NaiveCov, s::AbstractVector)
    nc.n += 1
    push!(nc.S, s)
end

function get_cov(nc::NaiveCov{T})::Matrix{T} where {T<:AbstractFloat}
    @assert nc.n >= 2 "Cannot get covariance with only one sample"
    return Statistics.cov(nc.S)
end

function reset!(nc::NaiveCov{T}) where {T<:AbstractFloat}
    nc.n = 0
    nc.S = Vector{Vector{T}}()
end

# Ref: https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/welford_covar_estimator.hpp
mutable struct WelfordCov{T<:AbstractFloat} <: CovEstimator{T}
    n :: Int
    μ :: Vector{T}
    M :: Matrix{T}
    δ :: Vector{T} # temporary
end

function WelfordCov(::Type{T}, d::Int) where {T}
    return WelfordCov(0, zeros(T, d), zeros(T, d, d), zeros(T, d))
end
WelfordCov(d::Int) = WelfordCov(Float64, d)

function reset!(wc::WelfordCov{T}) where {T<:AbstractFloat}
    wc.n = 0
    wc.μ .= zero(T)
    wc.M .= zero(T)
end

function add_sample!(wc::WelfordCov, s::AbstractVector)
    wc.n += 1
    @unpack δ, μ, n, M = wc
    δ .= s .- μ
    μ .+= δ ./ n
    M .+= (s .- μ) .* δ'
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/covar_adaptation.hpp
function get_cov(wc::WelfordCov{T}) where {T<:AbstractFloat}
    n, M = T(wc.n), wc.M
    @assert n >= 2 "Cannot get variance with only one sample"
    return (n / ((n + 5) * (n - 1))) .* M + T(1e-3) * (5 / (n + 5)) * LinearAlgebra.I
end

###
### Preconditioning matrix adaptors.
###

abstract type AbstractPreconditioner <: AbstractAdaptor end

finalize!(adaptor::T) where {T<:AbstractPreconditioner} = nothing

####
#### Preconditioning matrix adaption implementation.
####

# Unit
struct UnitPreconditioner{T} <: AbstractPreconditioner end
Base.show(io::IO, ::UnitPreconditioner) = print(io, "UnitPreconditioner")

UnitPreconditioner(::Type{T}=Float64) where {T} = UnitPreconditioner{T}()

string(::UnitPreconditioner) = "I"
reset!(::UnitPreconditioner) = nothing
getM⁻¹(dpc::UnitPreconditioner{T}) where {T<:AbstractFloat} = UniformScaling{T}(one(T))
adapt!(
    ::UnitPreconditioner,
    ::AbstractVector{<:Real},
    ::AbstractFloat,
    is_update::Bool=true
) = nothing

mutable struct DiagPreconditioner{
    T<:Real,
    AT<:AbstractVector{T},
    TEst<:VarEstimator{T}
} <: AbstractPreconditioner
    n_min   :: Int
    ve      :: TEst
    var     :: AT
end
Base.show(io::IO, ::DiagPreconditioner) = print(io, "DiagPreconditioner")

# Diagonal
DiagPreconditioner(d::Int, n_min::Int=10) = DiagPreconditioner(Float64, d, n_min)
function DiagPreconditioner(
    ::Type{T},
    d::Int;
    n_min::Int=10,
    var=Vector(ones(T, d))
) where {T}
    ve = WelfordVar(T, d)
    return DiagPreconditioner(n_min, ve, var)
end

string(dpc::DiagPreconditioner) = string(dpc.var)
reset!(dpc::DiagPreconditioner) = reset!(dpc.ve)
getM⁻¹(dpc::DiagPreconditioner) = dpc.var

function adapt!(
    dpc::DiagPreconditioner,
    θ::AbstractVector{T},
    α::AbstractFloat,
    is_update::Bool=true
) where {T<:Real}
    resize!(dpc, θ)
    add_sample!(dpc.ve, θ)
    if dpc.ve.n >= dpc.n_min && is_update
        # Can be made inplace
        dpc.var .= get_var(dpc.ve)
    end
end

# Dense
mutable struct DensePreconditioner{
    T<:AbstractFloat,
    TEst<:CovEstimator{T}
} <: AbstractPreconditioner
    n_min :: Int
    ce    :: TEst
    covar :: Matrix{T}
end
Base.show(io::IO, ::DensePreconditioner) = print(io, "DensePreconditioner")

DensePreconditioner(d::Int, n_min::Int=10) = DensePreconditioner(Float64, d, n_min)
function DensePreconditioner(
    ::Type{T},
    d::Int;
    n_min::Int=10,
    covar=LinearAlgebra.diagm(0 => ones(T, d))
) where {T}
    ce = WelfordCov(T, d)
    # TODO: take use of the line below when we have an interface to set which covariance estimator to use
    # ce = NaiveCov(T)
    return DensePreconditioner(n_min, ce, covar)
end

string(dpc::DensePreconditioner) = string(LinearAlgebra.diag(dpc.covar))
reset!(dpc::DensePreconditioner) = reset!(dpc.ce)
getM⁻¹(dpc::DensePreconditioner) = dpc.covar

function adapt!(
    dpc::DensePreconditioner,
    θ::AbstractVector{T},
    α::AbstractFloat,
    is_update::Bool=true
) where {T<:AbstractFloat}
    resize!(dpc, θ)
    add_sample!(dpc.ce, θ)
    if dpc.ce.n >= dpc.n_min && is_update
        # Can be made inplace
        dpc.covar .= get_cov(dpc.ce)
    end
end

# Resize pre-conditioner if necessary.
Base.resize!(
    pc::UnitPreconditioner,
    θ::AbstractVector{T}
) where {T<:Real} = nothing
function Base.resize!(
    dpc::DiagPreconditioner,
    θ::AbstractVector{T}
) where {T<:Real}
    if length(θ) != length(dpc.var)
        @assert dpc.ve.n == 0 "Cannot resize a var estimator when it contains samples."
        dpc.ve = WelfordVar(T, length(θ))
        dpc.var = ones(T, length(θ))
    end
end
function Base.resize!(
    dpc::DensePreconditioner,
    θ::AbstractVector{T}
) where {T<:Real}
    if length(θ) != size(dpc.covar,1)
        @assert dpc.ce.n == 0 "Cannot resize a var estimator when it contains samples."
        dpc.ce = WelfordCov(T, length(θ))
        dpc.covar = LinearAlgebra.diagm(0 => ones(T, length(θ)))
    end
end

####
#### Preconditioning mass matrix.
####

abstract type AbstractMetric end

_string_M⁻¹(mat::AbstractMatrix, n_chars::Int=32) = _string_M⁻¹(LinearAlgebra.diag(mat), n_chars)
function _string_M⁻¹(vec::AbstractVector, n_chars::Int=32)
    s_vec = string(vec)
    l = length(s_vec)
    s_dots = " ...]"
    n_diag_chars = n_chars - length(s_dots)
    return s_vec[1:min(n_diag_chars,end)] * (l > n_diag_chars ? s_dots : "")
end

struct UnitEuclideanMetric{T} <: AbstractMetric
    M⁻¹::UniformScaling{T}
    dim::Int
end

UnitEuclideanMetric(::Type{T}, dim::Int) where {T} = UnitEuclideanMetric{T}(UniformScaling{T}(one(T)), dim)
UnitEuclideanMetric(dim::Int) = UnitEuclideanMetric(Float64, dim)

renew(ue::UnitEuclideanMetric, M⁻¹) = UnitEuclideanMetric(M⁻¹, ue.dim)

Base.length(e::UnitEuclideanMetric) = e.dim
Base.show(io::IO, uem::UnitEuclideanMetric) = print(io, "UnitEuclideanMetric($(_string_M⁻¹(ones(uem.dim))))")

struct DiagEuclideanMetric{T,A<:AbstractVector{T}} <: AbstractMetric
    # Diagnal of the inverse of the mass matrix
    M⁻¹     ::  A
    # Sqare root of the inverse of the mass matrix
    sqrtM⁻¹ ::  A
    # Pre-allocation for intermediate variables
    _temp   ::  A
end

function DiagEuclideanMetric(M⁻¹::AbstractVector{T}) where {T<:Real}
    return DiagEuclideanMetric(M⁻¹, sqrt.(M⁻¹), Vector{T}(undef, size(M⁻¹, 1)))
end
DiagEuclideanMetric(::Type{T}, D::Int) where {T} = DiagEuclideanMetric(ones(T, D))
DiagEuclideanMetric(D::Int) = DiagEuclideanMetric(Float64, D)

renew(ue::DiagEuclideanMetric, M⁻¹) = DiagEuclideanMetric(M⁻¹)

Base.length(e::DiagEuclideanMetric) = size(e.M⁻¹, 1)
Base.show(io::IO, dem::DiagEuclideanMetric) = print(io, "DiagEuclideanMetric($(_string_M⁻¹(dem.M⁻¹)))")

function Base.getproperty(dem::DiagEuclideanMetric, d::Symbol)
    return d === :dim ? size(getfield(dem, :M⁻¹), 1) : getfield(dem, d)
end

struct DenseEuclideanMetric{
    T,
    AV<:AbstractVector{T},
    AM<:AbstractMatrix{T},
    TcholM⁻¹<:UpperTriangular{T},
} <: AbstractMetric
    # Inverse of the mass matrix
    M⁻¹::AM
    # U of the Cholesky decomposition of the mass matrix
    cholM⁻¹::TcholM⁻¹
    # Pre-allocation for intermediate variables
    _temp::AV
end

function DenseEuclideanMetric(M⁻¹::AbstractMatrix{T}) where {T<:Real}
    _temp = Vector{T}(undef, size(M⁻¹, 1))
    return DenseEuclideanMetric(M⁻¹, cholesky(Symmetric(M⁻¹)).U, _temp)
end
DenseEuclideanMetric(::Type{T}, D::Int) where {T} = DenseEuclideanMetric(Matrix{T}(I, D, D))
DenseEuclideanMetric(D::Int) = DenseEuclideanMetric(Float64, D)

renew(ue::DenseEuclideanMetric, M⁻¹) = DenseEuclideanMetric(M⁻¹)

Base.length(e::DenseEuclideanMetric) = size(e.M⁻¹, 1)
Base.show(io::IO, dem::DenseEuclideanMetric) = print(io, "DenseEuclideanMetric(diag=$(_string_M⁻¹(dem.M⁻¹)))")

function Base.getproperty(dem::DenseEuclideanMetric, d::Symbol)
    return d === :dim ? size(getfield(dem, :M⁻¹), 1) : getfield(dem, d)
end

# `rand` functions for `metric` types.
function Base.rand(
    rng::AbstractRNG,
    metric::UnitEuclideanMetric{T}
) where {T}
    r = randn(rng, T, metric.dim)
    return r
end

function Base.rand(
    rng::AbstractRNG,
    metric::DiagEuclideanMetric{T}
) where {T}
    r = randn(rng, T, metric.dim)
    r ./= metric.sqrtM⁻¹
    return r
end

function Base.rand(
    rng::AbstractRNG,
    metric::DenseEuclideanMetric{T}
) where {T}
    r = randn(rng, T, metric.dim)
    ldiv!(metric.cholM⁻¹, r)
    return r
end

Base.rand(metric::AbstractMetric) = rand(GLOBAL_RNG, metric)

####
#### Preconditioner constructors
####

Preconditioner(m::UnitEuclideanMetric{T}) where {T} = UnitPreconditioner(T)
Preconditioner(m::DiagEuclideanMetric{T}) where {T} = DiagPreconditioner(T, m.dim; var=m.M⁻¹)
Preconditioner(m::DenseEuclideanMetric{T}) where {T} = DensePreconditioner(T, m.dim; covar=m.M⁻¹)

Preconditioner(m::Type{TM}, dim::Int=2) where {TM<:AbstractMetric} = Preconditioner(Float64, m, dim)
Preconditioner(::Type{T}, ::Type{TM}, dim::Int=2) where {T<:AbstractFloat, TM<:AbstractMetric} = Preconditioner(TM(T, dim))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Mineos_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Mineos")
JLLWrappers.@generate_main_file("Mineos", UUID("77b29d01-de08-56f3-96c6-9f6f21a29421"))
end  # module Mineos_jll


text:
module TestProject

using DataFrames
using CSV
using Lale
using Test

function proj_test()

   strfeat= laleoperator("Project","lale")(columns=Dict("type"=>"string"))
   numfeat = laleoperator("Project","lale")(columns=Dict("type"=>"number"))

   OneHotEncoder = laleoperator("OneHotEncoder")
   ConcatFeatures = laleoperator("ConcatFeatures","lale")
   Tree = laleoperator("DecisionTreeClassifier")
   KNN = laleoperator("KNeighborsClassifier")
   RF =  laleoperator("RandomForestClassifier")
   PCA = laleoperator("PCA")
   Hyperopt = laleoperator("Hyperopt","lale")

   df = CSV.read(joinpath(dirname(pathof(Lale)),"../demo/old/credit.csv"),DataFrame)
   y = df[!,"class"] |> collect
   X = df[:,2:end]

   train_X,train_y,test_X,test_y = Lale.train_test_split(X, y,testprop=0.20)

   prep_to_numbers =
   ((strfeat >> OneHotEncoder(handle_unknown = "ignore")) & numfeat )>> ConcatFeatures
   planned_orig = prep_to_numbers >> ( Tree | KNN)
   lopt = LalePipeOptimizer(planned_orig,max_evals = 3,cv = 3)
   laletrained = fit(lopt,train_X,train_y)
   lalepred  = Lale.transform(laletrained,test_X)
   @test score(:accuracy,lalepred,test_y) > 0.50

   pca_tree_planned = prep_to_numbers >> PCA >> RF
   lopt = LalePipeOptimizer(pca_tree_planned,max_evals = 3,cv = 3)
   laletrained = fit(lopt,train_X,train_y)
   lalepred  = Lale.transform(laletrained,test_X)
   @test score(:accuracy,lalepred,test_y) > 0.50
end
@testset "Lale Projection" begin
   proj_test()
end

end