text:
using RemoteFiles
using Dates: DateTime
using Test

rm("image.png", force=true)
rm("tmp", force=true, recursive=true)

@testset "RemoteFiles" begin
    @testset "Backends" begin
        if RemoteFiles.CURL() in RemoteFiles.BACKENDS
            download(RemoteFiles.CURL(), "https://httpbin.org/image/png", "image.png")
            @test isfile("image.png")
            rm("image.png", force=true)
        end

        if RemoteFiles.Wget() in RemoteFiles.BACKENDS
            download(RemoteFiles.Wget(), "https://httpbin.org/image/png", "image.png")
            @test isfile("image.png")
            rm("image.png", force=true)
        end

        download(RemoteFiles.Http(), "https://httpbin.org/image/png", "image.png")
        @test isfile("image.png")
        rm("image.png", force=true)

        backends = copy(RemoteFiles.BACKENDS)
        override_backends(RemoteFiles.CURL())
        @test RemoteFiles.BACKENDS == RemoteFiles.AbstractBackend[RemoteFiles.CURL()]
        override_backends(RemoteFiles.Wget(), RemoteFiles.CURL())
        @test RemoteFiles.BACKENDS == RemoteFiles.AbstractBackend[RemoteFiles.Wget(), RemoteFiles.CURL()]
        reset_backends()
        @test RemoteFiles.BACKENDS == backends
    end
    @testset "RemoteFile" begin
        r = RemoteFile("https://httpbin.org/image/png")
        @test r.file == "png"
        r = RemoteFile("https://httpbin.org/image/png", file="image.png")
        @test r.file == "image.png"

        @test_logs((:info, r"Downloading file 'image.png' from 'https://httpbin.org/image/png'."),
                   (:info, r"File 'image.png' was successfully downloaded."),
                   match_mode=:any, download(r, verbose=true))
        @test isfile(r)
        rm(r, force=true)

        r = RemoteFile("https://httpbin.org/image/png", file="image.png", dir="tmp")
        download(r)
        @test isfile(r)
        rm(r, force=true)

        @test_throws ArgumentError RemoteFile("garbage")

        r = RemoteFile("https://garbage/garbage/garbage.garbage", wait=0, retries=0)
        @test_throws DownloadError download(r)

        r = RemoteFile("https://garbage/garbage/garbage.garbage", wait=0, retries=1, try_backends=false)
        @test_throws DownloadError download(r)

        # Fail early for unexpected errors
        r = RemoteFile("https://garbage/garbage/garbage.garbage", wait=0, retries=0, failed=:warn)
        @test_throws DownloadError download(r)

        r = RemoteFile("https://httpbin.org/image/png", file="image.png", updates=:never)
        download(r)
        c1 = lastupdate(r)
        @test_logs (:info, r"File 'image.png' is up-to-date.") download(r, verbose=true)
        c2 = lastupdate(r)
        @test c1 == c2
        rm(r, force=true)

        r = RemoteFile("https://httpbin.org/image/png", file="image.png", updates=:always)
        download(r)
        c1 = lastupdate(r)
        sleep(1)
        download(r)
        c2 = lastupdate(r)
        @test c1 == c2
        rm(r, force=true)

        r = RemoteFile("https://httpbin.org/image/png", file="image.png", updates=:always)
        download(r)
        r = RemoteFile("https://garbage/garbage/garbage.garbage", file="image.png",
                       wait=1, retries=1, failed=:warn, updates=:always)
        @test_logs((:warn, r"Download of 'https://garbage/garbage/garbage.garbage' failed after 1 retries. Local file was not updated."),
                   match_mode=:any, download(r, verbose=true))
        rm(r, force=true)

        @RemoteFile r "https://httpbin.org/image/png" file="image.png"
        download(r)
        @test isfile(r)
        rm(r, force=true)

        dir = "data"
        @RemoteFile r "https://httpbin.org/image/png" file="image.png" dir=dir
        r1 = @RemoteFile "https://httpbin.org/image/png" file="image.png" dir=dir
        download(r)
        download(r1)
        @test isdir(dir)
        @test isfile(r)
        @test isfile(r1)
        rm(dir, force=true, recursive=true)
    end
    @testset "RemoteFile backends" begin
        r = RemoteFile("https://httpbin.org/image/png", backends=[RemoteFiles.Http()], file="image.png")
        @test_logs((:info, r"Downloading file 'image.png' from 'https://httpbin.org/image/png'."),
                   (:info, r"File 'image.png' was successfully downloaded."),
                   match_mode=:any, download(r, verbose=true))
        @test isfile(r)
        rm(r, force=true)

        if RemoteFiles.CURL() in RemoteFiles.BACKENDS
            r = RemoteFile("https://httpbin.org/image/png", backends=[RemoteFiles.CURL()], file="image.png")
            @test_logs((:info, r"Downloading file 'image.png' from 'https://httpbin.org/image/png'."),
                       (:info, r"File 'image.png' was successfully downloaded."),
                       match_mode=:any, download(r, verbose=true))
            @test isfile(r)
            rm(r, force=true)
        end
        if RemoteFiles.Wget() in RemoteFiles.BACKENDS
            r = RemoteFile("https://httpbin.org/image/png", backends=[RemoteFiles.Wget()], file="image.png")
            @test_logs((:info, r"Downloading file 'image.png' from 'https://httpbin.org/image/png'."),
                       (:info, r"File 'image.png' was successfully downloaded."),
                       match_mode=:any, download(r, verbose=true))
            @test isfile(r)
            rm(r, force=true)
        end
    end
    @testset "RemoteFileSets" begin
        set = RemoteFileSet("Images",
                            file1=RemoteFile("https://httpbin.org/image/png", file="image1.png"),
                            file2=RemoteFile("https://httpbin.org/image/png", file="image2.png"),
                           )
        rm(set, force=true)
        download(set)
        @test isfile(set[:file1])
        @test isfile(set[:file2])
        @test isfile(set["file1"])
        @test isfile(set["file2"])
        @test isfile(set)
        rm(set)

        @RemoteFileSet set "Images" begin
            file1 = @RemoteFile "https://httpbin.org/image/png" file="image1.png"
            file2 = @RemoteFile "https://httpbin.org/image/png" file="image2.png"
        end
        download(set)
        @test isfile(set, :file1)
        @test isfile(set, :file2)
        @test all(map(isfile, paths(set)))
        @test isfile(set)
        rm(set, :file1, force=true)
        rm(set, :file2, force=true)

        dir = "data"
        @RemoteFileSet set "Images" begin
            file1 = @RemoteFile "https://httpbin.org/image/png" file="image1.png" dir=dir
        end
        download(set)
        @test isdir(dir)
        rm(dir, force=true, recursive=true)
    end
    @testset "Updates" begin
        @test RemoteFiles.samecontent(@__FILE__, @__FILE__) == true

        updates = :helllottatimes
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test_throws ErrorException RemoteFiles.isoutdated(last, now, updates)

        updates = :mondays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :tuesdays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 7)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :wednesdays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 2, 28)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :thursdays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :fridays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :saturdays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :sundays
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 2)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :yearly
        last = DateTime(2017, 2, 28)
        now = DateTime(2018, 3, 5)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :monthly
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 2, 20)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :weekly
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 6)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == false

        updates = :daily
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 3, 1)
        @test RemoteFiles.isoutdated(last, now, updates) == true
        last = DateTime(2017, 2, 28)
        now = DateTime(2017, 2, 28)
        @test RemoteFiles.isoutdated(last, now, updates) == false
    end
end


text:
"""
    erosion_one(img::AbstractArray, bgimg::AbstractArray; figtitle=" ", ref_img=nothing, ksize= 3, thrfun = Otsu(), nclose=1, col=RGBA(1.0,0.0,0.0,0.5), showPlot=true)

## Description
Function to process and plot a sand erosion case.
## Arguments
- img: image to extract pattern
- bgimg: image to subtract from
- figtitle: figure title
- ref_img: reference image (if "nothing", then no background image is used)
- ksize: size of the kernel (Kernel.gaussian)
- thrfun: threshold algorithm (default Otsu)
- nclose: number of times to apply dilatation followed by erosion (closing)
- col: color
- showPlot: if true returns a figure
returns plot with erosion pattern given two images
## Example
```jldoctest
julia> using VisualEFM

julia> fnames = "rpm" .* string.(collect(0:90:270)) .* ".jpg";

#(Not run)
julia> im = read_image.(fnames);

julia> erosion_one(im[3], im[1], figtitle="WT at 180 rpm", ref_img=im[1])
```
"""
function erosion_one(img::AbstractArray, bgimg::AbstractArray; figtitle=" ", ref_img=nothing, ksize= 3, thrfun = Otsu(),
                    nclose=1, col=RGBA(1.0,0.0,0.0,0.5), showPlot=true)

    img_blur = apply_gaussian(img, ksize=ksize)
    bgimg_blur = apply_gaussian(bgimg, ksize=ksize)

    img_erosion = getdiffpattern(img_blur, bgimg_blur, thrfun=thrfun, nclose=nclose)
    mask = Bool.(img_erosion)

    img_erosion = Float64.(img_erosion).*col

    if !isnothing(ref_img) # include reference image on the background
        ref_img = RGBA.(ref_img, 1.0)
        ref_img[mask] .= RGBA(0.0,0.0,0.0,0.0)
        img_erosion = ref_img + img_erosion
    end

    if showPlot
        fig = plot(img_erosion, axis=nothing, bordercolor="white")
        title!(string(figtitle))
        display(fig)
    end

    return [img_erosion, mask]
end

"""
    erosion_anim(sname::String, imgs::AbstractArray, bgimgs::AbstractArray; figtitle=" ", fps=1, ref_img=nothing, ksize= 3, thrfun = Otsu(), nclose=1, col=RGBA(1.0,0.0,0.0,0.5))

## Description
Function to create an animation of the erosion process.
## Arguments
- sname: name of the file to be saved (.gif)
- imgs: array of original images
- bgimg: array of images to subtract from
- figtitle: figure title
- fps: frames per second
- ref_img: reference image (if "nothing", then no background is used)
- ksize: size of the kernel (Kernel.gaussian)
- thrfun: threshold algorithm (default Otsu)
- nclose: number of times to apply dilatation followed by erosion (closing)
- col: color
saves gif animation
## Example
```jldoctest
julia> using VisualEFM

julia> fnames = "rpm" .* string.(collect(0:90:270)) .* ".jpg";

# (Not run)
julia> im = read_image.(fnames);

julia> erosion_anim("erosion_animation.gif", im, im[1], ref_img=im[1])
```
"""
function erosion_anim(sname::String, imgs::AbstractArray, bgimgs::AbstractArray; figtitle=" ",
                    fps=1, ref_img=nothing, ksize= 3, thrfun = Otsu(),
                    nclose=1, col=RGBA(1.0,0.0,0.0,0.5))

    n = length(imgs)
    if length(figtitle) == 1 # same figtitle for everything
        figtitle = repeat([figtitle], n)
    end

    if eltype(bgimgs) <: AbstractArray
        if length(bgimgs) != length(imgs)
            throw(DimensionMismatch("bgimgs must be an array with length equals to one or with
                                    the same length of imgs array"))
        end

        anim = @animate for i in 1:n
            erosion_one(imgs[i], bgimgs[i], figtitle=figtitle[i],
                        ref_img=ref_img, ksize=ksize, thrfun=thrfun,
                        nclose=nclose, col=col, showPlot=true)
        end

    else
        anim = @animate for i in 1:n
            erosion_one(imgs[i], bgimgs; figtitle=figtitle[i],
                            ref_img=ref_img, ksize=ksize,thrfun=thrfun,
                            nclose=nclose, col=col, showPlot=true)
        end
    end

    gif(anim, sname, fps=fps)

end

"""
    erosion_colormap(imgseq::AbstractArray, U::AbstractVector; figtitle=" ", cb_title=" ", ref_img=nothing, ksize= 3, thrfun = Otsu(), nclose=1, col=cfd, showPlot=true)

## Description
Function to create heatmap from the erosion images.
## Arguments
- imgseq: array with images in sequence
- U: erosion reference (array 1D of velocities)
- cb_title = colobar title (default="")
- ref_img: reference image (if "nothing", then imgseq[1] is used)
- ksize: size of the kernel (Kernel.gaussian)
- thrfun: threshold algorithm (default Otsu)
- nclose: number of times to apply dilatation followed by erosion (closing)
- col: color
- alpha: opacity parameter
- showPlot: true to return a plot
returns a color map plot
## Example
```jldoctest
julia> using VisualEFM

julia> fnames = "rpm" .* string.(collect(0:90:270)) .* ".jpg";

# (Not run)
julia> im = read_image.(fnames);

julia> erosion_colormap(im, [90,180,270], figtitle="Color Map", cb_title="Grayscale", ref_img=im[1])
```

"""
function erosion_colormap(imgseq::AbstractArray, U::AbstractVector; figtitle=" ", cb_title=" ",
                    ref_img=nothing, ksize= 3, thrfun = Otsu(), nclose=1, col=cfd, showPlot=true)

    n = length(imgseq) - 1
    U = float(U)
    if length(U) != n
        throw(DimensionMismatch("U length should a 1 be equal to the length of imgseq minus one"))
    end

    imgs_erosion =[]

    mapU = fill(0.0,size(imgseq[1]))
    mask = fill(false, size(imgseq[1]))

    for i in 1:n
        k = i + 1
        img_temp, mask_temp = erosion_one(imgseq[k], imgseq[i], figtitle = " ", ref_img=nothing,
                                ksize=ksize, thrfun=thrfun, nclose=nclose, showPlot=false)
        push!(imgs_erosion, img_temp)

        mask = mask .| mask_temp

        mask_temp[mapU .!= 0.0] .= false # ensures to remove intersections
        mapU = mapU + Float64.(mask_temp)*U[i]

    end

    id0 = mapU .== 0.0
    mapU[id0] .= NaN

    if showPlot
        if isnothing(ref_img)
            fig = heatmap(mapU[end:-1:1,:], color=col, title=string(figtitle),
                            colorbar_title=string(cb_title), aspect_ratio=1, axis=nothing,
                            bordercolor="white")
        else
            bg = copy(ref_img) # don't change original
            bg[mask] .= RGBA(0.0,0.0,0.0,0.0)
    	    fig = plot(bg)
            heatmap!(mapU, color=col, title=string(figtitle), colorbar_title=string(cb_title),
                    axis=nothing, bordercolor="white")
        end
        display(fig)
    end

    return [mapU, ref_img]
end

ext: jl


text:
export musescore_export, second2tick

using FileIO
# using MusicVisualizations: MUSESCORE, MUSESCORE_EXISTS, test_musescore

# Temporary fix until PR to MusicVisualizations is merged
const MUSESCORE = @static if Sys.iswindows() "MuseScore3" elseif Sys.islinux() "musescore" else "mscore" end
const MUSESCORE_EXISTS = [false]

function test_musescore()
    if !MUSESCORE_EXISTS[1]
        r = try
            r = run(`$(MUSESCORE) -v`)
        catch
            false
        end
        if r == false || ((typeof(r) == Base.Process) && r.exitcode != 0)
            throw(SystemError(
            """
            The command `$(MUSESCORE) -v` did not run, which probably means that
            MuseScore is not accessible from the command line.
			Please first install MuseScore
            on your computer and then add it to your PATH."""
            ))
        end
    end
    global MUSESCORE_EXISTS[1] = true
end

"""
    musescore_export(midifilename::String)
    musescore_export(midi::MIDIFile, cleanup=true)

Use MuseScore to convert the given midifile to MP3.
If cleanup is true, the temporary midifile is deleted.
"""
function musescore_export(midifilename::String)
    MUSESCORE_EXISTS[1] || test_musescore()

    name = midifilename[1:findfirst('.', midifilename)]
    audioname = name*"mp3"

    run(`$MUSESCORE $midifilename -o $audioname`)
    @info "Exported to $audioname"

    audioname
end

function musescore_export(midi::MIDIFile, cleanup=true)
    name = tempname()
    midiname = name*".mid"
    save(midiname, midi)
    musescore_export(midiname)
    cleanup && rm(midiname)
end

"""
    musescore_export(ns::NoteSequence, cleanup=true)

Use MuseScore to convert the given `NoteSequence` to MP3.
If cleanup is true, the temporary midifile is deleted.
"""
function musescore_export(ns::NoteSequence, cleanup=true)
    midi = midifile(ns)
    musescore_export(midi, cleanup)
end

"""
    second2tick(seconds::Real, tpq::Int=DEFAULT_TPQ,  qpm::Float64=DEFAULT_QPM)

Returns a MIDI tick corresponding to the given time in seconds,
quarter notes per minute and the amount of ticks per quarter note.
"""
function second2tick(seconds::Real, tpq::Int=DEFAULT_TPQ,  qpm::Float64=DEFAULT_QPM)
    round(seconds * 1e3 / ms_per_tick(tpq, qpm))
end


text:
using APES
using Base.Test

#
# Test signal
#

## Li & Stoica fig1
spectral_lines_f = [0.04, 0.05, 0.08, 0.09,    0.16:0.005:0.35]
spectral_lines_p = vec(repmat([pi/4], length(spectral_lines_f)))
spectral_lines_a = [10, 10, 10, 2,   vec(repmat([0.5], length(spectral_lines_f)-4, 1))]

## Simulation properties
N = 128
noise_mean = 0
noise_vari = 1

## Create example signal
function create_example(N, f, a, p)
    z = vec(zeros(N, 1))
    n = 1.0:N
    for i in 1:length(f)
        z = z + a[i] * sin(2 * pi * f[i] * n + p[i])
    end
    return z
end

# For passing in known cov matrix
function covariance_matrix{T <: FloatingPoint}(s::Array{T}; M::Int=int(length(s)/2))
    N = length(s)
    L = N - M + 1
    CovMatrix = zeros(M, M)
    for l = 0:N-M
        CovMatrix += s[l+1 : 1 : l+M-1+1] * s[l+1 : 1 : l+M-1+1]'
    end
    CovMatrix ./ L
end

# For reproducible results
srand(1)

#
# Test example from Li & Stoica
#

signal = create_example(N, spectral_lines_f, spectral_lines_a, spectral_lines_p)
noise  = sqrt(noise_vari) * randn(size(signal))
sn     = signal + noise
MSE    = mean((abs(apes(sn, spectral_lines_f)) .- spectral_lines_a).^2)

@test_approx_eq abs(apes(sn, spectral_lines_f)[1:4])  [10.150160442095107,9.871766839666593,10.02406973575788,2.131756870868149]
@test_approx_eq MSE  0.10794065421715185

#
# Print pretty figures
#

#=using Gadfly=#
#=using DSP=#
#=f = 0.0:0.001:0.5=#
#=a = abs(apes(sn, f))=#
#=a = layer(x=f, y=a, Geom.line,Theme(default_color=color("green")))=#
#=f = layer(x=freq(periodogram(sn, nfft=N*16)), y=sqrt(power(periodogram(sn, nfft=N*16))/(length(sn)/2)), Geom.line,Theme(default_color=color("red")))=#
#=p = plot(a, f, Guide.title("FFT (red) APES(green)"), Guide.xlabel("Frequency"), Guide.ylabel("Amplitude"))=#
#=draw(PNG("Eg.png", 28cm, 14cm), p)=#

#
# Test at lower SNR
#

signal  = create_example(N, spectral_lines_f, spectral_lines_a, spectral_lines_p)
noise   = sqrt(noise_vari) * randn(size(signal))
sn      = signal + noise * 50
MSE     = mean((abs(apes(sn, spectral_lines_f)) .- spectral_lines_a).^2)
MSE_w_Q = mean((abs(apes(sn, spectral_lines_f, Q = covariance_matrix(noise))) .- spectral_lines_a).^2)

@test_approx_eq MSE     92.84644647200092
@test_approx_eq MSE_w_Q 9.936179958764319


text:
"""
Base type of fractional differentiation algorithms, in FractionalCalculus.jl, all of the fractional derivative algorithms belong to ```FracDiffAlg```
"""
abstract type FracDiffAlg end

const FunctionAndNumber = Union{Function, Number}

include("Caputo.jl")
include("GL.jl")
include("Hadamard.jl")
include("Riesz.jl")
include("RL.jl")

include("SymbolicDiff.jl")

"""
    fracdiff(f, α, point, h, FracDiffAlg())

```fracdiff``` is the general function for computing fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x^2, 0.5, 1, 0.001, RLDiff_Approx())
1.5044908143658473
```

There are many algorithms can be used to computing fractional derivative, please refer to our manual for more details: http://scifracx.org/FractionalCalculus.jl/dev/Derivative/derivativeapi/.
"""
function fracdiff(f, α, point, h, ::FracDiffAlg)

end

function checkfracdiff(f, α, point, h)
    point > 0 ? nothing : DomainError("Only positive points")
end

## Macros for convenient computing.
"""
    @fracdiff(f, α, point)

Return the α-order derivative of **f** at specific point.

```julia-repl
julia> @fracdiff(x->x, 0.5, 1)
1.1283791670955188
```
"""
macro fracdiff(f, α, point)
    return :(fracdiff($f, $α, $point, 0.0001, RLDiffApprox()))
end

"""
    @semifracdiff(f, point)

Return the semi-derivative of **f** at spedific point.

```julia-repl
julia> @semifracdiff(x->x, 1)
1.1283791670955188
```
"""
macro semifracdiff(f::Union{Number, Function}, point)
    return :(fracdiff($f, 0.5, $point, 0.0001, RLDiffApprox()))
end


text:
# tests

@show ob = Object(Triangle(), (1,2), (3,4), (5,))
@show ob = Object(Triangle() ; center=(1,2))
@show ob = Object(Triangle(), (1,2), (3,4), (5,))
@show ob = Object(Triangle(), (1,2), (3,4), 5)
Object(Triangle(), (1,2), (3,4), 5., 6., nothing)
Object(Triangle(), (1,2), (3,4), (5.,), 6., nothing)

ob = Object(Triangle() ; center=(1,2))
@show ob
println("here")
rotate(ob, π)
throw(0)

ob = Object(Triangle() ; center=(1,2))
translate(ob, (2,3))
translate(ob, 2, 3)
scale(ob, (2,3))
scale(ob, 2)
throw(0)

#=
abstract type Shape2d <: AbstractShape end
abstract type Shape3d <: AbstractShape end
=#

struct Ellipse <: AbstractShape end
struct Ellipsoid <: AbstractShape end

Object(Ellipse(), (0,0), (1,2))
Object(Ellipse() ; center=(3,4) )
Object(Ellipsoid(), (0,0,0), (1,2,3))
Object(Ellipsoid() ; center=(2,3,4) )

Object(Ellipse(), (0,0), (1,2), π/2, 7, nothing)
Object(Ellipse(), (0,0), (1,2), π/2)
e2 = Object(Ellipse() ; value=7)
# e2 = Object(Ellipse(), 1, 2, 3, 4, π/2, 7)
e2 * 3
e2 / 3
3 * e2
eltype(e2)

translate(e2, -1, 1)
translate(e2, (-1,1))

scale(e2, (2,2))
scale(e2, 2)

rotate(e2, -π)

using Test

#=
=#

#=
struct tmp{D,T}
    a::Dims{D}
    b::NTuple{Da,T} where Da
end
t = @inferred tmp{2,Int}((1,2), (3,))
t = @inferred tmp{3,Int}((1,2,3), (3,4))
typeof(t)

struct Rect <: ShapeType end
struct Gauss <: ShapeType end

"""
    Ellipse(center::NTuple{2,Real}, radii::NTuple{2,Real}, ϕ::Real)
    Ellipse(x, y, rx, ry, ϕ)
    Ellipse(r) (circle of radius `r`)
    circle(x,y,r) (circle of radius `r` centered at `(x,y)`)
    circle((x,y), r) ditto
    circle(r) centered at origin
Construct `Ellipse` from parameters
"""
function Ellipse(center::NTuple{2,Real}, radii::NTuple{2,Real}, ϕ::Real)
	T = promote_type(eltype.(center)..., eltype.(radii)..., eltype(ϕ))
	Shape2d{Ellipse,T}(T.(center), T.(radii), T(ϕ))
end

Ellipse(x::Real, y::Real, rx::Real, ry::Real, ϕ::Real) =
	Ellipse((x, y), (rx, ry), ϕ)

circle(r::Real) = Ellipse((0, 0), (r, r), 0)
circle(x::Real, y::Real, r::Real) = Ellipse((x, y), (r, r), 0)
circle(center::NTuple{2,Real}, r::Real) = Ellipse(center, (r, r), 0)

"""
    Rect(center::NTuple{2,Real}, width::NTuple{2,Real}, ϕ::Real)
    Rect(x, y, rx, ry, ϕ)
    square(x,y,w) (square of width `w` centered at `(x,y)`)
    square((x,y), w) ditto
    square(w) centered at origin
Construct `Rect` from parameters.
"""
function Rect(center::NTuple{2,Real}, width::NTuple{2,Real}, ϕ::Real)
	T = promote_type(eltype.(center)..., eltype.(width)..., eltype(ϕ))
	Shape2d{Rect,T}(T.(center), T.(width), T(ϕ))
end

Rect(x::Real, y::Real, wx::Real, wy::Real, ϕ::Real) =
	Rect((x, y), (wx, wy), ϕ)

square(w::Real) = Rect((0, 0), (w, w), 0)
square(x::Real, y::Real, w::Real) = Rect((x, y), (w, w), 0)
square(center::NTuple{2,Real}, w::Real) = Rect(center, (w, w), 0)

function Gauss(center::NTuple{2,Real}, width::NTuple{2,Real}, ϕ::Real)
	T = promote_type(eltype.(center)..., eltype.(width)..., eltype(ϕ))
	Shape2d{Gauss,T}(T.(center), T.(width), T(ϕ))
end

Gauss(x::Real, y::Real, wx::Real, wy::Real, ϕ::Real) =
	Gauss((x, y), (wx, wy), ϕ)

Gauss(w::Real) = Gauss((0, 0), (w, w), 0)
Gauss(x::Real, y::Real, w::Real) = Gauss((x, y), (w, w), 0)
Gauss(center::NTuple{2,Real}, w::Real) = Gauss(center, (w, w), 0)

# methods

# convert
"""
   convert(shape, T::Type)
Convert `eltype` of `shape` to `T <: Real`
"""
Base.convert(shape::Shape2d{S,O}, T::Type) where {S <: ShapeType, O <: Real} =
	Shape2d{S,T}(T.(shape.center), T.(shape.width), T.(shape.ϕ))

"""
    coordfun(shape::Shape2d, x, y)
Put coordinates `(x,y)` in canonical axes associated with `shape`.
"""
function coordfun(shape::Shape2d{Rect}, x::AbstractArray, y::AbstractArray)
	size(x) == size(y) || throw("x,y size mismatch")
	x = (x .- shape.center[1]) ./ shape.width[1]
	y = (y .- shape.center[2]) ./ shape.width[2]
	c = cos(shape.ϕ)
	s = sin(shape.ϕ)
	return (x * c + y * s, -x * s + y * c) # todo check +/-
end
coordfun(shape::Shape2d{Rect}, x::Real, y::Real) = coordfun(shape, [x], [y])

"""
    imagefun(shape::Shape2d)
Return a function of `(x,y)` for displaying the shape as an image.
"""
function imagefun(shape::Shape2d{Rect})
	return (x,y) -> imagesquare.(coordfun(shape, x, y)...)
end

imagesquare(x,y) = (abs(x) < 1/2) & (abs(y) < 1/2)
imagecirc(x,y) = abs2(x) + abs2(y) < (1/2)^2
#imagegauss(x,y) = exp(-(abs2(x) + abs2(y)) / fwhmfactor) # todo

# tests
using Test

T = Float64
e = Shape2d{Ellipse,T}((1.,2.), (3., 4.), 5.)
e = Shape2d(Ellipse(), 1, 2, 3, 4., π)
e = Shape2d(Ellipse(), (1, 2), (3, 4.), π)
r = Shape2d(Rect(), (1, 2), (3, 4.), π)
@test rotate(r, -r.ϕ).ϕ == 0
e = Ellipse(1, 2, 3, 4., π)
c0 = circle(3)
r = Rect(1, 2, 3, 4., π)
s = square(1)
c1 = translate(c0, (2,4))
c2 = scale(c1, (2,4))
f = convert(e, BigFloat)
eltype(f) == BigFloat

rxy = imagefun(r)
=#


text:
using TinyCassette

using Test: @testset, @test, @test_throws

@testset "TinyCassette" begin
    include("test_if_statement.jl")
    include("test_subtype_matching.jl")
    include("test_vararg.jl")
    include("test_infer.jl")
    include("test_generated.jl")

    # this currently doesn't work: wait until
    # we can set world-age bounds on emitted
    # code info
    #include("test_recompilation.jl")
end


text:
using Base.Test

@testset "integration tests" begin
    include("integration/sequential.jl")
end


text:
module TestPipelines

using MLJ
using Test
using Statistics
import Random.seed!
seed!(1234)

@load KNNRegressor

NN = 7
X = MLJ.table(rand(NN, 3));
y = 2X.x1 - X.x2 + 0.05*rand(NN);
Xs = source(X); ys = source(y, kind=:target)

broadcast_mode(v) = mode.(v)

@testset "linear_learning_network" begin
    t = MLJ.table
    m = MLJ.matrix
    f = FeatureSelector()
    h = OneHotEncoder()
    k = KNNRegressor()
    u = UnivariateStandardizer()
    c = ConstantClassifier()

    models = [f, h]
    lin = MLJ.linear_learning_network(Xs, nothing, nothing, nothing,
                                      models...) |> MLJ.tree
    @test lin.operation == transform
    @test lin.model == h
    @test lin.arg1.operation == transform
    @test lin.arg1.model == f
    @test lin.arg1.arg1.source == Xs
    @test lin.arg1.train_arg1.source == Xs
    @test lin.train_arg1 == lin.arg1

    models = [f, k]
    lin = MLJ.linear_learning_network(Xs, ys, nothing, nothing,
                                      models...) |> MLJ.tree
    @test lin.operation == predict
    @test lin.model == k
    @test lin.arg1.operation == transform
    @test lin.arg1.model == f
    @test lin.arg1.arg1.source == Xs
    @test lin.arg1.train_arg1.source == Xs
    @test lin.train_arg1 == lin.arg1
    @test lin.train_arg2.source == ys

    models = [m, t]
    lin = MLJ.linear_learning_network(Xs, nothing, nothing, nothing,
                                      models...) |> MLJ.tree
    @test lin.operation == t
    @test lin.model == nothing
    @test lin.arg1.operation == m
    @test lin.arg1.model == nothing
    @test lin.arg1.arg1.source == Xs

    # with learned target transformation:
    models = [f, k]
    lin = MLJ.linear_learning_network(Xs, ys, u, nothing,
                                      models...) |> MLJ.tree
    @test lin.operation == inverse_transform
    @test lin.model == u
    @test lin.arg1.operation == predict
    @test lin.arg1.model == k
    @test lin.arg1.arg1.operation == transform
    @test lin.arg1.arg1.model == f
    @test lin.arg1.arg1.arg1.source == Xs
    @test lin.arg1.arg1.train_arg1.source == Xs
    @test lin.arg1.train_arg1 == lin.arg1.arg1
    @test lin.arg1.train_arg2.operation == transform
    @test lin.arg1.train_arg2.model == u
    @test lin.arg1.train_arg2.arg1.source == ys
    @test lin.arg1.train_arg2.train_arg1.source == ys
    @test lin.train_arg1.source == ys

    # with static target transformation:
    models = [f, k]
    lin = MLJ.linear_learning_network(Xs, ys, StaticTransformer(log),
                                      StaticTransformer(exp),
                                      models...) |> MLJ.tree
    @test lin.operation == transform
    @test lin.model.f == exp
    @test lin.arg1.operation == predict
    @test lin.arg1.model == k
    @test lin.arg1.arg1.operation == transform
    @test lin.arg1.arg1.model == f
    @test lin.arg1.arg1.arg1.source == Xs
    @test lin.arg1.arg1.train_arg1.source == Xs
    @test lin.arg1.train_arg1 == lin.arg1.arg1
    @test lin.arg1.train_arg2.operation == transform
    @test lin.arg1.train_arg2.model.f == log
    @test lin.arg1.train_arg2.arg1.source == ys

    # with supervised model not at end and static target transformation:
    models = [f, c, broadcast_mode]
    lin = MLJ.linear_learning_network(Xs, ys, StaticTransformer(log),
                                      StaticTransformer(exp),
                                      models...) |> MLJ.tree
    @test lin.operation == transform
    @test lin.model.f == exp
    @test lin.arg1.operation == broadcast_mode
    @test lin.arg1.model == nothing
    @test lin.arg1.arg1.operation == predict
    @test lin.arg1.arg1.model == c
    @test lin.arg1.arg1.arg1.operation == transform
    @test lin.arg1.arg1.arg1.model == f
    @test lin.arg1.arg1.arg1.arg1.source == Xs
    @test lin.arg1.arg1.arg1.train_arg1.source == Xs
    @test lin.arg1.arg1.train_arg1 == lin.arg1.arg1.arg1
    @test lin.arg1.arg1.train_arg2.operation == transform
    @test lin.arg1.arg1.train_arg2.model.f == log
    @test lin.arg1.arg1.train_arg2.arg1.source == ys

    # with supervised model not at end and with learned target transformation:
    models = [f, c, broadcast_mode]
    lin = MLJ.linear_learning_network(Xs, ys, u, nothing,
                                      models...) |> MLJ.tree
    @test lin.operation == inverse_transform
    @test lin.model == u
    @test lin.arg1.operation == broadcast_mode
    @test lin.arg1.model == nothing
    @test lin.arg1.arg1.operation == predict
    @test lin.arg1.arg1.model == c
    @test lin.arg1.arg1.arg1.operation == transform
    @test lin.arg1.arg1.arg1.model == f
    @test lin.arg1.arg1.arg1.arg1.source == Xs
    @test lin.arg1.arg1.arg1.train_arg1.source == Xs
    @test lin.arg1.arg1.train_arg1 == lin.arg1.arg1.arg1
    @test lin.arg1.arg1.train_arg2.operation == transform
    @test lin.arg1.arg1.train_arg2.model == u
    @test lin.arg1.arg1.train_arg2.arg1.source == ys
    @test lin.arg1.arg1.train_arg2.train_arg1.source == ys
    @test lin.train_arg1.source == ys

    # build a linear network for training:
    N = MLJ.linear_learning_network(Xs, ys, u, nothing, f, k)
    fM = machine(f, Xs)
    Xt = transform(fM, Xs)

    # build the same network by hand:
    uM = machine(u, ys)
    yt = transform(uM, ys)
    kM = machine(k, Xt, yt)
    zhat = predict(kM, Xt)
    N2 = inverse_transform(uM, zhat)

    # compare predictions
    fit!(N); fit!(N2)
    yhat =  N();
    @test yhat ≈ N2()
    k.K = 3; f.features = [:x3,]
    fit!(N); fit!(N2)
    @test !(yhat ≈ N())
    @test N() ≈ N2()
    global hand_built = N();

end

## PIPELINE_PREPROCESS

F = FeatureSelector()
H = OneHotEncoder()
K = KNNRegressor()
C = ConstantClassifier()
U = UnivariateStandardizer()

m = :(MLJ.matrix)
t = :(MLJ.table)
f = :(fea=F)
h = :(hot=H)
k = :(knn=K)
c = :(cnst=C)
e = :(target=exp)
l = :(inverse=log)
u = :(target=U)

function pipe(args...)
    ret = string("Pipe(", args[1])
    for x in args[2:end]
        ret *= string(", ", x)
    end
    return Meta.parse(ret*")")
end
MLJ.pipeline_preprocess(modl, ex) =
    MLJ.pipeline_preprocess(modl, ex, missing)

ex = pipe(f, m, t, h, k, u)
out = MLJ.pipeline_preprocess(TestPipelines, ex)
@test out[1] == :Pipe
@test out[2] == [:fea, :hot, :knn, :target]
@test eval.(out[3]) == [F, H, K, U]
@test eval.(out[4]) == [F, MLJ.matrix, MLJ.table, H, K]
@test eval.(out[5]) == U
@test eval.(out[6]) == nothing
@test out[7] == :DeterministicNetwork

ex =  pipe(f, m, t, h, k, e, l)
out = MLJ.pipeline_preprocess(TestPipelines, ex)
@test out[1] == :Pipe
@test out[2] == [:fea, :hot, :knn, :target, :inverse]
@test eval.(out[3])[1:3] == [F, H, K]
@test [eval.(out[3])[4:5]...]  isa AbstractVector{<:StaticTransformer}
@test eval.(out[4]) == [F, MLJ.matrix, MLJ.table, H, K]
@test eval.(out[5]).f == exp
@test eval.(out[6]).f == log
@test out[7] == :DeterministicNetwork

ex =  pipe(f, k)
out = MLJ.pipeline_preprocess(TestPipelines, ex)
@test out[1] == :Pipe
@test out[2] == [:fea, :knn]
@test eval.(out[3]) == [F, K]
@test eval.(out[4]) == [F, K]
@test eval.(out[5]) == nothing
@test eval.(out[6]) == nothing
@test out[7] == :DeterministicNetwork

ex =  pipe(f, c)
out = MLJ.pipeline_preprocess(TestPipelines, ex, :(is_probabilistic=true))
@test out[1] == :Pipe
@test out[2] == [:fea, :cnst]
@test eval.(out[3]) == [F, C]
@test eval.(out[4]) == [F, C]
@test eval.(out[5]) == nothing
@test eval.(out[6]) == nothing
@test out[7] == :ProbabilisticNetwork

ex =  pipe(f, h)
out = MLJ.pipeline_preprocess(TestPipelines, ex)
@test out[1] == :Pipe
@test out[2] == [:fea, :hot]
@test eval.(out[3]) == [F, H]
@test eval.(out[4]) == [F, H]
@test eval.(out[5]) == nothing
@test eval.(out[6]) == nothing
@test out[7] == :UnsupervisedNetwork

ex =  :(Pipe())
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex)

# target is a function but no inverse=...:
ex =  pipe(f, m, t, h, k, e)
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex)

# inverse but no target:
ex =  pipe(f, k, l)
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex)

# target specified but no component is supervised:
ex =  pipe(f, h, u)
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex)

# function as target but no inverse:
ex =  pipe(f, m, t, h, k, e)
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex)

# is_probabilistic=true declared but no supervised models
ex = pipe(f, m, t, h)
@test_throws ArgumentError MLJ.pipeline_preprocess(TestPipelines, ex,
                                                   :(is_probabilistic=true))

## SIMPLE SUPERVISED PIPELINE WITH TARGET TRANSFORM

# test a simple pipeline prediction agrees with prediction of
# hand-built learning network built earlier:
p = @pipeline(Pipe(sel=FeatureSelector(), knn=KNNRegressor(),
                   target=UnivariateStandardizer()))
p.knn.K = 3; p.sel.features = [:x3,]
mach = machine(p, X, y)
fit!(mach)
@test MLJ.tree(mach.fitresult).arg1.model.K == 3
MLJ.tree(mach.fitresult).arg1.arg1.model.features == [:x3, ]
@test predict(mach) ≈ hand_built

# test a simple probabilistic classifier pipeline:
X = MLJ.table(rand(7,3))
y = categorical(collect("ffmmfmf"))
Xs = source(X)
ys = source(y, kind=:target)
p = @pipeline(Pipe21(hot=OneHotEncoder(),
                    cnst=ConstantClassifier()),
              is_probabilistic=true)
mach = machine(p, X, y)
fit!(mach)
@test p isa ProbabilisticNetwork
pdf(predict(mach)[1], 'f') ≈ 4/7

# test a simple deterministic classifier pipeline:
X = MLJ.table(rand(7,3))
y = categorical(collect("ffmmfmf"))
Xs = source(X)
ys = source(y, kind=:target)
p = @pipeline(Piper3(hot=OneHotEncoder(), cnst=ConstantClassifier(),
                    broadcast_mode))
mach = machine(p, X, y)
fit!(mach)
@test predict(mach) == fill('f', 7)

# test a pipeline with static transformation of target:
NN = 100
X = (x1=rand(NN), x2=rand(NN), x3=categorical(rand("abc", NN)));
y = 1000*abs.(2X.x1 - X.x2 + 0.05*rand(NN))
# by hand:
Xs =source(X); ys = source(y, kind=:target);
hot = OneHotEncoder()
hot_=machine(hot, Xs)
W = transform(hot_, Xs)
sel = FeatureSelector(features=[:x1,:x3__a])
sel_ = machine(sel, W)
Wsmall = transform(sel_, W)
z = log(ys)
knn = KNNRegressor(K=4)
knn_ = machine(knn, Wsmall, z)
zhat = predict(knn_, Wsmall)
yhat = exp(zhat)
fit!(yhat)
pred1 = yhat()
# with pipeline:
p = @pipeline Pipe4(hot=OneHotEncoder(),
                    sel=FeatureSelector(),
                    knn=KNNRegressor(),
                    target=v->log.(v),
                    inverse=v->exp.(v))
p.sel.features = [:x1, :x3__a]
p.knn.K = 4
p_ = machine(p, X, y)
fit!(p_)
pred2 = predict(p_)
@test pred1 ≈ pred2

# and another:
X = (age =    [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
height = [67.0, 81.5, 55.6, 90.0, 61.1]
p = @pipeline Pipe9(X -> coerce(X, :age=>Continuous),
                    hot = OneHotEncoder(),
                    knn = KNNRegressor(K=3),
                    target = UnivariateStandardizer())
fit!(machine(p, X, height))

# ex = :(Pipe8(X -> coerce(Dict(:age=>Continuous), X),
#           hot = OneHotEncoder(),
#           knn = KNNRegressor(K=3),
#           target = UnivariateStandardizer()))

# MLJ.pipeline_preprocess(Main, ex, missing)
# eval.(lin[3])
# eval(eval.(lin[4])[1])
# eval(lin[5])
# eval(lin[6])
# lin[7]

# MLJ.pipeline_(Main, ex, :(is_probabilistic=missing))
# lin = fit!(machine(p, X, height))

end
true


text:
@testset "BSpline" begin
    include("constant.jl")
    include("linear.jl")
    include("quadratic.jl")
    include("cubic.jl")
    include("mixed.jl")
    include("multivalued.jl")
    include("non1.jl")
end


text:
using Test
include("../chroma.jl")
include("../dtw.jl")

stransf = Matrix{Complex{Cdouble}}(undef, 2049, 1)
stransf[:,1] = zeros(2049)
stransf[42] = 1000
# stranf represents a 440Hz frequency stft frame
@test chroma(stransf, 44100., Cuint(4096), 1., 1.) ==
	[0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]'

feat = Cdouble.([1 1 1; 2 2 2; 3 3 3; 4 4 4])
feat ./= sqrt.([1 2 3 4]*[1 2 3 4]')
# dtw implementation expects normalized inputs
@test dtw(feat, feat) == ([1, 2, 3], [1, 2, 3])

max_stars_repo_name: gabriel-ss/interpretation-switcher

max_issues_repo_name: gabriel-ss/interpretation-switcher

max_forks_repo_name: gabriel-ss/interpretation-switcher


text:
using CMAES, Optim, LinearAlgebra
using Test

f(x) = sum(abs2.(x))
D = 5

@testset "FullCovMatrix" begin

    T = CMAES.FullCovMatrix
    mfit = optimize(
       f,10*rand(D).-5,CMAES.CMA(T,5;dimension=D),
       Optim.Options(f_calls_limit=5000,store_trace=true,extended_trace=true),
    )
    @test Optim.minimum(mfit) < 1e-12
end

@testset "BD_CovMatrix" begin

    T = CMAES.BD_CovMatrix
    mfit = optimize(
       f,10*rand(D).-5,CMAES.CMA(T,5;dimension=D),
       Optim.Options(f_calls_limit=5000,store_trace=true,extended_trace=true),
    )
    @test Optim.minimum(mfit) < 1e-12
end

@testset "BD_CovMatrix+Diagonal" begin

    mfit = optimize(
       f,10*rand(D).-5,CMAES.CMA{Diagonal,CMAES.BD_CovMatrix}(9,5),
       Optim.Options(f_calls_limit=5000,store_trace=true,extended_trace=true),
    )
    @test Optim.minimum(mfit) < 1e-12
end


text:
#=
Non-equidistant Fourier Transform
Uso: julia NFT.jl "input.dat" f0 ff step "output.dat"
=#

I=im
len=length
str=string

#=
function SingleFT(t::Float64, y::Float64, f::Float64)
    """
    Single iteration of fourier transform
    """
    return y*exp(-2.0*pi*I*f*t)
end
=#

SFT=(t,y,f)->y*exp(-2.0*pi*I*f*t)

filename=ARGS[1]
f0=float(ARGS[2])
ff=float(ARGS[3])
step=float(ARGS[4])
outname=ARGS[5]

file=readlines(open(filename))
X=zeros(length(file))
Y=zeros(length(file))
for i in 1:length(file)
    X[i]+=float( split(file[i])[1] )
    Y[i]+=float( split(file[i])[2] )
end
Y=Y-mean(Y)

#=
transf=[]

for freq in f0:step:ff
    FT=(t,y)->SFT(t,y,freq)
    push!(transf,sum(map(FT,X,Y)))
end
=#

open(outname, "w") do outf
    for freq in f0:step:ff
        FT=(t,y)->SFT(t,y,freq)
        s = abs(sum(map(FT,X,Y)))
        write(outf, str(freq)*"\t"*str(s)*"\n" )
    end
end

max_stars_repo_name: santiagohenao/NFT

max_issues_repo_name: santiagohenao/NFT

max_forks_repo_name: santiagohenao/NFT


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module SubFoo
thissub = "Foo1.SubFoo"
end

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia


text:
using RandomBattles
using Test, BenchmarkTools

@testset "RandomBattles.jl" begin
    static_state1 = StaticState(["lanturn" "dragonair" "froslass" "azumarill" "sableye" "marowak_alolan"])
    dynamic_state1 = DynamicState(static_state1)
    static_state2 = StaticState(["azumarill" "sableye" "marowak_alolan" "lanturn" "dragonair" "froslass"])
    dynamic_state2 = DynamicState(static_state2)
    @test sizeof(static_state1) == 168
    @test isbits(static_state1)
    @test sizeof(dynamic_state1) == 20
    @test isbits(dynamic_state1)
    res = @benchmark play_battle($dynamic_state1, $static_state1)
    @test res.allocs == 0
    N = 100000
    @test sum(battle_scores(dynamic_state1, static_state1, N)) +
        sum(battle_scores(dynamic_state2, static_state2, N)) ≈ N atol=0.02*N
end


text:
macro import_winding_huge()
    :(
        using ComplexVisual:
            CV_WindingPainter,
            cv_clear_cache, CV_2DWindingColorbarPainter,
            CV_WindingColorBarCreateData, cv_create_winding_colorbar,
            cv_setup_winding_colorbar
    )
end

"""
stores informations about a found connected component.
"""
struct CV_ConnectedComponent
    point_inside   :: Tuple{Int32, Int32}  # a position inside
    mask           :: BitMatrix            # bit-mask: 1 if pixel is inside
    winding_number :: Int32                # computed winding number
end

# {{{ Flood filling algorithm (kind of hack to find connected component)
"""
Informations for flood filling (to find connected components).
"""
struct CV_FloodFill_info{planT}
    data          :: Matrix{UInt32}
    empty_value   :: UInt32
    plan          :: planT
end

"""
scan row `y` from `lx` to `rx` to find new fill seeds and append
this new sees to the fill `plan`.
"""
function cv_fill_scan(info::CV_FloodFill_info,
        lx::Int32, rx::Int32, y::Int32)  # {{{
    data, empty_value, plan = info.data, info.empty_value, info.plan

    added = false
    for x in lx:rx
        if data[x, y] != empty_value
            added = false
        elseif !added
            push!(plan, (x, y))
            added = true
        end
    end
    return nothing
end # }}}

"""
flood fill the region (connected component) starting at `(x,y)` with
the given `value`.
"""
function cv_flood_fill(info::CV_FloodFill_info,
        x::Int32, y::Int32, value::UInt32)  # {{{

    data, empty_value, plan = info.data, info.empty_value, info.plan
    X, Y = size(data)
    o = Int32(1)
    push!(plan, (x, y))
    while !isempty(plan)
        x, y = pop!(plan)
        data[x, y] != empty_value && continue
        data[x, y] = value
        lx, rx = x, x
        while lx > o && data[lx-o, y] == empty_value
            lx -= o
            data[lx ,y] = value
        end
        while rx < X && data[rx+o, y] == empty_value
            rx += o
            data[rx, y] = value
        end
        y < Y && cv_fill_scan(info, lx, rx, y+o)
        y > o && cv_fill_scan(info, lx, rx, y-o)
    end
    return nothing
end # }}}

"""
Uses `data` (typically a black and white image) to find the connected
components/regions. To mark the regions, entries in the same regions
in `data` get the same value (starting with 0x1).

returns the number of components found.

`data` is changed.
"""
function cv_mark_connected_components!(data::Matrix{UInt32},
        empty_value::UInt32)
    # Vector used to prevent dependency on DataStructures.jl
    plan = Vector{Tuple{Int32, Int32}}() # may be optimized by using a Deque
    info = CV_FloodFill_info(data, empty_value, plan)

    compcount = UInt32(0)

    while true
        next = findfirst(value -> value == empty_value, data)
        next === nothing && break
        compcount += 0x1

        cv_flood_fill(info, Int32(next[1]), Int32(next[2]), compcount)
    end
    return compcount
end
# }}}

"""
Objects used to compute winding numbers.

This is done with a kind of "hack".

* We use a `CV_LinePainter` to draw (with `ANTIALIAS_NONE`) the
  closed curve(s) in a "computational"/auxiliary `CV_Math2DCanvas`
* Then we use a flood fill algorithm to find and fill each connected region
  (connected component) with a unique value (we use the values `0x1`, `0x2`,
  etc.)
* We then compute/find for every connected component the maximal square
  inside this component. We use this as an heuristic to find a point
  inside the component which is "far" away from the curves/boundary.
* We take the center of each of such squares to use the "axis crossing method"
  to estimate/compute the winding number.
"""
mutable struct CV_WindingHelpers
    comp_canvas  :: CV_Math2DCanvas
    con_comps    :: Vector{CV_ConnectedComponent} # found connected components
    color_dict   :: Dict{Int32, UInt32}           # colors for winding numbers
end

const cv_windingpainter_line_painter_style = cv_black → cv_op_source →
        cv_linewidth(1) → cv_antialias(Cairo.ANTIALIAS_NONE)

"""
```
CV_WindingPainter{trafoT, linePainterT <: CV_StyledPainter}
    trafo                  trafoT
    hide_winding_numbers   Dict{Int32, Bool}
    styled_line_painter    linePainterT
    helpers                CV_WindingHelpers
    cache_flag             Bool
    cache                  CV_Math2DCanvasPainterCache
```

Painter which fills connected components by a color depending on the
winding number of a curve.
"""
struct CV_WindingPainter{trafoT,
        linePainterT <: CV_StyledPainter} <: CV_2DCanvasPainter  # {{{
    trafo                 :: trafoT
    hide_winding_numbers  :: Dict{Int32, Bool}    # default: false
    styled_line_painter   :: linePainterT
    helpers               :: CV_WindingHelpers
    cache_flag            :: Bool
    cache                 :: CV_Math2DCanvasPainterCache
end

"""
```
CV_WindingPainter(trafo, closed_curves, cache_flag=true)
    trafo
    closed_curves   CV_LineSegments,
    cache_flag      Bool
```
"""
function CV_WindingPainter(trafo::trafoT, closed_curves::CV_LineSegments,
        cache_flag::Bool=true) where {trafoT}  # {{{
    line_painter = CV_LinePainter(trafo, closed_curves, true)
    helpers = CV_WindingHelpers(
        CV_Math2DCanvas(-1.0+1.0im, 1.0-1.0im, 1),
        Vector{CV_ConnectedComponent}(),
        Dict{Int32, UInt32}())
    styled_painter = cv_windingpainter_line_painter_style ↦ line_painter
    return CV_WindingPainter{trafoT, typeof(styled_painter)}(
        trafo,
        Dict{Int32, Bool}(),
        styled_painter,
        helpers,
        cache_flag,
        CV_Math2DCanvasPainterCache())
end  # }}}

function cv_clear_cache(wp::CV_WindingPainter)
    if wp.cache_flag
        wp.cache.last_canvas = nothing
    end
    return nothing
end

"""
takes care that `wp.helpers.comp_canvas` is a clone (i.e. has the same
size and resolution) as `canvas`.

make sure all pixels have the value empty_value
"""
function cv_similar_comp_canvas(wp::CV_WindingPainter,
                                canvas::CV_Math2DCanvas, empty_value::UInt32) # {{{
    comp_canvas = wp.helpers.comp_canvas
    if comp_canvas.corner_ul != canvas.corner_ul         ||
            comp_canvas.corner_lr != canvas.corner_lr    ||
            comp_canvas.resolution != canvas.resolution
        # needs to change
        cv_destroy(comp_canvas)
        comp_canvas = CV_Math2DCanvas(
            canvas.corner_ul, canvas.corner_lr, canvas.resolution)
        wp.helpers.comp_canvas = comp_canvas
    end
    csurf = wp.helpers.comp_canvas.surface
    cdata = csurf.data
    Cairo.flush(csurf)
    cdata[:] .= empty_value
    Cairo.mark_dirty(csurf)
    return nothing
end # }}}

"""
computes connected components.

This is a completely discrete algorithm. A flood fill algorithm (on the
discretized pixels) is used to get a bitmask for every connected component.

see also `CV_WindingHelpers`.
"""
function cv_compute_conn_masks_for_canvas(
        wp::CV_WindingPainter, canvas::CV_Math2DCanvas, trafo) # {{{
    # Prepare computational canvas
    empty_value = 0xffffffff
    cv_similar_comp_canvas(wp, canvas, empty_value)

    # draw curves
    ccon = cv_create_context(wp.helpers.comp_canvas)
    cv_prepare(ccon, wp.styled_line_painter.style)
    cv_paint(ccon, wp.styled_line_painter.painter)
    cv_destroy(ccon)

    # find connected components with flood fill
    cdata = wp.helpers.comp_canvas.surface.data
    comp_count = cv_mark_connected_components!(cdata, empty_value)

    # get all bitmasks
    bitmasks = map(count -> cdata .== count, 0x1:comp_count)
    return bitmasks :: Vector{BitMatrix}
end # }}}

"""
finds largest square sub matrix inside `data` with ones-entries.

`use_mem` must be a matrix of the same size as `data`.

return value is `(x, y, len)`
"""
function cv_find_largest_square(data::BitMatrix, use_mem::Matrix{UInt32}) # {{{

    use_mem[:] .= 0x0
    use_mem[end, :] .= data[end, :]
    use_mem[:, end] .= data[:, end]
    for col in size(data, 2)-1:-1:1
        for row in size(data, 1)-1:-1:1
            use_mem[row, col] = data[row, col] ?
                0x1 + min(use_mem[row, col+1], use_mem[row+1, col],
                          use_mem[row+1, col+1]) : 0x0
        end
    end
    len, pos = findmax(use_mem)
    return (Int32(pos[1]), Int32(pos[2]), Int32(len))
end # }}}

function cv_get_winding_number(lp::CV_LinePainter, trafo, z::ComplexF64) # {{{
    wnr_twice = Int32(0)
    o, t = Int32(1), Int32(2)

    for segment in lp.segments
        z1, w1 = segment[1], trafo(segment[1]) - z
        for z2 in (segment[2:end]..., segment[1])
            w2 = trafo(z2) - z
            w1 == w2 && continue
            iw1, iw2 = imag(w1), imag(w2)

            # axis crossing test
            if iw1 > 0  &&  iw2 > 0
            elseif iw1 < 0  &&  iw2 < 0
            elseif iw1 == 0 == iw2
            elseif iw1 == 0
                if real(w1) ≥ 0
                    wnr_twice += (iw2 > 0) ? +o : -o    # +½ or -½
                end
            elseif iw2 == 0
                if real(w2) ≥ 0
                    wnr_twice += (iw1 > 0) ? -o : +o    # -½ or +½
                end
            else
                if (iw1*real(w2) - iw2*real(w1))/(iw1-iw2) ≥ 0
                    wnr_twice += (iw2 > 0) ? +t : -t  # +1 or -1
                end
            end

            z1, w1 = z2, w2
        end
    end

    return cv_half(wnr_twice)
end  # }}}

function cv_compute_connected_components_for_canvas(
        wp::CV_WindingPainter, canvas::CV_Math2DCanvas, trafo) # {{{
    empty!(wp.helpers.con_comps)  # Delete old data

    bitmasks = cv_compute_conn_masks_for_canvas(wp, canvas, trafo)

    # (mis-)use comp_canvas to compute largest sqaure in every bitmask
    wnr_min, wnr_max = typemax(Int32), typemin(Int32)
    comp_canvas = wp.helpers.comp_canvas
    cdata = comp_canvas.surface.data
    lp = wp.styled_line_painter.painter
    for bitmask in bitmasks
        sx, sy, len = cv_find_largest_square(bitmask, cdata)
        px, py = sx + cv_half(len), sy + cv_half(len)   # center of square
        x, y = cv_pixel2math(comp_canvas, px, py)
        wnr = cv_get_winding_number(lp, trafo, x + 1im*y)
        wnr_min, wnr_max = min(wnr_min, wnr), max(wnr_max, wnr)
        push!(wp.helpers.con_comps,
            CV_ConnectedComponent((px, py), bitmask, wnr))
    end
    return wnr_min, wnr_max
end # }}}

function cv_compute_color_dict(wp::CV_WindingPainter,
        wnr_min, wnr_max) # {{{
    color_dict = wp.helpers.color_dict
    empty!(color_dict)
    color_dict[0] = 0xffaaaaaa
    for wnr in wnr_min:wnr_max
        wnr == 0  &&  continue
        ratio = (wnr > 0) ? wnr/wnr_max : wnr/wnr_min # in (0, 1]
        res = round(UInt32, 255*0.9*sqrt(1-ratio))
        color_dict[wnr] = (wnr > 0) ?
            0xff000000 | 0x000000ff << 8   | res << 16 | res :
            0xff000000 | 0x000000ff << 16  | res <<  8 | res
    end
    return nothing
end # }}}

"""
```
cv_paint(cc, wp)
    cc     CV_2DCanvasContext
    wp     CV_WindingPainter
```
"""
function cv_paint(cc::CV_2DCanvasContext{canvasT}, wp::CV_WindingPainter
                  ) where {canvasT <: CV_Math2DCanvas}
    canvas, cache, trafo = cc.canvas, wp.cache, wp.trafo

    if wp.cache_flag && cv_is_in_cache(cache, canvas, nothing)
    else
        wnr_min, wnr_max = cv_compute_connected_components_for_canvas(
            wp, canvas, trafo)
        cv_compute_color_dict(wp, wnr_min, wnr_max)
        if wp.cache_flag
            cache.last_canvas = canvas
            cache.last_data = nothing
        end
    end

    surface = cc.canvas.surface
    Cairo.flush(surface)
    data = surface.data

    color_dict = wp.helpers.color_dict
    for con_comp in wp.helpers.con_comps
        wnr = con_comp.winding_number
        get(wp.hide_winding_numbers, wnr, false) && continue
        data[con_comp.mask] .= color_dict[wnr]
    end
    Cairo.mark_dirty(surface)

    return nothing
end

"""
A painter for drawing a colorbar for a `CV_WindingPainter`.

For each windingnumber `wnr` the rectangle `(wnr-0.5, start, 1.0, len)`
(in the math coordinate system) is filled with the color of the winding number.
"""
struct CV_2DWindingColorbarPainter <: CV_2DCanvasPainter # {{{
    winding_painter   :: CV_WindingPainter
    start             :: Float64
    len               :: Float64
    vertical          :: Bool
end

function cv_paint(cc::CV_2DCanvasContext,
                  wcb_painter::CV_2DWindingColorbarPainter)
    canvas, ctx = cc.canvas, cc.ctx

    # fill background (for unknown colors)
    set_operator(ctx, Cairo.OPERATOR_SOURCE)
    set_source_rgb(ctx, 1, 1, 1)
    rectangle(cc.ctx,
        real(canvas.corner_ul), imag(canvas.corner_lr),
        real(canvas.corner_lr) - real(canvas.corner_ul),
        imag(canvas.corner_ul) - imag(canvas.corner_lr))
    fill(ctx)

    # now every winding number
    color_dict = wcb_painter.winding_painter.helpers.color_dict
    hide_winding_numbers = wcb_painter.winding_painter.hide_winding_numbers
    for wnr in keys(color_dict)
        if get(hide_winding_numbers, wnr, false)
            set_source_rgb(ctx, 0, 0, 0)   # black for hide
        else
            value = color_dict[wnr]
            set_source_rgba(ctx,
                ((value & 0x00ff0000) >> 16)/255,
                ((value & 0x0000ff00) >>  8)/255,
                ((value & 0x000000ff)      )/255,
                ((value & 0xff000000) >> 24)/255)
        end
        if wcb_painter.vertical
            rectangle(ctx, wcb_painter.start, wnr - 0.5, wcb_painter.len, 1.0)
        else
            rectangle(ctx, wnr - 0.5, wcb_painter.start, 1.0, wcb_painter.len)
        end
        fill(ctx)
    end
    return nothing
end

# }}}

struct CV_WindingColorBarCreateData{ccsT, dcbT, contextT, rfuncT} <: CV_CreateData # {{{
    internal_slider     :: CV_SliderCreateData{ccsT, dcbT, contextT}
    colorbar_container  :: CV_SliderContainer{ccsT, dcbT}
    container_context   :: contextT
    react_func          :: rfuncT
    colorbar_painter    :: CV_2DWindingColorbarPainter
end

function cv_create_winding_colorbar(
        pixel_width::Integer, pixel_height::Integer,
        winding_painter::CV_WindingPainter,
        wnr_min::Integer, wnr_max::Integer,
        rulers::Union{NTuple{N, CV_Ruler}, Missing}=missing;
        attach::CV_AttachType=cv_south,
        decoraction_with_layout_and_position_callback=(inner_layout, pos) ->
            cv_border(inner_layout, pos, 1)) where {N} # {{{

    if rulers === missing
        rulers=(CV_Ruler(
            cv_format_ticks("%.0f", Float64.(wnr_min:1:wnr_max)...),
            CV_TickLabelAppearance()), )
    end

    slider_data = cv_create_hslider(pixel_width, pixel_height,
        wnr_min-0.5, wnr_max+0.5, rulers;
        attach, decoraction_with_layout_and_position_callback)

    react_func = z -> begin
        wnr = round(Int32, real(z))
        if haskey(winding_painter.helpers.color_dict, wnr)
            winding_painter.hide_winding_numbers[wnr] =
                !get(winding_painter.hide_winding_numbers, wnr, false)
            return CV_Response(; redraw_flag=true)
        else
            return nothing
        end
    end

    colorbar_painter = CV_2DWindingColorbarPainter(
        winding_painter,
        0.0, imag(slider_data.slider_container.can_slider.corner_ul), false)

    return CV_WindingColorBarCreateData(
        slider_data, slider_data.slider_container,
        slider_data.container_context, react_func, colorbar_painter)
end # }}}

function cv_setup_winding_colorbar(setup::CV_SceneSetupChain,
        colorbar_data::CV_WindingColorBarCreateData,
        cont_l::CV_2DLayoutPosition)
    return cv_setup_hslider(setup, colorbar_data.internal_slider,
        cont_l, colorbar_data.colorbar_painter,
        colorbar_data.react_func;
        react_to_actionpixel_update=false,
        react_to_statepixel_update=true)
end # }}}

# }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
type _section_dataset{T}
    this::T
end

function (_this::_section_dataset)(x::Int)
    0 <= x <= 9 || throw(ArgumentError("dataset number should be within 0 - 9"))

    while length(_this.this.datasets) < x+1
        push!(_this.this.datasets, DataSet())
    end

    _this.this.current_dataset = _this.this.datasets[x+1]
    _this.this
end

type _section_setlayout{T}
    this::T
end

function (_this::_section_setlayout)(layout::String)
    _this.this.layout = layout
end

type Section
    layout::String
    datasets::Vector{DataSet}

    current_dataset::DataSet

    dataset::_section_dataset{Section}

    Section() = begin
        this = new()

        this.layout = "auto"
        this.datasets = [DataSet()]
        this.current_dataset = this.datasets[1]

        this.dataset = _section_dataset{Section}(this)

        this
    end
end


text: println("Hello from Binder!")


text:
function throwError()::Nothing
    throw(ErrorException("Error occured"))
end

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
using PyPlot

function plotmodel(modelvect; ymaxG = nothing, ymaxJ = nothing)

    dt = 0.01
    # Create a time only dataset
    dϵ = timeline(t_start = 0.0, t_end = 10.0, step = dt)
    dσ = timeline(t_start = 0.0, t_end = 10.0, step = dt)

    # calculate a strain/stress data by appling a function of time (by defalut a unit step otherwise substitute hstep(amp = 2.))
    dϵ = strainfunction(dϵ, hstep())
    dσ = stressfunction(dσ, hstep()) - stressfunction(dσ, hstep(offset = 5.0))

    colplot = ["red", "#dc7633", "#229954", "#d4ac0d", "blue"]

    fig, ax = subplots(1,2, figsize=(15,7))

    for i = 1:1:length(modelvect)
        #Relaxation modulus
        dG = modelpredict(dϵ, modelvect[i])
        ax[1].plot(dG.t, dG.σ, color=colplot[i])
        ax[1].plot([-5*dt,0,dG.t[1]],[0.0, 0.0,dG.σ[1]], color=colplot[i])

        # Creep modulus
        dJ = modelpredict(dσ, modelvect[i,1])
        ax[2].plot(dJ.t, dJ.ϵ, color=colplot[i])
        ax[2].plot([-5*dt,0,dJ.t[1]],[0.0, 0.0,dJ.ϵ[1]], color=colplot[i])
    end

    # Plot settings

    ax[1].set_xlabel("Time", fontsize = 14);
    ax[1].set_ylabel("Stress", fontsize = 14);
    ax[1].tick_params("both", labelsize=12);
    ax[1].set_title("Relaxation response", fontsize = 14);
    isnothing(ymaxJ) ? ax[1].set_ylim(bottom = -dt, top = ymaxG) : ax[1].set_ylim(bottom = -dt)
    ax[1].tick_params("both", labelsize=12)
    ax[1].tick_params(axis="x", which="minor", bottom=true)
    ax[1].grid(b=true, alpha = 0.2, axis="y", which="major");
    ax[1].grid(b=true, alpha = 0.2, axis="x", which="major");

    ax[2].set_xlabel("Time", fontsize = 14);
    ax[2].set_ylabel("Strain", fontsize = 14);
    ax[2].tick_params("both", labelsize=12);
    ax[2].set_title("Creep response", fontsize = 14);
    isnothing(ymaxJ) ? ax[2].set_ylim(bottom = -dt, top = ymaxJ) : ax[2].set_ylim(bottom = -dt)
    ax[2].tick_params("both", labelsize=12)
    ax[2].tick_params(axis="x", which="minor", bottom=true)
    ax[2].grid(b=true, alpha = 0.2, axis="y", which="major");
    ax[2].grid(b=true, alpha = 0.2, axis="x", which="major");

    dω = frequencyspec()
    fig, ax = subplots(1,1, figsize=(5,5))
    for i = 1:1:length(modelvect)
        # Storage and Loss moduli
        d_Gp = dynamicmodelpredict(dω, modelvect[i])
        ax.loglog(d_Gp.ω, d_Gp.Gp, color=colplot[i], "-")
        ax.loglog(d_Gp.ω, d_Gp.Gpp, color=colplot[i], "--")
    end
    ax.set_xlabel("Frequency", fontsize = 10);
    ax.set_ylabel("Storage (—) and Loss (- -) moduli", fontsize = 10);
    ax.tick_params("both", labelsize=8);
    ax.grid(b=true, alpha = 0.2, axis="y", which="major");
    ax.grid(b=true, alpha = 0.2, axis="x", which="major");
    ax.set_title("Frequency response", fontsize = 10);

end


text:
__precompile__()

module DeferredFutures

using AutoHashEquals
using Distributed
using Serialization: AbstractSerializer, serialize_any, serialize_cycle, serialize_type

import Distributed: AbstractRemoteRef
import Serialization: serialize

export @defer, DeferredChannel, DeferredFuture, DeferredRemoteRef, reset!

"""
`DeferredRemoteRef` is the common supertype of `DeferredFuture` and `DeferredChannel` and is
the counterpart of `$AbstractRemoteRef`.
"""
abstract type DeferredRemoteRef <: AbstractRemoteRef end

@auto_hash_equals mutable struct DeferredFuture <: DeferredRemoteRef
    outer::RemoteChannel
end

"""
    DeferredFuture(pid::Integer=myid()) -> DeferredFuture

Create a `DeferredFuture` on process `pid`. The default `pid` is the current process.

Note that the data in the `DeferredFuture` will still be located wherever it was `put!`
from. The `pid` argument controlls where the outermost reference to that data is located.
"""
function DeferredFuture(pid::Integer=myid())
    ref = DeferredFuture(RemoteChannel(pid))
    finalizer(finalize_ref, ref)
    return ref
end

"""
    show(io::IO, ref::DeferredFuture)

Print a simplified string representation of the `DeferredFuture` with its RemoteChannel
parameters.
"""
function Base.show(io::IO, ref::DeferredFuture)
    rc = ref.outer
    print(io, "$(typeof(ref).name.name) at ($(rc.where),$(rc.whence),$(rc.id))")
end

"""
    serialize(s::AbstractSerializer, ref::DeferredFuture)

Serialize a DeferredFuture such that it can de deserialized by `deserialize` in a cluster.
"""
function serialize(s::AbstractSerializer, ref::DeferredFuture)
    serialize_cycle(s, ref) && return

    serialize_type(s, DeferredFuture, true)

    serialize_any(s, ref.outer)
end

@auto_hash_equals mutable struct DeferredChannel <: DeferredRemoteRef
    outer::RemoteChannel
    func::Function  # Channel generating function used for creating the `RemoteChannel`
end

"""
    DeferredChannel(pid::Integer=myid(), num::Integer=1; content::DataType=Any) -> DeferredChannel

Create a `DeferredChannel` with a reference to a remote channel of a specific size and type.
f() is a function that when executed on `pid` must return an implementation of an
`AbstractChannel`.

The default `pid` is the current process.
"""
function DeferredChannel(f::Function, pid::Integer=myid())
    ref = DeferredChannel(RemoteChannel(pid), f)
    finalizer(finalize_ref, ref)
    return ref
end

"""
    DeferredChannel(pid::Integer=myid(), num::Integer=1; content::DataType=Any) -> DeferredChannel

Create a `DeferredChannel`. The default `pid` is the current process. When initialized, the
`DeferredChannel` will reference a `Channel{content}(num)` on process `pid`.

Note that the data in the `DeferredChannel` will still be located wherever the first piece
of data was `put!` from. The `pid` argument controls where the outermost reference to that
data is located.
"""
function DeferredChannel(pid::Integer=myid(), num::Integer=1; content::DataType=Any)
    ref = DeferredChannel(RemoteChannel(pid), ()->Channel{content}(num))
    finalizer(finalize_ref, ref)
    return ref
end

"""
    show(io::IO, ref::DeferredChannel)

Print a simplified string representation of the `DeferredChannel` with its RemoteChannel
parameters and its function.
"""
function Base.show(io::IO, ref::DeferredChannel)
    rc = ref.outer
    print(
        io,
        "$(typeof(ref).name.name)($(ref.func)) at ($(rc.where),$(rc.whence),$(rc.id))"
    )
end

"""
    serialize(s::AbstractSerializer, ref::DeferredChannel)

Serialize a DeferredChannel such that it can de deserialized by `deserialize` in a cluster.
"""
function serialize(s::AbstractSerializer, ref::DeferredChannel)
    serialize_cycle(s, ref) && return

    serialize_type(s, DeferredChannel, true)
    serialize_any(s, ref.outer)
    serialize(s, ref.func)
end

"""
    finalize_ref(ref::DeferredRemoteRef)

This finalizer is attached to both `DeferredFuture` and `DeferredChannel` on construction
and finalizes the inner and outer `RemoteChannel`s.

For more information on finalizing remote references, see the Julia manual[^1].

[^1]: [Remote References and Distributed Garbage Collection](http://docs.julialang.org/en/latest/manual/parallel-computing.html#Remote-References-and-Distributed-Garbage-Collection-1)
"""
function finalize_ref(ref::DeferredRemoteRef)
    # finalizes as recommended in Julia docs:
    # http://docs.julialang.org/en/latest/manual/parallel-computing.html#Remote-References-and-Distributed-Garbage-Collection-1

    # check for ref.outer.where == 0 as the contained RemoteChannel may have already been
    # finalized
    if ref.outer.where > 0 && isready(ref.outer)
        inner = take!(ref.outer)

        finalize(inner)
    end

    finalize(ref.outer)

    return nothing
end

"""
    reset!{T<:DeferredRemoteRef}(ref::T) -> T

Removes any data from the `DeferredRemoteRef` and allows it to be reinitialized with data.

Returns the input `DeferredRemoteRef`.
"""
function reset!(ref::DeferredRemoteRef)
    if isready(ref.outer)
        inner = take!(ref.outer)

        # as recommended in Julia docs:
        # http://docs.julialang.org/en/latest/manual/parallel-computing.html#Remote-References-and-Distributed-Garbage-Collection-1
        finalize(inner)
    end

    return ref
end

"""
    put!(ref::DeferredFuture, v) -> DeferredFuture

Store a value to a `DeferredFuture`. `DeferredFuture`s, like `Future`s, are write-once
remote references. A `put!` on an already set `DeferredFuture` throws an `Exception`.
Returns its first argument.
"""
function Base.put!(ref::DeferredFuture, val)
    if !isready(ref.outer)
        inner = RemoteChannel()
        put!(ref.outer, inner)
        put!(fetch(ref.outer), val)

        return ref
    else
        throw(ErrorException("DeferredFuture can only be set once."))
    end
end

"""
    put!(rr::DeferredChannel, val) -> DeferredChannel

Store a value to the `DeferredChannel`. If the channel is full, blocks until space is
available. Returns its first argument.
"""
function Base.put!(ref::DeferredChannel, val)
    # On the first call to `put!` create the `RemoteChannel` and `put!` it in the `Future`
    if !isready(ref.outer)
        inner = RemoteChannel(ref.func)
        put!(ref.outer, inner)
    end

    # `fetch` the `RemoteChannel` and `put!` the value in there
    put!(fetch(ref.outer), val)

    return ref
end

"""
    isready(ref::DeferredRemoteRef) -> Bool

Determine whether a `DeferredRemoteRef` has a value stored to it. Note that this function
can cause race conditions, since by the time you receive its result it may no longer be
true.
"""
function Base.isready(ref::DeferredRemoteRef)
    isready(ref.outer) && isready(fetch(ref.outer))
end

"""
    fetch(ref::DeferredRemoteRef) -> Any

Wait for and get the value of a remote reference.
"""
function Base.fetch(ref::DeferredRemoteRef)
    fetch(fetch(ref.outer))
end

"""
    wait(ref::DeferredRemoteRef) -> DeferredRemoteRef

Block the current task until a value becomes available on the `DeferredRemoteRef`. Returns
its first argument.
"""
function Base.wait(ref::DeferredRemoteRef)
    wait(ref.outer)
    wait(fetch(ref.outer))
    return ref
end

# mimics the Future/RemoteChannel indexing behaviour in Base
Base.getindex(ref::DeferredRemoteRef, args...) = getindex(fetch(ref.outer), args...)

"""
    close(ref::DeferredChannel)

Closes a `DeferredChannel`. An exception is thrown by:
* `put!` on a closed `DeferredChannel`
* `take!` and `fetch` on an empty, closed `DeferredChannel`
"""
function Base.close(ref::DeferredChannel)
    if isready(ref.outer)
        inner = fetch(ref.outer)
        close(inner)
    else
        rc = RemoteChannel()
        close(rc)

        put!(ref.outer, rc)
    end

    return nothing
end

"""
    take!(ref::DeferredChannel, args...)

Fetch value(s) from a `DeferredChannel`, removing the value(s) in the processs. Note that
`take!` passes through `args...` to the innermost `AbstractChannel` and the default
`Channel` accepts no `args...`.
"""
Base.take!(ref::DeferredChannel, args...) = take!(fetch(ref.outer), args...)

"""
    @defer Future(...)
    @defer RemoteChannel(...)

`@defer` transforms a `Future` or `RemoteChannel` construction into a 'DeferredFuture' or
'DeferredChannel' construction.
"""
macro defer(ex::Expr)
    if ex.head != :call
        throw(AssertionError("Expected expression to be a function call, but got $(ex)."))
    end

    if ex.args[1] == :Future
        return Expr(:call, :DeferredFuture, ex.args[2:end]...)
    elseif ex.args[1] == :RemoteChannel
        return Expr(:call, :DeferredChannel, ex.args[2:end]...)
    else
        throw(AssertionError("Expected RemoteChannel or Future and got $(ex.args[1])."))
    end
end

end # module


text:
######################################################################################################################
# This file replicates the NICE model runs (with and without revenue recycling) presented in Budolfson et al. (2021),
# "Climate Action With Revenue Recycling Has Benefits For Poverty, Inequality, And Wellbeing," Nature Climate Change.
######################################################################################################################

# Activate the project for the paper and make sure all packages we need
# are installed.
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# Load required Julia packages.
using Interpolations, MimiFAIR13, NLopt

# Load NICE+recycling source code.
include("MimiNICE_recycle_time_varying.jl")

# ------------------------------------------------------------------------------------------------
# NICE + REVENUE RECYCLING PARAMETERS TO CHANGE
# ------------------------------------------------------------------------------------------------

# Pure rate of time preference.
ρ =  0.015

# Elasticity of marginal utility of consumption.
η =  1.5

# Income elasticity of climate damages (1 = proportional to income, -1 = inversely proportional to income).
damage_elasticity = 1.0

# Share of recycled carbon tax revenue that each region-quintile pair receives (row = region, column = quintile)
recycle_share = ones(12,5) .* 0.2

# Should the time-varying elasticity values only change across the range of GDP values from the studies?
# true = limit calculations to study gdp range, false = allow calculations for 0 to +Inf GDP.
bound_gdp_elasticity = false

# Quintile income distribution scenario (options = "constant", "lessInequality", "moreInequality", "SSP1", "SSP2", "SSP3", "SSP4", or "SSP5")
quintile_income_scenario = "constant"

# Do you also want to perform a reference case optimization run (with no revenue recycling)?
run_reference_case = true

# Name of folder to store your results in (a folder will be created with this name).
results_folder = "base_case"

# run and SAVE, or just run
save_results = true

# ------------------------------------------------------------------------------------------------
# CHOICES ABOUT YOUR ANALYSIS & OPTIMZATION
# ------------------------------------------------------------------------------------------------

# Number of 10-year timesteps to find optimal carbon tax for (after which model assumes full decarbonization).
n_objectives = 45

# Global optimization algorithm (:symbol) used for initial result. See options at http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms
global_opt_algorithm = :GN_DIRECT_L

# Local optimization algorithm (:symbol) used to "polish" global optimum solution.
local_opt_algorithm = :LN_SBPLX

# Maximum time in seconds to run global optimization (in case optimization does not converge).
global_stop_time = 600

# Maximum time in seconds to run local optimization (in case optimization does not converge).
local_stop_time = 300

# Relative tolerance criteria for global optimization convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
global_tolerance = 1e-8

# Relative tolerance criteria for global optimization convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
local_tolerance = 1e-12

# ------------------------------------------------------------------------------------------------
# RUN EVERYTHING & SAVE KEY RESULTS
# ------------------------------------------------------------------------------------------------

# ------------------------------
# Baseline Run
# ------------------------------
# Note: model parameters not controlled above are set in the "instantiate_model_in_interface.jl file"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------
# eta=2 and rho=0.1%
# ------------------------------
ρ =  0.001; η =  2
results_folder = "eta_2_rho_01"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------
# Damage Elasticity = 0.5
# ------------------------------
ρ =  0.015; η =  1.5
damage_elasticity = 0.5
results_folder = "damage_elasticity_05"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------
# Damage Elasticity = -1
# ------------------------------
damage_elasticity = -1.0
results_folder = "damage_elasticity_neg1"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# -------------------------------------
# Bounded GDP in Elasticity Regression
# -------------------------------------
damage_elasticity = 1.0
bound_gdp_elasticity = true
results_folder = "bounded_GDP"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# -----------------------------------------------------------
# Bottom Quintile Gets 40% of Revenue and Rest Get 15% Each
# -----------------------------------------------------------
bound_gdp_elasticity = false
recycle_share = [ones(12,1) .* 0.4 ones(12,4) * 0.15]
results_folder = "bottom_recycle_40"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))

# -------------------------------------------------
# No Tax Revenue Recycled to the Bottom Quintile
# -------------------------------------------------
recycle_share = [ones(12,1) .* 0.0 ones(12,4) * 0.25]
results_folder = "no_recycling_for_bottom_quintile"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))

# ------------------------------------------------
# All Revenue to the Bottom Two Quintiles
# ------------------------------------------------
recycle_share = [ones(12,2) .* 0.5 ones(12,3) * 0.0]
results_folder = "all_revenue_to_bottom_2_quintile"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))

# ------------------------------------------------
# 30% of Tax Revenue is the Administrative Cost
# ------------------------------------------------
recycle_share = ones(12,5) .* 0.2 .* 0.7
results_folder = "admin_burden_70_percent_tax"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))

# ---------------------
# No Climate Damages
# ---------------------
recycle_share = ones(12,5) .* 0.2
results_folder = "no_climate_damages"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :slrmultiplier, bau_model[:sealeveldamages,:slrmultiplier] .* 0)
update_param!(nice, :a1, bau_model[:damages,:a1] .* 0)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 0)
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :slrmultiplier, bau_model[:sealeveldamages,:slrmultiplier] .* 0)
update_param!(nice, :a1, bau_model[:damages,:a1] .* 0)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 0)
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Double Conventional Damages (Normal SLR Damages)
# ------------------------------------------------
results_folder = "double_climate_damages"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :a1, bau_model[:damages,:a1] .* 2)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 2)
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :a1, bau_model[:damages,:a1] .* 2)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 2)
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ---------------------------------------------------------------------------
# Double Conventional Damages (Normal SLR Damages) + Damage Elasticity = -1
# ---------------------------------------------------------------------------
damage_elasticity = -1.0
results_folder = "double_climate_damages_elasticity_neg_1"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :a1, bau_model[:damages,:a1] .* 2)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 2)
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :a1, bau_model[:damages,:a1] .* 2)
update_param!(nice, :a2, bau_model[:damages,:a2] .* 2)
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (less inequality)
# ------------------------------------------------
damage_elasticity = 1
quintile_income_scenario = "lessInequality"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (more inequality)
# ------------------------------------------------
quintile_income_scenario = "moreInequality"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (SSP1)
# ------------------------------------------------
quintile_income_scenario = "SSP1"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (SSP2)
# ------------------------------------------------
quintile_income_scenario = "SSP2"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (SSP3)
# ------------------------------------------------
quintile_income_scenario = "SSP3"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (SSP4)
# ------------------------------------------------
quintile_income_scenario = "SSP4"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Inequality Projections (SSP5)
# ------------------------------------------------
quintile_income_scenario = "SSP5"
results_folder = "time_varying_consumption_shares_"*quintile_income_scenario
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------------------
# Use More Basic Consumption Share Studies in Meta-Regression
# ------------------------------------------------------------
quintile_income_scenario = "constant"
results_folder = "expenditure_share_based_studies"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
meta_intercept, meta_slope = meta_regression(elasticity_studies[(elasticity_studies.Type .== "Direct") .| (elasticity_studies.Type .== "EMI"),:])
update_param!(nice, :elasticity_intercept, meta_intercept)
update_param!(nice, :elasticity_slope, meta_slope)
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
meta_intercept, meta_slope = meta_regression(elasticity_studies[(elasticity_studies.Type .== "Direct") .| (elasticity_studies.Type .== "EMI"),:])
update_param!(nice, :elasticity_intercept, meta_intercept)
update_param!(nice, :elasticity_slope, meta_slope)
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ----------------------------------------------------
# Use CGE and Input-Output Studies in Meta-Regression
# ----------------------------------------------------
results_folder = "CGE_IO_studies"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
meta_intercept, meta_slope = meta_regression(elasticity_studies[(elasticity_studies.Type .== "CGE") .| (elasticity_studies.Type .== "IO"),:])
update_param!(nice, :elasticity_intercept, meta_intercept)
update_param!(nice, :elasticity_slope, meta_slope)
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
meta_intercept, meta_slope = meta_regression(elasticity_studies[(elasticity_studies.Type .== "CGE") .| (elasticity_studies.Type .== "IO"),:])
update_param!(nice, :elasticity_intercept, meta_intercept)
update_param!(nice, :elasticity_slope, meta_slope)
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------
# Optimal Global Revenue Recycling
# ------------------------------------------------
results_folder = "optimal_equal_global_recycling"
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :global_recycle_share, ones(12))
include(joinpath("replication_helpers", "optimize_recycle.jl"))
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
update_param!(nice, :global_recycle_share, ones(12))
include(joinpath("replication_helpers", "optimize_reference.jl"))

# ------------------------------------------------------
# ------------------------------------------------------
# Add iterative optimzation with FAIRv1.3 Climate Model
# ------------------------------------------------------
# ------------------------------------------------------

# Reset model parameters to default values.
ρ = 0.015
η = 1.5
damage_elasticity = 1.0
recycle_share = ones(12,5) .* 0.2
bound_gdp_elasticity = false
quintile_income_scenario = "constant"

# Reset optimization parameters to default values
# Note: FAIR optimziation just uses a local optimization due to computational burden of iterative approach
run_reference_case = true
n_objectives = 45
n_fair_loops = 4
local_opt_algorithm = :LN_SBPLX
local_stop_time = 800
local_tolerance = 1e-11
results_folder = "FAIR_climate_model"

# Initialize BAU and base version of NICE with revenue recucling.
include(joinpath("replication_helpers", "instantiate_model_in_interface.jl"))
include(joinpath("replication_helpers", "optimize_with_FAIR.jl"))

# Reload NICE+recycling model.
include("MimiNICE_recycle_time_varying.jl")

# ------------------------------------------------------
# ------------------------------------------------------
# 2°C Recycling Runs (non-optimized)
# ------------------------------------------------------
# ------------------------------------------------------

# Get baseline instance of the model.
nice = create_nice_recycle()
run(nice)
backstopprices = nice[:emissions,:pbacktime] * 1000;

# Load 2°C carbon tax values and resulting mitiation rates.
twodegCP = DataFrame(load(joinpath(@__DIR__, "..", "data", "two_degree_carbon_tax.csv")))[!,1]
twodegMIU = Matrix(DataFrame(load(joinpath(@__DIR__, "..", "data", "two_degree_carbon_mitigation.csv"))))

# Set 2°C rates and prices and run model.
prices, rates = mu_from_tax(twodegCP[2:end,1], backstopprices);
update_param!(nice, :MIU, rates)
update_param!(nice, :global_carbon_tax, prices)
update_param!(nice, :min_study_gdp, 1e-10)
update_param!(nice, :max_study_gdp, +Inf)
run(nice)

# Get an instance of the BAU no-policy model.
bau_model = create_nice_recycle()
update_param!(bau_model, :MIU, rates * 0)
update_param!(bau_model, :min_study_gdp, 1e-10)
update_param!(bau_model, :max_study_gdp, +Inf)
run(bau_model)

# Save the results
output_directory = joinpath(@__DIR__, "..", "results", "two_degree_constraint")
mkpath(output_directory)
save_nice_recycle_results(nice, bau_model, prices, prices, output_directory, revenue_recycling=true)

# Save the no reycle runs.
update_param!(nice, :global_carbon_tax, prices .* 0)
run(nice)
save_nice_recycle_results(nice, bau_model, prices, prices, output_directory, revenue_recycling=false)

# Save the BAU runs.
output_directory = joinpath(@__DIR__, "..", "results", "bau_no_policy_at_all")
mkpath(output_directory)
save_nice_recycle_results(bau_model, bau_model, prices * 0, prices * 0, output_directory, revenue_recycling=false)

# ---------------------------------------------------------------------------------------
# Carry out runs with only CGE and IO studies, and L-shaped elasticity-income dependence
# ---------------------------------------------------------------------------------------
nice = create_nice_recycle()
update_param!(nice, :MIU, rates)
update_param!(nice, :global_carbon_tax, prices)

# First, perform a meta-regression based on study results to calculate elasticity vs. ln gdp per capita relationship.
meta_intercept, meta_slope = meta_regression(elasticity_studies[elasticity_studies.Type .!= "Direct",:])
update_param!(nice, :elasticity_intercept, 0.7)
update_param!(nice, :elasticity_slope, 0.0)
run(nice)

# Save recyling results.
output_directory = joinpath(@__DIR__, "..", "results", "two_degree_constraint_CGE_IO_studies_only")
mkpath(output_directory)
save_nice_recycle_results(nice, bau_model, prices, prices, output_directory, revenue_recycling=true)

# Save no-recycling results.
update_param!(nice, :global_carbon_tax, prices .* 0)
run(nice)
save_nice_recycle_results(nice, bau_model, prices, prices, output_directory, revenue_recycling=false)

println("All done!")


text:
# HSL_MA97 package-specific dependencies.

# Users should place hsl_ma97-x.y.z.tar.gz in ~/.julia/v0.j/HSL/deps/downloads
libhsl_ma97 = library_dependency("libhsl_ma97", depends=[libblas, liblapack, libmetis4])
# libhsl_ma97 = library_dependency("libhsl_ma97", depends=[libopenblas, libmetis4])

hsl_ma97_depsdir = BinDeps.depsdir(libhsl_ma97)
hsl_ma97_srcdir = joinpath(hsl_ma97_depsdir, "src")
hsl_ma97_src = joinpath(hsl_ma97_srcdir, "hsl_ma97-$hsl_ma97_version")
hsl_ma97_prefix = joinpath(hsl_ma97_depsdir, "usr")
hsl_ma97_libdir = joinpath(hsl_ma97_depsdir, "usr", "lib")

provides(SimpleBuild,
         (@build_steps begin
            BinDeps.ChecksumValidator(hsl_ma97_sha256, hsl_ma97_archive)
            CreateDirectory(hsl_ma97_srcdir)
            FileUnpacker(hsl_ma97_archive, hsl_ma97_srcdir, "hsl_ma97-$hsl_ma97_version")
            (@build_steps begin
              ChangeDirectory(hsl_ma97_src)
              `./configure F77=gfortran CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --prefix=$hsl_ma97_prefix --with-blas="$blas_lib" --with-lapack="$lapack_lib" --with-metis="$metis_libpath"`
              `make install`
              `gfortran -fPIC -shared -Wl,$all_load $hsl_ma97_libdir/libhsl_ma97.a $(split(blas_lib)) $(split(lapack_lib)) $(split(metis_libpath)) -lgomp -Wl,$noall_load -o $hsl_ma97_libdir/libhsl_ma97.$so`
            end)
          end), libhsl_ma97, os=:Unix)

hsl_modules[:libhsl_ma97] = :libhsl_ma97


text:
mpSim = MPsim()

attributes = Attribute.( ["atA", "atB", "atC"] )
setPossibleAttributeValues!( attributes[1], ["A1", "A2"] )
setInitialAttributeValues!( attributes[2], Dict( "B1" => 5, "B2" => 5 ) )
setPossibleAttributeValues!( attributes[3], ["C1", "C2", "C3"] )
setSimulationAttributes!( mpSim, attributes )

nodes = BaseNode.( ["A", "B", "C", "D", "E"] )
setNodeRequirements!( nodes[1], Dict( "atA" => "A1", "atC" => "C1" ) )
setNodeRequirements!( nodes[2], Dict( "atA" => "A2", "atC" => "C1" ) )
setNodeRequirements!( nodes[3], Dict( "atA" => "A1", "atC" => "C2" ) )
setNodeRequirements!( nodes[4], Dict( "atA" => "A2", "atC" => "C2" ) )
setNodeRequirements!( nodes[5], Dict( "atC" => "C3" ) )
setSimulationBaseNodes!( mpSim, nodes )

recruitments = Recruitment.( ["EW", "EW"] )
setRecruitmentTarget!.( recruitments, ["A", "B"] )
setRecruitmentFixed!.( recruitments, 15 )
setRecruitmentSchedule!.( recruitments, 12 )
setSimulationRecruitment!( mpSim, recruitments )

transitions = [Transition( "B+", "A", "C" ), Transition( "B+", "B", "D" ),
    Transition( "Adv", "C", "E" ), Transition( "Adv", "D", "E" ),
    Transition( "B-", "A" ), Transition( "B-", "B" )]
setTransitionSchedule!.( transitions, 12 )
setTransitionFluxLimits!.( transitions[1:4], [10, 10, 2, 2], [10, 10, 2, 2] )
addTransitionCondition!.( transitions,
    MPcondition.( "tenure", ==, [24, 24, 60, 60, 24, 24] ) )
setSimulationTransitions!( mpSim, transitions )

retirement = Retirement()
setRetirementSchedule!( retirement, 12 )
setRetirementCareerLength!( retirement, 120 )
setSimulationRetirement!( mpSim, retirement )

setSimulationLength!( mpSim, 600 )

attritions = Attrition.( ["A/B", "C/D", "E"] )
setAttritionPeriod!.( attritions, 12 )
setAttritionRate!.( attritions, [0.2, 0.1, 0.05] )
addSimulationAttrition!( mpSim, attritions... )