text:
using Test: @test, @testset, @test_logs, @test_throws

using HTTP: HTTP, Form, Request, Response
using JSON: JSON

using BrokenRecord: BrokenRecord, FORMAT, configure!, playback

const BR = BrokenRecord
const url = "https://httpbin.org"

@testset "BrokenRecord.jl" begin
    mktempdir() do dir
        path = joinpath(dir, "test1.json")
        @test playback(() -> 1, path) == 1
        @test isfile(path)
        @test JSON.parsefile(path) == Dict("responses" => [], "format" => string(FORMAT))

        path = joinpath(dir, "test2.json")
        resp1 = playback(() -> HTTP.get("$url/get"), path)
        @test isfile(path)
        resp2 = playback(() -> HTTP.get("$url/get"), path)
        @test resp1.body == resp2.body

        path = joinpath(dir, "test3.json")
        playback(() -> HTTP.post("$url/anything"; body="hi"), path)
        # Wrong body.
        @test_throws Exception playback(() -> HTTP.post("$url/anything"; body="bye"), path)
        # Wrong method.
        @test_throws Exception playback(() -> HTTP.put("$url/anything"; body="hi"), path)
        # Wrong headers.
        @test_throws Exception playback(path) do
            HTTP.post("$url/anything"; body="hi", headers=["foo" => "bar"])
        end
        # Wrong path.
        @test_throws Exception playback(() -> HTTP.post("$url/post"; body="hi"), path)
        # Wrong query.
        @test_throws Exception playback(path) do
            HTTP.post("$url/anything"; query=Dict("foo" => "bar"), body="hi")
        end

        path = joinpath(dir, "test4.json")
        resp1 = playback(() -> HTTP.get("$url/get"), path)
        resp2 = playback(path; ignore_headers=["foo"]) do
            HTTP.get("$url/get"; headers=["foo" => "bar"])
        end
        @test resp1.body == resp2.body
        resp3 = playback(path; ignore_query=["foo"]) do
            HTTP.get("$url/get"; query=Dict("foo" => "bar"))
        end
        @test resp1.body == resp3.body

        path = joinpath(dir, "test5.json")
        resp1 = playback(() -> HTTP.get("$url"), path)
        resp2 = playback(() -> HTTP.get("$url"), path)
        @test resp1.body == resp2.body

        path = joinpath(dir, "test6.json")
        playback(path) do
            HTTP.get("$url/get")
            HTTP.post("$url/post")
        end
        @test_throws Exception playback(() -> HTTP.get("$url/get"), path)

        path = joinpath(dir, "test7.json")
        playback(path) do
            resp = HTTP.get("$url/get")
            empty!(resp.body)
        end
        @test !isempty(playback(() -> HTTP.get("$url/get").body, path))
    end

    mktempdir() do dir
        playback(() -> (), joinpath(dir, "foo bar.json"))
        @test !isfile(joinpath(dir, "foo bar.json"))
        @test isfile(joinpath(dir, "foo_bar.json"))
    end

    mktempdir() do dir
        configure!(; path=dir, ignore_headers=["foo"], ignore_query=["bar"])
        path = "test.json"
        resp1 = playback(() -> HTTP.get("$url/get"), path)
        @test !isfile("test.json")
        @test isfile(joinpath(dir, "test.json"))
        resp2 = playback(path) do
            HTTP.get("$url/get"; headers=["foo" => "bar"], query=Dict("bar" => "baz"))
        end
        @test resp1.body == resp2.body

        path = joinpath(dir, "test2.json")
        resp1 = playback(path) do
            HTTP.get("$url/response-headers"; query=Dict("foo" => "bar"))
        end
        @test HTTP.header(resp1, "foo") == "bar"
        resp2 = playback(path) do
            HTTP.get("$url/response-headers"; query=Dict("foo" => "bar"))
        end
        @test isempty(HTTP.header(resp2, "foo"))
    end

    mktempdir() do dir
        path = joinpath(dir, "test.json")
        playback(path) do
            open(@__FILE__) do f
                HTTP.post("$url/post"; body=Form(Dict(:file => f)))
            end
        end
        playback(path) do
            open(@__FILE__) do f
                @test_logs (:warn, r"streamed") HTTP.post("$url/post"; body=Form(Dict(:file => f)))
            end
        end
    end

    @testset "Storage backends" begin
        @testset "Backend identification" begin
            @test BR.get_storage("foo", "yml") == (BR.YAMLStorage, "foo.yml")
            @test BR.get_storage("foo.yml", "yml") == (BR.YAMLStorage, "foo.yml")
            @test BR.get_storage("foo.yml", "json") == (BR.YAMLStorage, "foo.yml")
            @test_throws Exception BR.get_storage("foo.hello", "yml")
        end

        @testset "Converting Response to/from Dict" begin
            resp = Response(200, ["foo" => "bar"]; body="baz")
            resp.request = Request("GET", "/", ["a" => "b"], "ok"; version=v"1.1")
            @test BR.resp_to_dict(resp) == Dict(
                "status" => 200,
                "headers" => ["foo" => "bar"],
                "body" => "baz",
                "version" => "1.1.0",
                "request" => Dict(
                    "method" => "GET",
                    "headers" => ["a" => "b"],
                    "target" => "/",
                    "body" => "ok",
                    "txcount" => 0,
                    "version" => "1.1.0",
                ),
            )
        end

        exts = Dict(v => k for (k, v) in BR.EXTENSIONS)
        mktempdir() do dir
            @testset "$(nameof(S))" for S in keys(exts)
                ext = exts[S]
                path = joinpath(dir, "test1.$ext")
                resp1 = playback(() -> HTTP.get("$url/get"), path)
                resp2 = playback(() -> HTTP.get("$url/get"), path)
                @test resp1.body == resp2.body
            end
        end
    end
end


text:
"""
   portfolio_data()

Returns a tuple `(R, Rtest)` with data for the portfolio optimization
example in section 17.1.3.
`R` is a 2000 x 20 matrix with daily returns over a period of 2000 days.
The first 19 columns are returns for 19 stocks;
the last column is for a risk-free asset.
`Rtest` is a 500 x 20 matrix with daily returns over a different period
of 500 days.
"""
function portfolio_data()
    pth = splitdir(pathof(VMLS))[1]
    prices = readdlm(joinpath(pth, "portfolio_data.txt"), comments = true);
    p_changes = prices[2:end,:] - prices[1:end-1,:];
    returns = p_changes ./ prices[1:end-1,:];
    R = returns[1:2000,:];
    Rtest = returns[2001:2500,:];
    return R, Rtest
end

ext: jl

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
module MaxHelpingHandPersistentMusicFadeTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/PersistentMusicFadeTrigger" PersistentMusicFadeTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight;
	direction::String="leftToRight", fadeA::Number=0.0, fadeB::Number=1.0, parameter::String="")

const placements = Ahorn.PlacementDict(
    "Persistent Music Fade (max480's Helping Hand)" => Ahorn.EntityPlacement(
        PersistentMusicFadeTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::PersistentMusicFadeTrigger)
    return Dict{String, Any}(
        "direction" => Maple.music_fade_trigger_directions
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module MountainCarExperiment

import Random

using GVFHordes
using Curiosity
using MinimalRLCore
using SparseArrays
using ProgressMeter

import Flux: Descent
const SRCU = Curiosity.SRCreationUtils

const MCU = Curiosity.MountainCarUtils

# VERY IMPORTANT THAT LEARNED POLICY NAMES MATCHES ORDER OF EVAL SET!!!!!
default_args() =
    Dict(
        "steps" => 10,
        "seed" => 1,

        #Tile coding params used by Rich textbook for mountain car
        "num_tilings" => 16,
        "num_tiles" => 2,
        "behaviour_num_tilings" => 8,
        "behaviour_num_tiles" => 2,
        "behaviour_reward_projector" => "tilecoding",
        "behaviour_eta" => 3^-1,

        "learned_policy" => true,
        "learned_policy_names" => ["Wall","Goal"],
        "behaviour_learner" => "Q",
        "behaviour_update" => "ESARSA",
        "behaviour_opt" => "Descent",
        # "behaviour_rew" => "env",
        "behaviour_gamma" => 0.99,
        "behaviour_lambda" => 0.95,
        "behaviour_w_init" => 0.0,

        # "eta" => 0.5,
        "alpha_init" => 0.1,

        "intrinsic_reward" =>"weight_change",
        "behaviour_trace" => "ReplacingTraces",
        "use_external_reward" => true,
        "exploration_strategy" => "epsilon_greedy",
        "exploration_param" => 1.0,
        "random_first_action" => false,

        "demon_learner" => "SR",
        "demon_update" => "TB",
        "demon_opt" => "Descent",
        "demon_lambda" => 0.9,
        "demon_rep" => "ideal",
        "demon_eta" => 3^-1,
        "demon_num_tilings" => 4,
        "demon_num_tiles" => 4,

        "exploring_starts"=>true,
        "save_dir" => "MountainCarExperiment",
        "logger_keys" => [LoggerKey.EPISODE_LENGTH, LoggerKey.MC_UNIFORM_ERROR, LoggerKey.MC_START_ERROR],
    )

function construct_agent(parsed)
    observation_size = 2
    action_space = 3

    if "eta" in keys(parsed)
        prefixes = ["behaviour","demon"]
        @show "demon_eta being overwritten"
        for prefix in prefixes
            parsed[join([prefix, "eta"], "_")] = parsed["eta"]
        end
    end
    if "alpha_init" in keys(parsed)
        prefixes = ["behaviour","demon"]
        for prefix in prefixes
            parsed[join([prefix, "alpha_init"], "_")] = parsed["alpha_init"]
        end
    end
    if parsed["demon_opt"] == "Descent"
        if "demon_eta" in keys(parsed)
            parsed["demon_eta"] = parsed["demon_eta"] / parsed["num_tilings"]
        end
    end
    if parsed["behaviour_opt"] == "Descent"
        if "behaviour_eta" in keys(parsed)
            parsed["behaviour_eta"] = parsed["behaviour_eta"] / parsed["num_tilings"]
        end
    end

    behaviour_learner = parsed["behaviour_learner"]
    behaviour_lu = parsed["behaviour_update"]

    behaviour_gamma = parsed["behaviour_gamma"]
    behaviour_trace = parsed["behaviour_trace"]
    intrinsic_reward_type = parsed["intrinsic_reward"]
    use_external_reward = parsed["use_external_reward"]

    fc = Curiosity.FeatureSubset(
        Curiosity.SparseTileCoder(parsed["num_tilings"], parsed["num_tiles"], 2),
        1:2)

    demon_reward_features = if parsed["demon_rep"] == "ideal"
        # Curiosity.ActionValueFeatureProjector(
            MCU.IdealDemonFeatures()
            # action_space)
    elseif parsed["demon_rep"] == "tilecoding"
        Curiosity.FeatureProjector(
                        Curiosity.SparseTileCoder(parsed["demon_num_tilings"], parsed["demon_num_tiles"], 2), true)
    else
        throw(ArgumentError("Invalid demon rep for SR"))
    end

    behaviour_reward_features =
    if parsed["behaviour_reward_projector"] == "ideal"
        Curiosity.ActionValueFeatureProjector(
            MCU.IdealDemonFeatures(),
            action_space)
    elseif parsed["behaviour_reward_projector"] == "tilecoding"
        Curiosity.ActionValueFeatureProjector(
                        Curiosity.FeatureProjector(
                        Curiosity.SparseTileCoder(parsed["behaviour_num_tilings"], parsed["behaviour_num_tiles"], 2),
                        false),
                    action_space)
    elseif parsed["behaviour_reward_projector"] == "ideal_state"
        MCU.IdealDemonFeatures()
    else
        throw(ArgumentError("SAD"))
    end

    feature_size = size(fc)
    demon_feature_size = size(demon_reward_features)
    behaviour_feature_size = size(behaviour_reward_features)

    demons = MCU.create_demons(parsed,
                         demon_reward_features)

    demon_learner = Curiosity.get_linear_learner(parsed,
                                                 feature_size,
                                                 action_space,
                                                 demons,
                                                 "demon",
                                                 demon_reward_features)

    exploration_strategy = Curiosity.get_exploration_strategy(parsed, 1:action_space)

    behaviour_num_tasks = 1
    num_SFs = 2
    num_demons = if parsed["behaviour_learner"] ∈ ["GPI"]
        # num_SFs * size(behaviour_reward_projector) * action_space + behaviour_num_tasks
        num_SFs * size(behaviour_reward_features) + behaviour_num_tasks

    elseif parsed["behaviour_learner"] ∈ ["Q"]
        behaviour_num_tasks
    elseif parsed["behaviour_learner"] == "RoundRobin"
        0
    else
        throw(ArgumentError("Hacky thing not working, yay!"))
    end
    behaviour_learner = Curiosity.get_linear_learner(parsed,
                                                     size(fc),
                                                     action_space,
                                                     num_demons,
                                                     behaviour_num_tasks,
                                                     "behaviour",
                                                     behaviour_reward_features)

    behaviour_demons = if behaviour_learner isa GPI
        @assert !(behaviour_reward_features isa Nothing)
        bh_gvf = MCU.make_behaviour_gvf(behaviour_learner,
                                          0.0,
                                          fc,
                                          exploration_strategy)
        pred_horde = GVFHordes.Horde([bh_gvf])
        policies = MCU.get_policies(parsed)

        gvfs = [MCU.steps_to_wall_gvf(policies[1]), MCU.steps_to_goal_gvf(policies[2])]

        SF_policies = [gvfs[1].policy, gvfs[2].policy]
        SF_discounts = [gvfs[1].discount, gvfs[2].discount]
        num_SFs = length(SF_policies)
        SF_horde = SRCU.create_SF_horde(SF_policies, SF_discounts, behaviour_reward_features, 1:action_space)
        Curiosity.GVFSRHordes.SRHorde(pred_horde, SF_horde, num_SFs, behaviour_reward_features)
    elseif behaviour_learner isa QLearner
        bh_gvf = MCU.make_behaviour_gvf(behaviour_learner,
                                          parsed["behaviour_gamma"],
                                          fc,
                                          exploration_strategy)
        GVFHordes.Horde([bh_gvf])
    else
        throw(ArgumentError("goes with which horde? " ))
    end

    random_first_action = parsed["random_first_action"]
    Agent(demons,
          feature_size,
          behaviour_learner,
          behaviour_demons,
          behaviour_gamma,
          demon_learner,
          observation_size,
          action_space,
          intrinsic_reward_type,
          fc,
          use_external_reward,
          exploration_strategy,
          random_first_action)
end

function main_experiment(parsed=default_args(); progress=false, working=false)

    num_steps = parsed["steps"]

    Random.seed!(parsed["seed"])

    normalized = true
    env = MountainCar(0.0,0.0,normalized)

    agent = construct_agent(parsed)

    logger_init_dict = Dict(
        LoggerInitKey.TOTAL_STEPS => num_steps,
        LoggerInitKey.INTERVAL => 1000,
    )

    Curiosity.experiment_wrapper(parsed, logger_init_dict, working) do parsed, logger
        eps = 1
        max_num_steps = num_steps
        steps = Int[]
        logger_start!(logger, env, agent)

        prg_bar = ProgressMeter.Progress(num_steps, "Step: ")
        while sum(steps) < max_num_steps
            is_terminal = false

            max_episode_steps = min(max_num_steps - sum(steps), 5000)
            tr, stp =
                run_episode!(env, agent, max_episode_steps) do (s, a, s_next, r, t)
                    logger_step!(logger, env, agent, s, a, s_next, r, t)
                    if progress
                        next!(prg_bar)
                    end
                end
            logger_episode_end!(logger)

            push!(steps, stp)
            eps += 1
        end
        logger
    end
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using StandardizedPredictors
using StatsModels
using StatsBase
using Statistics
using Test

string_mime(mime, x) = sprint(show, mime, x)

@testset "StandardizedPredictors.jl" begin
    include("centering.jl")
    include("scaling.jl")
    include("zscoring.jl")
end


text:
using Distributions
using Base.Test

v = 7.0
S = eye(2)
S[1, 2] = S[2, 1] = 0.5

W = Wishart(v,S)
IW = InverseWishart(v,S)

for d in [W,IW]
    @test size(d) == size(rand(d))
    @test length(d) == length(rand(d))
    @test typeof(d)(params(d)...) == d
    @test partype(d) == Float64

    ### conversions
    D = typeof(d).name.primary
    @test typeof(convert(D{Float32}, d)) == typeof(D(Float32(v), Array{Float32}(S)))

    @test typeof(convert(D{Float32}, params(d)...)) == typeof(D(Float32(v), Array{Float32}(S)))
end

@test partype(Wishart(7, eye(Float32, 2))) == Float32
@test partype(InverseWishart(7, eye(Float32, 2))) == Float32

@test_approx_eq_eps mean(rand(W,100000)) mean(W) 0.1
@test_approx_eq_eps mean(rand(IW,100000)) mean(IW) 0.1

v = 3.0

@test_approx_eq_eps pdf(Wishart(v,S), S) 0.04507168 1e-8
@test_approx_eq_eps pdf(Wishart(v,S), inv(S)) 0.01327698 1e-8
@test_approx_eq_eps pdf(Wishart(v,inv(S)),S) 0.0148086 1e-8
@test_approx_eq_eps pdf(Wishart(v,inv(S)),inv(S)) 0.01901462 1e-8

@test_approx_eq logpdf(Wishart(v,S), S) log(pdf(Wishart(v,S), S))
@test_approx_eq logpdf(Wishart(v,S), inv(S)) log(pdf(Wishart(v,S), inv(S)))
@test_approx_eq logpdf(Wishart(v,inv(S)), S) log(pdf(Wishart(v,inv(S)), S))
@test_approx_eq logpdf(Wishart(v,inv(S)), inv(S)) log(pdf(Wishart(v,inv(S)), inv(S)))

@test_approx_eq_eps pdf(InverseWishart(v,S), S) 0.04507168 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,S), inv(S)) 0.006247377 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,inv(S)),S)  0.03147137 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,inv(S)),inv(S)) 0.01901462 1e-8

@test_approx_eq logpdf(InverseWishart(v,S), S) log(pdf(InverseWishart(v,S), S))
@test_approx_eq logpdf(InverseWishart(v,S), inv(S)) log(pdf(InverseWishart(v,S), inv(S)))
@test_approx_eq logpdf(InverseWishart(v,inv(S)), S) log(pdf(InverseWishart(v,inv(S)), S))
@test_approx_eq logpdf(InverseWishart(v,inv(S)), inv(S)) log(pdf(InverseWishart(v,inv(S)), inv(S)))


text:
######## annualized return ###################

# returns a single value in simple terms
function annualized_return{T}(ta::TimeArray{T,1})
    exp(sum(diff(log(ta.values)))) ^ (252/length(ta)) - 1
end

function annualized_return(fa::Array{Float64, 1})
    exp(sum(diff(log(fa)))) ^ (252/size(fa,1)) - 1
end

######## equity curve ########################

function equity{T}(ta::TimeArray{T,1}; prices=false)
    prices ?
    TimeArray(ta.timestamp, [1.0, expm1(cumsum(diff(log(ta.values)))) + 1], ["equity"]) :
    TimeArray(ta.timestamp, [expm1(cumsum(ta.values)) + 1], ["equity"])
end


text:
struct Involute{T} <: JointType{T}
	i::SVector{3, T}
	j::SVector{3, T}
	k::SVector{3, T}
	ρ::T
	ρ_inv::T

	function Involute{T}(b_axis::AbstractVector,n_axis::AbstractVector,ρ::T) where {T}
		j, k = map(axis -> normalize(SVector{3}(axis)), (n_axis, b_axis))
        @assert isapprox(k ⋅ j, 0; atol = 100 * eps(T))
        @assert ρ > 0
        new{T}(j × k, j, k, ρ, inv(ρ))
    end
end

@propagate_inbounds function joint_transform(jt::Involute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D, q::AbstractVector)
	rot = RotMatrix(AngleAxis(jt.ρ_inv*q[1],jt.k[1],jt.k[2],jt.k[3],false))
	trans = rot*(jt.ρ*jt.j+q[1]*jt.i)
	Transform3D(frame_after,frame_before,rot,trans)
end

@propagate_inbounds function joint_twist(jt::Involute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
	angular = jt.ρ_inv*jt.k*v[1]
	linear = jt.ρ_inv*q[1]*jt.j*v[1]
	Twist(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_spacial_acceleration(jt::Involute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector, vd::AbstractVector)
	S = promote_eltype(jt, q, v, vd)
	angular = jt.ρ_inv*jt.k*vd[1]
	linear = (jt.ρ_inv*q[1]*jt.j)*vd[1] + (jt.j-jt.ρ_inv*q[1]*jt.i)*jt.ρ_inv*v[1]^2
	SpatialAcceleration{S}(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function bias_acceleration(jt::Involute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(jt, q, v)
    linear = zero(SVector{3,S})
    angular = (jt.j-jt.ρ_inv*q[1]*jt.i)*jt.ρ_inv*v[1]^2
    SpatialAcceleration{S}(frame_after, frame_before, frame_after, angular,linear)
end

@propagate_inbounds function joint_torque!(τ::AbstractVector, jt::Involute, q::AbstractVector, joint_wrench::Wrench)
	R = RotMatrix(AngleAxis(-jt.ρ_inv*q[1],jt.k[1],jt.k[2],jt.k[3],false))
	τ[1] = jt.ρ_inv*q[1]*jt.j'*linear(joint_wrench) + ρ_inv*k'*angular(joint_wrench)
	nothing
end

@inline function motion_subspace(jt::Involute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    S = promote_eltype(jt, q)
	angular = hcat(zero(SMatrix{3,2,S}),jt.k)
	linear = hcat(zero(SVector{3,S}),jt.j,zero(SVector{3,S}))
	GeometricJacobian(frame_after, frame_before, frame_after, angular, linear)
end

@inline function constraint_wrench_subspace(jt::Involute, joint_transform::Transform3D)
    S = promote_eltype(jt, joint_transform)
    angular = hcat(jt.i, jt.j, zero(SVector{3, S}))
    linear = hcat(jt.i, zero(SVector{3, S}), jt.k )
    WrenchMatrix(joint_transform.from, angular, linear)
end


text:
import Base.Dates: trunc, DatePeriod, TimePeriod
import Base.Dates: firstdayofweek, lastdayofweek, firstdayofmonth, lastdayofmonth,
    firstdayofyear, lastdayofyear, firstdayofquarter, lastdayofquarter

# Truncation
# TODO: Just utilize floor code for truncation?
function trunc{P<:DatePeriod}(zdt::ZonedDateTime, ::Type{P})
    ZonedDateTime(trunc(localtime(zdt), P), timezone(zdt))
end
function trunc{P<:TimePeriod}(zdt::ZonedDateTime, ::Type{P})
    local_dt = trunc(localtime(zdt), P)
    utc_dt = local_dt - zdt.zone.offset
    ZonedDateTime(utc_dt, timezone(zdt); from_utc=true)
end
trunc(zdt::ZonedDateTime, ::Type{Millisecond}) = zdt

# Adjusters
for prefix in ("firstdayof", "lastdayof"), suffix in ("week", "month", "year", "quarter")
    func = Symbol(prefix * suffix)
    @eval begin
        $func(dt::ZonedDateTime) = ZonedDateTime($func(localtime(dt)), dt.timezone)
    end
end


text:
# Taken from SnoopCompile
function snoop_vanilla(filename, path)
    code_object = """
    while !eof(STDIN)
        eval(Main, deserialize(STDIN))
    end
    """
    julia_cmd = build_julia_cmd(
        get_backup!(false, nothing), nothing, nothing, nothing, nothing, nothing,
        nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing
    )
    @show julia_cmd
    in, io = open(`$julia_cmd --eval $code_object`, "w", STDOUT)
    serialize(in, quote
        import SnoopCompile
    end)
    # Now that the new process knows about SnoopCompile, it can
    # expand the macro in this next expression
    serialize(in, quote
          SnoopCompile.@snoop1 $filename include($(escape_string(path)))
    end)
    close(in)
    wait(io)
    println("done.")
    nothing
end

function snoop(path, compilationfile, csv)
    snoop_vanilla(abspath(csv), abspath(path))
    data = SnoopCompile.read(csv)
    pc = SnoopCompile.parcel(reverse!(data[2]))
    delims = r"([\{\} \n\(\),])_([\{\} \n\(\),])"
    tmp_mod = eval(:(module $(gensym()) end))
    open(compilationfile, "w") do io
        println(io, "Sys.__init__()")
        println(io, "Base.early_init()")
        for (k, v) in pc
            k == :unknown && continue
            try
                eval(tmp_mod, :(import $k))
                println(io, "import $k")
                info("import $k")
            catch e
                warn("Module not found: $k")
            end
        end
        for (k, v) in pc
            for ln in v
                # replace `_` for free parameters, which print out a warning otherwise
                ln = replace(ln, delims, s"\1XXX\2")
                # only print out valid lines
                # TODO figure out the actual problems and why snoop compile emits invalid code
                try
                    parse(ln) # parse to make sure expression is parsing without error
                    # wrap in try catch to catch problematic code emitted by SnoopCompile
                    # without interupting the whole precompilation
                    # (usually, SnoopCompile emits 1% erroring statements)
                    println(io, "try\n    ", ln, "\nend")
                catch e
                    warn("Not emitted because code couldn't parse: ", ln)
                end
            end
        end
    end
end

function static_library_snoop()
    for (k, v) in pc
        for ln in v
            # replace `_` for free parameters, which print out a warning otherwise
            ln = replace(ln, delims, s"\1XXX\2")
            # only print out valid lines
            # TODO figure out the actual problems and why snoop compile emits invalid code
            try
                parse(ln) # parse to make sure expression is parsing without error
                # wrap in try catch to catch problematic code emitted by SnoopCompile
                # without interupting the whole precompilation
                # (usually, SnoopCompile emits 1% erroring statements)
                println(io, "try\n    ", ln, "\nend")
            catch e
                warn("Not emitted because code couldn't parse: ", ln)
            end
        end
    end
end

"""
    snoop_userimg(userimg, packages::Tuple{String, String}...)

    Traces all function calls in packages and writes out `precompile` statements into the file `userimg`
"""
function snoop_userimg(userimg, packages::Tuple{String, String}...)
    snooped_precompiles = map(packages) do package_snoopfile
        package, snoopfile = package_snoopfile
        abs_package_path = if ispath(package)
            normpath(abspath(package))
        else
            Pkg.dir(package)
        end
        file2snoop = normpath(abspath(joinpath(abs_package_path, snoopfile)))
        package = package_folder(get_root_dir(abs_package_path))
        isdir(package) || mkpath(package)
        precompile_file = joinpath(package, "precompile.jl")
        snoop(
            file2snoop,
            precompile_file,
            joinpath(package, "snooped.csv")
        )
        precompile_file
    end
    open(userimg, "w") do io
        for path in snooped_precompiles
            write(io, open(read, path))
            println(io)
        end
    end
    userimg
end


text:
using StableArrays
using LinearAlgebra

A = rand(10, 10)

# non stable

@show tr(A^500) # Inf in most cases

# stable

B = stabilize!(A)

@show tr(B^500) # Not Inf


text:
using Distributions, Plots; pyplot()

dist = TriangularDist(4,6,5)
N = 10^6
data = rand(dist,N)
yData=(data .- 5).^2

println("Mean: ", mean(yData), " Variance: ", var(data))

p1 = histogram(data, xlabel="x", bins=80, normed=true, ylims=(0,1.1))
p2 = histogram(yData, xlabel="y", bins=80, normed=true, ylims=(0,15))
plot(p1, p2, ylabel="Proportion", size=(800, 400), legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
function sum_vectors!(x, y, z)
n = length(x)
for i = 1:n
x[i] = y[i] + z[i]
end
end

max_stars_repo_path: Module 2/Chapter06/SIMD.jl

max_issues_repo_path: Module 2/Chapter06/SIMD.jl

max_forks_repo_path: Module 2/Chapter06/SIMD.jl


text:
using Documenter, ComputacionCientificaEnJulia

DocMeta.setdocmeta!(ComputacionCientificaEnJulia, :DocTestSetup,
                   :(using ComputacionCientificaEnJulia); recursive=true)

makedocs(
    format = Documenter.HTML(prettyurls=haskey(ENV, "GITHUB_ACTIONS"),
                             assets = ["assets/aligned.css"]),
    sitename = "Computación Científica en Julia",
    doctest = false,
    strict = false,
    pages = [
        "Generalidades" => "index.md",
        "Funciones" => Any["Definición" => "funciones/funciones.md",
                           "Ejemplos" => "funciones/ejemplos.md"]
            ]
)

deploydocs(
    repo = "github.com/mforets/computacion-cientifica-en-julia.git",
    push_preview=false,
)

max_stars_repo_name: mforets/computacion-cientifica-en-julia

max_issues_repo_name: mforets/computacion-cientifica-en-julia

max_forks_repo_name: mforets/computacion-cientifica-en-julia


text:
module StateSpaceFilters

using LinearAlgebra

abstract type AbstractStateSpaceModel end

struct PDMatChol{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
    dim::Int                    # matrix dimension
    chol::Cholesky{T,S}         # Cholesky factorization of mat
end

struct StateSpaceModel
    dynamics
    measurement
end

struct StationaryKalmanFilter
    model
    F
    G
    v
    H
end

struct InformationFilter end

struct SquareRootFilter end

struct FilterResults
    states_predicted
    states_filtered
    innovations
end

function run_filter(filter, y, x0)

    N = length(y)
    x_pred = Array{typeof(x0)}(undef, N)
    x_filtered = Array{typeof(x0)}(undef, N)
    e = Array{typeof(y)}(undef, N)

    x_pred[1] = x0
    for n = 1:N
        model_obs = get_observation_model(filter, n)
        x_filt[n], e[n] = measurement_update(model_obs, y[n], x_pred[n])

        if n != N
            model_dynamics = get_dynamics_model(filter, n)
            x_pred[n+1] = time_update(model_dynamics, x_filt[n])
        end
    end
    return  FilterResults(
        states_predicted,
        states_filtered,
        innovations
    )
end

end # module


text:
#=

This file specifies how to truncate the AerosolOptics struct, given the truncation type

=#

"""
    $(FUNCTIONNAME)(mod::δBGE, aero::AerosolOptics))

Returns the truncated aerosol optical properties as [`AerosolOptics`](@ref)
- `mod` a [`δBGE`](@ref) struct that defines the truncation order (new length of greek parameters) and exclusion angle
- `aero` a [`AerosolOptics`](@ref) set of aerosol optical properties that is to be truncated
"""
function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}; reportFit=false, err_β=nothing,err_ϵ=nothing,err_γ=nothing) where {FT}
    @unpack greek_coefs, ω̃, k = aero
    @unpack α, β, γ, δ, ϵ, ζ = greek_coefs
    @unpack l_max, Δ_angle =  mod

    # Obtain Gauss-Legendre quadrature points and weights for phase function:
    μ, w_μ = gausslegendre(length(β));

    # Reconstruct phase matrix elements:
    scattering_matrix, P, P² = reconstruct_phase(greek_coefs, μ; returnLeg=true)

    @unpack f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ = scattering_matrix

    # Find elements that exclude the peak (if wanted!)
    iμ = findall(x -> x < cosd(Δ_angle), μ)

    # Prefactor for P2:
    fac = zeros(FT,l_max);
    for l = 2:l_max - 1
        fac[l + 1] = sqrt(FT(1) / ( ( l - FT(1)) * l * (l + FT(1)) * (l + FT(2)) ));
    end

    # Create subsets (for Ax=y weighted least-squares fits):
    y₁₁ = view(f₁₁, iμ)
    y₁₂ = view(f₁₂, iμ)
    y₃₄ = view(f₃₄, iμ)
    A   = view(P, iμ, 1:l_max)
    B   = fac[3:end]' .* view(P², iμ, 3:l_max)

    # Weights (also avoid division by 0)
    minY = zeros(length(iμ)) .+ FT(1e-8);
    W₁₁ = Diagonal(w_μ[iμ] ./ max(abs.(y₁₁), minY));
    W₁₂ = Diagonal(w_μ[iμ] ./ max(abs.(y₁₂), minY));
    W₃₄ = Diagonal(w_μ[iμ] ./ max(abs.(y₃₄), minY));
    # W₁₂ = Diagonal(w_μ[iμ]);
    # W₃₄ = Diagonal(w_μ[iμ]);
    # Julia backslash operator for least squares (just like Matlab);
    cl = ((W₁₁ * A) \ (W₁₁ * y₁₁))   # B in δ-BGR (β)
    γᵗ = similar(cl); γᵗ[1:2] .=0
    ϵᵗ = similar(cl); ϵᵗ[1:2] .=0
    γᵗ[3:end] = ((W₁₂ * B) \ (W₁₂ * y₁₂))   # G in δ-BGE (γ)
    ϵᵗ[3:end] = ((W₃₄ * B) \ (W₃₄ * y₃₄))   # E in δ-BGE (ϵ)

    if reportFit
        println("Errors in δ-BGE fits:")
        mod_y = convert.(FT, A * cl)
        mod_γ = convert.(FT, B * γᵗ[3:end])
        mod_ϵ = convert.(FT, B * ϵᵗ[3:end])
        # push!(err_β, StatsBase.rmsd(W₁₁ * mod_y, W₁₁ * y₁₁; normalize=true))
        # push!(err_γ, StatsBase.rmsd(W₁₂ * mod_γ, W₁₂ * y₁₂; normalize=true))
        # push!(err_ϵ, StatsBase.rmsd(W₃₄ * mod_ϵ, W₃₄ * y₃₄; normalize=true))
        @show StatsBase.rmsd(mod_y, y₁₁; normalize=true)
        @show StatsBase.rmsd(mod_γ, y₁₂; normalize=true)
        @show StatsBase.rmsd(mod_ϵ, y₃₄; normalize=true)
    end

    # Integrate truncated function for later renormalization (here: fraction that IS still scattered):
    c₀ = FT(cl[1]) # ( w_μ' * (P[:,1:l_max] * cl) ) / 2
    # @show c₀, cl[1]
    # Compute truncated greek coefficients:
    βᵗ = cl / c₀                                    # Eq. 38a, B in δ-BGR (β)
    δᵗ = (δ[1:l_max] .- (β[1:l_max] .- cl)) / c₀    # Eq. 38b, derived from β
    αᵗ = (α[1:l_max] .- (β[1:l_max] .- cl)) / c₀    # Eq. 38c, derived from β
    ζᵗ = (ζ[1:l_max] .- (β[1:l_max] .- cl)) / c₀    # Eq. 38d, derived from β

    # Adjust scattering and extinction cross section!
    greek_coefs = GreekCoefs(αᵗ, βᵗ, γᵗ, δᵗ, ϵᵗ, ζᵗ  )

    # C_sca  = (ω̃ * k);
    # C_scaᵗ = C_sca * c₀;
    # C_ext  = k - (C_sca - C_scaᵗ);
    @show typeof(ω̃), typeof(k),typeof(c₀)
    # return AerosolOptics(greek_coefs = greek_coefs, ω̃=C_scaᵗ / C_ext, k=C_ext, fᵗ = 1-c₀)
    return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=k, fᵗ=(FT(1) - c₀))
end


text:
"""
    AOGCM1D(ndays=60)

Atmosphere-Ocean, coupled, one-dimensional column model.

```
outputs,parameters=AOGCM1D(10)
p = dirname(pathof(AirSeaFluxes))
include(joinpath(p,"recipes_plots.jl"))
p1,p2,p3=plot_final(outputs,parameters)
display(p3)
```
"""
function AOGCM1D(ndays=60)

cpl=1;
moist=1;

dt=60; #seconds
#ndays=60; # duration of simulation in days

n=Int(round(ndays*24*60*60/dt)); # number of time steps

H=50; # atmospheric layer thickness
dz=1; # oceanic layer thickness
vlevo=50; # number of oceanic vertical levels
vleva=50; # number of atmospheric vertical levels
rhow0=1029; # kg/m^3
cp=4184; #J/K/kg
cpa=1005; #J/K/kg
sb = 5.670e-8;
d2k=273.15;
Av=2500000;
gamma_blk=0.01;
cvapor_fac=640380;
cvapor_exp=5107.4;
cvapor_exp_ice=5897.8;
saltsat=0.980;
p0=101300; #reference pressure
gravity_mks=9.81;
humid_fac=0.606;
tref=25; #reference temperature
sref=25;
SO=sref; #salinity not implemented
alpha=2e-4; #expansion coefficients for temperature
beta=7.4e-4; #expansion coefficients for salt
emissivity=0.97; # ocean emissivity
kappa=1e-5; # m^2/kg air absorption coefficient
Rgas=287.05;
karman=0.4;

ZA=H*[1:vleva;];
ZO=zeros(vlevo+1);
ZO[2:vlevo+1]=-[1:vlevo;]*dz;
kT0=ones(vlevo)*1e-4; # SST vertical mixing
kT=kT0; # SST vertical mixing

Aimp=0.5; # Atmospheric implicitness
Oimp=0.5; # Oceanic implicitness

R=0.62; g1=0.6; g2=20;

f=R*exp.(ZO/g1)+(1-R)*exp.(ZO/g2);
f[ZO.<-200].=0;
f=-diff(f,dims=1);

albedo=0.06;
lat=30;
lon=0;
t=[1:n;]*dt/3600;
jstart=70;# julian day, borial summer
jadd=round.((t.-t[1])/24);
jd=jstart.+jadd;

AOI=angle_of_incidence.(lat,lon,jd,t); AOI[AOI.<0].=0;
SW=(1362*(1-albedo).*f'.*AOI);
#mean(SW(:,1))
#SW=(repmat((sin(2*pi/n*(1:n)*ndays+pi/2)+1)/2*1362*(1-albedo),vlev,1).*
#repmat(f',1,n))'.*cos(lat/180*pi);

TO=zeros(n+1,vlevo);
UO=zeros(n+1,vlevo);
SST=zeros(n+1);

UA=zeros(n+1,vleva);
qA=zeros(n+1,vleva);
THETA=zeros(n+1,vleva);
rhoa=zeros(n+1,vleva);
TA=zeros(n+1,vleva)

UOst=zeros(vlevo); TOst=zeros(vlevo);
UAst=zeros(vleva); qAst=zeros(vleva); THETAst=zeros(vleva);

############################################################################
# Initial conditions
############################################################################
#U0=ones(n,1)*5;
#THETA0=(sin(2*pi/n*(1:n)*ndays/5-pi/2)+21);
TA[1,:]=25 .-0.007.*H.*[1:vleva;];
P=zeros(vleva);
P[1]=101300;
P[2:vleva]=P[1].*exp.(-cumsum(H./(Rgas.*(TA[1,2:end].+d2k)./gravity_mks)));
rhoa[1,:]=P./(Rgas.*(TA[1,:].+d2k));
THETA[1,:]=(TA[1,:].+d2k).*(P[1]./P).^(2/7).-d2k;
THETA[1,1:20].=THETA[1,1];
THETA0=ones(n,1)*THETA[1,end].+1;#+gamma_blk*H+0.01;#+0.51;#+SW*kappa/cpa*dt;
T0=ones(n+1,vlevo).*26;#repmat((16+(vlev-1:-1:0)/vlev*8),n+1,1);
TO[:,:].=T0;
TO[:,1]=T0[:,1]; TO[1,:]=T0[1,1].+15 .*([vlevo:-1:1;]/vlevo.-1);
SST.=TO[:,1];
ssq = saltsat*cvapor_fac*exp(-cvapor_exp./(TA[1]+d2k))./rhoa[1,1];
RH=0.8;
q0=0.5*ssq;
#Tv=(THETA(i-1,:)+d2k).*(1+humid_fac*q(i-1));
qA[1,:]=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA[1,:].+d2k))./rhoa[1,:].*RH;
UA[1,:]=(0.2/karman).*(log.((ZA.-H/2)./1e-4));
#U(1,:)=U(1,end);
U0=(0.2/karman).*(log((ZA[end]+H/2)/1e-4)).*ones(1,n+1);
atemp=THETA[1,1]+d2k; aqh=qA[1,1];
speed=abs(UA[1,1]); sst=SST[1];
tmp0=bulkformulae(atemp,aqh,speed,sst,H/2,H/2,H/2,H/2,rhoa[1,1]);
huol0=tmp0["huol"]
TAup=270;
#U0=(sin(2*pi/n*(1:n)*ndays/5+2*pi/n*3600/dt*sf)*8+U(1,end));
#UO(:)=0.;
#U0(3600*24*5/dt:end)=U(1,end); one pertubation
############################################################################

LH=zeros(1,n); SH=zeros(1,n); TAU=zeros(1,n); TAUO=zeros(1,n);
LW=zeros(1,n); E=zeros(1,n); ssq=zeros(1,n);
xi=zeros(1,n); ce=zeros(1,n); ch=zeros(1,n);
ustar=zeros(1,n); qstar=zeros(1,n); tstar=zeros(1,n);
psimh=zeros(1,n); psixh=zeros(1,n); Ri=zeros(n,vleva);
rd=zeros(1,n); rh=zeros(1,n); re=zeros(1,n);
Qnet=zeros(1,n);
KAm=zeros(n,vleva);KAt=zeros(n,vleva);
KOm=zeros(n,vlevo);KOt=zeros(n,vlevo);

for i=2:n+1
    if rem(i,1440*60/dt)==0
        println(string("Day  ",round(i/(1440*60/dt),digits=5)))
    end

    ########################################################################
    # surface
    ########################################################################
    tmp=bulkformulae(TA[i-1,1]+d2k,qA[i-1,1],abs(UA[i-1,1]-UO[i-1,1]),SST[i-1,1],
        H/2,H/2,H/2,10,rhoa[i-1,1]);

    hl=tmp["hl"]
    hs=tmp["hs"]
    ce[i-1]=tmp["ce"]
    ch[i-1]=tmp["ch"]
    ssq[i-1]=tmp["ssq"]
    xi[i-1]=tmp["huol"]
    rd[i-1]=tmp["rd"]
    rh[i-1]=tmp["rh"]
    re[i-1]=tmp["re"]
    ustar[i-1]=tmp["ustar"]
    qstar[i-1]=tmp["qstar"]
    tstar[i-1]=tmp["tstar"]
    psimh[i-1]=tmp["psimh"]
    psixh[i-1]=tmp["psixh"]

    LH[i-1]=-hl;
    SH[i-1]=-hs;
    E[i-1]=-hl/Av;
    ϵa=0.725+0.17*log10(sum(qA[1,:]/1000*1.2*H)*100)
    LW[i-1]=emissivity*sb*((SST[i-1]+d2k).^4 .- ϵa.*(TAup).^4);
    TAU[i-1]=rhoa[i-1,1].*ustar[i-1].^2 .*sign(ustar[i-1]);
    Qnet[i-1]=SW[i-1,1]-LH[i-1]-SH[i-1]-LW[i-1];
    ########################################################################
    # Atmospheric vertical diffusion
    ########################################################################

    Tv=(THETA[i-1,:].+d2k).*(1 .+humid_fac*qA[i-1,:]);
    KAm[i-1,:], KAt[i-1,:]= holtslag(Tv,qA[i-1,1],UA[i-1,:],
                 UA[i-1,:].*0,LH[i-1],SH[i-1],ustar[i-1],[0; ZA]);
    kam=KAm[i-1,:];
    kat=KAt[i-1,:];

    ########################################################################
    # Oceanic vertical diffusion
    ########################################################################

    rhow=rhow0.*(1 .-alpha.*(TO[i-1,:].-tref).+beta.*(SO.-sref));
    B=gravity_mks.*(rhow0.-rhow)./rhow0;
    Uo=(vlevo:-1:1)/vlevo/10; Vo=zeros(size(Uo)); FWflux=zeros(size(Qnet));
    KOm[i-1,:],KOt[i-1,:],_=
        kpp(TO[i-1,:],TO[i-1,:].*0,UO[i,:],Vo,rhow,Qnet[i-1]
             ,0,ustar[i-1],ZO,lat);
    kot=KOt[i-1,:];
    kom=KOm[i-1,:];

    ########################################################################
    # Atmospheric explicit time stepping
    ########################################################################

    TAUS=sign(UA[i-1,1]-UO[i-1,1]);

    UAst[1]=UA[i-1,1]-TAU[i-1]*TAUS/(H*rhoa[i-1,1])*dt-(1-Aimp)*(
        kam[1].*(UA[i-1,1]-UA[i-1,2])/H^2*dt);
    UAst[end]=UA[i-1,end]-(1-Aimp)*(
        kam[end]*(UA[i-1,end]-U0[i-1])/H^2*dt);
    UAst[2:end-1]=UA[i-1,2:end-1]+(1-Aimp)*(
        kam[1:end-2].*(UA[i-1,1:end-2]-UA[i-1,2:end-1])/H^2*dt-
        kam[2:end-1].*(UA[i-1,2:end-1]-UA[i-1,3:end  ])/H^2*dt);

    qAst[1]=qA[i-1,1]-E[i-1]/(H*rhow[1])*dt-
        (1-Aimp)*(kat[1]*(qA[i-1,1]-qA[i-1,2])/H^2*dt);
    qAst[end]=qA[i-1,end]-(1-Aimp)*(kat[end]*(qA[i-1,end]-q0)/H^2*dt);
    qAst[2:end-1]=qA[i-1,2:end-1]+(1-Aimp)*(
        kat[1:end-2].*(qA[i-1,1:end-2]-qA[i-1,2:end-1])/H^2*dt-
        kat[2:end-1].*(qA[i-1,2:end-1]-qA[i-1,3:end  ])/H^2*dt);

    THETAst[1]=THETA[i-1,1]+(SH[i-1])/cpa/(H*rhoa[i-1,1])*dt-
         (1-Aimp)*(kat[1]*(THETA[i-1,1]-THETA[i-1,2])/H^2*dt);
    THETAst[end]=THETA[i-1,end]-
         (1-Aimp)*(kat[end]*(THETA[i-1,end]-THETA0[i-1])/H^2*dt);
    THETAst[2:end-1]=THETA[i-1,2:end-1]+(1-Aimp)*(
         kat[1:end-2].*(THETA[i-1,1:end-2]-THETA[i-1,2:end-1])/H^2*dt-
         kat[2:end-1].*(THETA[i-1,2:end-1]-THETA[i-1,3:end  ])/H^2*dt);

    ########################################################################
    # Oceanic explicit time stepping
    ########################################################################

    if cpl==1

        TOst[1]=TO[i-1,1]+(SW[i-1,1]-LH[i-1]-SH[i-1]-LW[i-1])/cp/(dz*rhow[1])*dt-
            (1-Oimp)*(kot[1]*(TO[i-1,1]-TO[i-1,2])/dz^2*dt);
        TOst[end]=TO[i-1,end]+SW[i-1,end]/cp/(dz*rhow[end])*dt+
            (1-Oimp)*(kot[end-1].*(TO[i-1,end-1]-TO[i-1,end])/dz^2*dt);
        TOst[2:end-1]=TO[i-1,2:end-1]+(1-Oimp)*(
            kot[1:end-2].*(TO[i-1,1:end-2]-TO[i-1,2:end-1])/dz^2*dt-
            kot[2:end-1].*(TO[i-1,2:end-1]-TO[i-1,3:end  ])/dz^2*dt)+
            SW[i-1,2:end-1]./cp./(dz*rhow[2:end-1])*dt;

        UOst[1]=UO[i-1,1]+
             (TAU[i-1])*TAUS/(rhow[1])/dz*dt-
             (1-Oimp)*(kom[1]*(UO[i-1,1]-UO[i-1,2])/dz^2*dt);
        UOst[end]=UO[i-1,end]+
             (1-Oimp)*(kom[end-1].*(UO[i-1,end-1]-UO[i-1,end])/dz^2*dt);
        UOst[2:end-1]=UO[i-1,2:end-1]+(1-Oimp)*(
             kom[1:end-2].*(UO[i-1,1:end-2]-UO[i-1,2:end-1])/dz^2*dt-
             kom[2:end-1].*(UO[i-1,2:end-1]-UO[i-1,3:end  ])/dz^2*dt);

    end

    if Aimp==0
        UA[i,:]=UAst; qA[i,:]=qAst; THETA[i,:]=THETAst;
    else
        du=-Aimp*kam[1:end-1].*dt./H^2; dl=-Aimp*kam[2:end].*dt./H^2;
        d=(1 .+Aimp.*dt./H^2 .*([kam[1:end-1]; 0]+[0; kam[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        UA[i,:]=A*UAst;
        du=-Aimp*kat[1:end-1].*dt./H^2; dl=-Aimp*kat[2:end].*dt./H^2;
        d=(1 .+Aimp.*dt./H^2 .*([kat[1:end-1]; 0]+[0; kat[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        qA[i,:]=A*qAst; THETA[i,:]=A*THETAst;
    end

    if Oimp==0
        UO[i,:]=UOst; TO[i,:]=TOst;
    else
        du=-Oimp*kom[1:end-1].*dt./dz^2; dl=-Oimp*kom[2:end].*dt./dz^2;
        d=(1 .+Oimp.*dt./dz^2 .*([kom[1:end-1]; 0]+[0; kom[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        UO[i,:]=A*UOst;
        du=-Oimp*kot[1:end-1].*dt./dz^2; dl=-Oimp*kot[2:end].*dt./dz^2;
        d=(1 .+Oimp.*dt./dz^2 .*([kot[1:end-1]; 0]+[0; kot[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        TO[i,:]=A*TOst;
    end
    SST[i]=TO[i,1]

    TA[i,:]=(THETA[i,:].+d2k).*(P./P[1]).^(2/7).-d2k;
    rhoa[i,:]=P./(Rgas.*(TA[i,:].+d2k));

    if moist==1
        qsat=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA[i,:].+d2k))./rhoa[i-1,:];
        dq=(qA[i,:].-qsat).*((qA[i,:].-qsat).>0);
        dq=min.(dq,qA[i,:]);
        dTHETA=Av.*dq./cp; #[J/kg][j-1 kg K]
        #println([i maximum(dTHETA) maximum(dq)])
        ind=qA[i,:].>qsat;
        qA[i,ind]=qA[i,ind]-dq[ind];
        #THETA[i,:]=THETA[i,: ]+dTHETA;
    end

end

outputs=Dict("U0" => U0,"UA" => UA,"THETA" => THETA,"SST" => SST,"TO" => TO,
"xi" => xi,"ustar" => ustar,"psimh" => psimh,"psixh" => psixh,
"rd" => rd,"rh" => rh,"re" => re,"SW" => SW,"LW" => LW,"SH" => SH,"LH" => LH)

parameters=Dict("dt"=>dt, "n"=>n, "ndays"=>ndays, "cp"=>cp, "dz"=>dz, "rhow0"=>rhow0)

return outputs,parameters

end


text:
#SigmoidLoss

Pkg.add("Vega")

Plots.plot(SigmoidLoss(), c=:olive, leg=:topright, linewidth=2, annotation=(1, 1.5, L"L\left( a \right) = 1 - tanh(a)", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
fast_hashed(v::AbstractVector) = isbitstype(eltype(v)) ? PooledArray(v) : v

function indices_iterator(cols)
    isempty(cols) && return Ref(Colon())
    grouping_sa = StructArray(map(refarray∘fast_hashed, cols))
    gp = GroupPerm(grouping_sa)
    return (sortperm(gp)[rg] for rg in gp)
end

splitapply(x::AbstractArray{<:Pair}) = x

splitapply(le::Entry) = splitapply(identity, le)

# TODO: make sure to also pool categorical positional mappings here
function splitapply(f, le::Entry)
    labels, mappings = map(getlabel, le.mappings), map(getvalue, le.mappings)
    axs = Broadcast.combine_axes(mappings.positional..., values(mappings.named)...)
    iter = (m for (k, m) in mappings.named if m isa AbstractVector && !iscontinuous(m))
    grouping_cols = Tuple(iter)
    list = Entry[]
    foreach(indices_iterator(grouping_cols)) do idxs
        for c in CartesianIndices(tail(axs))
            submappings = map(labels, mappings) do label, v
                I = ntuple(ndims(v)) do n
                    i = n == 1 ? idxs : c[n-1]
                    return adjust_index(axs[n], axes(v, n), i)
                end
                return Labeled(label, view(v, I...))
            end
            discrete, continuous = separate!(submappings)
            new_entries = maybewrap(f(Entry(le.plottype, continuous, le.attributes)))
            for new_entry in maybewrap(new_entries)
                push!(list, recombine!(discrete, new_entry))
            end
        end
    end
    return list
end


text:
# Automatically generated using Clang.jl
const IMPLOT_AUTO = Cint(-1)
const IMPLOT_AUTO_COL = ImVec4(0,0,0,-1)

const ImPlotContext = Ptr{Cvoid}

"Options for plots."
@cenum ImPlotFlags::Cint begin
    ImPlotFlags_None        = 0         # Default
    ImPlotFlags_NoTitle     = 1 << 0
    ImPlotFlags_NoLegend    = 1 << 1
    ImPlotFlags_NoMenus     = 1 << 2
    ImPlotFlags_NoBoxSelect = 1 << 3
    ImPlotFlags_NoMousePos  = 1 << 4
    ImPlotFlags_NoHighlight = 1 << 5
    ImPlotFlags_NoChild     = 1 << 6
    ImPlotFlags_Equal       = 1 << 7
    ImPlotFlags_YAxis2      = 1 << 8
    ImPlotFlags_YAxis3      = 1 << 9
    ImPlotFlags_Query       = 1 << 10
    ImPlotFlags_Crosshairs  = 1 << 11
    ImPlotFlags_AntiAliased = 1 << 12
    ImPlotFlags_CanvasOnly  = (1 << 0) | (1 << 1) | (1 << 2) | (1 << 3) | (1 << 4)
end

"Options for plot axes (X and Y)"
@cenum ImPlotAxisFlags::Cint begin
    ImPlotAxisFlags_None            = 0
    ImPlotAxisFlags_NoGridLines     = 1 << 0
    ImPlotAxisFlags_NoTickMarks     = 1 << 1
    ImPlotAxisFlags_NoTickLabels    = 1 << 2
    ImPlotAxisFlags_LogScale        = 1 << 3
    ImPlotAxisFlags_Time            = 1 << 4
    ImPlotAxisFlags_Invert          = 1 << 5
    ImPlotAxisFlags_LockMin         = 1 << 6
    ImPlotAxisFlags_LockMax         = 1 << 7
    ImPlotAxisFlags_Lock            = (1 << 6) | (1 << 7)
    ImPlotAxisFlags_NoDecorations   = (1 << 0) | (1 << 1) | (1 << 2)
end

"Plot styling colors"
@cenum ImPlotCol::Cint begin
    # Item styling colors
    ImPlotCol_Line              = 0
    ImPlotCol_Fill              = 1
    ImPlotCol_MarkerOutline     = 2
    ImPlotCol_MarkerFill        = 3
    ImPlotCol_ErrorBar          = 4
    # Plot styling colors
    ImPlotCol_FrameBg           = 5
    ImPlotCol_PlotBg            = 6
    ImPlotCol_PlotBorder        = 7
    ImPlotCol_LegendBg          = 8
    ImPlotCol_LegendBorder      = 9
    ImPlotCol_LegendText        = 10
    ImPlotCol_TitleText         = 11
    ImPlotCol_InlayText         = 12
    ImPlotCol_XAxis             = 13
    ImPlotCol_XAxisGrid         = 14
    ImPlotCol_YAxis             = 15
    ImPlotCol_YAxisGrid         = 16
    ImPlotCol_YAxis2            = 17
    ImPlotCol_YAxisGrid2        = 18
    ImPlotCol_YAxis3            = 19
    ImPlotCol_YAxisGrid3        = 20
    ImPlotCol_Selection         = 21
    ImPlotCol_Query             = 22
    ImPlotCol_Crosshairs        = 23
    ImPlotCol_COUNT             = 24
end

"Plot styling variables"
@cenum ImPlotStyleVar::Cint begin
    # Item styling variables
    ImPlotStyleVar_LineWeight           = 0   # Float32,  plot item line weight in pixels
    ImPlotStyleVar_Marker               = 1   # Int32,    marker specification
    ImPlotStyleVar_MarkerSize           = 2   # Float32,  marker size in pixels (roughly the marker's "radius")
    ImPlotStyleVar_MarkerWeight         = 3   # Float32,  plot outline weight of markers in pixels
    ImPlotStyleVar_FillAlpha            = 4   # Float32,  alpha modifier applied to all plot item fills
    ImPlotStyleVar_ErrorBarSize         = 5   # Float32,  error bar whisker width in pixels
    ImPlotStyleVar_ErrorBarWeight       = 6   # Float32,  error bar whisker weight in pixels
    ImPlotStyleVar_DigitalBitHeight     = 7   # Float32,  digital channels bit height (at 1) in pixels
    ImPlotStyleVar_DigitalBitGap        = 8   # Float32,  digital channels bit padding gap in pixels
    # plot styling variables
    ImPlotStyleVar_PlotBorderSize       = 9   # Float32,  thickness of border around plot area
    ImPlotStyleVar_MinorAlpha           = 10   # Float32,  alpha multiplier applied to minor axis grid lines
    ImPlotStyleVar_MajorTickLen         = 11   # ImVec2, major tick lengths for X and Y axes
    ImPlotStyleVar_MinorTickLen         = 12   # ImVec2, minor tick lengths for X and Y axes
    ImPlotStyleVar_MajorTickSize        = 13   # ImVec2, line thickness of major ticks
    ImPlotStyleVar_MinorTickSize        = 14   # ImVec2, line thickness of minor ticks
    ImPlotStyleVar_MajorGridSize        = 15   # ImVec2, line thickness of major grid lines
    ImPlotStyleVar_MinorGridSize        = 16   # ImVec2, line thickness of minor grid lines
    ImPlotStyleVar_PlotPadding          = 17   # ImVec2, padding between widget frame and plot area,
                                               # labels, or outside legends (i.e. main padding)
    ImPlotStyleVar_LabelPadding         = 18   # ImVec2, padding between axes labels, tick labels, and plot edge
    ImPlotStyleVar_LegendPadding        = 19   # ImVec2, legend padding from plot edges
    ImPlotStyleVar_LegendInnerPadding   = 20   # ImVec2, legend inner padding from legend edges
    ImPlotStyleVar_LegendSpacing        = 21   # ImVec2, spacing between legend entries
    ImPlotStyleVar_MousePosPadding      = 22   # ImVec2, padding between plot edge and interior info text
    ImPlotStyleVar_AnnotationPadding    = 23   # ImVec2, text padding around annotation labels
    ImPlotStyleVar_FitPadding           = 24       # ImVec2, additional fit padding as a percentage
                                               # of the fit extents (e.g. ImVec2(0.1f,0.1f) adds
                                               # 10% to the fit extents of X and Y)
    ImPlotStyleVar_PlotDefaultSize      = 25   # ImVec2, default size used when ImVec2(0,0) is passed to BeginPlot
    ImPlotStyleVar_PlotMinSize          = 26   # ImVec2, minimum size plot frame can be when shrunk
    ImPlotStyleVar_COUNT                = 27
end

"Marker specifications"
@cenum ImPlotMarker::Cint begin
    ImPlotMarker_None       = -1    # no marker
    ImPlotMarker_Circle     = 0     # a circle marker
    ImPlotMarker_Square     = 1     # a square maker
    ImPlotMarker_Diamond    = 2     # a diamond marker
    ImPlotMarker_Up         = 3     # an upward-pointing triangle marker
    ImPlotMarker_Down       = 4     # an downward-pointing triangle marker
    ImPlotMarker_Left       = 5     # an leftward-pointing triangle marker
    ImPlotMarker_Right      = 6     # an rightward-pointing triangle marker
    ImPlotMarker_Cross      = 7     # a cross marker (not fillable)
    ImPlotMarker_Plus       = 8     # a plus marker (not fillable)
    ImPlotMarker_Asterisk   = 9     # a asterisk marker (not fillable)
    ImPlotMarker_COUNT      = 10
end

"Built-in colormaps"
@cenum ImPlotColormap::Cint begin
    ImPlotColormap_Default  =  0    # ImPlot default colormap         (n=10)
    ImPlotColormap_Deep     =  1    # a.k.a. seaborn deep             (n=10)
    ImPlotColormap_Dark     =  2    # a.k.a. matplotlib "Set1"        (n=9)
    ImPlotColormap_Pastel   =  3    # a.k.a. matplotlib "Pastel1"     (n=9)
    ImPlotColormap_Paired   =  4    # a.k.a. matplotlib "Paired"      (n=12)
    ImPlotColormap_Viridis  =  5    # a.k.a. matplotlib "viridis"     (n=11)
    ImPlotColormap_Plasma   =  6    # a.k.a. matplotlib "plasma"      (n=11)
    ImPlotColormap_Hot      =  7    # a.k.a. matplotlib/MATLAB "hot"  (n=11)
    ImPlotColormap_Cool     =  8    # a.k.a. matplotlib/MATLAB "cool" (n=11)
    ImPlotColormap_Pink     =  9    # a.k.a. matplotlib/MATLAB "pink" (n=11)
    ImPlotColormap_Jet      =  10   # a.k.a. MATLAB "jet"             (n=11)
    ImPlotColormap_COUNT    =  11
end

"Used to position items on a plot (e.g. legends, labels, etc.)"
@cenum ImPlotLocation::Cint begin
    ImPlotLocation_Center    = 0                   # center-center
    ImPlotLocation_North     = 1 << 0              # top-center
    ImPlotLocation_South     = 1 << 1              # bottom-center
    ImPlotLocation_West      = 1 << 2              # center-left
    ImPlotLocation_East      = 1 << 3              # center-right
    ImPlotLocation_NorthWest = (1 << 0) | (1 << 2) # top-left
    ImPlotLocation_NorthEast = (1 << 0) | (1 << 3) # top-right
    ImPlotLocation_SouthWest = (1 << 1) | (1 << 2) # bottom-left
    ImPlotLocation_SouthEast = (1 << 1) | (1 << 3) # bottom-right
end

@cenum ImPlotOrientation::Cint begin
    ImPlotOrientation_Horizontal = 0
    ImPlotOrientation_Vertical   = 1
end

@cenum ImPlotYAxis::Cint begin
    ImPlotYAxis_1   = 0 # left (default
    ImPlotYaxis_2   = 1 # first on right side
    ImPlotYAxis_3   = 2 # second on right side
end

struct ImPlotPoint
    x::Cdouble
    y::Cdouble
end

struct ImPlotRange
    Min::Cdouble
    Max::Cdouble
end

struct ImPlotLimits
    X::ImPlotRange
    Y::ImPlotRange
end

struct ImPlotStyle
    LineWeight::Cfloat
    Marker::ImPlotMarker
    MarkerSize::Cfloat
    MarkerWeight::Cfloat
    FillAlpha::Cfloat
    ErrorBarSize::Cfloat
    ErrorBarWeight::Cfloat
    DigitalBitHeight::Cfloat
    DigitalBitGap::Cfloat
    # Plot styling variables
    PlotBorderSize::Cfloat
    MinorAlpha::Cfloat
    MajorTickLen::ImVec2
    MinorTickLen::ImVec2
    MajorTickSize::ImVec2
    MinorTickSize::ImVec2
    MajorGridSize::ImVec2
    MinorGridSize::ImVec2
    PlotPadding::ImVec2
    LabelPadding::ImVec2
    LegendPadding::ImVec2
    LegendInnerPadding::ImVec2
    LegendSpacing::ImVec2
    MousePosPadding::ImVec2
    AnnotationPadding::ImVec2
    FitPadding::ImVec2
    PlotDefaultSize::ImVec2
    PlotMinSize::ImVec2
    Colors::NTuple{Int(ImPlotCol_COUNT), ImVec4}
    AntiAliasedLines::Bool
    UseLocalTime::Bool
    UseISO8601::Bool
    Use24HourClock::Bool
end

struct ImPlotInputMap
    PanButton::ImGuiMouseButton
    PanMod::ImGuiKeyModFlags
    FitButton::ImGuiMouseButton
    ContextMenuButton::ImGuiMouseButton
    BoxSelectButton::ImGuiMouseButton
    BoxSelectMod::ImGuiKeyModFlags
    BoxSelectCancelButton::ImGuiMouseButton
    QueryButton::ImGuiMouseButton
    QueryMod::ImGuiKeyModFlags
    QueryToggleMod::ImGuiKeyModFlags
    HorizontalMod::ImGuiKeyModFlags
    VerticalMod::ImGuiKeyModFlags
end


text:
import Dolo

model = Dolo.yaml_import("rbc_iid.yaml")

# define
include("tmp_module.jl")
import temp

# get grid for endogenous
gg = model.options.grid
gg.a
gg.orders             # this is imported from the yaml file
grid = CartesianGrid(gg.a, gg.b, gg.orders) # temporary compatibility

# now look at exogenous
# everything below should work for any type of process (normal, markov chain, ar1, ...)

ptype = "MC"
if ptype == "IID"
    process = MvNormal(0.01)
elseif ptype == "MC"
    values = [ -0.01  0.01 ]'
    transitions = [
        0.9 0.1;
        0.1 0.9
    ]
    process = DiscreteMarkovProcess(transitions, values)
end

dprocess = discretize(process)

model.calibration.flat[:beta]

model.calibration[:parameters]
#beta=model.calibration.flat[:beta]
β = 0.5

drv_answer = value_iteration(model, dprocess, process, grid, β, true, 300)

println(drv_answer)
using Gadfly
kvec = linspace(8,12,100)

drv_answer.grid
yvec = [drv_answer(1,[k])[1] for k in kvec]
plot(x=kvec, y=yvec, Geom.line)


text:
module HyperDualMatrixTools

using HyperDualNumbers, LinearAlgebra, SparseArrays, SuiteSparse

# Personal add-on: Allow for multiplication by ε of a sparse matrix
import Base.*
"""
    *(d::Hyper, M::SparseMatrixCSC)

Exactly multiplies a sparse matrix `M` by the hyperdual number `d`.
This overload allows, e.g., `ε₁ * M` to be non-zero.
(This is to correct the effect of, e.g., the `ε₁ == 0` feature in HyperDualNumbers.jl,
which silently removes non-zeros from sparse matrices when multiplied by `ε₁`.)
"""
function *(h::Hyper, M::SparseMatrixCSC)
    i, j, v = findnz(M)
    m, n = size(M)
    return sparse(i, j, h .* v, m, n)
end
export *

"""
    HyperDualFactors

Container type to work efficiently with backslash on hyperdual-valued sparse matrices.

`factorize(M)` will create an instance containing
- `Af = factorize(realpart.(M))` — the factors of the real part
- `B = ε₁part.(M)` — the ``\\varepsilon_1`` part
- `C = ε₂part.(M)` — the ``\\varepsilon_2`` part
- `D = ε₁ε₂part.(M)` — the ``\\varepsilon_1\\varepsilon_2`` part
for a hyperdual-valued matrix `M`.

This is because only the factors of the real part are needed when solving a linear system of the type ``M x = b`` for a hyperdual-valued matrix ``M = A + \\varepsilon_1 B + \\varepsilon_2 C + \\varepsilon_1 \\varepsilon_2 D``.
In fact, the inverse of ``M`` is given by
``M^{-1} = (I - \\varepsilon_1 A^{-1} B - \\varepsilon_2 A^{-1} C - \\varepsilon_1\\varepsilon_2 A^{-1} (D - B A^{-1} C - C A^{-1} B)) A^{-1}``.
"""
mutable struct HyperDualFactors{TAf,T}
    Af::TAf # the factors of the real part
    B ::T   # the ε₁ part
    C ::T   # the ε₂ part
    D ::T   # the ε₁ε₂ part
end
export HyperDualFactors

# Factorization functions
for f in (:lu, :qr, :cholesky, :factorize)
    @eval begin
        import LinearAlgebra: $f
        """
        $($f)(M::SparseMatrixCSC{<:Hyper,<:Int})

        Invokes `$($f)` on just the real part of `M` and stores it along with the dual parts into a `HyperDualFactors` object.
        """
        $f(M::SparseMatrixCSC{<:Hyper,<:Int}) = HyperDualFactors($f(realpart.(M)), ε₁part.(M), ε₂part.(M), ε₁ε₂part.(M))
        """
        $($f)(M::Array{<:Hyper,2})

        Invokes `$($f)` on just the real part of `M` and stores it along with the dual parts into a `HyperDualFactors` object.
        """
        $f(M::Array{<:Hyper,2}) = HyperDualFactors($f(realpart.(M)), ε₁part.(M), ε₂part.(M), ε₁ε₂part.(M))
        export $f
    end
end

# In-place factorization for the case where the real part is already stored
import LinearAlgebra: factorize
function factorize(Mf::HyperDualFactors, M; update_factors = false)
    Mf.B = ε₁part.(M)
    Mf.C = ε₂part.(M)
    Mf.D = ε₁ε₂part.(M)
    if update_factors
        Mf.Af = factorize(realpart.(M))
    end
    return Mf
end
export factorize

# Adjoint and transpose definitions for `HyperDualFactors`
for f in (:adjoint, :transpose)
    @eval begin
        import Base: $f
        """
        $($f)(M::HyperDualFactors)

        Invokes `$($f)` on `M.Af`, `M.B`, `M.C`, and `M.D` and returns them into a new `HyperDualFactors` object.
        """
        $f(M::HyperDualFactors) = HyperDualFactors($f(M.Af), $f(M.B), $f(M.C), $f(M.D))
        export $f
    end
end

import Base.\
"""
    \\(M::HyperDualFactors, y::AbstractVecOrMat{Float64})

Backsubstitution for `HyperDualFactors`.
See `HyperDualFactors` for details.
"""
function \(M::HyperDualFactors, a::AbstractVecOrMat{Float64})
    A, B, C, D = M.Af, M.B, M.C, M.D
    A⁻¹a = A \ a
    _A⁻¹BA⁻¹a = A \ (-B * A⁻¹a)
    _A⁻¹CA⁻¹a = A \ (-C * A⁻¹a)
    return A⁻¹a + ε₁ * _A⁻¹BA⁻¹a + ε₂ * _A⁻¹CA⁻¹a +
        ε₁ε₂ * (A \ (-D * A⁻¹a - C * _A⁻¹BA⁻¹a - B * _A⁻¹CA⁻¹a))
end

"""
    \\(M::HyperDualFactors, y::AbstractVecOrMat{Hyper256})

Backsubstitution for `HyperDualFactors`.
See `HyperDualFactors` for details.
"""
function \(M::HyperDualFactors, y::AbstractVecOrMat{Hyper256})
    a, b, c, d = realpart.(y), ε₁part.(y), ε₂part.(y), ε₁ε₂part.(y)
    A, B, C, D = M.Af, M.B, M.C, M.D
    A⁻¹a = A \ a
    A⁻¹b_A⁻¹BA⁻¹a = A \ (b - B * A⁻¹a)
    A⁻¹c_A⁻¹CA⁻¹a = A \ (c - C * A⁻¹a)
    return A⁻¹a + ε₁ * A⁻¹b_A⁻¹BA⁻¹a + ε₂ * A⁻¹c_A⁻¹CA⁻¹a +
        ε₁ε₂ * (A \ (d - D * A⁻¹a - C * A⁻¹b_A⁻¹BA⁻¹a - B * A⁻¹c_A⁻¹CA⁻¹a))
end

"""
    \\(Af::Factorization{Float64}, y::AbstractVecOrMat{Hyper256})

Backsubstitution for HyperDual-valued RHS.
"""
function \(Af::Factorization{Float64}, y::AbstractVecOrMat{Hyper256})
    return (Af \ realpart.(y)) + ε₁ * (Af \ ε₁part.(y)) + ε₂ * (Af \ ε₂part.(y)) + ε₁ε₂ * (Af \ ε₁ε₂part.(y))
end

"""
    \\(M::AbstractArray{<:Hyper,2}, y::AbstractVecOrMat)

Backslash (factorization and backsubstitution) for Dual-valued matrix `M`.
"""
\(M::SparseMatrixCSC{<:Hyper,<:Int}, y::AbstractVecOrMat) = factorize(M) \ y
\(M::Array{<:Hyper,2}, y::AbstractVecOrMat) = factorize(M) \ y
export \

import Base.isapprox
function isapprox(x::AbstractVecOrMat{Hyper256}, y::AbstractVecOrMat{Hyper256})
    bigx = [realpart.(x) ε₁part.(x) ε₂part.(x) ε₁ε₂part.(x)]
    bigy = [realpart.(y) ε₁part.(y) ε₂part.(y) ε₁ε₂part.(y)]
    return isapprox(bigx, bigy)
end
isapprox(x::AbstractVecOrMat, y::AbstractVecOrMat{Hyper256}) = isapprox(hyper.(x), y)
isapprox(x::AbstractVecOrMat{Hyper256}, y::AbstractVecOrMat) = isapprox(x, hyper.(y))
function isapprox(x::Hyper256, y::Hyper256)
    bigx = [realpart(x) ε₁part(x) ε₂part(x) ε₁ε₂part(x)]
    bigy = [realpart(y) ε₁part(y) ε₂part(y) ε₁ε₂part(y)]
    return isapprox(bigx, bigy)
end
isapprox(x::Float64, y::Hyper256) = isapprox(hyper(x), y)
isapprox(x::Hyper256, y::Float64) = isapprox(x, hyper(y))

end # module


text:
using Documenter
using AtiyahBott

makedocs(
    sitename = "AtiyahBott",
    format = Documenter.HTML(),
    modules = [AtiyahBott]
)

deploydocs(
    repo = "github.com/mgemath/AtiyahBott.jl.git",
    versions = nothing,
    #=,
    target = "build",
    push_preview = true,=#
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays, ArrayLayouts, SparseArrays
import LazyArrays: CachedArray, colsupport, rowsupport, LazyArrayStyle, broadcasted,
            PaddedLayout, ApplyLayout, BroadcastLayout, AddArray, LazyLayout
import ArrayLayouts: OnesLayout

@testset "Lazy MemoryLayout" begin
    @testset "ApplyArray" begin
        A = [1.0 2; 3 4]
        @test eltype(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))) == Float64
        @test @inferred(MemoryLayout(typeof(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))))) ==
            ApplyLayout{typeof(+)}()
    end

    @testset "BroadcastArray" begin
        A = [1.0 2; 3 4]

        @test @inferred(MemoryLayout(typeof(BroadcastArray(+, A, Fill(0, (2, 2)), Zeros(2, 2))))) ==
            BroadcastLayout{typeof(+)}()

        @test MemoryLayout(typeof(Diagonal(BroadcastArray(exp,randn(5))))) == DiagonalLayout{LazyLayout}()
    end

    @testset "Vcat" begin
        @test @inferred(MemoryLayout(typeof(Vcat(Ones(10),Zeros(10))))) == PaddedLayout{OnesLayout}()
        @test @inferred(MemoryLayout(typeof(Vcat([1.],Zeros(10))))) == PaddedLayout{DenseColumnMajor}()
    end

    @testset "adjtrans/symherm" begin
        A = ApplyArray(+, randn(2,2), randn(2,2))
        B = ApplyArray(+, randn(2,2), im*randn(2,2))
        @test MemoryLayout(Symmetric(A)) isa SymmetricLayout{LazyLayout}
        @test MemoryLayout(Symmetric(B)) isa SymmetricLayout{LazyLayout}
        @test MemoryLayout(Hermitian(A)) isa SymmetricLayout{LazyLayout}
        @test MemoryLayout(Hermitian(B)) isa HermitianLayout{LazyLayout}
        @test MemoryLayout(A') isa LazyLayout
        @test MemoryLayout(B') isa LazyLayout
        @test MemoryLayout(transpose(A)) isa LazyLayout
        @test MemoryLayout(transpose(B)) isa LazyLayout
    end
end

include("applytests.jl")
include("multests.jl")
include("ldivtests.jl")
include("addtests.jl")
include("setoptests.jl")
include("macrotests.jl")
include("lazymultests.jl")
include("concattests.jl")
include("paddedtests.jl")
include("broadcasttests.jl")
include("cachetests.jl")

@testset "Kron"  begin
    A = [1,2,3]
    B = [4,5,6,7]

    @test Array(@inferred(Kron(A))) == A == copyto!(similar(A), Kron(A))
    K = @inferred(Kron(A,B))
    @test size(K) == (12,)
    @test size(K,1) == 12
    @test size(K,2) == 1
    @test axes(K) == (Base.OneTo(12),)
    @test [K[k] for k=1:length(K)] == Array(K) == kron(A,B)

    A = randn(3)
    K = @inferred(Kron(A,B))
    @test K isa Kron{Float64}
    @test all(K.args .=== (A,B))
    @test [K[k] for k=1:length(K)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B) == copyto!(similar(K), K)

    C = [7,8,9,10,11]
    K = Kron(A,B,C)
    @test [K[k] for k=1:length(K)] == Array(K) == kron(A,B,C) == copyto!(similar(K), K)

    @testset "Rectangular Factors: Compatibility with Matrix functions" begin
        A = randn(3,2)
        B = randn(4,6)
        K, k = Kron(A,B), kron(A,B)
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == k
        @test det(K) == 0  # kronecker of rectangular factors
        @test logdet(K) == -Inf
        @test logabsdet(K) == (-Inf, 0)
        @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
        @test tr(K) ≈ tr(k)

        @test K' == Kron(A', B') == k'
        @test transpose(K) == Kron(transpose(A), transpose(B)) == transpose(k)
        @test pinv(K) == Kron(pinv(A), pinv(B)) ≈ pinv(k)
        @test_throws SingularException inv(K)

        K, k = Kron(A,B'), kron(A,B')
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B')) == k
        @test_throws DimensionMismatch det(K)
        @test_throws DimensionMismatch tr(K)

        K, k = Kron(A',B), kron(A',B)
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B)) == k
        @test_throws DimensionMismatch det(K)
        @test_throws DimensionMismatch tr(K)

        K, k = Kron(A',B'), kron(A',B')
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B')) == k
        @test det(K) == 0  # kronecker of rectangular factors
        @test logdet(K) == -Inf
        @test logabsdet(K) == (-Inf, 0)
        @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
        @test tr(K) ≈ tr(k)
    end

    @testset "Square Factors: Compatibility with Matrix functions" begin
        # need A,B,C to be invertible
        A, B, C = 0, 0, 0
        while det(A) == 0
            A = randn(3,3)
        end
        while det(B) == 0
            B = randn(6,6)
        end
        while det(C) == 0
            C = randn(2,2)
        end
        K, k = Kron(A,B,C), kron(A,B,C)
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B,C)) == k

        @test K' == Kron(A', B', C') == k'
        @test transpose(K) == Kron(transpose(A), transpose(B), transpose(C)) == transpose(k)
        @test pinv(K) == Kron(pinv(A), pinv(B), pinv(C)) ≈ pinv(k)
        @test inv(K) == Kron(inv(A), inv(B), inv(C)) ≈ inv(k)

        @test logdet(K) ≈ logdet(k)
        @test all(logabsdet(K) .≈ logabsdet(k))
        @test det(K) ≈ det(k)
        @test diag(K) ≈ diag(k)
        @test tr(K) ≈ tr(k)
    end

    @testset "Zero Factor" begin
        A = zeros(3,3)
        B = randn(4,4)
        K, k = Kron(A, B), kron(A, B)
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == k == zeros(12, 12)
        @test det(K) ≈ det(k) ≈ 0
        @test isinf(logdet(K))
        @test logdet(K) ≈ logdet(k) == -Inf
        @test all(logabsdet(K) .≈ logabsdet(k) .≈ (-Inf, 0))

        c = randn(10)
        @test_throws DimensionMismatch K * c
        C = randn(10, 12)
        @test_throws DimensionMismatch K * C
        c = randn(12)
        @test (K * c) == zeros(12)
    end

    @testset "Identity Factor" begin
        A = I(3)
        B = randn(4,4)
        K, k = Kron(A, B), kron(A, B)
        @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == k
        @test det(K) ≈ det(k) ≈ det(B)^size(A,1)

        c = randn(12)
        @test (K * c) ≈ (k * c)
    end

    A = randn(3,2)
    B = randn(4,6)
    K = Kron(A,B)
    C = similar(K)
    copyto!(C, K)
    @test C == kron(A,B)
    @test_throws DimensionMismatch randn(2,2) .= K

    A = rand(Int,3,2)
    K = Kron(A,B)
    @test K isa Kron{Float64}
    @test all(K.args .=== (A,B))
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B)

    K = @inferred(Kron{Float64}(Eye{Float64}(1), zeros(4)))
    @test Array(K) == zeros(4,1)

    @testset "Applied bug" begin
        A = randn(5,5)
        K = Kron(A,A)
        k = applied(kron,A,A)
        @test K[1,1] == k[1,1] == A[1,1]^2
        x = randn(5)
        K = Kron(A,x)
        k = applied(kron,A,x)
        @test K[1,1] == k[1,1] == A[1,1]*x[1]
        K = Kron(x,x)
        k = applied(kron,x,x)
        @test K[1] == k[1] == x[1]^2
        K = Kron(A)
        k = applied(kron,A)
        @test K[1,1] == k[1,1] == A[1,1]
        K = Kron(x)
        k = applied(kron,x)
        @test K[1] == k[1] == x[1]
    end

    @testset "triple vector" begin
        x = randn(5)
        y = randn(6)
        z = randn(4)
        K = Kron(x,y,z)
        @test K[1] == K[1,1] == x[1]y[1]z[1]
        @test K == kron(x,y,z)
    end

    @testset "2-factor kron-mul" begin
        A, B = randn(4, 4), randn(3, 2)
        K, k = Kron(A, B), kron(A, B)
        x = randn(size(k, 2))
        X = randn(size(k, 2), 7)

        res_vec = K * x
        @test size(res_vec, 1) == size(K, 1) == size(k, 1)
        @test res_vec ≈ (k * x)

        res_mat = K * X
        @test size(res_mat, 1) == size(K, 1) == size(k, 1)
        @test res_mat ≈ (k * X)
    end

    @testset "3-factor kron-mul" begin
        A, B, C = randn(4, 4), randn(3, 2), randn(5, 6)
        K, k = Kron(A, B, C), kron(A, B, C)
        x = randn(size(k, 2))
        X = randn(size(k, 2), 8)

        res_vec = K * x
        @test size(res_vec, 1) == size(K, 1) == size(k, 1)
        @test res_vec ≈ (k * x)

        res_mat = K * X
        @test size(res_mat, 1) == size(K, 1) == size(k, 1)
        @test res_mat ≈ (k * X)
    end

    @testset "3-factor all-square kron-mul" begin
        A, B, C = randn(4, 4), randn(3, 3), randn(5, 5)
        K, k = Kron(A, B, C), kron(A, B, C)
        x = randn(size(k, 2))
        X = randn(size(k, 2), 8)

        res_vec = K * x
        @test size(res_vec, 1) == size(K, 1) == size(k, 1)
        @test res_vec ≈ (k * x)

        res_mat = K * X
        @test size(res_mat, 1) == size(K, 1) == size(k, 1)
        @test res_mat ≈ (k * X)
    end

    @testset "3-factor sparse kron-mul" begin
        A, B, C = sprandn(4, 4, 0.2), sprandn(3, 3, 0.2), sprandn(5, 5, 0.2)
        K, k = Kron(A, B, C), kron(A, B, C)
        x = randn(size(k, 2))
        X = randn(size(k, 2), 8)

        res_vec = K * x
        @test size(res_vec, 1) == size(K, 1) == size(k, 1)
        @test res_vec ≈ (k * x)

        res_mat = K * X
        @test size(res_mat, 1) == size(K, 1) == size(k, 1)
        @test res_mat ≈ (k * X)
    end

    @testset "3-factor mixed-sparsity kron-mul" begin
        A, B, C = sprandn(4, 4, 0.2), randn(3, 3), sprandn(5, 5, 0.2)
        K, k = Kron(A, B, C), kron(A, B, C)
        x = randn(size(k, 2))
        X = randn(size(k, 2), 8)

        res_vec = K * x
        @test size(res_vec, 1) == size(K, 1) == size(k, 1)
        @test res_vec ≈ (k * x)

        res_mat = K * X
        @test size(res_mat, 1) == size(K, 1) == size(k, 1)
        @test res_mat ≈ (k * X)
    end

    @testset "kron-by-kron mul" begin
        A, B = randn(3, 4), rand(4, 6)
        C, D = randn(4, 5), rand(6, 4)

        K1, K2 = Kron(A, B), Kron(C, D)
        res = K1 * K2
        @test res ≈ (kron(A, B) * kron(C, D))
        @test res isa Kron{eltype(A), 2}
    end

    @testset "kron-by-kron mul (unaligned)" begin
        A, B = randn(3, 8), rand(4, 3)
        C, D = randn(4, 5), rand(6, 4)

        K1, K2 = Kron(A, B), Kron(C, D)
        res = K1 * K2
        @test res ≈ (kron(A, B) * kron(C, D))
        @test res isa Matrix{eltype(A)}
    end

end

@testset "Diff and Cumsum" begin
    x = Vcat([3,4], [1,1,1,1,1], 3)
    y = @inferred(cumsum(x))
    @test y isa Vcat
    @test y == cumsum(Vector(x)) == Cumsum(x)
    @test @inferred(diff(x)) == diff(Vector(x)) == Diff(x)
    @test diff(x) isa Vcat
    @test Cumsum(x) == Cumsum(x)

    @test sum(x) == sum(Vector(x)) == last(y)
    @test cumsum(Vcat(4)) === Vcat(4)
    @test diff(Vcat(4)) === Vcat{Int}()

    A = randn(3,4)
    @test Cumsum(A; dims=1) == cumsum(A; dims=1)
    @test Cumsum(A; dims=2) == cumsum(A; dims=2)
    @test Diff(A; dims=1) == diff(A; dims=1)
    @test Diff(A; dims=2) == diff(A; dims=2)

    @test Cumsum(A; dims=1) == Cumsum(A; dims=1)

    a = Vcat([1,2,3], Fill(2,100_000_000))
    @test @inferred(cumsum(a))[end] == sum(a) == 200000006
    a = Vcat(2, Fill(2,100_000_000))
    @test @inferred(cumsum(a))[end] == sum(a) == 200000002

    @testset "empty" begin
        @test cumsum(Vcat(Int[], 1:5)) == cumsum(1:5)
        a = Vcat(Int[], [1,2,3])
        @test @inferred(cumsum(a)) == cumsum(Vector(a))
        a = Vcat(1, [1,2], [4,5,6])
        @test @inferred(cumsum(a)) == cumsum(Vector(a))
        a = Vcat(1, Int[], [4,5,6])
        @test cumsum(a) == cumsum(Vector(a))
        a = Vcat(1, Int[], Int[], [4,5,6])
        @test cumsum(a) == cumsum(Vector(a))
    end

    @testset "lazy cumsum" begin
        c = BroadcastArray(exp, 1:10)
        @test cumsum(c) == Cumsum(BroadcastArray(exp, 1:10))
        @test cumsum(BroadcastArray(exp, 1:10)) isa typeof(Cumsum(BroadcastArray(exp, 1:10)))
        @test cumsum(ApplyArray(+, 1:10)) == Cumsum(ApplyArray(+, 1:10))
        @test cumsum(ApplyArray(+, 1:10)) isa typeof(Cumsum(ApplyArray(+, 1:10)))

        @test copyto!(similar(c), cumsum(c)) == cumsum(Vector(c))
    end

    @testset "Cumprod" begin
        a = Accumulate(*, 1:5)
        @test IndexStyle(typeof(a)) == IndexLinear()
        @test a == cumprod(1:5)
        v = BroadcastArray(+, 1, BroadcastArray(^, 1:10_000_000, -2.0))
        a = accumulate(*, v)
        @test a isa Accumulate
        @test copy(a) isa Accumulate
        @test copy(a) == a
        @test a[end] ≈ prod(1 .+ (1:10_000_000).^(-2.0))
        @test LazyArrays.AccumulateAbstractVector(*, 1:5) == Accumulate(*, 1:5)
        @test LazyArrays.AccumulateAbstractVector(*, 1:5) isa LazyArrays.AccumulateAbstractVector
    end
end

@testset "col/rowsupport" begin
    A = randn(5,6)
    @test rowsupport(A,1) === Base.OneTo(6)
    @test colsupport(A,1) === Base.OneTo(5)
    D = Diagonal(randn(5))
    @test rowsupport(D,3) === colsupport(D,3) === 3
    Z = Zeros(5)
    @test rowsupport(Z,1) === colsupport(Z,1) === 1:0
    @test_broken cache(D)
    C = cache(Array,D);
    @test colsupport(C,2) === 2:2
    @test @inferred(colsupport(C,1)) === 1:1
    @test colsupport(cache(Zeros(5,5)),1) == 1:0
    C = cache(Zeros(5));
    @test colsupport(C,1) == 1:0
    C[3] = 1
    @test colsupport(C,1) == 1:3

    LazyArrays.zero!(C)
    @test colsupport(C,1) == 1:3
    @test C == zeros(5)

    @testset "convexunion" begin
        # bug from BandedMartrices.jl
        @test LazyArrays.convexunion(7:10,9:8) == LazyArrays.convexunion(9:8,7:10) == 7:10
        @test LazyArrays.convexunion(1,5) == LazyArrays.convexunion(1:3,5) == LazyArrays.convexunion(5,1:3) == 1:5
    end
end

@testset "triu/tril" begin
    A = ApplyArray(triu,randn(2,2))
    @test A isa ApplyArray{Float64}
    @test A[2,1] == 0
    @test A[1,1] == A.args[1][1,1]
    @test A == triu(A.args[1])
    A = ApplyArray(tril,randn(2,2))
    @test A isa ApplyArray{Float64}
    @test A[1,2] == 0
    @test A[1,1] == A.args[1][1,1]
    @test A == tril(A.args[1])
    A = ApplyArray(triu,randn(2,2),1)
    @test A[1,1] == 0
    @test A[1,2] == A.args[1][1,2]
    @test A isa ApplyArray{Float64}
    @test A == triu(A.args[1],1)
    A = ApplyArray(tril,randn(2,2),-1)
    @test A isa ApplyArray{Float64}
    @test A[1,1] == 0
    @test A[2,1] == A.args[1][2,1]
    @test A == tril(A.args[1],-1)

    A = ApplyMatrix(exp,randn(2,2))
    @test triu(A) isa ApplyMatrix{Float64,typeof(triu)}
    @test tril(A) isa ApplyMatrix{Float64,typeof(tril)}
    @test triu(A,1) isa ApplyMatrix{Float64,typeof(triu)}
    @test tril(A,1) isa ApplyMatrix{Float64,typeof(tril)}
end

@testset "BroadcastArray" begin
    bc = broadcasted(exp,[1,2,3])
    v = BroadcastArray(exp, [1,2,3])
    @test BroadcastArray(bc) == BroadcastVector(bc) == BroadcastVector{Float64,typeof(exp),typeof(bc.args)}(bc) ==
        v == BroadcastVector(exp, [1,2,3]) == exp.([1,2,3])

    Base.IndexStyle(typeof(BroadcastVector(exp, [1,2,3]))) == IndexLinear()

    bc = broadcasted(exp,[1 2; 3 4])
    M = BroadcastArray(exp, [1 2; 3 4])
    @test BroadcastArray(bc) == BroadcastMatrix(bc) == BroadcastMatrix{Float64,typeof(exp),typeof(bc.args)}(bc) ==
        M == BroadcastMatrix(BroadcastMatrix(bc)) == BroadcastMatrix(exp,[1 2; 3 4]) == exp.([1 2; 3 4])

    @test exp.(v') isa BroadcastMatrix
    @test exp.(transpose(v)) isa BroadcastMatrix
    @test exp.(M') isa BroadcastMatrix
    @test exp.(transpose(M)) isa BroadcastMatrix

    bc = BroadcastArray(broadcasted(+, 1:10, broadcasted(sin, 1:10)))
    @test bc[1:10] == (1:10) .+ sin.(1:10)

    bc = BroadcastArray(broadcasted(+,1:10,broadcasted(+,1,2)))
    @test bc.args[2] == 3
end

@testset "padded columns" begin
    A = randn(5,5)
    U = UpperTriangular(A)
    v = view(U,:,3)
    @test MemoryLayout(v) isa PaddedLayout{DenseColumnMajor}
    @test layout_getindex(v,1:4) == U[1:4,3]
    @test layout_getindex(v,1:4) isa Vcat

    L = LowerTriangular(A)
    w = view(L,3,:)
    @test MemoryLayout(w) isa PaddedLayout{ArrayLayouts.StridedLayout}
    @test layout_getindex(w,1:4) == L[3,1:4]
    @test layout_getindex(w,1:4) isa Vcat
end


text:
module Stratum

import JSON, UUIDs

export StratumEndpoint, send_notification, send_request, send_success_response, send_error_response

include("core.jl")
include("proxy.jl")

end


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

JLLWrappers.@generate_main_file_header("LibModbus")
JLLWrappers.@generate_main_file("LibModbus", UUID("f28bf708-275a-52ef-9ea6-4821c17573f9"))
end  # module LibModbus_jll


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigInt()
    l = Exo1SplitBinion(Exo1Real(a, b))
    r = Exo1SplitBinion(a, 0, b, 0)
    l == r
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo1SplitBinion(a))
end

@test begin
    a = 1
    !isreal(Exo1SplitBinion(a, a, a, a))
end

@test begin
    z = random(Exo1SplitBinion{Int})
    z == +(z)
end

@test begin
    a = 1
    real(Exo1SplitBinion(a)) == a
end

@test begin
    a = 1
    b = 2
    c = 3
    l = unreal(Exo1SplitBinion(0, a, b, c))
    r = [a, b, c]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1SplitBinion(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, s: 2, a: 3, sa: 4]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo1SplitBinion(SplitBinion(a, b))
    r = Exo1SplitBinion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo1SplitBinion(a, b, c)
    r = Exo1SplitBinion(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigFloat})
    y = random(Exo1SplitBinion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1SplitBinion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1SplitBinion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = associator(x, y, z)
    r = zero(Exo1SplitBinion{BigInt})
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo1SplitBinion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{Rational{BigInt}})
    y = random(Exo1SplitBinion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{Rational{BigInt}})
    y = random(Exo1SplitBinion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo1SplitBinion{Rational{BigInt}})
    y = random(Exo1SplitBinion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1SplitBinion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(SplitBinion{BigInt})
    x = random(Exo1SplitBinion{BigInt})
    l = a * x
    r = x * a
    l != r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    z = random(Exo1SplitBinion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) ||quadrance(x) == big(0)
end

@test begin
    x = random(Exo1SplitBinion{BigInt})
    y = random(Exo1SplitBinion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo1SplitBinion(0, 0.0))
end

@test_throws ErrorException begin
    inv(Exo1SplitBinion(0, 0, 1, 1))
end

@test_throws ErrorException begin
    inv(Exo1SplitBinion(1, 1, 2, 1))
end

@test_throws ErrorException begin
    random(Exo1SplitBinion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo1SplitBinion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1SplitBinion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1SplitBinion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo1SplitBinion{BigInt})
    z == z + zero(Exo1SplitBinion{BigInt})
end

@test begin
    z = random(Exo1SplitBinion{BigInt})
    z == z * one(Exo1SplitBinion{BigInt})
end

@test begin
    z = random(Exo1SplitBinion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo1SplitBinion{BigInt})
    z == z * one(z)
end


text:
using QuantumStatistics

function integrand(config)
	X = config.var[1]
	Y = config.var[2]
	if config.curr == 1
		x = X[1]
		y = Y[1]
		return exp(-(x^2 + y^2))
	elseif config.curr == 2
		x1, x2 = X[1], X[2]
		y1, y2 = Y[1], Y[2]
		return exp(-(x1^2 + x2^2 + y1^2 + y2^2))
	else
		throw("not possible!")
	end
end

function measure(config)
    factor = 1.0 / config.reweight[config.curr]
    weight = integrand(config)
    config.observable[config.curr] += weight / abs(weight) * factor
end

X = MonteCarlo.Tau(1.0, 0.2)
Y = MonteCarlo.Tau(2.0, 0.2)
dof = [[1, 1], [2, 2]] # 1 integral, two variables
config = MonteCarlo.Configuration(1e6, (X, Y), dof, [0.0, 0.0])
avg, err = MonteCarlo.sample(config, integrand, measure; Nblock=32)

if isnothing(avg) == false
	println("$(avg[1]) +- $(err[1])")
	println("$(avg[2]) +- $(err[2])")
end


text:
function solve_last_corner!(cost::CostParameters, mp::RLSModel, τ)
    c = cost.Cmin
    P1 = 0.
    vN1 = (mp.β*Φ(0, -K(c, cost)))/(1 - mp.β)
    vI1 = vN1 - K(c, cost)
    vN2 = (mp.β*Φ(0, -K(c, cost)))/(1 - mp.β)
    vI2 = vN2 - K(c, cost)
    P2 = 0.
    τ[end].EQs[end, end, 1] = EQ(P1, vN1, vI1, P2, vN2, vI2)
end

function solve_last_edge!(cost::CostParameters, mp::RLSModel, τ)
    nC = cost.nC
    c = cost.Cmin

    # if at (c1,c,c) edge
    for ic1 = 1:nC-1
        c1 = cost.C[ic1]
        vI1 =  - K(0., cost) + mp.β*Φ(τ[end].EQs[end, end, 1].vN1,τ[end].EQs[end, end, 1].vI1)
        vN1 = 0. # Since vI1 in the corner is negative!
        P1 = 0.
        vN2 = (r2(c1, c) + mp.β*(P1*0+(1-P1)*Φ(0,-K(c, cost))))/(1-mp.β*(1-P1))
        vI2 = vN2 - K(c, cost)
        P2 = 0.

        τ[end].EQs[ic1, end, 1] = EQ(P1, vN1, vI1, P2, vN2, vI2)
    end
    for ic2 = 1:nC-1
        c2 = cost.C[ic2]
        vI2 = r2(c, c2) - K(c, cost) + mp.β*Φ(τ[end].EQs[end, end, 1].vN2,τ[end].EQs[end, end, 1].vI2)
        vN2 = 0.
        P2 = 0.
        vN1 = (r1(c, c2) + mp.β*(P2*0+(1-P2)*Φ(0, -K(c, cost))))/(1 - mp.β*(1-P2))
        vI1 = vN1-K(c, cost)
        P1 = 0.

        τ[end].EQs[end, ic2, 1] = EQ(P1, vN1, vI1, P2, vN2, vI2)
    end
end

function solve_last_interior!(cost::CostParameters, mp::RLSModel, ic1::Int64, ic2::Int64, τ, ess)
            find_interior_last!(cost, mp, ic1, ic2, ess.nC, cost.C[ess.nC], τ, ess)
end

function solve_last_interior!(cost::CostParameters, mp::RLSModel, τ, ess)
    for ic2 in 1:ess.nC-1
        for ic1 in 1:ess.nC-1
            find_interior_last!(cost, mp, ic1, ic2, ess.nC, cost.C[ess.nC], τ, ess)
        end
    end
end

function find_interior_last!(cost::CostParameters, mp::RLSModel, ic1, ic2, ic, c, τ, ESS)
    pi = cost.π[ic]

    ## First player 1 -> leads to P2 candidates
    a = r1(cost.C[ic1], cost.C[ic2]) - K(c, cost) + mp.β*H1(ic, ic2, pi, τ, ESS)
    b = mp.β*(H1(ic, ic, pi, τ, ESS)-H1(ic, ic2, pi, τ, ESS))
    d = r1(cost.C[ic1],cost.C[ic2])
    e = mp.β*H1(ic1, ic, pi, τ, ESS)

    pa = - mp.β * pi * b
    pb = mp.β*H1(ic1, ic, pi, τ, ESS) + (mp.β*pi-1)*b - mp.β*pi*a
    pc = r1(cost.C[ic1],cost.C[ic2]) + (mp.β*pi-1) * a

    pstar2 = quad(pa, pb, pc)

    ## Then player 2 -> leads to P1 candidates
    A = r2(cost.C[ic1], cost.C[ic2]) - K(c, cost) + mp.β*H2(ic1, ic, pi, τ, ESS)
    B = mp.β*(H2(ic, ic, pi, τ, ESS)-H2(ic1, ic, pi, τ, ESS))
    D = r2(cost.C[ic1],cost.C[ic2])
    E = mp.β*H2(ic, ic2, pi, τ, ESS)

    qa = - mp.β * pi * B
    qb = mp.β*H2(ic, ic2, pi, τ, ESS) + (mp.β*pi-1) * B - mp.β*pi*A
    qc = r2(cost.C[ic1],cost.C[ic2]) + (mp.β*pi-1) * A
    pstar1 = quad(qa, qb, qc)

    count = 0
    for i = 1:length(pstar1)
        for j = 1:length(pstar2)
            if i in 1:2 && j in 1:2
                exP1 = pc + pb * pstar2[j] + pa * pstar2[j]^2 < 0. ? 1. : 0.
                exP2 = qc + qb * pstar1[i] + qa * pstar1[i]^2 < 0. ? 1. : 0.

                if abs(exP1 - pstar1[i]) < 1e-8 && abs(exP2-pstar2[j]) < 1e-8
                    count += 1

                    @inbounds τ[end].EQs[ic1, ic2, count] = EQ(pstar1[i],
                                                               (d + e*pstar2[j] + mp.β*pi*(1-pstar2[j])*(a+b*pstar2[j]))*pstar1[i]+(1-pstar1[i])*(d+e*pstar2[j])/(1-mp.β*pi*(1-pstar2[j])),
                                                               a + b*pstar2[j],
                                                               pstar2[j],
                                                               (D + E*pstar1[i] + mp.β*pi*(1-pstar1[i])*(A+B*pstar1[i]))*pstar2[j]+(1-pstar2[j])*(D+E*pstar1[i])/(1-mp.β*pi*(1-pstar1[i])),
                                                               A + B*pstar1[i])
                end
            elseif i > 2 && j > 2 && pstar1[i] >= 0 && pstar2[j] >= 0 && pstar1[i] <= 1 && pstar2[j] <= 1
                count += 1
                @inbounds τ[end].EQs[ic1, ic2, count] = EQ(pstar1[i],
                                                           a + b * pstar2[j],
                                                           a + b * pstar2[j],
                                                           pstar2[j],
                                                           A + B * pstar1[i],
                                                           A + B * pstar1[i])
            end
        end
    end
    ESS.bases[1 + div(ESS.nC*(ESS.nC+1)*(2*ESS.nC+1),6) - div(ic*(ic+1)*(2*ic+1),6) + 2*(ic - 1) + sub2ind((ic-1,ic-1),ic1,ic2)]=count
end


text:
__precompile__()

module DistributedArrays

using Distributed
using Serialization
using LinearAlgebra
using Statistics

import Base: +, -, *, div, mod, rem, &, |, xor
import Base.Callable
import LinearAlgebra: axpy!, dot, norm, mul!

import Primes
import Primes: factor

# DArray exports
export DArray, SubDArray, SubOrDArray, @DArray
export dzeros, dones, dfill, drand, drandn, distribute, localpart, localindices, ppeval

# non-array distributed data
export ddata, gather

# immediate release of localparts
export close, d_closeall

include("darray.jl")
include("core.jl")
include("serialize.jl")
include("broadcast.jl")
include("mapreduce.jl")
include("linalg.jl")
include("sort.jl")

include("spmd.jl")
export SPMD

end # module