text:
# -*- coding: utf-8 -*-
# + {}
module Transformations
import Base: ∘, show, convert, promote, one, zero, inv, *, ^, -

using SemanticModels
using SemanticModels.ExprModels.Parsers

export Transformation, ConcatTransformation, Product, Pow

postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

abstract type Transformation end

struct ConcatTransformation
    seq::Vector{Transformation}
end

promote(a::ConcatTransformation, b::Transformation) = (a, ConcatTransformation([b]))
convert(::ConcatTransformation, t::Transformation) = ConcatTransformation([t])

function ∘(f::ConcatTransformation, g::ConcatTransformation)
    append!(g.seq, f.seq)
    return g
end

function (f::ConcatTransformation)(m::AbstractModel)
    return foldl((m, t)->t(m), [m; f.seq])
end

"""    Pow{T}

Pow{T} represents raising terms in an equation to a power. The type Pow{Int} forms
a group that is isomorphic to the Integers (Z, +, 0). The group operation is ∘ to
match the notation that they are compositions of functions that act on models.

This is an example of encoding meaning into the type system. We are saying
"treat this number `inc` as a function on models, but also remember that instances
of this type form a group isomorphic to the integers.
"""
struct Pow{T} <: Transformation
    inc::T
end

# this is the group operation for the Pow transformation group
∘(p::Pow, q::Pow) = Pow(p.inc + q.inc)

# we implement these algebra functions to make Pow{Int} feel more like Int.
one(::Type{Pow}) = Pow(1)
zero(::Type{Pow}) = Pow(0)
inv(p::Pow) = Pow(-p.inc)

"""    ^(p::Transformation, n::Int)

the universal definition of ^ as repeated application of ∘.
"""
^(p::Transformation, n::Int) = begin
    n >= 0 || error("cannot compute negative powers in monoid/group")
    n == 0 && return zero(p)
    s = p
    for i in 1:n-1
        s = p∘s
    end
    return s
end

function show(io::IO, p::Pow)
    write(io, "Pow($(p.inc))")
end

struct Product{T} <: Transformation
    dims::T
end

function (f::Product)(m::AbstractModel)
    return foldl((m, t)->t(m), [m; f.dims])
end

function show(io::IO, p::Product)
    write(io, "$(p.dims)")
end

∘(p::Transformation, args...) = ∘(promote(p, args...)...)
∘(p::Product, q::Product) = Product(p.dims .∘ q.dims)
promote(p::Product{T}, t::T) where T= (p, Product(t))

# TODO make tests
Product((Pow(1),Pow(2))) ∘ (Pow(1),Pow(2))

end
# -


text:
function Sum()
    return Sum(())
end

function Sum(arg0::Sum)
    return Sum((Sum,), arg0)
end

function aggregate(obj::Sum, arg0::Sum)
    return jcall(obj, "aggregate", void, (Sum,), arg0)
end

function clear(obj::Sum)
    return jcall(obj, "clear", void, ())
end

function copy(obj::Sum)
    return jcall(obj, "copy", Sum, ())
end

function evaluate(obj::Sum, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end

function evaluate(obj::Sum, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::Sum)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::Sum)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::Sum, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end

ext: jl


text:
"""
Package: ForecastPlots

    ccf(x1::{AbstractVector,DataFrame},
        x2::{AbstractVector,DataFrame};
        type = "cor",
        lag = Integer(ceil(10*log10(length(x1)))),
        alpha = (0.95,0.99);
        plot = true)

Compute the cross-correlation or cros-covariance of two univariate series with an optional plot

The results are normalized to preserve homoscedasticity. The distribution used to normalize the data is an approximation of a Fisher Transformation via a Normal Distribution. There is a plot recipe for the returned object, if the type is `cor` the plot will also show confidence intervals for the given alpha values.

If, for a given integer `k`, `x2` repeats `x1` values such that x1[t] = x2[t+k] for all `i` then high correlation value will be placed *at the right from the center* in the results. That is, this convention will be represented in the plots as `x1_t = x2_{t+k} -> _____0__k__` meaning x2 behavior can be predicted by x1 in k units.

# Arguments
- `x1`: Vector or uni-dimensional DataFrame of data.
- `x2`: Vector or uni-dimensional DataFrame of data.
- `type`: Valid values are "cor" for correlation (default) and "cov" for convariance.
- `lag`: Maximum number of lags.
- `alpha`: A tuple with two thresholds (t1,t2) with t1 <= t2 to plot confidence intervals. The default values are 0.95 and 0.99.
- `plot`: When true displays a plot with the results.

# Returns

Cros-correlation vector with an optional plot

# Examples
```julia-repl
x1 = rand(100);
x2 = circshift(x1,6);
ccf(x1, x2; type="cor");
ccf(x1, x2; type="cor", plot = false)
41-element Vector{Float64}:
  0.09861519494432992
 -0.04384418688776631
 -0.1900182513825246
  [...]
```
"""
function ccf(x1::AbstractVector{<:A},
             x2::AbstractVector{<:B};
             type::String = "cor",
             lag::Integer = Integer(ceil(10*log10(length(x1)))),
             alpha::Tuple{AbstractFloat,AbstractFloat} = (0.95,0.99),
             plot::Bool = true,
             kw...) where {A<:Real, B<:Real}

    N = length(x1)
    @assert N == length(x2) "Vectors should be of equal size"
    @assert N >= 4 "Vectors should have at least four values"
    @assert isnothing(findfirst(isnan,x1)) "No missing values allowed"
    @assert isnothing(findfirst(isnan,x2)) "No missing values allowed"
    @assert isnothing(findfirst(ismissing,x1)) "No missing values allowed"
    @assert isnothing(findfirst(ismissing,x2)) "No missing values allowed"
    @assert type in  ["cor","cov"] "The options for `type` are `cor` and `cov`"
    @assert 1 <= lag <= N-4
    @assert length(alpha) == 2
    @assert 0.0 < alpha[1] < alpha[2] < 1.0

    # Promote to same Float type
    x1,x2,_ = Base.promote(x1,x2,[Base.promote_op(/,A,B)(1)])

    # auto-ccf
    auto = x1 == x2

    ft = (type == "cor" ? Statistics.cor : Statistics.cov)

    x = []

    for i in 1:lag
        push!(x,ft(x1[i+1:N],x2[1:(N-i)]))
    end

    if auto
        kx = collect(N-1:-1:N-lag)
    else
        x = vcat(ft(x1,x2),x) # adding central value (lag=0)
        x = reverse(x)
        for i in 1:lag
            push!(x,ft(x2[i+1:N],x1[1:(N-i)]))
        end
        kx = vcat(collect(N-lag:N-1), N, collect(N-1:-1:N-lag))
    end

    # ccf normalized with the standard error for the
    # Normal distsribution of an approximation for
    # the Fisher transformation
    function fse(v::Real)::AbstractFloat
        1/sqrt(v-3)
    end
    k = fse.(kx)/fse(N)
    cres = x ./ k

    call = "ccf("* (auto ? "x" : "x1, x2") *
        "; type=\""*type*
        "\", lag="*string(lag)*
        ", alpha="*string(alpha)*")"

    # Confidence Intervals
    a1 = alpha[1]
    a2 = alpha[2]
    z1 = Distributions.quantile(Normal(), a1 + (1-a1)/2)
    ci1 = z1*fse(N)
    z2 = Distributions.quantile(Normal(), a2 + (1-a2)/2)
    ci2 = z2*fse(N)
    ci = (ci1,ci2)

    plot && display(cf_gr((cres, N=N, type, lag, alpha, ci, auto); kw...))

    cres

end

function cf_gr(args; kw...)

    # plot configuration

    titlefont = font(8, "Courier")
    xtickfont = font(3, "Courier")
    ytickfont = font(3, "Courier")
    yguidefont = font(5, "Courier")
    xguidefont = font(5, "Courier")

    N = args.N
    if split(args.type,'_')[1] == "pacf"
        type = (args.type) == "pacf_step" ? "Stepwise PACF" : "Real PACF"
        type = (args.type) == "pacf_stepwise-real" ? "Stepwise and Real PACF   |s|r|" : type
        auto_prefix = ""
    else
        type = (args.type == "cor") ? "Correlation" : "Covariance"
        auto_prefix = args.auto ? "Auto-" : "Cross-"
    end

    yguide = auto_prefix*type
    xguide = "Lag"

    # Center ticks
    lr = length(args.cres)
    gap = div(lr,11)
    mp = div(lr+1,2)
    rs = mp:gap:lr
    ls = mp-gap:-gap:1

    if args.auto
        xticks = 1:gap:lr
    else
        xt = vcat(reverse(collect(ls)),collect(rs))
        xticks = (xt, xt .- xt[div(end+1,2)])
    end

    plot(collect(1:size(args.cres)[1]+1) .- ((ndims(args.cres) == 2) ? .1 : 0),
         ndims(args.cres) == 2 ? args.cres[:,1] : args.cres;
         seriestype = :sticks, linewidth = 100/(args.lag+1),
         legend = :none, kw...)

    if ndims(args.cres) == 2

        plot!(collect(1:size(args.cres)[1]+1) .+ .1,
              args.cres[:,2];
              seriestype = :sticks, linewidth = 100/(args.lag+1))

    end

    if (args.type == "cor") | (split(args.type,'_')[1] == "pacf")
        lg = args.lag
        a1 = args.alpha[1]; c1 = args.ci[1]
        a2 = args.alpha[2]; c2 = args.ci[2]
        dc = c2-c1

        plot!([-c1,c1]; seriestype = :hline, seriescolor = :black, linealpha = 0.5, linestyle = :dash)

        ax = args.auto ? lg : 2*lg + 1
        annotations =
            [(ax,c1+dc/8,string("CI %",a1*100),font(3, "Courier")),
             (ax,c2+dc/8,string("CI %",a2*100),font(3, "Courier"))]

        plot!([-c2,c2]; seriestype = :hline, seriescolor = :black,
              linealpha = 0.5, linestyle = :dot, annotations)

    end

end


text:
include("testio.jl")
include("testevaluate.jl")
#include("testsgm.jl")

include("testseg.jl")
include("testaff.jl")


text:
using StorageMirrorServer: timeout_call, TimeoutException, with_cache_dir

@testset "timeout_call" begin
    timeout_err = TimeoutException(1)
    msg = @capture_err try
        throw(timeout_err)
    catch err
        @warn err
    end
    @test occursin(string(timeout_err), msg)

    try
        response = timeout_call(2) do
            sleep(1)
            return true
        end
        @test response == true
    catch err
        @test false # should never reach here
    end

    try
        response = timeout_call(1) do
            sleep(2)
            return true
        end
        @test false # should never reach here
    catch err
        @test err isa TimeoutException && err.timeout == 1
    end

    try
        response = timeout_call(2) do
            sleep(1)
            error("reached here")
        end
        @test false # should never reach here
    catch err
        @test err isa ErrorException && err.msg == "reached here"
    end

    try
        response = timeout_call(1) do
            sleep(2)
            @test false # should never reach here
        end
        @test false # should never reach here
    catch err
        @test err isa TimeoutException && err.timeout == 1
    end

    start_time = now()
    Threads.@threads for i in 1:8
        try
            response = timeout_call(5) do
                sleep(1)
                return true
            end
            response == true || error("should never reach here")
        catch err
            error("should never reach here")
        end
    end
    # each iteration should end in about 1 - 1.5s
    @test 8_000 < (now() - start_time).value * Threads.nthreads() < 1.5 * 8_000

    start_time = now()
    Threads.@threads for i in 1:8
        try
            response = timeout_call(1) do
                sleep(5)
                error("should never reach here")
            end
            error("should never reach here")
        catch err
            if !(err isa TimeoutException && err.timeout == 1)
                rethrow(err)
            end
        end
    end
    # each iteration should end in about 1 - 1.5s
    @test 8_000 < (now() - start_time).value * Threads.nthreads() < 1.5 * 8_000
end

@testset "with_cache_dir" begin
    @testset "case 1" begin
        # case 1: the first function call succeeds

        cache_dir = tempname()
        logfile = joinpath(cache_dir, "log.txt")

        @test !isdir(cache_dir) && !isfile(logfile)
        with_cache_dir(cache_dir) do
            mkpath(cache_dir)
            open(logfile, "w") do io
                print(io, "This should be called.")
            end
        end
        @test isdir(cache_dir) && isfile(logfile)
        isfile(logfile) && @test read(logfile, String) == "This should be called."

        with_cache_dir(cache_dir) do
            open(logfile, "w") do io
                print(io, "This doesn't get called")
            end
            @test false # mark a fail case if this get called
        end
        @test isfile(logfile)
        isfile(logfile) && @test read(logfile, String) == "This should be called."
    end

    @testset "case 2" begin
        # case 2: the first function get called but returns a false

        cache_dir = tempname()
        logfile = joinpath(cache_dir, "log.txt")

        @test !isdir(cache_dir) && !isfile(logfile)
        with_cache_dir(cache_dir) do
            mkpath(cache_dir)
            open(logfile, "w") do io
                print(io, "This get called.")
            end
            return false
        end
        @test !isdir(cache_dir) && !isfile(logfile)

        with_cache_dir(cache_dir) do
            mkpath(cache_dir)
            open(logfile, "w") do io
                print(io, "This get called again.")
            end
        end
        @test isfile(logfile)
        isfile(logfile) && @test read(logfile, String) == "This get called again."
    end

    @testset "case 3" begin
        # case 3: the first get called but an exception is thrown

        cache_dir = tempname()
        logfile = joinpath(cache_dir, "log.txt")

        @test !isdir(cache_dir) && !isfile(logfile)
        try
            with_cache_dir(cache_dir) do
                mkpath(cache_dir)
                open(logfile, "w") do io
                    print(io, "This get called.")
                end
                error("some error")
            end
        catch err
            @test err isa ErrorException && err.msg == "some error"
        end
        @test !isdir(cache_dir) && !isfile(logfile)

        with_cache_dir(cache_dir) do
            mkpath(cache_dir)
            open(logfile, "w") do io
                print(io, "This get called again.")
            end
        end
        @test isfile(logfile)
        isfile(logfile) && @test read(logfile, String) == "This get called again."
    end

    @testset "case 4" begin
        # case 4: cache isn't valid anymore
        cache_dir = tempname()
        logfile = joinpath(cache_dir, "log.txt")

        is_valid_cache(cache_dir) = "log.txt" in readdir(cache_dir)

        @test !isdir(cache_dir)
        with_cache_dir(cache_dir; by=is_valid_cache) do
            mkpath(cache_dir)
            @test !is_valid_cache(cache_dir)
        end
        @test !isdir(cache_dir)

        with_cache_dir(cache_dir; by=is_valid_cache) do
            mkpath(cache_dir)
            open(logfile, "w") do io
                print(io, "This get called.")
            end
        end
        @test isfile(logfile)
        isfile(logfile) && @test read(logfile, String) == "This get called."
    end

end


text:
#We  are using this to make the Cd v. Velocity Plots
#import necessary libraries

using Plots

function air_Density_Calculation(h)
    temp = 59 - .00356 * h;         #h = altitude(ft), temp = Fereinheidt.

    rho = 2116 * [(temp + 459.7)/ 518.6]^5.256;
    #Density as a function of tempurature. Only of temp below 36152 ft.

    return rho;

end
function root_Mean_Chord(c_r, c_t)
    c_mac = (2/3) * (c_r + c_t + ((c_r*c_t)/(c_r+c_t)));
    return c_mac
end

function form_Factor_Calc(t_c, c_mac, Delta)
    #requires tip chord, root mean Chord, and sweep angle (Delta)
    #returns value of k

    Z = 2*cos(Delta)
    k = 1 + Z*(t_c/c_mac) + 100*(t_c/c_mac)^4

end
#Function Needs: Speed, Altitude, Surface Area, Tip Chord, Root Mean Chord, Sweep Angle, and mu.
function parasitic_Drag_Calculation(V_inf, q_inf, S_wet, tip2chord, c_mac, Delta, mu)
    #S_Wet = 2*(1+ 0.2*(tip2chord)) * S_wing;  #S_Wing passed in vector.
    #q_inf = (1/2) * p_inf * (V_inf^2);
    k = form_Factor_Calc(t, c, Delta);
    #Pass root chord and tip chord lengths, Only needed for Cessna, so don't do it here.
    RE = (p_inf * V_inf * c)/(mu);
    C_f = 1.328/sqrt(RE);           #Laminar flow, we won't worry about turbulent flow.

    #Use atmospheric caluculator online to find specific values for mu, T, etc.

    #mu is passed on vector, c is the mean aerodynamic chord of the aircraft.
    #mu is the dynamic viscocity of the fluid.

    #fineness ratio, fr, is the ratio of an airplanes length to its width.c

    #D_p = k * C_f * q_inf * S_Wet
    #q_inf = (1/2) * p_inf * (V_inf)^2
    #RE = (p*V*c) / (mu)
    #   p = fluid density
    #   mu = fluid viscocity

    D_p = k * C_f * q_inf * S_Wet;  #So far only need k and C_f

    return D_p;
end
function induced_Drag_Calculation(L, v_inf, h, )
    p_inf

    D_i = L^2 / (q_inf * pi * b^2 * e_inv)
end
function main()
    #We need 2 values, induced drag and parasitic drag.
    #If plane cruises above Mach 0.3 (>230 mph), include Compressiblity Drag

    #Choose 2 planes whose values you mimic. One should be glider esq. (high aspect ratio) and the other more like a fighter jet (low aspect ratio)

    #Cessna 172R
        # Aspect Ratio: 7.32 (according to wikipedia)
        # Length: 27 ft, 2 in
        # Wing Span: 36 ft, 1 in.
        # Surface Area: 174 ft^2, S_wing
        # Aspect Ratio: 4.24153... (according to my calculations)
        # Cruising Speed: 140 mph (top 188 mph)
        # Stalling speed: 54 mph
        # Loaded Wight: 2450 lb
        # Cruising Alt: 6000 ft., mu = 3.6162e-7   lb-s / ft^2
        # Overall Length: 26 ft, 11 in.

    #Vought V-173, Flying Pancake (oldmachinepress.com)
        # Span: 23 ft, 4 in.
        # Surface Area: 427 ft^2, S_wing
        # Aspect Ratio = 1.275....
        # Max Speed: 138 mph
        # Cruising Speed: 75 mph
        # Stalling Speed:
        # Cruising Altitude: 28000 ft., mu = 3.1510e-7  lb-s / ft^2
        # Root mean chord: SurfaceArea / Span

        # Loaded Weight = 3050 lb

        #Calculate the root mean chord, c, before calling drag functions.

        rho = air_Density_Calculation(h);
        S_Wet = 2*(1+ 0.2*(tip2chord)) * S_wing;  #S_Wing passed in vector.
        q_inf = (1/2) * p_inf * (V_inf^2);

        mu = 3.1510e-7; #28000 ft.
        # mu = 3.6162e-7; #6000 ft.

    #Parasitic Drag:
        D_p = parasitic_Drag_Calculation(velocity, rho, S_Wet, tip2chord, c_mac, delta, mu);
        #D_p = k * C_f * q_inf * S_wet
        #k = Pressure drag constant
        #C_f = Skin Friction drag constant
        #q_inf = Wind speed over wing
        #S_wet = 2*(1+ 0.2*(t/c)) * S_wing
        # c_mean is the root mean chord.
        # t_c is the thickness to chord ratio, as a decimal-percent

        #Function Needs: Speed, Altitude, Surface Area, Tip Chord, Root Mean Chord, Sweep Angle, and mu.

    #Induced Drag
        df = 1;                         #Diameter of the fuselage.
        e_inv = 0.98 *(1-(2*(d_f/b)^2))
        D_i = induced_Drag_Calcuation();
        #D_i = L^2 / (q_inf * pi * b^2 * e_inv)
        #C_di = C_l^2 / (pi * AR * e_inv)

        #e_inv(Invicide span efficiency) = 0.98 * [1-2*(d_f/b)^2]
        #d_f = diameter of the fuselage
        #L is lift, which is just equal to the weightof the aircraft.

        #D = D_p + D_i

        #Velocity range (0, cruising speed), Do this to start
        #Find the point at which the two drags equalize, this is the minimum
        # value we want
end


text:
const microclimoz_files = [
    "RH120cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92496.1",
    "RH1cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92758.1",
    "RH1cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92759.1",
    "RH1cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92760.1",
    "RH1cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92757.1",
    "SNOWDEP_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92675.1",
    "SNOWDEP_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92676.1",
    "SNOWDEP_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92677.1",
    "SNOWDEP_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92678.1",
    "SOLR.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92487.1",
    "TSKY_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92512.1",
    "TSKY_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92514.1",
    "TSKY_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92515.1",
    "TSKY_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92513.1",
    "Tair120cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92526.1",
    "Tair1cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92502.1",
    "Tair1cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92524.1",
    "Tair1cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92525.1",
    "Tair1cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92511.1",
    "V120cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92495.1",
    "V1cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92497.1",
    "ZEN.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92486.1",
    "ausborder.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92754.1",
    "base_data.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92748.1",
    "hum0cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92681.1",
    "hum0cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92682.1",
    "hum0cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92683.1",
    "hum0cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92684.1",
    "hum100cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92715.1",
    "hum100cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92716.1",
    "hum100cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92717.1",
    "hum100cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92718.1",
    "hum10cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92695.1",
    "hum10cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92696.1",
    "hum10cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92697.1",
    "hum10cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92698.1",
    "hum15cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92699.1",
    "hum15cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92700.1",
    "hum15cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92701.1",
    "hum15cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92703.1",
    "hum2.5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92687.1",
    "hum2.5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92688.1",
    "hum2.5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92689.1",
    "hum2.5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92690.1",
    "hum200cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92680.1",
    "hum20cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92705.1",
    "hum20cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92704.1",
    "hum20cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92706.1",
    "hum20cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92702.1",
    "hum30cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92707.1",
    "hum30cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92708.1",
    "hum30cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92709.1",
    "hum30cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92710.1",
    "hum50cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92711.1",
    "hum50cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92712.1",
    "hum50cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92713.1",
    "hum50cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92714.1",
    "hum5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92691.1",
    "hum5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92692.1",
    "hum5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92693.1",
    "hum5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92694.1",
    "moist0cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92625.1",
    "moist0cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92626.1",
    "moist0cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92636.1",
    "moist0cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92637.1",
    "moist100cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92670.1",
    "moist100cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92671.1",
    "moist100cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92672.1",
    "moist100cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92673.1",
    "moist10cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92650.1",
    "moist10cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92651.1",
    "moist10cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92652.1",
    "moist10cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92653.1",
    "moist15cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92654.1",
    "moist15cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92655.1",
    "moist15cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92656.1",
    "moist15cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92657.1",
    "moist2.5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92638.1",
    "moist2.5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92639.1",
    "moist2.5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92640.1",
    "moist2.5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92641.1",
    "moist200cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92619.1",
    "moist20cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92658.1",
    "moist20cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92659.1",
    "moist20cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92660.1",
    "moist20cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92661.1",
    "moist30cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92662.1",
    "moist30cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92663.1",
    "moist30cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92664.1",
    "moist30cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92665.1",
    "moist50cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92666.1",
    "moist50cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92667.1",
    "moist50cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92668.1",
    "moist50cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92669.1",
    "moist5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92644.1",
    "moist5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92645.1",
    "moist5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92648.1",
    "moist5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92649.1",
    "pot0cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92574.1",
    "pot0cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92575.1",
    "pot0cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92576.1",
    "pot0cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92577.1",
    "pot100cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92610.1",
    "pot100cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92628.1",
    "pot100cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92629.1",
    "pot100cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92634.1",
    "pot10cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92570.1",
    "pot10cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92571.1",
    "pot10cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92572.1",
    "pot10cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92573.1",
    "pot15cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92584.1",
    "pot15cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92585.1",
    "pot15cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92588.1",
    "pot15cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92589.1",
    "pot2.5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92578.1",
    "pot2.5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92579.1",
    "pot2.5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92580.1",
    "pot2.5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92581.1",
    "pot200cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92679.1",
    "pot20cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92592.1",
    "pot20cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92593.1",
    "pot20cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92596.1",
    "pot20cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92597.1",
    "pot30cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92600.1",
    "pot30cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92601.1",
    "pot30cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92602.1",
    "pot30cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92603.1",
    "pot50cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92604.1",
    "pot50cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92605.1",
    "pot50cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92606.1",
    "pot50cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92608.1",
    "pot5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92509.1",
    "pot5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92527.1",
    "pot5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92528.1",
    "pot5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92529.1",
    "soil0cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92503.1",
    "soil0cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92518.1",
    "soil0cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92516.1",
    "soil0cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92523.1",
    "soil100cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92562.1",
    "soil100cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92564.1",
    "soil100cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92565.1",
    "soil100cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92566.1",
    "soil10cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92538.1",
    "soil10cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92539.1",
    "soil10cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92540.1",
    "soil10cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92541.1",
    "soil2.5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92530.1",
    "soil2.5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92531.1",
    "soil2.5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92532.1",
    "soil2.5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92533.1",
    "soil200cm.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92568.1",
    "soil20cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92546.1",
    "soil20cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92547.1",
    "soil20cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92548.1",
    "soil20cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92549.1",
    "soil30cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92550.1",
    "soil30cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92551.1",
    "soil30cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92552.1",
    "soil30cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92553.1",
    "soil50cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92554.1",
    "soil50cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92555.1",
    "soil50cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92556.1",
    "soil50cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92557.1",
    "soil5cm_0pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92534.1",
    "soil5cm_50pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92535.1",
    "soil5cm_70pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92536.1",
    "soil5cm_90pctShade.zip" => "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/knb.92537.1",
]

"""
    download_microclim(dir::String; overwrite=false)

Download the microclim-oz dataset to `dir`.

WARNING: This may take many hours, and requires 380 GB of free space when unzipped.

Passing `overwrite=true` will force downloading and unzipping of all files.
Otherwise only missing files are downloaded and unzipped.
"""
function download_microclim(dir::String; overwrite=false)
    for (filename, url) in microclimoz_files

        # Download zip file
        filepath = joinpath(dir, filename)
        if !isfile(filepath) || overwrite
            println("Downloading: ", filepath, "\n    from: ", url)
            download(url, filepath)
        else
            println("Skipping ", filename, ", file exists")
        end

        # Unzip
        needsunzip = true
        unzipdir = if '_' in filename
            joinpath(dir, match(r"([^_]*).*\.zip", filename).captures[1])
        else
            dir
        end
        out = Pipe()
        zipinfo = run(pipeline(`zipinfo -1 $filepath`; stdout=out), wait=false)
        close(out.in)
        if success(zipinfo)
            zipfiles = map(s -> joinpath(unzipdir, s), split(String(read(out)), "\n"; keepempty=false))
            needsunzip = !all(isfile.(zipfiles))
        else
            error("`zipinfo` for $filepath failed: the file may be corrupted. Delete and run this method again")
        end
        needsunzip && run(`unzip -o $filepath -d $unzipdir`)
    end
end


text:
# TODO: have a doc for ops here.
#"@knet function dot(w,x) is matrix multiplication."
#"@knet function input() fetches the next network input."
# ### mul2 element-wise multiplication:

abstract Op

# Each op must provide the following:
# back_reads_y (tosave)
# back_reads_x (tosave)
# ninputs (netcomp1)
# infersize (used by netinit)
# canoverwrite (used by initforw)
# ysize, loss (not used any more?)

forw(l::Op, y, x...; o...) = error("$(typeof(l)) has not implemented forw")
back(l::Op, dy, dx...; o...)   = error("$(typeof(l)) has not implemented back")
loss(l::Op, dy, y; o...)   = error("$(typeof(l)) has not implemented loss")

ninputs(l::Op)          = error("$(typeof(l)) has not implemented ninputs")
ysize(l::Op, x...)      = error("$(typeof(l)) has not implemented ysize")
canoverwrite(l::Op)       = error("$(typeof(l)) has not implemented canoverwrite")
back_reads_x(l::Op)     = error("$(typeof(l)) has not implemented back_reads_x")
back_reads_y(l::Op)     = error("$(typeof(l)) has not implemented back_reads_y")
# Base.eltype(::Op) = nothing

function Base.isequal(a::Op,b::Op)
    typeof(a)==typeof(b) || return false
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end

### DEAD CODE
# Each Op implements some common methods, stubs are given below.
# forw takes input x and returns output y, possibly setting some state.
# back takes dy, the loss gradient wrt y, calculates loss gradient wrt
# parameters and optionally returns dx, the loss gradient wrt x.
# Some layers overwrite their inputs.


text:
##########
##########
# Definition of functions used to look at species rewiring in relevant paper
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-05-10
##########
##########

# set up =======================================================================

# set packages
using DataFrames, Distributions, Random, StatsBase, Documenter, Colors
using LightGraphs, GraphPlot, GraphRecipes, Plots

include("./src/network_creation.jl")

# create matrices for 1000 iterations of each setting ==========================

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
# p7.m - accuracy of periodic spectral differentiation

# Compute derivatives for various values of N:
Nmax = 50; allN = 6:2:Nmax; E = zeros(4,length(allN));
for N = 6:2:Nmax
    h = 2*pi/N; x = h*(1:N);
    column = [0; @. .5*(-1)^(1:N-1)*cot((1:N-1)*h/2)];
    D = toeplitz(column,column[[1;N:-1:2]]);
    v = @. abs(sin(x))^3;                     # 3rd deriv in BV
    vprime = @. 3*sin(x)*cos(x)*abs(sin(x));
    j = round(Int,N/2-2);
    E[1,j] = norm(D*v-vprime,Inf);
    v = @. exp(-sin(x/2)^(-2));               # C-infinity
    vprime = @. .5*v*sin(x)/sin(x/2)^4;
    E[2,j] = norm(D*v-vprime,Inf);
    v = @. 1/(1+sin(x/2)^2);                 # analytic in a strip
    vprime = @. -sin(x/2)*cos(x/2)*v^2;
    E[3,j] = norm(D*v-vprime,Inf);
    v = sin.(10*x); vprime = 10*cos.(10*x);   # band-limited
    E[4,j] = norm(D*v-vprime,Inf);
end

# Plot results:
titles = [L"|\sin(x)|^3",L"\exp(-\sin^{-2}(x/2))",L"1/(1+\sin^2(x/2))",L"\sin(10x)"]; clf();
for iplot = 1:4
    subplot(2,2,iplot);
    semilogy(allN,E[iplot,:],".-",markersize=6);
    axis([0,Nmax,1e-16,1e3]); grid(true);
    xticks(0:10:Nmax); yticks(10.0.^(-15:5:0));
    title(titles[iplot]);
    iplot > 2 ? xlabel("N") : nothing;
    iplot%2 > 0? ylabel("error") : nothing;
end


text:
using FaultTolerantControl
using Test
using LinearAlgebra
using Transducers

@testset "costs" begin
    Q_p = Matrix(I, 3, 3)
    Q_ω = Matrix(I, 3, 3)
    F_p = Matrix(I, 3, 3)
    F_ω = Matrix(I, 3, 3)
    cf = PositionAngularVelocityCostFunctional(Q_p, Q_ω, F_p, F_ω)
    # data
    t0, tf = 0.0, 1.0
    Δt = tf - t0
    ts = t0:0.01:tf
    e_p_nom = ones(3)
    e_ω_nom = ones(3)
    e_ps = ts |> Map(t -> e_p_nom) |> collect
    e_ωs = ts |> Map(t -> e_ω_nom) |> collect
    J = cost(cf, ts, e_ps, e_ωs)
    @test J ≈ e_p_nom'*Q_p*e_p_nom*Δt + e_ω_nom'*Q_ω*e_ω_nom + e_ps[end]'*F_p*e_ps[end] + e_ωs[end]'*F_ω*e_ωs[end]
end


text:
# #calculate expected value
# tst = Float64[]
#
# for i in 1:length(res)
#     push!(txn, res[i].value)
#     push!(tst, txn.value.quantity * txn.value.basis)
# end
#
# eq = 100_000.00
# for t in tst
#     eq -= t
# end
#
# const Eval = eq # 96214.00

facts("AssetTransaction type constructor works") do

    context("Good things work") do
        @fact 1 => 1
    end

    context("Bad things don't work") do
        @fact_throws AssetTransaction("foo")
    end
end


text:
module Construct

greet() = print("Hello World!")

include("mars_weather.jl")
export weather_past_week

end # module


text:
module LotkaVolterra2dSingularSPARK

    const Δt = 0.1
    const nt = 100000

    const PLOT_DIR = "figures"
    const SYMP_DIR = "symplecticity"

    using GeometricIntegrators
    using GeometricProblems.LotkaVolterra2dSingular
    using GeometricProblems.LotkaVolterra2dPlots

    set_config(:nls_atol, 1E-14)
    set_config(:nls_rtol, 1E-14)
    set_config(:nls_atol_break, 1E6)
    set_config(:nls_rtol_break, 1E6)
    set_config(:nls_stol_break, 1E6)
    set_config(:nls_nmax, 100);

    include("common.jl")
    include("tableau_lists.jl")

    function plot(sol, stages, equ, dir, file, fig_suff)
        if !isdir(dir)
            mkdir(dir)
        end

        plot_lotka_volterra_2d(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * fig_suff)

        plot_lotka_volterra_2d_solution(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_solution" * fig_suff)

        plot_lotka_volterra_2d_traces(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_traces" * fig_suff)

        _plot(sol, stages, equ, dir, file, fig_suff)
    end

end


text:
using PyCallJLD
using Test
using PyCall, JLD

tmp = joinpath(tempdir(), "temp.jld")
const deque = pyimport("collections")["deque"]

try
    obj = [deque([1,2,3]), deque([4,5,6])]
    save(tmp, "x", obj)
    @test load(tmp, "x") == obj

    save(tmp, "x", PyObject("fé"))
    @test load(tmp, "x") == PyObject("fé")

    obj = PyObject(rand(UInt8, 1000))
    save(tmp, "x", obj)
    @test load(tmp, "x") == obj
finally
    rm(tmp, force=true)
end


text:
import Random
using AbstractAlgebra

function test_params_isgb(rng, nvariables, exps, nterms,
                        npolys, grounds, coeffssize, orderings)

    without(x, k) = x[1:end .!= k]

    for n in nvariables
        for e in exps
            for nt in nterms
                for np in npolys
                    for gr in grounds
                        for ord in orderings
                            for csz in coeffssize
                                set = Groebner.generate_set(
                                    n, e, nt, np, csz, rng, gr, ord
                                )
                                isempty(set) && continue
                                gb = Groebner.groebner(set)

                                @test Groebner.isgroebner(gb)
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "isgroebner random stress tests" begin

    rng = Random.MersenneTwister(5)

    # :degrevlex finite case tests
    nvariables = [2, 3, 5]
    exps       = [1:2, 2:4, 4:5]
    nterms     = [1:1, 1:2, 4:5]
    npolys     = [1:1, 1:3, 4:5]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:degrevlex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for isgroebner"
    test_params_isgb(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :lex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:lex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for isgroebner"
    test_params_isgb(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :deglex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:deglex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for isgroebner"
    test_params_isgb(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

end


text:
"""
    solve

Solve equation.
"""
function solve end

include("solve_advection.jl")
include("solve_diffusion.jl")
include("solve_burgers.jl")


text:
using RobustAndOptimalControl, Random, Statistics, MonteCarloMeasurements, ControlSystems
Random.seed!(0)
unsafe_comparisons()
u1 = 2*δ()
u2 = 5*δ()
u3 = 0.2*δ()
nx, nu, ny = 10, 5, 5
A = randn(nx,nx) .+ 0 .* δ.()
A[1,1] = 4 + u1
A[2, 3] = 9 + u2
B = randn(nx, nu) .+ 0 .* δ.()
B[4, 1] = 2 + u3

sys = ss(A, B, randn(ny,nx), 0)
l, res = find_lft(sys, [u1, u2, u3])

sysh = ss(l,sys)

@test pmean.(sysh.A) ≈ pmean.(sys.A) atol=1e-5
@test pstd.(sysh.A) ≈ pstd.(sys.A) rtol=1e-5
@test pmean.(sysh.B) ≈ pmean.(sys.B) atol=1e-5
@test pstd.(sysh.B) ≈ pstd.(sys.B) rtol=1e-5
@test pmean.(sysh.C) ≈ pmean.(sys.C) atol=1e-5
@test pstd.(sysh.C) ≈ pstd.(sys.C) rtol=1e-5
@test pmean.(sysh.D) ≈ pmean.(sys.D) atol=1e-5
@test pstd.(sysh.D) ≈ pstd.(sys.D) rtol=1e-5

## TODO: add example from p. 173 Zhou

δm = δr(32)
δC = δr(32)
δk = δr(32)

m̄ = 2
c̄ = 1
k̄ = 10

m = m̄*(1 + 0.1*δm)
c = c̄*(1 + 0.2*δC)
k = k̄*(1 + 0.3*δk)

M = [
    0 1 0 0 0 0
    -k̄/m̄ -c̄/m̄ 1/m̄ -1/m̄ -1/m̄ -0.1/m̄
    0.3k̄ 0 0 0 0 0
    0 0.2c̄ 0 0 0 0
    -k̄ -c̄ 1 -1 -1 -0.1
]

A = [
    0 1
    -k/m -c/m
]
B = [
    0
    1/m
]

sys = ss(A, B, I(2), 0)
using Optim, ComponentArrays
# Optim.similar_axis(x::ComponentArray, n) = x .* zeros(length(x), n)
unsafe_comparisons()
l, res = find_lft(sys, [δm, δC, δk])
# l, res = find_lft(sys, 3)
Mh = Matrix(l)
norm(M - Mh) / norm(M)
sysh = ss(l,sys)

@test pmean.(sysh.A) ≈ pmean.(sys.A) rtol=1e-2
@test pstd.(sysh.A) ≈ pstd.(sys.A) rtol=1e-2
@test pmean.(sysh.B) ≈ pmean.(sys.B) rtol=1e-2
@test pstd.(sysh.B) ≈ pstd.(sys.B) rtol=1e-2


text:
using Base.Dates  # just using Dates in versions > 0.6
println("Time zero (the epoch) is $(unix2datetime(0)).")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("../../src/SPICE.jl")
using .SPICE

# ---- State-change matrix
const ν = [[1 0 0 0 0 0 0];[0 0 0 0 0 0 0];[0 1 1 0 0 0 0];[1 0 0 0 0 0 0];[0 0 0 0 1 1 0];[0 0 0 0 0 0 1];[0 0 0 1 0 0 0];[0 0 1 0 0 0 0]]-
[[0 0 0 0 0 0 0];[1 0 0 0 0 0 0];[1 1 0 0 0 0 0];[0 0 1 0 0 0 0];[0 0 1 1 0 0 0];[0 0 0 0 1 0 0];[0 0 0 0 0 1 1];[0 0 0 0 0 0 0]]

# ---- Initial state-space
x0 = [500, 4, 110, 300, 2, 20, 90]

# ---- Highest orders of reaction for tau-leaping
const hor = [2,2,2,2,1,2,2]

# ---- Bounds on the initial search space for parameters
# ---- Real parameters are [.38, .04, .082, .12, .021, .1,.005, 13.21]
const bounds = [[1e-6,1e-6,1e-6,1e-6,1e-6,1e-6,1e-6,1e-6], [10.0,10.0,10.0,10.0,10.0,10.0,10.0,100.0]]

# ---- Hazard functions
function F(p)
	H(p)
	p.a .*= p.θ
end

function H(p)
	p.a[1] = 1.0
    p.a[2] = p.x[1]
    p.a[3] = p.x[1] * p.x[2]
    p.a[4] = p.x[3]
    p.a[5] = p.x[3] * p.x[4]
    p.a[6] = p.x[5]
    p.a[7] = p.x[6] * p.x[7]
    p.a[8] = 1.0
end
# ---- Configure system in SPICE
system = System(Model(x0,F,H,ν,bounds,hor=hor, obsname=[:R,:L,:RL,:G,:Ga,:Gbg,:Gd]), "./data/5-sets",routine=CEM(ssa=:Direct, nElite = 10, nRepeat = 1, nSamples=1000, maxIter=250, mSamples=20000, shoot=false, splitting=false, sampling=:log))

estimate(system, 100, "yeast")

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
Base.@kwdef mutable struct SimParam
    # Reynolds number
    Re::Real = 250.
    # In-flow [m2/s]
    Q::Real = 1.
    # CFL
    cfl::Real = 0.5
    # Relaxation coefficient
    alpha::Real = 1.97
    # Ramp time for bump
    dtdτ::Real = 1.
    # Tolerance for Gauss Seidel
    tol::Real = 1e-3
    # Ramp for khi
    ramp::Int = 200
end
getuc(pa::SimParam, gb::GridBox) = 0.75*pa.Q/gb.D
getum(pa::SimParam, gb::GridBox) = pa.Q / gb.D
getnu(pa::SimParam, gb::GridBox) = getuc(pa, gb) * gb.D / pa.Re

mutable struct Simulation
    u::UGrid
    v::VGrid
    p::PGrid
    params::SimParam
    dt::Real
    canal::CanalConfig
    logger::AbstractLogger
    step::Int
end

function Simulation(params::SimParam, gridbox::GridBox, canal::CanalConfig, logger = SimpleLogger())
    dt = params.cfl * gridbox.dx / getuc(params, gridbox)
    Simulation(
        UGrid(gridbox),
        VGrid(gridbox),
        PGrid(gridbox),
        params,
        dt,
        canal,
        logger,
        0
        )
end

function Base.copy(sim::Simulation)
    Simulation(
        copy(sim.u),
        copy(sim.v),
        copy(sim.p),
        sim.params,
        sim.dt,
        sim.canal,
        sim.logger,
        sim.step
    )
end


text:
#=
A backend for POVRay

WARNING: Install POVRay and then install http://www.ignorancia.org/index.php/technical/lightsys/
so that realistic skies can be used

=#
export POVRay,
       povray,
       povray_material,
       default_povray_material

# We will use MIME types to encode for POVRay
const MIMEPOVRay = MIME"text/povray"

# First, basic types
show(io::IO, ::MIMEPOVRay, s::String) =
  show(io, s)
show(io::IO, ::MIMEPOVRay, s::Symbol) =
  print(io, s)
show(io::IO, ::MIMEPOVRay, r::Real) =
  show(io, r)
show(io::IO, mime::MIMEPOVRay, v::Vector) =
  begin
    write(io, "<")
    for (i, e) in enumerate(v)
      if i > 1
        write(io, ", ")
      end
      show(io, mime, e)
    end
    write(io, ">")
    nothing
  end
show(io::IO, ::MIMEPOVRay, p::Union{Loc,Vec}) =
  # swap y with z to make it consistent with POVRay coordinate system
  let p = in_world(p)
    print(io, "<$(p.x), $(p.z), $(p.y)>")
  end
show(io::IO, ::MIMEPOVRay, c::RGB) =
  print(io, "color rgb <$(Float64(red(c))), $(Float64(green(c))), $(Float64(blue(c)))>")
show(io::IO, ::MIMEPOVRay, c::RGBA) =
  print(io, "color rgbt <$(Float64(red(c))), $(Float64(green(c))), $(Float64(blue(c))), $(Float64(1-alpha(c)))>")

write_povray_object(f::Function, io::IO, type, material, args...) =
  let mime = MIMEPOVRay()
    write(io, "$(type) {")
    for (i, arg) in enumerate(args)
      if i == 1
        write(io, "\n  ")
      elseif i > 1
        write(io, ", ")
      end
      show(io, mime, arg)
    end
    write(io, '\n')
    if ! isnothing(material)
      write_povray_material(io, material)
    end
    let res = f()
      write(io, "}\n")
      res
    end
  end

write_povray_object(io::IO, type, material, args...) =
  write_povray_object(io, type, material, args...) do
  end

write_povray_value(io::IO, value::Any) =
  show(io, MIMEPOVRay(), value)

write_povray_param(io::IO, name::String, value::Any) =
  begin
    write(io, "  ", name, " ")
    write_povray_value(io, value)
    write(io, '\n')
  end

write_povray_material(io::IO, material) =
  begin
    write(io, "  ")
    write_povray_value(io, material)
    write(io, '\n')
  end

write_povray_matrix(io::IO, p::Loc) =
  let t = (translated_cs(p.cs, p.x, p.y, p.z).transform)
    write_povray_param(io, "matrix",
      [t[1,1], t[3,1], t[2,1],
       t[1,2], t[3,2], t[2,2],
       t[1,3], t[3,3], t[2,3],
       t[1,4], t[3,4], t[2,4]])
  end

write_povray_camera(io::IO, camera, target, lens) =
  write_povray_object(io, "camera", nothing) do
    write_povray_param(io, "location", camera)
    write_povray_param(io, "look_at", target)
    write_povray_param(io, "right", Symbol("x*image_width/image_height"))
    write_povray_param(io, "angle", view_angles(lens)[1])
  end

write_povray_pointlight(io::IO, location, color) =
  write_povray_object(io, "light_source", nothing, location, color)

#=
struct POVRayPigment
  color::RGB
end

show(io::IO, mime::MIMEPOVRay, p::POVRayPigment) =
  write_povray_object(io, mime, "pigment") do
    show(io, mime, p.color)
  end

#=
NORMAL:
  normal { [NORMAL_IDENTIFIER] [NORMAL_TYPE] [NORMAL_MODIFIER...] }
NORMAL_TYPE:
  PATTERN_TYPE Amount |
  bump_map { BITMAP_TYPE "bitmap.ext" [BUMP_MAP_MODS...]}
NORMAL_MODIFIER:
  PATTERN_MODIFIER | NORMAL_LIST | normal_map { NORMAL_MAP_BODY } |
  slope_map{ SLOPE_MAP_BODY } | bump_size Amount |
  no_bump_scale Bool | accuracy Float
=#

struct POVRayNormal
  type::String
  amount::Real
end

povray_normal(; bumps) =
  POVRayNormal("bumps", bumps)

show(io::IO, mime::MIMEPOVRay, n::POVRayNormal) =
  write_povray_object(io, mime, "normal") do
    write_povray_param(io, mime, n.type, n.amount)
  end

#=
FINISH:
  finish { [FINISH_IDENTIFIER] [FINISH_ITEMS...] }
FINISH_ITEMS:
  ambient COLOR | diffuse [albedo] Amount [, Amount] | emission COLOR |
  brilliance Amount | phong [albedo] Amount | phong_size Amount | specular [albedo] Amount |
  roughness Amount | metallic [Amount] | reflection COLOR |
  crand Amount | conserve_energy BOOL_ON_OFF |
  reflection { Color_Reflecting_Min [REFLECTION_ITEMS...] } |
  subsurface { translucency COLOR } |
  irid { Irid_Amount [IRID_ITEMS...] }
REFLECTION_ITEMS:
  COLOR_REFLECTION_MAX | fresnel BOOL_ON_OFF |
  falloff FLOAT_FALLOFF | exponent FLOAT_EXPONENT |
  metallic FLOAT_METALLIC
IRID_ITEMS:
  thickness Amount | turbulence Amount
=#

struct POVRayFinish
  ambient::Union{RGB,Nothing}
  diffuse::Union{Real,Nothing}
  emission::Union{RGB,Nothing}
  brilliance::Union{Real,Nothing}
  phong::Union{Real,Nothing}
  phong_size::Union{Real,Nothing}
  specular::Union{Real,Nothing}
  roughness::Union{Real,Nothing}
  metallic::Union{Real,Nothing}
  reflection::Union{RGB,Nothing}
  crand::Union{Real,Nothing}
  conserve_energy::Union{Bool,Nothing}
  subsurface::Any
  irid::Any
end

povray_finish(; ambient::Union{RGB,Nothing}=nothing,
                diffuse::Union{Real,Nothing}=nothing,
                emission::Union{RGB,Nothing}=nothing,
                brilliance::Union{Real,Nothing}=nothing,
                phong::Union{Real,Nothing}=nothing,
                phong_size::Union{Real,Nothing}=nothing,
                specular::Union{Real,Nothing}=nothing,
                roughness::Union{Real,Nothing}=nothing,
                metallic::Union{Real,Nothing}=nothing,
                reflection::Union{RGB,Nothing}=nothing,
                crand::Union{Real,Nothing}=nothing,
                conserve_energy::Union{Bool,Nothing}=nothing,
                subsurface::Any=nothing,
                irid::Any=nothing) =
  POVRayFinish(ambient,
               diffuse,
               emission,
               brilliance,
               phong,
               phong_size,
               specular,
               roughness,
               metallic,
               reflection,
               crand,
               conserve_energy,
               subsurface,
               irid)

show(io::IO, mime::MIMEPOVRay, f::POVRayFinish) =
  write_povray_object(io, mime, "finish") do
    for name in fieldnames(typeof(f))
      let v = getfield(f, name)
        if ! isnothing(v)
          write_povray_param(io, mime, string(name), v)
        end
      end
    end
  end

struct POVRayTexture
  pygment::Union{POVRayPigment,Nothing}
  normal::Union{POVRayNormal,Nothing}
  finish::Union{POVRayFinish,Nothing}
end

struct POVRayLibraryTexture
  name::String
end

show(io::IO, mime::MIMEPOVRay, m::POVRayLibraryTexture) =
  write_povray_object(io, mime, "texture") do
    show(io, mime, m.name)
    #show(io, mime, m.interior)
  end

show(io::IO, mime::MIMEPOVRay, t::POVRayTexture) =
  write_povray_object(io, mime, "texture") do
    show(io, mime, t.pygment)
    show(io, mime, t.normal)
    show(io, mime, t.finish)
  end

struct POVRayInterior
  pygment::Union{POVRayPigment,Nothing}
  normal::Union{POVRayNormal,Nothing}
  finish::Union{POVRayFinish,Nothing}
end

struct POVRayMaterial
  name::String
  texture::POVRayTexture
  interior::Union{POVRayInterior,Nothing}
end
=#

abstract type POVRayMaterial end

struct POVRayDefinition <: POVRayMaterial
  name::String
  kind::String
  description::String
end

convert_to_povray_identifier(name) = replace(name, " " => "_")

write_povray_definition(io::IO, d::POVRayDefinition) =
  write(io, "#declare $(convert_to_povray_identifier(d.name)) =\n  $(d.kind) $(d.description)\n")

show(io::IO, ::MIMEPOVRay, d::POVRayDefinition) =
  write(io, "$(d.kind) { $(convert_to_povray_identifier(d.name)) }")

struct POVRayInclude <: POVRayMaterial
  filename::AbstractString
  kind::String
  name::String
end

write_povray_definition(io::IO, m::POVRayInclude) =
  write(io, "#include \"$(m.filename)\"\n")

show(io::IO, ::MIMEPOVRay, m::POVRayInclude) =
  write(io, "$(m.kind) { $(m.name) }")

export povray_definition, povray_include
const povray_definition = POVRayDefinition
const povray_include = POVRayInclude

const povray_concrete =
  povray_definition("Concrete", "texture", """{
   pigment {
     granite turbulence 1.5 color_map {
       [0  .25 color White color Gray75] [.25  .5 color White color Gray75]
       [.5 .75 color White color Gray75] [.75 1.1 color White color Gray75]}}
  finish {
    ambient 0.2 diffuse 0.3 crand 0.03 reflection 0 }
  normal {
    dents .5 scale .5 }}""")

const povray_neutral =
  povray_definition("Material", "texture",
    "{ pigment { color rgb 1 } finish { reflection 0 ambient 0 }}")

const default_povray_material = Parameter{POVRayMaterial}(povray_neutral)

povray_material(name::String;
                gray::Real=0.3,
                red::Real=gray, green::Real=gray, blue::Real=gray,
                specularity=0, roughness=0,
                transmissivity=nothing, transmitted_specular=nothing) =
  povray_definition(name, "texture", """{
  pigment { rgb <$(Float64(red)),$(Float64(green)),$(Float64(blue))> }
  finish { specular $(specularity) roughness $(roughness) }
}""")

#=
IMAGE_MAP:
  pigment {
    image_map {
      [BITMAP_TYPE] "bitmap[.ext]" [gamma GAMMA] [premultiplied BOOL]
      [IMAGE_MAP_MODS...]
      }
  [PIGMENT_MODFIERS...]
  }
 IMAGE_MAP:
  pigment {
   image_map {
     FUNCTION_IMAGE
     }
  [PIGMENT_MODFIERS...]
  }
 BITMAP_TYPE:
   exr | gif | hdr | iff | jpeg | pgm | png | ppm | sys | tga | tiff
 IMAGE_MAP_MODS:
   map_type Type | once | interpolate Type |
   filter Palette, Amount | filter all Amount |
   transmit Palette, Amount | transmit all Amount
 FUNCTION_IMAGE:
   function I_WIDTH, I_HEIGHT { FUNCTION_IMAGE_BODY }
 FUNCTION_IMAGE_BODY:
   PIGMENT | FN_FLOAT | pattern { PATTERN [PATTERN_MODIFIERS] }
=#
povray_image_map_material(name::String; image_map, map_type=0) =
  povray_definition(name, "texture", """{
  pigment { image_map $(image_map) map_type $(map_type) }
  }""")

####################################################
# Sky models

povray_realistic_sky_string(turbidity, inner) = """
#version 3.7;
#include "colors.inc"
#include "CIE.inc"
#include "lightsys.inc"
#include "lightsys_constants.inc"
#include "sunpos.inc"
#declare Current_Turbidity = $(turbidity);
//To solve a bug in CIE_Skylight.in
#local fiLuminous=finish{ambient 0 emission 1 diffuse 0 specular 0 phong 0 reflection 0 crand 0 irid{0}}
#include "CIE_Skylight.inc"
global_settings {
  assumed_gamma 1.0
  radiosity {
  }
}
#default {finish {ambient 0 diffuse 1}}
light_source{
  $(inner)
  Light_Color(SunColor,5)
  translate SolarPosition
}
"""
# This one seems to work better. Still 3.5 though.

#=
povray_realistic_sky_string(turbidity, inner) = """
#version 3.5;
#include "colors.inc"
#include "CIE.inc"
#include "lightsys.inc"
#include "lightsys_constants.inc"
#include "sunpos.inc"
global_settings {
  assumed_gamma 1.0
  radiosity {
  }
}
#default {finish {ambient 0 diffuse 1}}
CIE_ColorSystemWhitepoint(Beta_ColSys, Daylight2Whitepoint(Kt_Daylight_Film))
CIE_GamutMapping(off)
#declare Lightsys_Brightness = 1.0;
#declare Lightsys_Filter = <1,1,1>;
#declare Al=38;    // sun altitude
#declare Az=100;   // sun rotation
#declare North=-z;
#declare DomeSize=1e5;
#declare Current_Turbidity = 5.0;
#declare Intensity_Mult = 0.7;
#include "CIE_Skylight.inc"
light_source{ 0
  Light_Color(SunColor,5)
  translate SolarPosition
}
"""
=#

povray_realistic_sky_string(altitude, azimuth, turbidity, withsun) =
  povray_realistic_sky_string(
    turbidity,
    "vrotate(<0,0,1000000000>,<-$(altitude),$(azimuth+180),0>)")

povray_realistic_sky_string(date, latitude, longitude, meridian, turbidity, withsun) =
  povray_realistic_sky_string(
    turbidity,
    """
#local xpto = SunPos($(year(date)), $(month(date)), $(day(date)), $(hour(date)), $(minute(date)), $(meridian), $(latitude), $(longitude));
vrotate(<0,0,1000000000>,<-Al,Az,0>)
""")

############################################
# Ground models

povray_ground_string(level, c) =
  "plane { y, $(cz(level)) pigment { rgb <$(Float64(red(c))), $(Float64(green(c))), $(Float64(blue(c)))> } }\n"
####################################################
# Clay models
povray_clay_settings_string() =
"""
#version 3.7;
global_settings {
  assumed_gamma 1.0
  radiosity {
      pretrace_start 64/image_width  //0.04
      pretrace_end 1/image_width     //0.002
      count 1000
      nearest_count 10
      error_bound 0.1
      recursion_limit 1
      low_error_factor 0.2 //0.7
      gray_threshold 0
      minimum_reuse 0.001  //0.01
      brightness 1.2
      adc_bailout 0.01/2
  }
}
sky_sphere {
    pigment {
      color rgb 1
    }
  }
"""

####################################################
abstract type POVRayKey end
const POVRayId = Int
const POVRayRef = GenericRef{POVRayKey, POVRayId}
const POVRayNativeRef = NativeRef{POVRayKey, POVRayId}
const POVRayUnionRef = UnionRef{POVRayKey, POVRayId}
const POVRaySubtractionRef = SubtractionRef{POVRayKey, POVRayId}

mutable struct POVRayBackend{K,T} <: LazyBackend{K,T}
  shapes::Shapes
  shape_material::Dict{Shape,POVRayMaterial}
  materials::Dict{POVRayMaterial,POVRayMaterial}
  sky::String
  ground::String
  buffer::LazyParameter{IOBuffer}
  camera::Loc
  target::Loc
  lens::Real
  sun_altitude::Real
  sun_azimuth::Real
end

const POVRay = POVRayBackend{POVRayKey, POVRayId}
# Traits
has_boolean_ops(::Type{POVRay}) = HasBooleanOps{false}()

# In POVRay, everytime we save a shape, we attach the default_povray_material
save_shape!(b::POVRay, s::Shape) =
  begin
    prepend!(b.shapes, [s])
    b.shape_material[s] = default_povray_material()
    s
  end

#=
The POVRay backend cannot realize shapes immediately, only when requested.
=#

void_ref(b::POVRay) = POVRayNativeRef(-1)

const povray =
  POVRay(Shape[],
         Dict{Shape,POVRayMaterial}(),
         Dict{POVRayMaterial,POVRayMaterial}(),
         povray_realistic_sky_string(DateTime(2020, 9, 21, 10, 0, 0), 39, 9, 0, 5, true),
         povray_ground_string(z(0), rgb(0.8,0.8,0.8)),
         LazyParameter(IOBuffer, IOBuffer),
         xyz(10,10,10),
         xyz(0,0,0),
         35,
         90,
         0)

buffer(b::POVRay) = b.buffer()
get_material(b::POVRay, key) = get!(b.materials, key, key)
get_material(b::POVRay, s::Shape) = get_material(b, get(b.shape_material, s, default_povray_material()))

#

set_view(camera::Loc, target::Loc, lens::Real, aperture::Real, b::POVRay) =
  begin
    b.camera = camera
    b.target = target
    b.lens = lens
  end

get_view(b::POVRay) =
  b.camera, b.target, b.lens

###################################

set_sun(altitude, azimuth, b::POVRay) =
  begin
    b.altitude = altitude
    b.azimuth = azimuth
  end

backend_realistic_sky(b::POVRay, date, latitude, longitude, meridian, turbidity, withsun) =
  b.sky = povray_realistic_sky_string(date, latitude, longitude, meridian, turbidity, withsun)

backend_realistic_sky(b::POVRay, altitude, azimuth, turbidity, withsun) =
  b.sky = povray_realistic_sky_string(altitude, azimuth, turbidity, withsun)

backend_ground(b::POVRay, level::Loc, color::RGB) =
  b.ground = povray_ground_string(level, color)

#
delete_all_shapes(b::POVRay) =
  begin
    (empty!(b.shapes); empty!(b.materials); empty!(b.shape_material); nothing)
  end

backend_delete_shapes(b::POVRay, shapes::Shapes) =
  begin
    b.shapes = filter(s->isnothing(findfirst(s1->s1===s, shapes)), b.shapes)
    for s in shapes delete!(b.shape_material, s) end
  end

realize(b::POVRay, s::Sphere) =
  let mat = get_material(b, s)
    write_povray_object(buffer(b), "sphere", mat, s.center, s.radius)
    void_ref(b)
  end

realize(b::POVRay, s::Torus) =
  let buf = buffer(b)
    write_povray_object(buf, "torus", get_material(b, s), s.re, s.ri) do
      write_povray_matrix(buf, s.center)
    end
    void_ref(b)
  end

#=
realize(b::ACAD, s::Cuboid) =
  ACADIrregularPyramidFrustum(connection(b), [s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3])
realize(b::ACAD, s::RegularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b),
                                regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                                regular_polygon_vertices(s.edges, add_z(s.cb, s.h), s.rt, s.angle, s.inscribed))
realize(b::ACAD, s::RegularPyramid) =
  ACADIrregularPyramid(connection(b),
                          regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                          add_z(s.cb, s.h))
realize(b::ACAD, s::IrregularPyramid) =
  ACADIrregularPyramid(connection(b), s.bs, s.t)
realize(b::ACAD, s::RegularPrism) =
  let ps = regular_polygon_vertices(s.edges, s.cb, s.r, s.angle, s.inscribed)
    ACADIrregularPyramidFrustum(connection(b),
                                   ps,
                                   map(p -> add_z(p, s.h), ps))
  end
realize(b::ACAD, s::IrregularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b), s.bs, s.ts)

realize(b::ACAD, s::IrregularPrism) =
  write_povray_object(buffer(b), "prism", get_material(b, s), 0, norm(s.v)) do
    s.bs,
                              map(p -> (p + s.v), s.bs))
## FIXME: deal with the rotation angle
realize(b::ACAD, s::RightCuboid) =
  ACADCenteredBox(connection(b), s.cb, s.width, s.height, s.h)
=#

realize(b::POVRay, s::Box) =
  let buf = buffer(b),
      mat = get_material(b, s)
    write_povray_object(buf, "box", mat, [0,0,0], [s.dx, s.dy, s.dz]) do
      write_povray_matrix(buf, s.c)
    end
    void_ref(b)
  end

realize(b::POVRay, s::Cone) =
  let buf = buffer(b),
      bot = s.cb,
      top = s.cb + vz(s.h, s.cb.cs),
      mat = get_material(b, s)
    write_povray_object(buf, "cone", mat, bot, s.r, top, 0)
    void_ref(b)
  end

realize(b::POVRay, s::ConeFrustum) =
  let buf = buffer(b),
      bot = s.cb,
      top = s.cb + vz(s.h, s.cb.cs),
      mat = get_material(b, s)
    write_povray_object(buf, "cone", mat, bot, s.rb, top, s.rt)
    void_ref(b)
  end

realize(b::POVRay, s::Cylinder) =
  let buf = buffer(b),
      mat = get_material(b, s)
    write_povray_object(buf, "cylinder", mat, [0,0,0], [0,0,s.h], s.r) do
      write_povray_matrix(buf, s.cb)
    end
    void_ref(b)
  end

write_povray_mesh(buf::IO, mat, points, closed_u, closed_v, smooth_u, smooth_v) =
  let si = size(points, 1),
      sj = size(points, 2),
      pts = points,
      vcs = add_z.(points, 1) .- points,
      idxs = quad_grid_indexes(si, sj, closed_u, closed_v)
#    let ps = reshape(permutedims(pts), :)
#      delete_all_shapes(autocad)
#      for tr in idxs
#        surface_polygon(ps[map(x->x+1,tr)], backend=autocad)
#      end
#    end
    write_povray_object(buf, "mesh2", mat) do
      write_povray_object(buf, "vertex_vectors", nothing, si*sj, reshape(permutedims(pts), :)...)
      # Must understand how to handle smoothness along one direction
      if smooth_u && smooth_v
        write_povray_object(buf, "normal_vectors", nothing, si*sj, reshape(permutedims(vcs), :)...)
      end
      write_povray_object(buf, "face_indices", nothing, length(idxs), idxs...)
      # Must understand how to handle smoothness along one direction
      #write_povray_object(buf, "normal_indices", nothing, length(idxs), idxs...)
    end
  end

realize(b::POVRay, s::SurfaceGrid) =
  let buf = buffer(b),
      mat = get_material(b, s)
    write_povray_mesh(
      buf,
      mat,
      convert(AbstractMatrix{<:Loc}, map_division(identity, s, size(s.points,1)-1, size(s.points,2)-1)),
      s.closed_u, s.closed_v,
      s.smooth_u, s.smooth_v)
    void_ref(b)
  end

realize(b::POVRay, s::SweepPath) =
  let vertices = in_world.(path_vertices(s.profile)),
      frames = map_division(identity, s.path, 20),
      buf = buffer(b),
      mat = get_material(b, s)
    write_povray_mesh(
      buf,
      mat,
      [xyz(cx(p), cy(p), cz(p), frame.cs) for p in vertices, frame in frames],
      is_closed_path(s.profile),
      is_closed_path(s.path),
      is_smooth_path(s.profile),
      is_smooth_path(s.path))
    void_ref(b)
  end

# HACK: JUST FOR TESTING
realize(b::POVRay, s::Thicken) =
  realize(b, s.shape)

realize(b::POVRay, s::EmptyShape) = void_ref(b)
realize(b::POVRay, s::UniversalShape) = void_ref(b)
#=
realize(b::POVRay, s::Move) =
    let r = map_ref(s.shape) do r
                POVRayMove(connection(b), r, s.v)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::POVRay, s::Scale) =
    let r = map_ref(s.shape) do r
                POVRayScale(connection(b), r, s.p, s.s)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::POVRay, s::Rotate) =
    let r = map_ref(s.shape) do r
                POVRayRotate(connection(b), r, s.p, s.v, s.angle)
                r
            end
        mark_deleted(s.shape)
        r
    end

=#

realize(b::POVRay, s::UnionShape) =
  let shapes = filter(! is_empty_shape, s.shapes)
    length(shapes) == 1 ?
      (ref(shapes[1]); delete_shape(shapes[1])) :
      write_povray_object(buffer(b), "union", get_material(b, s)) do
        for ss in shapes
          ref(ss)
          delete_shape(ss)
        end
      end
    void_ref(b)
  end

realize(b::POVRay, s::IntersectionShape) =
  write_povray_object(buffer(b), "intersection", get_material(b, s)) do
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end

realize(b::POVRay, s::SubtractionShape3D) =
  write_povray_object(buffer(b), "difference", get_material(b, s)) do
    ref(s.shape)
    delete_shape(s.shape)
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end

# BIM

realize_box(b::POVRay, mat, p, dx, dy, dz) =
  let buf = buffer(b),
      bot = p,
      top = add_xyz(p, dx, dy, dz)
    write_povray_object(buf, "box", mat, bot, top)
    void_ref(b)
  end

realize_prism(b::POVRay, top, bot, side, path::PathSet, h::Real) =
  # PathSets require a different approach
  let buf = buffer(b),
      v = planar_path_normal(path)*h,
      bot_vss = map(path_vertices, path.paths),
      top_vss = map(path_vertices, translate(path, v).paths)
    write_povray_polygons(buf, bot, map(reverse, bot_vss))
    write_povray_polygons(buf, top, top_vss)
    for (bot_vs, top_vs) in zip(bot_vss, top_vss)
      for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
        write_povray_polygon(buf, side, vs)
      end
    end
    void_ref(b)
  end

realize_pyramid_frustum(b::POVRay, top, bot, side, bot_vs::Locs, top_vs::Locs, closed=true) =
  let buf = buffer(b)
    if closed
      write_povray_polygon(buf, bot, reverse(bot_vs))
      write_povray_polygon(buf, top, top_vs)
    end
    for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
      write_povray_polygon(buf, side, vs)
    end
    void_ref(b)
  end

realize_polygon(b::POVRay, mat, path::PathSet, acw=true) =
  acw ?
    write_povray_polygons(buffer(b), mat, map(path_vertices, path.paths)) :
    write_povray_polygons(buffer(b), mat, map(reverse ∘ path_vertices, path.paths))

realize_polygon(b::POVRay, mat, vs::Locs, acw=true) =
  let buf = buffer(b)
    write_povray_polygon(buf, mat, acw ? vs : reverse(vs))
  end

write_povray_polygon(io::IO, mat, vs) =
  write_povray_object(io, "polygon", mat, length(vs)+1, vs..., vs[1])

write_povray_polygons(io::IO, mat, vss) =
  write_povray_object(io, "polygon", mat,
    mapreduce(length, +, vss) + length(vss),
    mapreduce(vs->[vs..., vs[1]], vcat, vss)...)

# Polygons with holes need PathSets in POVRay

subtract_paths(b::POVRay, c_r_w_path::PathSet, c_l_w_path::PathSet, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path.paths..., c_r_op_path),
  path_set(c_l_w_path.paths..., c_l_op_path)

subtract_paths(b::POVRay, c_r_w_path::Path, c_l_w_path::Path, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path, c_r_op_path),
  path_set(c_l_w_path, c_l_op_path)

#=
POVRay families need to know the different kinds of materials
that go on each surface.
In some cases it might be the same material, but in others, such
as slabs, outside walls, etc, we will have different materials.
=#

const POVRayMaterialFamily = BackendMaterialFamily{POVRayMaterial}
povray_material_family(mat::POVRayMaterial) =
  POVRayMaterialFamily(mat)

const POVRaySlabFamily = BackendSlabFamily{POVRayMaterial}
povray_slab_family(top::POVRayMaterial, bot::POVRayMaterial=top, side::POVRayMaterial=bot) =
  POVRaySlabFamily(top, bot, side)

const POVRayRoofFamily = BackendRoofFamily{POVRayMaterial}
povray_roof_family(top::POVRayMaterial, bot::POVRayMaterial=top, side::POVRayMaterial=bot) =
  POVRayRoofFamily(top, bot, side)

const POVRayWallFamily = BackendWallFamily{POVRayMaterial}
povray_wall_family(right::POVRayMaterial, left::POVRayMaterial=right) =
  POVRayWallFamily(right, left)

export povray_material_family,
       povray_slab_family,
       povray_roof_family,
       povray_wall_family,
       default_povray_material

povray_stone = povray_include("stones2.inc", "texture", "T_Stone35")
povray_metal = povray_include("textures.inc", "texture", "Chrome_Metal")
#povray_wood = povray_include("textures.inc", "texture", "DMFWood1")
povray_wood = povray_include("woods.inc", "texture", "T_Wood10")
povray_glass = povray_include("textures.inc", "material", "M_Glass")

export povray_stone, povray_metal, povray_wood, povray_glass
set_backend_family(default_wall_family(), povray,
  povray_wall_family(povray_material("InteriorWall70", gray=0.7)))
set_backend_family(default_slab_family(), povray,
  povray_slab_family(povray_material("GenericFloor20", gray=0.2), povray_material("GenericCeiling80", gray=0.8)))
set_backend_family(default_roof_family(), povray,
  povray_roof_family(povray_material("GenericFloor20", gray=0.2), povray_material("GenericCeiling30", gray=0.3)))
set_backend_family(default_beam_family(), povray, povray_material_family(povray_metal))
set_backend_family(default_column_family(), povray, povray_material_family(povray_metal))
set_backend_family(default_door_family(), povray, povray_material_family(povray_wood))
set_backend_family(default_panel_family(), povray, povray_material_family(povray_glass))
set_backend_family(default_table_family(), povray, povray_material_family(povray_wood))
set_backend_family(default_chair_family(), povray, povray_material_family(povray_wood))
set_backend_family(default_table_chair_family(), povray, povray_material_family(povray_wood))
set_backend_family(default_truss_node_family(), povray, povray_material_family(povray_metal))
set_backend_family(default_truss_bar_family(), povray, povray_material_family(povray_metal))

use_family_in_layer(b::POVRay) = true

# Layers

current_layer(b::POVRay) =
  default_povray_material()

current_layer(layer, b::POVRay) =
  default_povray_material(layer)

backend_create_layer(b::POVRay, name::String, active::Bool, color::RGB) =
  begin
    @assert active
    povray_material(name, red=red(color), green=green(color), blue=blue(color))
  end

#=
create_ground_plane(shapes, material=default_povray_ground_material()) =
  if shapes == []
    error("No shapes selected for analysis. Use add-povray-shape!.")
  else
    let (p0, p1) = bounding_box(union(shapes)),
        (center, ratio) = (quad_center(p0, p1, p2, p3),
                  distance(p0, p4)/distance(p0, p2));
     ratio == 0 ?
      error("Couldn"t compute height. Use add-povray-shape!.") :
      let pts = map(p -> intermediate_loc(center, p, ratio*10), [p0, p1, p2, p3]);
         create_surface_layer(pts, 0, ground_layer(), material)
        end
       end
  end
        w = max(floor_extra_factor()*distance(p0, p1), floor_extra_width())
        with(current_layer,floor_layer()) do
          box(xyz(min(p0.x, p1.x)-w, min(p0.y, p1.y)-w, p0.z-1-floor_distance()),
              xyz(max(p0.x, p1.x)+w, max(p0.y, p1.y)+w, p0.z-0-floor_distance()))
        end
      end
    end

=#

#=
#FIXME define the family parameters for beams
realize(b::POVRay, s::Beam) =
    ref(right_cuboid(s.p0, 0.2, 0.2, s.p1, 0))
=#
realize(b::POVRay, s::Union{Door, Window}) =
  void_ref(b)

####################################################

export_to_povray(path::String, b::POVRay=current_backend()) =
  let buf = b.buffer()
    # First pass, to fill material dictionary
    take!(buf)
    # We cannot do this because the array might be updated during the iteration
    for s in b.shapes
      mark_deleted(s)
    end
    i = 1
    while i <= length(b.shapes)
      ref(b.shapes[i])
      i += 1
    end
    open(path, "w") do out
      # write the sky
      write(out, b.sky)
      # write the ground
      write(out, b.ground)
      # write materials
      for (k,v) in b.materials
        write_povray_definition(out, k)
      end
      # write the objects
      write(out, String(take!(buf)))
      # write the view
      write_povray_camera(out, b.camera, b.target, b.lens)
    end
  end

#=
Simulations need to be done on a temporary folder, so that we can have multiple
simulations running at the same time.
=#

povray_simulation_path() = joinpath(string(@__DIR__), "FOO.pov")
  #=
  let (path, io) = mktemp(mktempdir(tempdir(), prefix="POVRay_"))
    close(io)
    path
  end
 =#
export povray_folder
const povray_folder = Parameter("C:/Program Files/POV-Ray/v3.7/bin/")

povray_cmd(cmd::AbstractString="pvengine64") = povray_folder() * cmd

##########################################
render_view(path::String, b::POVRay) =
  let povpath = path_replace_suffix(path, ".pov")
    @info povpath
    export_to_povray(povpath)
    film_active() ?
      run(`$(povray_cmd()) +A +HR Width=$(render_width()) Height=$(render_height()) -D /EXIT /RENDER $(povpath)`, wait=true) :
      run(`$(povray_cmd()) +A +HR Width=$(render_width()) Height=$(render_height()) /RENDER $(povpath)`, wait=false)
  end

export clay_model
clay_model(level::Loc=z(0), b::POVRay=povray) =
  begin
    b.sky = povray_clay_settings_string()
    b.ground = "plane { y, $(cz(level)) texture{ pigment { color rgb 3 } finish { reflection 0 ambient 0 }}}\n"
  end


text:
type OptimizationState
    state::Any
    f_state::Float64
    iteration::Int64
    metadata::Dict
end
function OptimizationState(s::Any, f::Float64, i::Integer)
    OptimizationState(s, f, i, Dict())
end

type OptimizationTrace
    states::Vector{OptimizationState}
end
OptimizationTrace() = OptimizationTrace(Array(OptimizationState, 0))

type OptimizationResults
    method::String
    initial_x::Vector{Float64}
    minimum::Vector{Float64}
    f_minimum::Float64
    iterations::Int64
    converged::Bool
    trace::OptimizationTrace
    f_calls::Int64
    g_calls::Int64
end

immutable DifferentiableFunction
    f::Function
    g!::Function
    fg!::Function
end

immutable TwiceDifferentiableFunction
    f::Function
    g!::Function
    fg!::Function
    h!::Function
end

function show(io::IO, o_trace::OptimizationState)
    print(io, "State of Optimization Algorithm\n")
    print(io, " * Iteration: $(o_trace.iteration)\n")
    print(io, " * State: $(o_trace.state)\n")
    print(io, " * f(State): $(o_trace.f_state)\n")
    print(io, " * Additional Information:\n")
    for (key, value) in o_trace.metadata
        @printf io "   * %s: %s\n" key value
    end
end
function push!(tr::OptimizationTrace, s::OptimizationState)
    push!(tr.states, s)
end
function getindex(tr::OptimizationTrace, i::Integer)
    getindex(tr.states, i)
end
function setindex!(tr::OptimizationTrace,
                   s::OptimizationState,
                   i::Integer)
    setindex!(tr.states, s, i)
end
function show(io::IO, o_trace::OptimizationTrace)
    for state in o_trace.states
        show(io, state)
    end
end

function show(io::IO, results::OptimizationResults)
    print(io, "Results of Optimization Algorithm\n")
    print(io, " * Algorithm: $(results.method)\n")
    print(io, " * Starting Point: $(results.initial_x)\n")
    print(io, " * Minimum: $(results.minimum)\n")
    print(io, " * Value of Function at Minimum: $(results.f_minimum)\n")
    print(io, " * Iterations: $(results.iterations)\n")
    print(io, " * Self-Reported Convergence: $(results.converged)\n")
    print(io, " * Objective Function Calls: $(results.f_calls)\n")
    print(io, " * Gradient Call: $(results.g_calls)")
end

function DifferentiableFunction(f::Function)
    function g!(x::Vector, storage::Vector)
        Calculus.finite_difference!(f, x, storage, :central)
        return
    end
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    return DifferentiableFunction(f, g!, fg!)
end

function DifferentiableFunction(f::Function, g!::Function)
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    return DifferentiableFunction(f, g!, fg!)
end

function TwiceDifferentiableFunction(f::Function)
    function g!(x::Vector, storage::Vector)
        Calculus.finite_difference!(f, x, storage, :central)
        return
    end
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    function h!(x::Vector, storage::Matrix)
        Calculus.finite_difference_hessian!(f, x, storage)
        return
    end
    return TwiceDifferentiableFunction(f, g!, fg!, h!)
end

function TwiceDifferentiableFunction(f::Function,
                                     g!::Function,
                                     h!::Function)
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    return TwiceDifferentiableFunction(f, g!, fg!, h!)
end


text:
export id,∂id

"identity function"
function id(rᵢ)
  rₒ=rᵢ
end
function ∂id(rᵢ,rₒ)
  ∂ij=Diagm(ones(rₒ))
end


text:
#=
The conversion functions for each field's allowed units to units into
SI units.
=#
const unitMappings = Dict{Symbol, Dict{String, Function}}(
    :pressure => Dict{String, Function}(
        "kPa"     => (value::Float64) -> value,
        "bar"     => (value::Float64) -> value * 100.0,
        "inHg"    => (value::Float64) -> value * 3.38639,
        "kgm/s2"  => (value::Float64) -> value,
        "mbar"    => (value::Float64) -> value * 0.1,
        "mmHg"    => (value::Float64) -> value * 0.133322,
        "N/m2"    => (value::Float64) -> value * 0.001,
        "Pa"      => (value::Float64) -> value * 0.001,
        "Pascals" => (value::Float64) -> value * 0.001
    ),
    :relativeHumidity => Dict{String, Function}(
        "percent" => (value::Float64) -> value
    ),
    :solarRadiation => Dict{String, Function}(
        "W/m2" => (value::Float64) -> value
    ),
    :temperature => Dict{String, Function}(
        "C"          => (value::Float64) -> value,
        "Celsius"    => (value::Float64) -> value,
        "F"          => (value::Float64) -> (value - 32.0) * (5.0 / 9.0),
        "Fahrenheit" => (value::Float64) -> (value - 32.0) * (5.0 / 9.0),
        "K"          => (value::Float64) -> value - 273.0,
        "Kelvin"     => (value::Float64) -> value - 273.0
    ),
    :windDirection => Dict{String, Function}(
        "degrees" => (value::Float64) -> value
    ),
    :windSpeed => Dict{String, Function}(
        "m/s"   => (value::Float64) -> value,
        "Knots" => (value::Float64) -> value * 0.514444,
        "k/hr"  => (value::Float64) -> value * 1000.0 / (60 * 60),
        "mph"   => (value::Float64) -> (value / 0.62137) * 1000.0 / (60 * 60)
    )
)

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer
using RadiativeTransfer.Scattering
using CUDA
using Test
# using TensorOperations
using KernelAbstractions
using KernelAbstractions.Extras
using StaticArrays
using BenchmarkTools
using TensorOperations
using LinearAlgebra
using NNlib

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

# Needs some warning if memory is getting too large !
FT = Float32
n = 32
nSpec = 20000

# Create CPU Matrices:
r⁻⁺_ = (randn(FT, n, n, nSpec));
t⁺⁺_ = (randn(FT, n, n, nSpec));
r⁺⁻_ = (randn(FT, n, n, nSpec));
t⁻⁻_ = (randn(FT, n, n, nSpec));
R⁻⁺_ = (randn(FT, n, n, nSpec));
T⁺⁺_ = (randn(FT, n, n, nSpec));
R⁺⁻_ = (randn(FT, n, n, nSpec));
T⁻⁻_ = (randn(FT, n, n, nSpec));

# And move to GPU as CuArray
r⁻⁺ = CuArray(r⁻⁺_);
t⁺⁺ = CuArray(t⁺⁺_);
r⁺⁻ = CuArray(r⁺⁻_);
t⁻⁻ = CuArray(t⁻⁻_);
R⁻⁺ = CuArray(R⁻⁺_);
T⁺⁺ = CuArray(T⁺⁺_);
R⁺⁻ = CuArray(R⁺⁻_);
T⁻⁻ = CuArray(T⁻⁻_);

# This is still a bit stupid but noy sure how to make this better...
I_static_ = Diagonal{FT}(ones(n))
I_static  = Diagonal(CuArray(I_static_));

# This should come from poltype later:
D_IQUV  = FT[1, 1, -1, -1]

# D = Diagonal(repeat(pol_type.D, size(qp_μ)[1]))
# Needs to be using Architecture framework in rtm main code and pol_type.D as above:
# Pol numbers:
nP = 4
D_ = (zeros(FT, n, n, 1));
_arr = (repeat(D_IQUV, n ÷ nP));
for i = 1:n
    D_[i,i,1] = _arr[i];
end
D = CuArray(D_)

function rt_interaction!(R⁻⁺::AbstractArray{FT,3},
                         T⁺⁺::AbstractArray{FT,3},
                         R⁺⁻::AbstractArray{FT,3},
                         T⁻⁻::AbstractArray{FT,3},
                         r⁻⁺::AbstractArray{FT,3},
                         t⁺⁺::AbstractArray{FT,3},
                         r⁺⁻::AbstractArray{FT,3},
                         t⁻⁻::AbstractArray{FT,3},
                         I_static::AbstractArray) where {FT}
    aux1 = similar(R⁻⁺)
    aux2 = similar(R⁻⁺)
    aux3 = similar(R⁻⁺)
    # Compute M1 = (I - R⁺⁻ * r⁻⁺) \ T⁺⁺
    aux1[:] = I_static .- R⁺⁻ ⊠ r⁻⁺
    batch_solve!(aux2, aux1, T⁺⁺)

    # Compute t_R⁻⁺ = R⁻⁺ + T⁻⁻ * r⁻⁺ * M1
    aux1[:] = r⁻⁺ ⊠ aux2  # r⁻⁺ * M1
    aux3[:] = T⁻⁻ ⊠ aux1  #
    R⁻⁺[:]  = R⁻⁺ + aux3

    # T⁺⁺ = t⁺⁺ * M1
    T⁺⁺[:] = t⁺⁺ ⊠ aux2

    # Repeating for mirror-reflected directions
    # Compute M1 = (I - r⁻⁺ * R⁺⁻) \ t⁻⁻
    aux1[:] = I_static .- r⁻⁺ ⊠ R⁺⁻
    batch_solve!(aux2, aux1, t⁻⁻)

    # t_R⁺⁻ = r⁺⁻ + t⁺⁺ * R⁺⁻ * M1
    aux3[:] = R⁺⁻ ⊠ aux2
    aux1[:] = t⁺⁺ ⊠ aux3
    R⁺⁻[:]  = r⁺⁻ + aux1
    T⁻⁻[:]  = T⁺⁺ ⊠ aux2

    # synchronize()
end

function rt_doubling!(ndoubl::Int,
                      r⁻⁺::AbstractArray{FT,3},
                      t⁺⁺::AbstractArray{FT,3},
                      r⁺⁻::AbstractArray{FT,3},
                      t⁻⁻::AbstractArray{FT,3},
                      D::AbstractArray{FT,3},
                      I_static::AbstractArray) where {FT}
    # # ToDo: Important output doubling applied to elemental layer, using same variables r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺ (can be renamed to t⁺⁺, etc)
    # Need to check with paper nomenclature. This is basically eqs. 23-28 in vSmartMOM but using simplifications in eq. 29-32)

    aux1 = similar(t⁺⁺)
    aux2 = similar(t⁺⁺)
    aux3 = similar(t⁺⁺)

    for n = 1:ndoubl

        # M1 = (I - r⁻⁺ * r⁻⁺) \ t⁺⁺
        aux1[:] = I_static .- r⁻⁺ ⊠ r⁻⁺     # (I - r⁻⁺ * r⁻⁺)
        batch_solve!(aux2, aux1, t⁺⁺)   # M1 = (I - r⁻⁺ * r⁻⁺) \ t⁺⁺

        # r⁻⁺[:] = r⁻⁺ + t⁺⁺ * r⁻⁺ * M1
        aux1[:] = r⁻⁺ ⊠ aux2               # r⁻⁺ * M1
        aux3[:] = t⁺⁺ ⊠ aux1               # t⁺⁺ * r⁻⁺ * M1
        r⁻⁺[:]  = r⁻⁺ + aux3               # r⁻⁺[:] = r⁻⁺ + t⁺⁺ * r⁻⁺ * M1

        # t⁺⁺[:] = t⁺⁺ * M1
        aux1[:] = t⁺⁺ ⊠ aux2           # t⁺⁺ * M1
        t⁺⁺[:]  = aux1                   # t⁺⁺[:] = t⁺⁺ * M1
    end

    # After doubling, revert D(DR)->R, where D = Diagonal{1,1,-1,-1}
    r⁻⁺[:] = D ⊠ r⁻⁺ ⊠ D
    # Using r⁺⁻ = Dr⁻⁺D
    r⁺⁻[:] = D ⊠ r⁻⁺ ⊠ D
    # Using t⁻⁻ = Dt⁺⁺D
    t⁻⁻[:] = D ⊠ t⁺⁺ ⊠ D
    return nothing
end

function rt_elemental!(pol_type, dτ, ϖ, Z⁺⁺, Z⁻⁺, m,
                       ndoubl::Int, qp_μ, wt_μ,
                       r⁻⁺::AbstractArray{FT,3},
                       t⁺⁺::AbstractArray{FT,3},
                       r⁺⁻::AbstractArray{FT,3},
                       t⁻⁻::AbstractArray{FT,3},
                       D::AbstractArray{FT,3},
                       I_static::AbstractArray) where {FT}

    aux1 = similar(Z⁻⁺)

    qp_μN = repeat(qp_μ, 3)
    wt_μN= repeat(wt_μ, 3)
    # println(typeof(qp_μN))
    #  reduce(vcat, (fill.(qp_μ,[pol_type.n])))

    # wt_μ4 = reduce(vcat, (fill.(wt_μ,[pol_type.n])))
    Nquadn = length(qp_μN)

    wct = m == 0 ? 0.50 * ϖ * wt_μN  : 0.25 * ϖ * wt_μN

    d_qp = r⁺⁻ isa CuArray{Float32,3} ? CuArray(Diagonal(1 ./ Array(qp_μN))) : Diagonal(1 ./ Array(qp_μN))
    # d_wct = r⁺⁻ isa CuArray{Float32,3} ? CuArray(Diagonal(1 ./ Array(qp_μN))) : Diagonal(1 ./ Array(qp_μN))
    d_wct = r⁺⁻ isa CuArray{Float32,3} ? CuArray(Diagonal(Array(wct))) : Diagonal(Array(wct))

    aux1 = d_qp ⊠ Z⁻⁺
    r⁻⁺ = aux1 ⊠ (d_wct * dτ)

    t⁺⁺ = I_static .- (d_qp ⊠ ((I_static .- Z⁺⁺ ⊠ d_wct) * dτ))

    if ndoubl < 1
        r⁺⁻ = D ⊠ r⁻⁺ ⊠ D
        t⁻⁻ = D ⊠ t⁺⁺ ⊠ D
    else
        r⁻⁺ = D ⊠ r⁻⁺
    end
end
function batch_solve!(X::CuArray{FT,3}, A::CuArray{FT,3}, B::CuArray{FT,3}) where {FT}
    temp = similar(A)
    pivot, info   = CUBLAS.getrf_strided_batched!(A, true);
    synchronize()
    getri_strided_batched!(A, temp, pivot); # inv stored in aux1
    NNlib.batched_mul!(X, temp, B)
    synchronize()
end

# batch Matrix inversion for CPU
# function batch_solve!(X::CuArray{FT,3}, A::CuArray{FT,3}, B::CuArray{FT,3}) where {FT}
#    temp = similar(A)
#    pivot, info   = CUBLAS.getrf_strided_batched!(A, true);
#    synchronize()
#    getri_strided_batched!(A, temp, pivot); # inv stored in aux1
#    NNlib.batched_mul!(X, temp, B)
#    synchronize()
    # getri_strided_batched(A, temp, pivot) # inv stored in aux1
    # X[:] = temp ⊠ B  # inv(I - R⁺⁻ * r⁻⁺) * T⁺⁺
    # synchronize()
# end

# batch Matrix inversion for CPU
function batch_solve!(X::AbstractArray{FT,3}, A::AbstractArray{FT,3}, B::AbstractArray{FT,3}) where {FT}
    for i = 1:size(A, 3)
        @views ldiv!(X[:,:,i], qr!(A[:,:,i]), B[:,:,i])
    end
end

# CUDA has no strided batched getri, but we can at least avoid constructing costly views (copied this over from gertf)
function getri_strided_batched(A::AbstractArray{Float32,3}, C::AbstractArray{Float32,3}, pivotArray::CuMatrix{Cint})
    m, n = size(A, 1), size(A, 2)
    if m != n
        throw(DimensionMismatch("All matrices must be square!"))
    end
    n = size(A, 1)
    ldc = max(1, stride(C, 2))
    lda = max(1, stride(A, 2))
    info = CUDA.zeros(Cint, size(A, 3))
    Cptrs = CUBLAS.unsafe_strided_batch(C)
    Aptrs = CUBLAS.unsafe_strided_batch(A)
    CUBLAS.cublasSgetriBatched(CUBLAS.handle(), n, Aptrs, lda, pivotArray, Cptrs, ldc, info, size(A, 3))
    return nothing
end

# CUDA has no strided batched getri, but we can at least avoid constructing costly views (copied this over from gertf)
function getri_strided_batched(A::AbstractArray{Float64,3}, C::AbstractArray{Float64,3}, pivotArray::CuMatrix{Cint})
    m, n = size(A, 1), size(A, 2)
    if m != n
        throw(DimensionMismatch("All matrices must be square!"))
    end
    n = size(A, 1)
    ldc = max(1, stride(C, 2))
    lda = max(1, stride(A, 2))
    info = CUDA.zeros(Cint, size(A, 3))
    # Cptrs = CUBLAS.unsafe_strided_batch(C)
    Cptrs = CUBLAS.unsafe_strided_batch(C)
    Aptrs = CUBLAS.unsafe_strided_batch(A)
    CUBLAS.cublasDgetriBatched(CUBLAS.handle(), n, Aptrs, lda, pivotArray, Cptrs, ldc, info, size(A, 3))
    return nothing
end

@show n, nSpec
function run_rt_interaction(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
    rt_interaction!(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
    synchronize()
end

function run_rt_doubling(nDoubling, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)
    rt_doubling!(nDoubling, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)
    synchronize()
end

function run_rt_elemental(pol_type, dτ, ϖ, Z⁺⁺, Z⁻⁺, m,
                          ndoubl, qp_μ, wt_μ,
                          r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)

    rt_elemental!(pol_type, dτ, ϖ, Z⁺⁺, Z⁻⁺, m,
                  ndoubl, qp_μ, wt_μ,
                  r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)
    synchronize()
end

# Testing RT Interaction time:
println("RT Interaction GPU time:")
@btime run_rt_interaction(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
println("RT Interaction CPU time (uses multi-threading in LAPACK/BLAS!):")
@btime run_rt_interaction(R⁻⁺_, T⁺⁺_, R⁺⁻_, T⁻⁻_, r⁻⁺_, t⁺⁺_, r⁺⁻_, t⁻⁻_, I_static_)

@testset "GPU-CPU consistency" begin
    for (matGPU, matCPU) in ((R⁻⁺, R⁻⁺_),
                             (T⁺⁺, T⁺⁺_),
                             (R⁺⁻, R⁺⁻_),
                             (T⁻⁻, T⁻⁻_),
                             (r⁻⁺, r⁻⁺_),
                             (t⁺⁺, t⁺⁺_),
                             (r⁺⁻, r⁺⁻_),
                             (t⁻⁻, t⁻⁻_))
        @test Array(matGPU) ≈ matCPU
    end
end

nDoubling = 4
@show nDoubling
# Testing Doubling:
println("RT Doubling GPU time:")
@btime run_rt_doubling(nDoubling, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)
println("RT Doubling CPU time:")
@btime run_rt_doubling(nDoubling, r⁻⁺_, t⁺⁺_, r⁺⁻_, t⁻⁻_, D_, I_static_)

@testset "GPU-CPU consistency" begin
    for (matGPU, matCPU) in ((R⁻⁺, R⁻⁺_),
                                (T⁺⁺, T⁺⁺_),
                                (R⁺⁻, R⁺⁻_),
                                (T⁻⁻, T⁻⁻_),
                                (r⁻⁺, r⁻⁺_),
                                (t⁺⁺, t⁺⁺_),
                                (r⁺⁻, r⁺⁻_),
                                (t⁻⁻, t⁻⁻_))
        @test Array(matGPU) ≈ matCPU
    end
end

pol_type = Stokes_IQU{FT}()

m = 0
ndoubl = 0
weight = m == 0 ? 0.5 : 1.0
dτ = 2.458321474082143e-7
ϖ = 1.0
Z⁺⁺ = randn(57, 57, 3)
Z⁻⁺ = randn(57, 57, 3)
vza = [60., 45., 30., 15., 0., 15., 30., 45., 60.]
vaz = [180., 180., 180., 180., 0., 0., 0., 0., 0.]
sza = 60.
Ltrunc = 14
Nquad, qp_μ, wt_μ = rt_set_streams(vSmartMOM.RadauQuad(), Ltrunc, sza, vza);
I_static_ = Diagonal{FT}(ones(57))
I_static  = Diagonal(CuArray(I_static_));

n = 57
D_ = (zeros(FT, n, n, 1));
_arr = (repeat(D_IQUV, n ÷ nP + 1));
for i = 1:n
    D_[i,i,1] = _arr[i];
end
D = CuArray(D_)

println("RT Elemental GPU time:")
@btime run_rt_elemental(pol_type, dτ, ϖ, CuArray(Z⁺⁺), CuArray(Z⁻⁺), m,
                        ndoubl, CuArray(qp_μ), CuArray(wt_μ),
                        r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, D, I_static)

println("RT Elemental CPU time:")
@btime run_rt_elemental(pol_type, dτ, ϖ, Z⁺⁺, Z⁻⁺, m,
                        ndoubl, qp_μ, wt_μ,
                        r⁻⁺_, t⁺⁺_, r⁺⁻_, t⁻⁻_, D_, I_static_)

@testset "GPU-CPU consistency" begin
    for (matGPU, matCPU) in ((R⁻⁺, R⁻⁺_),
                                (T⁺⁺, T⁺⁺_),
                                (R⁺⁻, R⁺⁻_),
                                (T⁻⁻, T⁻⁻_),
                                (r⁻⁺, r⁻⁺_),
                                (t⁺⁺, t⁺⁺_),
                                (r⁺⁻, r⁺⁻_),
                                (t⁻⁻, t⁻⁻_))
        @test Array(matGPU) ≈ matCPU
    end
end

# @inline function unsafe_strided_batch_fake(strided::CuArray{T}, batchsize::Int) where {T}
#    # batchsize = last(size(strided))
#    stride = prod(size(strided)[1:end - 1])
#    ptrs = [pointer(strided, stride + 1) for i in 1:batchsize]
#    return CuArray(ptrs)
# end