text:
test("Add path to clipper") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = Clip()

    @test add_path!(c, path, PolyTypeSubject, true) == false

    push!(path, IntPoint(1, 1))

    @test add_path!(c, path, PolyTypeSubject, true) == true
end

test("Union") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(1, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(1, 0))
    push!(path2, IntPoint(1, 1))
    push!(path2, IntPoint(2, 1))
    push!(path2, IntPoint(2, 0))

    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)
    add_path!(c, path2, PolyTypeSubject, true)

    result, polys = execute(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

    @test result == true
    @test polys[1][1] == Clipper.IntPoint(0, 0)
    @test polys[1][2] == Clipper.IntPoint(2, 0)
    @test polys[1][3] == Clipper.IntPoint(2, 1)
    @test polys[1][4] == Clipper.IntPoint(0, 1)

    result, pt = execute_pt(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)
    @test result == true
    @test string(pt) == "Top-level PolyNode with 1 immediate children."
    @test length(children(pt)) === 1

    pn = children(pt)[1]

    @test length(children(pn)) == 0
    @test contour(pn)[1] == Clipper.IntPoint(0, 0)
    @test contour(pn)[2] == Clipper.IntPoint(2, 0)
    @test contour(pn)[3] == Clipper.IntPoint(2, 1)
    @test contour(pn)[4] == Clipper.IntPoint(0, 1)
end

test("Difference") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 10))
    push!(path1, IntPoint(10, 10))
    push!(path1, IntPoint(10, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(4, 0))
    push!(path2, IntPoint(4, 10))
    push!(path2, IntPoint(6, 10))
    push!(path2, IntPoint(6, 0))

    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)
    add_path!(c, path2, PolyTypeClip, true)

    result, polys = execute(c, ClipTypeDifference, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

    @test result == true
    @test polys[1][1] == Clipper.IntPoint(10, 10)
    @test polys[1][2] == Clipper.IntPoint(6, 10)
    @test polys[1][3] == Clipper.IntPoint(6, 0)
    @test polys[1][4] == Clipper.IntPoint(10, 0)

    @test polys[2][1] == Clipper.IntPoint(0, 10)
    @test polys[2][2] == Clipper.IntPoint(0, 0)
    @test polys[2][3] == Clipper.IntPoint(4, 0)
    @test polys[2][4] == Clipper.IntPoint(4, 10)

    result, pt = execute_pt(c, ClipTypeDifference, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)
    @test result == true
    @test string(pt) == "Top-level PolyNode with 2 immediate children."
    @test length(children(pt)) === 2

    pn1, pn2 = (children(pt)...,)
    @test length(children(pn1)) == 0
    @test length(children(pn2)) == 0

    @test contour(pn1)[1] == Clipper.IntPoint(10, 10)
    @test contour(pn1)[2] == Clipper.IntPoint(6, 10)
    @test contour(pn1)[3] == Clipper.IntPoint(6, 0)
    @test contour(pn1)[4] == Clipper.IntPoint(10, 0)

    @test contour(pn2)[1] == Clipper.IntPoint(0, 10)
    @test contour(pn2)[2] == Clipper.IntPoint(0, 0)
    @test contour(pn2)[3] == Clipper.IntPoint(4, 0)
    @test contour(pn2)[4] == Clipper.IntPoint(4, 10)
end

test("GetBounds") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(1, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(1, 0))
    push!(path2, IntPoint(1, 1))
    push!(path2, IntPoint(2, 1))
    push!(path2, IntPoint(2, 0))

    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)
    add_path!(c, path2, PolyTypeSubject, true)

    rect = get_bounds(c)

    @test rect.left == 0
    @test rect.top == 0
    @test rect.right == 2
    @test rect.bottom == 1
end

test("Clear") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 10))
    push!(path1, IntPoint(10, 10))
    push!(path1, IntPoint(10, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(4, 0))
    push!(path2, IntPoint(4, 10))
    push!(path2, IntPoint(6, 10))
    push!(path2, IntPoint(6, 0))

    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)
    add_path!(c, path2, PolyTypeSubject, true)

    Clipper.clear!(c)

    rect = get_bounds(c)

    @test rect.left == 0
    @test rect.top == 0
    @test rect.right == 0
    @test rect.bottom == 0
end

test("AddPaths") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(1, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(1, 0))
    push!(path2, IntPoint(1, 1))
    push!(path2, IntPoint(2, 1))
    push!(path2, IntPoint(2, 0))

    paths = Vector{IntPoint}[path1, path2]

    c = Clip()
    add_paths!(c, paths, PolyTypeSubject, true)

    result, polys = execute(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

    @test result == true
    @test polys[1][1] == Clipper.IntPoint(0, 0)
    @test polys[1][2] == Clipper.IntPoint(2, 0)
    @test polys[1][3] == Clipper.IntPoint(2, 1)
    @test polys[1][4] == Clipper.IntPoint(0, 1)
end

test("Orientation") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(1, 0))

    @test orientation(path1) == false

    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(1, 0))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(0, 0))

    @test orientation(path1) == true
end

test("Area") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))

    @test area(path1) == -0.5
end

test("Point in polygon") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(4, 0))
    push!(path1, IntPoint(0, 4))

    # Returns -1 if on boundary (check all the boundary points)
    @test pointinpolygon(IntPoint(0, 0), path1) == -1
    @test pointinpolygon(IntPoint(1, 0), path1) == -1
    @test pointinpolygon(IntPoint(2, 0), path1) == -1
    @test pointinpolygon(IntPoint(3, 0), path1) == -1
    @test pointinpolygon(IntPoint(4, 0), path1) == -1
    @test pointinpolygon(IntPoint(3, 1), path1) == -1
    @test pointinpolygon(IntPoint(2, 2), path1) == -1
    @test pointinpolygon(IntPoint(1, 3), path1) == -1
    @test pointinpolygon(IntPoint(0, 4), path1) == -1
    @test pointinpolygon(IntPoint(0, 3), path1) == -1
    @test pointinpolygon(IntPoint(0, 2), path1) == -1
    @test pointinpolygon(IntPoint(0, 1), path1) == -1

    # Returns 1 if inside (check all the interior points)
    @test pointinpolygon(IntPoint(1, 1), path1) == 1
    @test pointinpolygon(IntPoint(2, 1), path1) == 1
    @test pointinpolygon(IntPoint(1, 2), path1) == 1

    # Returns 0 if outside (check a few places outside)
    @test pointinpolygon(IntPoint(10, 10), path1) == 0
    @test pointinpolygon(IntPoint(-1, -1), path1) == 0
end

struct IntPoint2
    X::Int64
    Y::Int64
end
Base.convert(::Type{IntPoint2}, x::IntPoint) = IntPoint2(x.X, x.Y)
Base.convert(::Type{IntPoint}, x::IntPoint2) = IntPoint(x.X, x.Y)

test("PolyTrees / PolyNodes") do
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(8, 8))
    push!(path1, IntPoint(0, 8))
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(8, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(1, 1))
    push!(path2, IntPoint(1, 7))
    push!(path2, IntPoint(7, 7))
    push!(path2, IntPoint(7, 1))

    path3 = Vector{IntPoint}()
    push!(path3, IntPoint(6, 6))
    push!(path3, IntPoint(2, 6))
    push!(path3, IntPoint(2, 2))
    push!(path3, IntPoint(6, 2))

    paths = Vector{IntPoint}[path1, path2, path3]

    c = Clip()
    add_paths!(c, paths, PolyTypeSubject, true)

    result, pt = execute_pt(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

    # test expected PolyTree structure
    @test result == true
    @test isa(pt, PolyNode{IntPoint})
    @test parent(pt) === pt     # in the wrapper we set the parent of top level to itself
    @test length(children(pt)) === 1
    @test contour(pt) == IntPoint[]     # top level has no contour
    @test length(children(pt)) == 1

    pn1 = children(pt)[1]
    @test !ishole(pn1)
    @test !isopen(pn1)
    @test contour(pn1) == path1
    @test length(children(pn1)) === 1
    @test parent(pn1) === pt
    @test length(children(pn1)) == 1
    @test contour(pn1)[1] == Clipper.IntPoint(8, 8)
    @test contour(pn1)[2] == Clipper.IntPoint(0, 8)
    @test contour(pn1)[3] == Clipper.IntPoint(0, 0)
    @test contour(pn1)[4] == Clipper.IntPoint(8, 0)

    pn2 = children(pn1)[1]
    @test ishole(pn2)
    @test !isopen(pn2)
    @test contour(pn2) == path2
    @test length(children(pn2)) === 1
    @test parent(pn2) === pn1
    @test length(children(pn2)) == 1
    @test contour(pn2)[1] == Clipper.IntPoint(1, 1)
    @test contour(pn2)[2] == Clipper.IntPoint(1, 7)
    @test contour(pn2)[3] == Clipper.IntPoint(7, 7)
    @test contour(pn2)[4] == Clipper.IntPoint(7, 1)

    pn3 = children(pn2)[1]
    @test !ishole(pn3)
    @test !isopen(pn3)
    @test contour(pn3) == path3
    @test isempty(children(pn3))
    @test parent(pn3) === pn2
    @test length(children(pn3)) == 0
    @test contour(pn3)[1] == Clipper.IntPoint(6, 6)
    @test contour(pn3)[2] == Clipper.IntPoint(2, 6)
    @test contour(pn3)[3] == Clipper.IntPoint(2, 2)
    @test contour(pn3)[4] == Clipper.IntPoint(6, 2)

    # Test that we can preserve the tree structure when converting between types.
    pt2 = convert(PolyNode{IntPoint2}, pt)
    pt3 = convert(PolyNode{IntPoint}, pt2)

    @test result == true
    @test isa(pt3, PolyNode{IntPoint})
    @test parent(pt3) === pt3     # in the wrapper we set the parent of top level to itself
    @test length(children(pt3)) === 1
    @test contour(pt3) == IntPoint[]     # top level has no contour

    pn1 = children(pt3)[1]
    @test !ishole(pn1)
    @test !isopen(pn1)
    @test contour(pn1) == path1
    @test length(children(pn1)) === 1
    @test parent(pn1) === pt3

    pn2 = children(pn1)[1]
    @test ishole(pn2)
    @test !isopen(pn2)
    @test contour(pn2) == path2
    @test length(children(pn2)) === 1
    @test parent(pn2) === pn1

    pn3 = children(pn2)[1]
    @test !ishole(pn3)
    @test !isopen(pn3)
    @test contour(pn3) == path3
    @test isempty(children(pn3))
    @test parent(pn3) === pn2

    # Only convert top-level PolyNodes (i.e. PolyTrees)
    @test_throws ErrorException convert(PolyNode{IntPoint2}, pn3)
end

test("IntPoint from floating point") do
    value1 = 58.4309
    value2 = 9.39547
    value3 = 10.393
    value4 = 132.45
    value5 = 0.0045

    test1 = IntPoint(value1, value2, 2, 4)
    @test IntPoint(5843, 940) == test1

    test2 = IntPoint(value3, value4, 3, 5)
    @test IntPoint(1039, 13245) == test2

    test3 = IntPoint(value1, value4, 6, 2)
    @test IntPoint(0, 0) == test3

    test4 = IntPoint(value2, value3, 1, 4)
    @test IntPoint(9395, 10393) == test4

    test5 = IntPoint(value4, value5, -2, 5)
    @test IntPoint(1324500000, 45000) == test5

    test6 = IntPoint(value4, value5, 3, 5)
    @test IntPoint(13245, 0) == test6
end

test("IntPoint to floating point") do
    point = IntPoint(13245, 340)

    @test tofloat(point, 3, 5) == (132.45, 3.4)
    @test tofloat(point, -4, 5) == (0.000013245, 0.00000034)
    @test tofloat(point, 3, 4) == (1324.5, 34)
    @test tofloat(point, -4, 4) == (0.00013245, 0.0000034)
end

test("Minkowski sum") do
	tri = [IntPoint(0,0), IntPoint(10,0), IntPoint(0,20)]
	@test minkowski_sum(tri,tri) ==
		[[IntPoint(0,0), IntPoint(20,0), IntPoint(0,40)]]
end


text:
module CloudFormation

using AWSCore
using AWSCore.Services: cloudformation
using EzXML
using MbedTLS: MbedException
using Memento
using Mocking
using OrderedCollections: OrderedDict
using XMLDict

const logger = getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(MyModule)`
# NOTE: If this line is not included then the precompiled `MyModule.logger` won't be registered at runtime.
__init__() = Memento.register(logger)

export raw_stack_description, stack_output, stack_description

describe_stacks(config; kwargs...) = cloudformation(config, "DescribeStacks"; kwargs...)

const _NRETRIES = 5

# for APIs we don't want to hammer
function cautious_delays(; kwargs...)
    return ExponentialBackOff(; n=_NRETRIES, first_delay=5, max_delay=300, kwargs...)
end

function minimal_delays(; kwargs...)
    return ExponentialBackOff(; n=_NRETRIES, first_delay=0.1, max_delay=60, kwargs...)
end

"""
    raw_stack_description(stack_name::AbstractString) -> String

Returns the description for the specified stack. Can optionally pass in the aws `config`
as a keyword argument.
"""
function raw_stack_description(
    stack_name::AbstractString;
    config::AWSConfig=aws_config()
)
    function retry_cond(s, e)
        if e isa AWSException
            if 500 <= e.cause.status <= 504
                debug(logger, "CloudFormation request encountered $(e.code); retrying")
                return (s, true)
            elseif e.cause.status == 429 || (e.cause.status == 400 && e.code == "Throttling")
                debug(logger, "CloudFormation request encountered $(e.code); retrying")
                return (s, true)
            end
        elseif e isa MbedException
            debug(logger, "CloudFormation request encountered $e; retrying")
            return (s, true)
        end

        return (s, false)
    end

    f = retry(delays=cautious_delays(; jitter=0.2), check=retry_cond) do
        aws_raw_config = deepcopy(config)
        aws_raw_config[:return_raw] = true
        @mock describe_stacks(aws_raw_config, StackName=stack_name)
    end

    response = String(f())

    return response
end

"""
    stack_output(stack_name::AbstractString) -> OrderedDict

The stack's OutputKey and OutputValue values as a dictionary. Can pass in the aws `config`
as a keyword argument.
"""
function stack_output(stack_name::AbstractString; config::AWSConfig=aws_config())
    outputs = OrderedDict{String, String}()
    description = raw_stack_description(stack_name; config=config)

    xml = root(parsexml(description))
    ns = ["ns" => namespace(xml)]
    output_elements = findall("//ns:Outputs/ns:member", xml, ns)

    for el in output_elements
        key = nodecontent(findfirst("ns:OutputKey", el, ns))
        val = nodecontent(findfirst("ns:OutputValue", el, ns))
        outputs[key] = val
    end

    return outputs
end

function output_pair(item::AbstractDict)
    key = item["OutputKey"]::String
    value = if isa(item["OutputValue"], String)
        item["OutputValue"]::String
    elseif isa(item["OutputValue"], AbstractDict) && isempty(item["OutputValue"])
        ""
    else
        throw(ArgumentError("Unhandled output value: $(repr(item["OutputValue"]))"))
    end

    return key => value
end

# BEGIN AWSTools.Cloudformation 0.8.1 deprecations

function stack_description(
    stack_name::AbstractString;
    config::AWSConfig=aws_config()
)
    dep_msg = """
        `stack_description(::AbstractString; ::AWSConfig)` is deprecated and will be removed.
        Please use `raw_stack_description(::AbstractString; ::AWSConfig)` instead and handle XML parsing in the calling function.
        We recommend using EzXML.
        """
    Base.depwarn(dep_msg, :stack_description)

    response = xml_dict(raw_stack_description(stack_name))

    return response["DescribeStacksResponse"]["DescribeStacksResult"]["Stacks"]["member"]
end

# END AWSTools.Cloudformation 0.8.1 deprecations

end  # CloudFormation

ext: jl


text:
# Tesing internals for elliptic curve parsing with verificatum
using Test

import ShuffleProofs: marshal_s_Gq, unmarshal, decode, marshal, unmarshal_publickey, marshal_publickey, marshal_privatekey, unmarshal_privatekey, ElGamal, load_verificatum_proposition, load_verificatum_proof, load_verificatum_simulator, ro_prefix, verify
import CryptoGroups: spec

s_Gq = "com.verificatum.arithm.ECqPGroup(P-256)::00000000020100000020636f6d2e766572696669636174756d2e61726974686d2e4543715047726f75700100000005502d323536"

tree = decode(split(s_Gq, "::")[2])
g = unmarshal(tree)

tree′ = marshal(g)

@test tree == tree′

### Encoding and decoding public key

BASE_DIR = "$(@__DIR__)/../../ref/P256/"
#auxsid = "default"
#NIZKP = "$(@__DIR__)/../../ref/demo/dir/nizkp/$auxsid/"

PUBLIC_KEY = "$BASE_DIR/publicKey"

pk_tree = decode(read(PUBLIC_KEY))
pk, g = unmarshal_publickey(pk_tree)

# So that things can be properly understood
@test unmarshal_publickey(marshal_publickey(pk, g)) == (pk, g)
@test pk_tree == marshal_publickey(pk, g)

sk = BigInt(123234)

@test unmarshal_privatekey(marshal_privatekey(g, sk)) == (sk, g)

### Decoding of ciphertexts

CIPHERTEXT_FILE = "$BASE_DIR/ciphertexts"

G = typeof(g)

bytes = read(CIPHERTEXT_FILE)
tree = decode(bytes)
𝐚, 𝐛 = convert(Tuple{Vector{G}, Vector{G}}, tree)

ciphertexts = ElGamal{G}(𝐚, 𝐛)

### Loading of the whole Proposition

proposition = load_verificatum_proposition(BASE_DIR, "default")

proof = load_verificatum_proof(BASE_DIR * "/dir/nizkp/default/proofs/", g)

simulator = load_verificatum_simulator(BASE_DIR)

@test ro_prefix(simulator.verifier) == hex2bytes("355806458d6cd42655a52be242705c8e824584ccdb6b1c016cad36c591413de4")

@test verify(simulator)


text:
using DataStructures
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "DictionarySorting" begin
    forward = OrderedDict(zip('a':'z', 26:-1:1))
    rev = OrderedDict(zip(reverse('a':'z'), 1:26))

    d = copy(rev)
    @test d == rev
    @test sort(d) == forward
    sort!(d)
    @test d == forward
    @test sort(d; rev=true) == rev
    sort!(d; rev=true)
    @test d == rev
    @test sort(d; byvalue=true) == rev
    @test sort(d; byvalue=true, rev=true) == forward
    sort!(d; byvalue=true, rev=true)
    @test d == forward
    @test sort(d; byvalue=true) == rev
    @test sort(d; byvalue=true, rev=true) == forward
    sort!(d; byvalue=true)
    @test d == rev

    unordered = Dict(zip('a':'z', 26:-1:1))
    @test sort(unordered) == forward
    @test sort(unordered; rev=true) == rev
    @test sort(unordered; byvalue=true) == rev
    @test sort(unordered; byvalue=true, rev=true) == forward
end


text:
using PyPlot

using TransitionalMCMC

# Prior Bounds
lb, ub  = -5, 5

# Prior Density and sampler
priorDen(x) = pdf(Uniform(lb, ub), x[1,:]) .* pdf(Uniform(lb, ub), x[2,:])
priorRnd(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

# Log Likelihood
function logLik(x)
    sleep(0.005)
    return -1 .* ((x[1,:].^2 .+ x[2,:] .- 11).^2 .+ (x[1,:] .+ x[2,:].^2 .- 7).^2)
end

Nsamples = 200

samps, acc = tmcmc(logLik, priorDen, priorRnd, Nsamples, 5, 2)

plt.scatter(samps[:,1], samps[:,2])


text:
####
# Some Closed-Form Formulae and Helpers
####
#returns max coeff_dev for dist S, M
comp_Sc(S, M) = (S + M - 1)/M
max_dev(S, M) = M * (S - 1) / (S + M - 1)
function max_cv(S, M)
	Sc = comp_Sc(S, M)
	M * sqrt((Sc - 1))
end

#Returns the minimum B
function min_gm(S, M, mu)
	Sc = comp_Sc(S, M)
	out = (Sc - 1)/Sc * log(1 - M) + 1/Sc * log(M * (S - 1) + 1)
	mu * exp(out)
end

#helper function for a robust (indefinite) quadratic
#ax^2 + bx +c >= 0 for all x in [l, u]
#adds everything to model m
function _addRQ!(m, a, b, c, l, u)
    y1 = @variable(m)  #theta
    @constraint(m, y1 >= 0)
    y = @variable(m)    #t
    @constraint(m, y >= 0)
    t2 = @variable(m)
    @constraint(m, 4 * y1 * y >= t2^2)
    @constraint(m, t2 == b - (l + u) * (y1 - a))
    @constraint(m, y1 - a >= 0)
    @constraint(m, c + l * u * (y1 - a) - y >= 0)
end

#conjecture that maximal CV occurs for uniform distribution among unimodals
#Not formally proven
function max_std_cv_unimodal_guess(mu, S, M, isSym)
    #convert to standardized distribution
    c = mu * (1 - M)
    Sc = vopp.comp_Sc(S, M)

    #adjust for symmetry if needed
    if isSym
        if Sc > 2
            Sc = 2
        end
        if Sc < 2
            lb = 2 - Sc
        else
            lb = 0
        end
    end

    return 2 / sqrt(12) * (Sc - lb) / (Sc + lb)
end

#for the standardized distribution
#document indexes 0 to N
function geom_price_ladder(S, delta)
    N = ceil(Int, 1 + log(S/delta) / log(1 + delta))
    ps = zeros(N + 1)
    ps[0 + 1] = 0.
    ps[1 + 1] = delta
    ps[N + 1] = S
    for i = 2:N-1
        ps[i + 1] = ps[i - 1 + 1] * (1 + delta)
    end
    ps
end

vopp_ub_scale(S, M) = -lambertw(-M / exp(1) / (S + M - 1))

#assumes standardized
function tight_dist_ub_scale(x, S)
    alpha = 1/vopp_ub_scale(S, 1)
    if 0 <= x <= alpha
        return 1.
    elseif x <= S
        return alpha / x
    else
        return 0.
    end
    return -1.
end

#Approximately computes the minimum and maximal valid q for distn
#assumes it is standardized
function min_max_q_IC_uni(S, mode, v0, supp_size=100)

    function G(p, m, t)
        if p > max(m, t)
            return 0.
        elseif p < min(m, t)
            return 1.
        elseif m == t
            return m >= p ? 1. : 0.
        end
        return (max(m, t) - p) / abs(m-t)
    end

    m = Model(solver=GurobiSolver(OutputFlag=false))

    t_grid = range(0, stop=S, length=supp_size)

    @variable(m, ws[1:supp_size] >= 0)
    @constraint(m, sum(ws) == 1)
    @constraint(m, sum(ws[i] * (mode + t_grid[i])/2 for i = 1:supp_size) == 1)
    @objective(m, Min, sum(ws[i] * G(v0, mode, t_grid[i]) for i = 1:supp_size))

    status = solve(m)
    @assert status == :Optimal
    minq = getobjectivevalue(m)

    @objective(m, Max, sum(ws[i] * G(v0, mode, t_grid[i]) for i = 1:supp_size))

    status = solve(m)
    @assert status == :Optimal
    maxq = getobjectivevalue(m)
    return minq, maxq
end

function min_max_q_IC_uni(S, M, mu, mode, phat, supp_size=100)
    c = mu * (1 - M)
    Sc = vopp.comp_Sc(S, M)
    mode_c = (mode - c) / (mu - c)
    v0 = (phat - 1)/ M + 1
    return min_max_q_IC_uni(Sc, mode_c, v0, supp_size)

end


text:
"""
    Shape{T1, T2, T3, T4, T5, T6, T7}

A polyhedral shape model of an asteroid.

# Fields

- `num_node` : Number of nodes
- `num_face` : Number of faces
- `nodes`    : 1-D array of node positions
- `faces`    : 1-D array of vertex indices of faces

- `facets`   : 1-D array of surface facets (`Facet`)

- `AREA`     : Surface area
- `VOLUME`   : Volume
- `COF`      : Center-of-figure
- `MOI`      : Moment of inertia tensor

– `force`    : Thermal recoil force at body-fixed frame (Yarkovsky effect)
- `torque`   : Thermal recoil torque at body-fixed frame (YORP effect)
- `Tz⁺`      : Pre-allocated vector for update of temperature profile on each facets (`Facet`)
"""
struct Shape{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10}
    num_node::T1
    num_face::T1
    nodes   ::T2
    faces   ::T3

    facets::T4

    AREA  ::T5
    VOLUME::T5
    COF   ::T6
    MOI   ::T7

    force ::T8
    torque::T9
    Tz⁺   ::T10
end

function Base.show(io::IO, shape::Shape)
    # println(io, "Shape model")
    println("Shape model")
    println("-----------")

    println("Nodes             : ", shape.num_node)
    println("Faces             : ", shape.num_face)
    println("Surface area      : ", shape.AREA)
    println("Volume            : ", shape.VOLUME)
    println("Equivalent radius : ", equivalent_radius(shape))
    println("Center-of-Figure  : ", shape.COF)
    println("Inertia tensor    : ")
    println("    | Ixx Ixy Ixz |   ", shape.MOI[1, :])
    println("    | Iyx Iyy Iyz | = ", shape.MOI[2, :])
    println("    | Izx Izy Izz |   ", shape.MOI[3, :])
end

function Shape(shapepath; scale=1, find_visible_facets=false, save_shape=false)

    ext = splitext(shapepath)[2]

    if ext == ".obj"
        nodes, faces = loadobj(shapepath; scale=scale, static=true, message=false)

        num_node = length(nodes)
        num_face = length(faces)

        facets = getfacets(nodes, faces)
        find_visible_facets && findVisibleFacets!(facets)

        AREA   = sum(facets.area)
        VOLUME = getvolume(facets)
        COF    = center_of_figure(facets)
        MOI    = moment_of_inertia(facets)

        force  = zeros(3)
        torque = zeros(3)
        Tz⁺    = similar(facets[begin].Tz)

        shape = Shape(
                num_node, num_face, nodes, faces, facets,
                AREA, VOLUME, COF, MOI,
                force, torque, Tz⁺
        )
        save_shape && save(splitext(shapepath)[1] * ".jld2", Dict("shape" => shape))

    elseif ext == ".jld2"
        shape = load(shapepath)["shape"]
    else
        return ArgumentError("Give a filepath of *.obj or *.jld2.")
    end

    return shape
end

equivalent_radius(VOLUME) = (3VOLUME/4π)^(1/3)
equivalent_radius(shape::Shape) = equivalent_radius(shape.VOLUME)

findVisibleFacets!(shape::Shape) = findVisibleFacets!(shape.facets)
isIlluminated(obs::Facet, r̂☉, shape::Shape) = isIlluminated(obs, r̂☉, shape.facets)
isIlluminated(r̂☉, shape::Shape) = [isIlluminated(obs, r̂☉, shape) for obs in shape.facets]

surface_temperature(shape) = [facet.Tz[begin] for facet in shape.facets]

################################################################
#                      Shape properites
################################################################

"""
    getvolume(facets) -> VOLUME

Calculate volume of a polyhedral
"""
getvolume(facets) = sum(((facet.A × facet.B) ⋅ facet.C) / 6 for facet in facets)

"""
    center_of_figure(facets) -> COF

Calculate center-of-figure position of a polyhedral
"""
function center_of_figure(facets)
    VOLUME = getvolume(facets)
    COF = zeros(3)

    for facet in facets
        @unpack A, B, C = facet
        volume = ((A × B) ⋅ C) / 6  # volume of pyramid element O-A-B-C
        center = (A + B + C) / 4    # center of pyramid element O-A-B-C
        COF += volume * center
    end
    COF / VOLUME
end

"""
    moment_of_inertia(facets) -> MOI

Calculate moment of inertia tensor of a polyhedron
"""
function moment_of_inertia(facets)
    MOI = zeros(3, 3)

    for facet in facets
        # v1, v2, v3 = m.vs

        # I[3, 3] += (v1[1]*v1[1] + v1[1]*v2[1] + v2[1]*v2[1] + v1[1]*v3[1] + v2[1]*v3[1] + v3[1]*v3[1] + v1[2]*v1[2] + v1[2]*v2[2] + v2[2]*v2[2] + v1[2]*v3[2] + v2[2]*v3[2] + v3[2]*v3[2]) / 60
    end
    MOI
end

################################################################
#              3D visualization of polyhedral shape
################################################################

"""
Vector of vector to 2-D matrix
"""
function VectorVector2Matrix(v)
    m = Matrix{eltype(v[end])}(undef, length(v), 3)
    for i in eachindex(v)
        m[i, :] .= v[i]
    end
    m
end

"""
    face2node(nodes, faces, data) -> node_data

Convert face-based data to node-based data

# Arguments
- `nodes` : 2-D matrix of nodes
- `faces` : 2-D matrix of faces
- `data`  : Face-based data

# Return
- `node_data` : Node-based data
"""
function face2node(nodes, faces, data)
    node_data = Vector{eltype(data)}(undef, size(nodes, 1))

    for node_index in eachindex(node_data)
        face_indices = findall(faces.==node_index)
        if length(face_indices) == 0
            node_data[node_index] = NaN
            continue
        end
        node_data[node_index] = mean(data[face_index[1]] for face_index in face_indices)
    end
    node_data
end

"""
    showshape(shape; data=nothing)

"""
function draw(shape::Shape; data=nothing, r̂☉=[1,0,0.])
    nodes = VectorVector2Matrix(shape.nodes)
    faces = VectorVector2Matrix(shape.faces)

    set_theme!(backgroundcolor=:black)

    if data == nothing
        color = :gray
    elseif data == :radius
        color = [norm(v) for v in eachrow(nodes)]
    elseif data == :temperature
        surf_temps = surface_temperature(shape)
        color = face2node(nodes, faces, surf_temps)
    elseif data == :illumination
        illumination = Float64.(isIlluminated(r̂☉, shape))
        color = face2node(nodes, faces, illumination)
    else
        color = face2node(nodes, faces, data)
    end

    scene = mesh(nodes, faces, color=color)
    display(scene)
end

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
## Packages
using Plots

## Files
include("VFI.jl")
include("LQ_Vaughan.jl")

################################### Running three methods ########################
md = Model();
kss,hss,lss,css = getSS(md);

## Construct capital grid
nk = 1000; # number of capital grid
kmin = 0.5*kss;
kmax = 1.5*kss;
kGrid = LinRange(kmin,kmax,nk);

V_VFI,kpol_VFI,hpol_VFI,cpol_VFI = VFI(md,nk,kGrid)

A_tld,B_tld,Q_tld,R,W = getM(md,kss,hss);
P_LQ,F_LQ = LQ(A_tld,B_tld,R,Q_tld);
Pv,Fv = Vaughan(A_tld,B_tld,Q_tld,R,W);

################################### Plotting ########################
## VFI
plot(kGrid,V_VFI[:,1])
plot!(kGrid,V_VFI[:,2])

plot(kGrid,kGrid,label = "", color = :black, linestyle = :dash)
plot!(kGrid,kpol_VFI[:,1],label = "low",legend = :topleft)
plot!(kGrid,kpol_VFI[:,3],label = "SS",legend = :topleft)
plot!(kGrid,kpol_VFI[:,5],label = "high",legend = :topleft)

plot(kGrid,hpol[:,1],label = "low",legend = :topleft)
plot!(kGrid,hpol[:,3],label = "high",legend = :topleft)

plot(kGrid,cpol[:,1],label = "low",legend = :topleft)
plot!(kGrid,cpol[:,3],label = "high",legend = :topleft)

## LQ and Vaughan
# Low shock
pol_L_LQ = zeros(2, 1, nk);
pol_L_v = zeros(2, 1, nk);
for i = 1:nk
    pol_L_LQ[:,:,i] = -F_LQ*[kGrid[i]; -0.5; 1];
    pol_L_v[:,:,i] = -Fv*[kGrid[i]; -0.5; 1];
end

kpol_L_LQ = zeros(1,nk);
hpol_L_LQ = zeros(1,nk);
kpol_L_v = zeros(1,nk);
hpol_L_v = zeros(1,nk);

for i = 1:nk
    kpol_L_LQ[1,i] = pol_L_LQ[1,:,i][1];
    hpol_L_LQ[1,i] = pol_L_LQ[2,:,i][1];
    kpol_L_v[1,i] = pol_L_v[1,:,i][1];
    hpol_L_v[1,i] = pol_L_v[2,:,i][1];
end

# High shock
pol_H_LQ = zeros(2, 1, nk);
pol_H_v = zeros(2, 1, nk);

for i = 1:nk
    pol_H_LQ[:,:,i] = -F_LQ*[kGrid[i]; 0.5; 1];
    pol_H_v[:,:,i] = -Fv*[kGrid[i]; 0.5; 1];
end

kpol_H_LQ = zeros(1,nk);
hpol_H_LQ = zeros(1,nk);
kpol_H_v = zeros(1,nk);
hpol_H_v = zeros(1,nk);

for i = 1:nk
    kpol_H_LQ[1,i] = pol_H_LQ[1,:,i][1];
    hpol_H_LQ[1,i] = pol_H_LQ[2,:,i][1];
    kpol_H_v[1,i] = pol_H_v[1,:,i][1];
    hpol_H_v[1,i] = pol_H_v[2,:,i][1];
end

# SS shock
pol_SS_LQ = zeros(2, 1, nk);
pol_SS_v = zeros(2, 1, nk);

for i = 1:nk
    pol_SS_LQ[:,:,i] = -F_LQ*[kGrid[i]; 0.0; 1];
    pol_SS_v[:,:,i] = -Fv*[kGrid[i]; 0.0; 1];
end

kpol_SS_LQ = zeros(1,nk);
hpol_SS_LQ = zeros(1,nk);
kpol_SS_v = zeros(1,nk);
hpol_SS_v = zeros(1,nk);

for i = 1:nk
    kpol_SS_LQ[1,i] = pol_SS_LQ[1,:,i][1];
    hpol_SS_LQ[1,i] = pol_SS_LQ[2,:,i][1];
    kpol_SS_v[1,i] = pol_SS_v[1,:,i][1];
    hpol_SS_v[1,i] = pol_SS_v[2,:,i][1];
end

################### Plots ##################################
# capital policy LQ
plot(kGrid,kGrid,label = "45 degree line", color = :black, linestyle = :dash,
title = "Optimal capital policy function (LQ)",xlabel = "k_{t}", ylabel = "k_{t+1}")
plot!(kGrid,vec(kpol_L_LQ),label = "Low_LQ",legend=:topleft)
plot!(kGrid,vec(kpol_H_LQ),label = "High_LQ")
plot!(kGrid,vec(kpol_SS_LQ),label = "SS_LQ")
savefig("figs/kpol_LQ.png")

# labor policy LQ
plot(kGrid,vec(hpol_L_LQ),label = "Low",legend=:topright,
title = "Optimal labor policy function (LQ)",xlabel = "k_{t}", ylabel = "h_{t}")
plot!(kGrid,vec(hpol_H_LQ),label = "High")
plot!(kGrid,vec(hpol_SS_LQ),label = "SS")
savefig("figs/hpol_LQ.png")

# capital policy Vaughan
plot(kGrid,kGrid,label = "45 degree line", color = :black, linestyle = :dash,
title = "Optimal capital policy function (Vaughan)",xlabel = "k_{t}", ylabel = "k_{t+1}")
plot!(kGrid,vec(kpol_L_v),label = "Low_Vaughan",legend=:topleft)
plot!(kGrid,vec(kpol_H_v),label = "High_Vaughan")
plot!(kGrid,vec(kpol_SS_v),label = "SS_Vaughan")
savefig("figs/kpol_Vaughan.png")

# labor policy Vaughan
plot(kGrid,vec(hpol_L_v),label = "Low",legend=:topright,
title = "Optimal labor policy function (Vaughan)",xlabel = "k_{t}", ylabel = "h_{t}")
plot!(kGrid,vec(hpol_H_v),label = "High")
plot!(kGrid,vec(hpol_SS_v),label = "SS")
savefig("figs/hpol_Vaughan.png")


text:
export stage, cached_stage, compute, debug_compute, free!, cleanup

###### Scheduler #######

compute(x; options=nothing) = compute(Context(), x; options=options)
compute(ctx, c::Chunk; options=nothing) = c

collect(ctx::Context, t::Thunk; options=nothing) =
    collect(ctx, compute(ctx, t; options=options); options=options)
collect(d::Union{Chunk,Thunk}; options=nothing) =
    collect(Context(), d; options=options)

abstract type Computation end

compute(ctx, c::Computation; options=nothing) =
    compute(ctx, stage(ctx, c); options=options)
collect(c::Computation; options=nothing) =
    collect(Context(), c; options=options)

"""
    compute(ctx::Context, d::Thunk; options=nothing) -> Chunk

Compute a Thunk - creates the DAG, assigns ranks to nodes for tie breaking and
runs the scheduler with the specified options. Returns a Chunk which references
the result.
"""
function compute(ctx::Context, d::Thunk; options=nothing)
    if !(:scheduler in keys(PLUGINS))
        PLUGINS[:scheduler] = get_type(PLUGIN_CONFIGS[:scheduler])
    end
    scheduler = PLUGINS[:scheduler]
    (scheduler).compute_dag(ctx, d; options=options)
end

function debug_compute(ctx::Context, args...; profile=false, options=nothing)
    @time res = compute(ctx, args...; options=options)
    get_logs!(ctx.log_sink), res
end

function debug_compute(arg; profile=false, options=nothing)
    ctx = Context()
    dbgctx = Context(procs(ctx), LocalEventLog(), profile)
    debug_compute(dbgctx, arg; options=options)
end

Base.@deprecate gather(ctx, x) collect(ctx, x)
Base.@deprecate gather(x) collect(x)

cleanup() = cleanup(Context())
function cleanup(ctx::Context)
    if :scheduler in keys(PLUGINS)
        scheduler = PLUGINS[:scheduler]
        (scheduler).cleanup(ctx)
        delete!(PLUGINS, :scheduler)
    end
    nothing
end

function get_type(s::String)
    local T
    for t in split(s, ".")
        t = Symbol(t)
        if !@isdefined(T)
            T = Base.require(@__MODULE__, t)
        else
            T = Core.eval(T, t)
        end
    end
    T
end

##### Dag utilities #####

"""
    dependents(node::Thunk, deps=Dict{Thunk, Set{Thunk}}()) -> Dict{Thunk, Set{Thunk}}

Find the set of direct dependents for each task.
"""
function dependents(node::Thunk, deps=Dict{Thunk, Set{Thunk}}())
    if !haskey(deps, node)
        deps[node] = Set{Thunk}()
    end
    for inp = inputs(node)
        if isa(inp, Thunk)
            s::Set{Thunk} = get!(()->Set{Thunk}(), deps, inp)
            push!(s, node)
            dependents(inp, deps)
        end
    end
    deps
end

"""
    noffspring(dpents::Dict{Thunk, Set{Thunk}}) -> Dict{Thunk, Int}

Recursively find the number of tasks dependent on each task in the DAG.
Takes a Dict as returned by [`dependents`](@ref).
"""
function noffspring(dpents::Dict{Thunk, Set{Thunk}})
    Dict(node => noffspring(node, dpents) for node in keys(dpents))
end

function noffspring(n, dpents)
    if haskey(dpents, n)
        ds = dpents[n]
        reduce(+, (noffspring(d, dpents) for d in ds), init = length(ds))
    else
        0
    end
end

"""
    order(node::Thunk, ndeps) -> Dict{Thunk,Int}

Given a root node of the DAG, calculates a total order for tie-breaking.

  * Root node gets score 1,
  * rest of the nodes are explored in DFS fashion but chunks
    of each node are explored in order of `noffspring`,
    i.e. total number of tasks depending on the result of the said node.

Args:
- node: root node
- ndeps: result of [`noffspring`](@ref)
"""
function order(node::Thunk, ndeps)
    function recur(nodes, s)
        for n in nodes
            output[n] = s += 1
            parents = collect(Iterators.filter(istask, inputs(n)))
            s = recur(sort!(parents, by=k->get(ndeps,k,0)), s)
        end
        return s
    end
    output = Dict{Thunk,Int}()
    recur([node], 0)
    return output
end


text:
# JUDI precon make sure it runs
model, model0, dm = setup_model()
q, src_geometry, rec_geometry, info = setup_geom(model)

ftol = 5f-5

nsrc = q.nsrc
w = judiWeights(randn(Float32, model0.n))
w_out = similar(w)

data = Array{Array}(undef, nsrc)
for j=1:nsrc
    data[j] = randn(Float32, rec_geometry.nt[j], length(rec_geometry.xloc[j]))
end

dobs = judiVector(rec_geometry, data)
dobs_out = similar(dobs)

F = judiFilter(rec_geometry, .002, .030)
Md = judiMarineTopmute2D(0, rec_geometry)
D = judiDepthScaling(model)
M = judiTopmute(model.n, 20, 1)

mul!(dobs_out, F, dobs)
@test isapprox(dobs_out, F*dobs; rtol=ftol)
mul!(dobs_out, F', dobs)
@test isapprox(dobs_out, F'*dobs; rtol=ftol)
mul!(dobs_out, Md, dobs)
@test isapprox(dobs_out, Md*dobs; rtol=ftol)
mul!(dobs_out, Md', dobs)
@test isapprox(dobs_out, Md'*dobs; rtol=ftol)

mul!(w_out, D, w)
@test isapprox(w_out, D*w; rtol=ftol)
@test isapprox(w_out.weights[1][end, end]/w.weights[1][end, end], sqrt((model.n[2]-1)*model.d[2]); rtol=ftol)
mul!(w_out, D', w)
@test isapprox(w_out, D'*w; rtol=ftol)
@test isapprox(w_out.weights[1][end, end]/w.weights[1][end, end], sqrt((model.n[2]-1)*model.d[2]); rtol=ftol)

mul!(w_out, M, w)
@test isapprox(w_out, M*w; rtol=ftol)
@test isapprox(norm(w_out.weights[1][:, 1:19]), 0f0; rtol=ftol)
mul!(w_out, M', w)
@test isapprox(w_out, M'*w; rtol=ftol)

# test the output of depth scaling and topmute operators, and test if they are out-of-place

dobs1 = deepcopy(dobs)
for Op in [F, F', Md , Md']
    m = Op*dobs
    # Test that dobs wasn't modified
    @test isapprox(dobs,dobs1,rtol=eps())
    # Test that it did compute something
    @test m != dobs
end

w1 = deepcopy(w)

for Op in [D, D']
    m = Op*w
    # Test that w wasn't modified
    @test isapprox(w,w1,rtol=eps())

    w_expect = deepcopy(w)
    for j = 1:model.n[2]
        w_expect.weights[1][:,j] = w.weights[1][:,j] * Float32(sqrt(model.d[2]*(j-1)))
    end
    @test isapprox(w_expect,m)
end

for Op in [M , M']
    m = Op*w
    # Test that w wasn't modified
    @test isapprox(w,w1,rtol=eps())

    @test all(isapprox.(m.weights[1][:,1:18], 0))
    @test isapprox(m.weights[1][:,21:end],w.weights[1][:,21:end])
end

n = (100,100,100)
d = (10f0,10f0,10f0)
o = (0.,0.,0.)
m = 0.25*ones(Float32,n)
model3D = Model(n,d,o,m)

D3 = judiDepthScaling(model3D)

dm = randn(Float32,prod(n))
dm1 = deepcopy(dm)
for Op in [D3, D3']
    opt_out = Op*dm
    # Test that dm wasn't modified
    @test dm1 == dm

    dm_expect = zeros(Float32,model3D.n)
    for j = 1:model3D.n[3]
        dm_expect[:,:,j] = reshape(dm,model3D.n)[:,:,j] * Float32(sqrt(model3D.d[3]*(j-1)))
    end
    @test isapprox(vec(dm_expect),opt_out)
end

M3 = judiTopmute(model3D.n, 20, 1)

for Op in [M3, M3']
    opt_out = Op*dm
    # Test that dm wasn't modified
    @test dm1 == dm

    @test all(isapprox.(reshape(opt_out,model3D.n)[:,:,1:18], 0))
    @test isapprox(reshape(opt_out,model3D.n)[:,:,21:end],reshape(dm,model3D.n)[:,:,21:end])
end

# test find_water_bottom in 3D

dm3D = zeros(Float32,10,10,10)
dm3D[:,:,4:end] .= 1f0
@test find_water_bottom(dm3D) == 4*ones(Integer,10,10)


text:
using Knet
import Gym
import AutoGrad
import Random
using Statistics

mutable struct History
    nS::Int
    nA::Int
    γ::Float64
    states::Vector{Float64}
    actions::Vector{Int}
    rewards::Vector{Float64}
end
History(nS, nA, γ) = History(nS, nA, γ, zeros(0),zeros(Int, 0),zeros(0))

function discount(rewards, γ)
    R = similar(rewards)
    R[end] = rewards[end]
    for k = length(rewards)-1:-1:1
        R[k] = γ * R[k+1] + rewards[k]
    end
    # return R
    return (R .- mean(R)) ./ (std(R) + 1e-10) #speeds up training a lot
end

function initweights(hidden, xsize=4, ysize=2)
    w = []
    x = xsize
    for y in [hidden...]
        push!(w, xavier(y, x))
        push!(w, zeros(y, 1))
        x = y
    end
    push!(w, xavier(ysize, x)) # Prob actions
    push!(w, zeros(ysize, 1))
    push!(w, xavier(1, x)) # Value function
    push!(w, zeros(1, 1))

    return w
end

function predict(w, x)
    for i = 1:2:length(w)-4
        x = relu.(w[i] * x .+ w[i+1])
    end
    prob_act = w[end-3] * x .+ w[end-2]
    value = w[end-1] * x .+ w[end]
    return prob_act, value
end

function sample_action(probs)
    @assert size(probs, 2) == 1
    cprobs = cumsum(probs, dims=1)
    sampled = cprobs .> rand()
    return mapslices(argmax, sampled, dims=1)[1]
end

function loss(w, history)
    nS, nA = history.nS, history.nA
    M = length(history.states)÷nS
    states = reshape(history.states, nS, M)
    R = discount(history.rewards, history.γ)

    p, V = predict(w, states)
    V = vec(V)
    A = R .- V   # advantage
    inds = history.actions + nA*(0:M-1)
    lp = logsoftmax(p, dims=1)[inds] # lp is a vector

    return -mean(lp .* AutoGrad.value(A)) + L2Reg(A)
end

L2Reg(x) = mean(x .* x)

function main(;
    hidden = [100], # width inner layers
    lr = 1e-2,
    γ = 0.99, #discount rate
    episodes = 500,
    render = true,
    seed = -1,
    infotime = 50)

    env = Gym.GymEnv("CartPole-v1")
    seed > 0 && (Random.seed!(seed); Gym.seed!(env, seed))
    nS, nA = 4, 2
    w = initweights(hidden, nS, nA)
    opt = [Adam(lr=lr) for _=1:length(w)]

    avgreward = 0
    for episode=1:episodes
        state = Gym.reset!(env)
        episode_rewards = 0
        history = History(nS, nA, γ)
        for t=1:10000
            p, V = predict(w, state)
            p = softmax(p, dims=1)
            action = sample_action(p)

            next_state, reward, done, info = Gym.step!(env, action-1)
            append!(history.states, state)
            push!(history.actions, action)
            push!(history.rewards, reward)
            state = next_state
            episode_rewards += reward

            episode % infotime == 0 && render && Gym.render(env)
            done && break
        end

        avgreward = 0.1 * episode_rewards + avgreward * 0.9
        if episode % infotime == 0
            println("(episode:$episode, avgreward:$avgreward)")
            Gym.close!(env)
        end

        dw = grad(loss)(w, history)
        update!(w, dw, opt)
    end

    return w
end


text:
using GLVisualize.ComposeBackend, Gadfly, Distributions
using Colors, GLVisualize

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end

const static_example = true

sds = [1, 1/2, 1/4, 1/8, 1/16, 1/32]
n = 10
ys = [mean(rand(Distributions.Normal(0, sd), n)) for sd in sds]
ymins = ys .- (1.96 * sds / sqrt(n))
ymaxs = ys .+ (1.96 * sds / sqrt(n))

p = plot(x=1:length(sds), y=ys, ymin=ymins, ymax=ymaxs,
     Geom.point, Geom.errorbar)

draw(composebackend, p)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# flowpipe interface
include("AbstractFlowpipe.jl")

# concrete implementations
include("Flowpipe.jl")
include("ShiftedFlowpipe.jl")
include("MappedFlowpipe.jl")
include("HybridFlowpipe.jl")
include("MixedFlowpipe.jl")
include("PartitionedFlowpipe.jl")


text:
mutable struct BioAFMResult
    Posteriors::Array{<:AbstractFloat,1}
    MostPlausiblePDB::Int64
    MostPlausibleQuate::Int64
end

BioAFMResult() = BioAFMResult([], 0, 0)

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
using EmpiricalModeDecomposition
using Base.Test

# Testing Script for EMD module
using Base.Test

samples = 50000
t = linspace(1,100,samples)
modei = [sin.(5*i.*t) for i in 5:-1:1]

modesi = [sum(modei[1:N]) for N in 2:5]

imfsi = [IMF(modesi[N-1], N=N) for N in 2:5]

for (n,imfs) in enumerate(imfsi)
    errs = [sum(abs.(imfs[:,N] .- modei[N]))/sum(abs.(modei[N])) for N in 1:(n+1)]
    println("Average Mode Recovery Error for $(n+1) Modes is: $errs")
end

println("""
    Assessing performance on simulated data:
        5 sin waves at freqs 25,20,15,10,5 Hz.
        Samples: $samples
        Noise: None
    """
)

for (n,imfs) in enumerate(imfsi)
    errs = [sum(abs.(imfs[:,N] .- modei[N]))/sum(abs.(modei[N])) for N in 1:(n+1)]
    println("Mode Recovery Error for $(n+1) Modes is:\n\tAverage $errs")

    edge = floor(Int,samples/5)
    errs = [sum(abs.(imfs[1:edge,N] .- modei[N][1:edge]))/sum(abs.(modei[N][1:edge])) for N in 1:(n+1)]
    println("\tEdge Effects: $errs")

    edge = floor(Int,samples/5)
    errs = [sum(abs.(imfs[edge:(end-edge),N] .- modei[N][edge:(end-edge)])
            )/sum(abs.(modei[N][edge:(end-edge)])) for N in 1:(n+1)]
    println("\tErrors without Edge Effects: $errs\n\n")
end

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition


text:
#Demo 1: Simple sine wave plot
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import NumericIO: UEXPONENT_SI

#==Input
===============================================================================#
#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
tcycle = 1e-9
npts = 2000 #Relatively large dataset
ncycles = 4
tmax = tcycle*ncycles
Δ = tmax/(npts-1)
t = collect(0:Δ:tmax)
y = sin.(t*(ncycles*2pi/tmax))

t[end] = t[1] #x-values are no longer ordered

#Lower resolution t-vector:
t_lres = collect(0:(tmax/10):tmax)

#==Generate plot
===============================================================================#
plot = InspectDR.transientplot(:lin, title="Sample Plot (π𝜋α𝛼β𝛽γ𝛾Ω𝛺℧ϕ𝜑𝜙ω𝜔f𝑓ϵ𝜀λ𝜆,⇒x∊ℝ)")
	InspectDR.overwritefont!(plot.layout, fontname="monospace", fontscale=1.5)
#TODO: change how layout[] works. Not obvious that you can't just set .bold property of a font
tf = plot.layout[:font_title]; tf.bold = false #Must read first
plot.layout[:font_title] = tf #Then assign to flag property changed
plot.layout[:enable_legend] = true
plot.layout[:halloc_legend] = 150
plot.layout[:enable_timestamp] = true
plot.layout[:format_xtick] = InspectDR.TickLabelStyle(UEXPONENT_SI) #Use SI notation on x-axis

graph = plot.strips[1]
graph.grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=true)

style = :dashdot #solid/dashdot/...
wfrm = add(plot, t, y .+1, id="sin(2πt)+1")
	wfrm.line = line(color=blue, width=1, style=style)
wfrm = add(plot, t, y .-1, id="sin(2πt)-1")
	wfrm.line = line(color=red, width=5, style=style)
wfrm = add(plot, t_lres, -2 .+4*(t_lres./tmax), id="-2+4t/tmax")
	wfrm.line = line(color=blue, width=3, style=:dash)
	wfrm.glyph = glyph(shape=:*, size=10)

a = plot.annotation
	a.xlabel = "Time (s)"
	a.ylabels = ["Signal Voltage (V)"]

let wfrm #HIDEWARN_0.7
#Show if uses f1 acceleration:
for wfrm in plot.data
	id = wfrm.id
	f1accel = isa(wfrm.ds, InspectDR.IDataset{true})
	@show id, f1accel
end
end

gplot = display(InspectDR.GtkDisplay(), plot)

InspectDR.write_png("export_plot.png", plot)
InspectDR.write_svg("export_plot.svg", plot)
InspectDR.write_eps("export_plot.eps", plot)
InspectDR.write_pdf("export_plot.pdf", plot)

:DONE


text:
function Fraction(arg0::jdouble)
    return Fraction((jdouble,), arg0)
end

function Fraction(arg0::jdouble, arg1::jdouble, arg2::jint)
    return Fraction((jdouble, jdouble, jint), arg0, arg1, arg2)
end

function Fraction(arg0::jdouble, arg1::jint)
    return Fraction((jdouble, jint), arg0, arg1)
end

function Fraction(arg0::jint)
    return Fraction((jint,), arg0)
end

function Fraction(arg0::jint, arg1::jint)
    return Fraction((jint, jint), arg0, arg1)
end

function abs(obj::Fraction)
    return jcall(obj, "abs", Fraction, ())
end

function add(obj::Fraction, arg0::Fraction)
    return jcall(obj, "add", Fraction, (Fraction,), arg0)
end

function add(obj::Fraction, arg0::jint)
    return jcall(obj, "add", Fraction, (jint,), arg0)
end

function compare_to(obj::Fraction, arg0::Fraction)
    return jcall(obj, "compareTo", jint, (Fraction,), arg0)
end

function divide(obj::Fraction, arg0::Fraction)
    return jcall(obj, "divide", Fraction, (Fraction,), arg0)
end

function divide(obj::Fraction, arg0::jint)
    return jcall(obj, "divide", Fraction, (jint,), arg0)
end

function double_value(obj::Fraction)
    return jcall(obj, "doubleValue", jdouble, ())
end

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

function float_value(obj::Fraction)
    return jcall(obj, "floatValue", jfloat, ())
end

function get_denominator(obj::Fraction)
    return jcall(obj, "getDenominator", jint, ())
end

function get_field(obj::Fraction)
    return jcall(obj, "getField", FractionField, ())
end

function get_numerator(obj::Fraction)
    return jcall(obj, "getNumerator", jint, ())
end

function get_reduced_fraction(::Type{Fraction}, arg0::jint, arg1::jint)
    return jcall(Fraction, "getReducedFraction", Fraction, (jint, jint), arg0, arg1)
end

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

function int_value(obj::Fraction)
    return jcall(obj, "intValue", jint, ())
end

function long_value(obj::Fraction)
    return jcall(obj, "longValue", jlong, ())
end

function multiply(obj::Fraction, arg0::Fraction)
    return jcall(obj, "multiply", Fraction, (Fraction,), arg0)
end

function multiply(obj::Fraction, arg0::jint)
    return jcall(obj, "multiply", Fraction, (jint,), arg0)
end

function negate(obj::Fraction)
    return jcall(obj, "negate", Fraction, ())
end

function percentage_value(obj::Fraction)
    return jcall(obj, "percentageValue", jdouble, ())
end

function reciprocal(obj::Fraction)
    return jcall(obj, "reciprocal", Fraction, ())
end

function signum(obj::Fraction)
    return jcall(obj, "signum", jint, ())
end

function subtract(obj::Fraction, arg0::Fraction)
    return jcall(obj, "subtract", Fraction, (Fraction,), arg0)
end

function subtract(obj::Fraction, arg0::jint)
    return jcall(obj, "subtract", Fraction, (jint,), arg0)
end

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


text:
using MinFEM

function solve_elasticity(file_name::String)

  mesh = import_mesh(file_name)

  L = asmElasticity(mesh,1.0,1.0)
  Mb = asmBoundaryMassMatrix(mesh, Edges=mesh.Boundaries[1001].Edges, qdim=2)

  f(x) = [x[1]; x[2]]
  s = evaluateMeshFunction(mesh, f, region=mesh.Boundaries[1001].Nodes, qdim=2)

  boundary = mesh.Boundaries[1003].Nodes

  pde = PDESystem(A=L, b=Mb*s, bc=zeros(2*mesh.nnodes), DI=boundary, qdim=2)
  solve(pde)

  return true

end

@test solve_elasticity("test_v4.msh")


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct SummaryDescription <: ProtoType
    type_hint::AbstractString
    SummaryDescription(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SummaryDescription

mutable struct HistogramProto <: ProtoType
    min::Float64
    max::Float64
    num::Float64
    sum::Float64
    sum_squares::Float64
    bucket_limit::Base.Vector{Float64}
    bucket::Base.Vector{Float64}
    HistogramProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HistogramProto
const __pack_HistogramProto = Symbol[:bucket_limit,:bucket]
meta(t::Type{HistogramProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_HistogramProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct SummaryMetadata_PluginData <: ProtoType
    plugin_name::AbstractString
    content::Array{UInt8,1}
    SummaryMetadata_PluginData(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SummaryMetadata_PluginData

mutable struct SummaryMetadata <: ProtoType
    plugin_data::SummaryMetadata_PluginData
    display_name::AbstractString
    summary_description::AbstractString
    SummaryMetadata(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SummaryMetadata

mutable struct Summary_Image <: ProtoType
    height::Int32
    width::Int32
    colorspace::Int32
    encoded_image_string::Array{UInt8,1}
    Summary_Image(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Summary_Image

mutable struct Summary_Audio <: ProtoType
    sample_rate::Float32
    num_channels::Int64
    length_frames::Int64
    encoded_audio_string::Array{UInt8,1}
    content_type::AbstractString
    Summary_Audio(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Summary_Audio

mutable struct Summary_Value <: ProtoType
    node_name::AbstractString
    tag::AbstractString
    metadata::SummaryMetadata
    simple_value::Float32
    obsolete_old_style_histogram::Array{UInt8,1}
    image::Summary_Image
    histo::HistogramProto
    audio::Summary_Audio
    tensor::TensorProto
    Summary_Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Summary_Value
const __fnum_Summary_Value = Int[7,1,9,2,3,4,5,6,8]
const __oneofs_Summary_Value = Int[0,0,0,1,1,1,1,1,1]
const __oneof_names_Summary_Value = [Symbol("value")]
meta(t::Type{Summary_Value}) = meta(t, ProtoBuf.DEF_REQ, __fnum_Summary_Value, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_Summary_Value, __oneof_names_Summary_Value, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Summary <: ProtoType
    value::Base.Vector{Summary_Value}
    Summary(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Summary

export SummaryDescription, HistogramProto, SummaryMetadata_PluginData, SummaryMetadata, Summary_Image, Summary_Audio, Summary_Value, Summary


text:
#
# TODO: fix normalization to use getter functions from MathematicalSystems
#
#export add_dimension,
#       next_set,
#       iscontinuoussystem,
#       ishybridsystem

import Base: *
import LazySets.constrained_dimensions

# continuous systems that are handled by this library
iscontinuoussystem(T::Type{<:AbstractSystem}) = false
iscontinuoussystem(T::Type{<:LCS}) = true
iscontinuoussystem(T::Type{<:CLCS}) = true
iscontinuoussystem(T::Type{<:CLCCS}) = true
iscontinuoussystem(T::Type{<:ACS}) = true
iscontinuoussystem(T::Type{<:CACCS}) = true
iscontinuoussystem(T::Type{<:CACS}) = true
iscontinuoussystem(T::Type{<:BBCS}) = true
iscontinuoussystem(T::Type{<:CBBCS}) = true
iscontinuoussystem(T::Type{<:CBBCCS}) = true
iscontinuoussystem(T::Type{<:SOLCS}) = true
iscontinuoussystem(T::Type{<:SOACS}) = true
iscontinuoussystem(T::Type{<:SOCLCCS}) = true
iscontinuoussystem(T::Type{<:SOCACCS}) = true

# hybrid systems that are handled by this library
ishybridsystem(T::Type{<:AbstractSystem}) = false
ishybridsystem(T::Type{<:HybridSystem}) = true

# systems that are second order
function is_second_order(::ST) where {ST<:AbstractContinuousSystem}
    ST <: SOLCS || ST <: SOACS || ST <: SOCLCCS || ST <: SOCACCS
end

#export LCS, LDS, CLCS, CLDS, CLCCS, CLCDS, CACCS, CACDS, CACS, CADS, IVP, BBCS,
#       CBBCS, CBBCCS

*(M::AbstractMatrix, input::ConstantInput) =  ConstantInput(M * input.U)

# convenience functions
next_set(inputs::ConstantInput) = collect(nextinput(inputs, 1))[1]
next_set(inputs::AbstractInput, state::Int64) = collect(nextinput(inputs, state))[1]

"""
    add_dimension(A::AbstractMatrix, m=1)

Append one or more zero rows and columns to a matrix.

### Input

- `A` -- matrix
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_test
julia> A = [0.4 0.25; 0.46 -0.67]
2×2 Array{Float64,2}:
 0.4    0.25
 0.46  -0.67

julia> add_dimension(A)
3×3 Array{Float64,2}:
 0.4    0.25  0.0
 0.46  -0.67  0.0
 0.0    0.0   0.0
```
To append more than one zero row-column, use the second argument `m`:

```jldoctest add_dimension_test
julia> add_dimension(A, 2)
4×4 Array{Float64,2}:
 0.4    0.25  0.0  0.0
 0.46  -0.67  0.0  0.0
 0.0    0.0   0.0  0.0
 0.0    0.0   0.0  0.0
```
"""
function add_dimension(A::AbstractMatrix, m=1)
    n = size(A, 1)
    return vcat(hcat(A, zeros(n, m)), zeros(m, n+m))
end

"""
    add_dimension(X::LazySet, m=1)

Adds an extra dimension to a LazySet through a Cartesian product.

### Input

- `X` -- a lazy set
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_set
julia> X = BallInf(ones(9), 0.5);

julia> dim(X)
9

julia> Xext = add_dimension(X);

julia> dim(Xext)
10

julia> X = ZeroSet(4);

julia> dim(add_dimension(X))
5

julia> typeof(X)
ZeroSet{Float64}
```

More than one dimension can be added passing the second argument:

```jldoctest add_dimension_set
julia> Xext = add_dimension(BallInf(zeros(10), 0.1), 4);

julia> dim(Xext)
14
```

### Notes

In the special case that the given set is a zero set, instead of cartesian product
a new zero set with extended dimensions is returned.
"""
function add_dimension(X::LazySet, m=1)
    return X * ZeroSet(m)
end

function add_dimension(X::ZeroSet, m=1)
    return ZeroSet(dim(X)+m)
end

"""
    add_dimension(cs, m=1)

Adds an extra dimension to a continuous system.

### Input

- `cs` -- continuous system
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_cont_sys
julia> using MathematicalSystems, SparseArrays

julia> A = sprandn(3, 3, 0.5);

julia> X0 = BallInf(zeros(3), 1.0);

julia> s = InitialValueProblem(LinearContinuousSystem(A), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4
```

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = ConstantInput(Ball2(ones(3), 0.1));

julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4

julia> dim(next_set(inputset(sext)))
4
```

Extending a system with a varying input set:

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = VaryingInput([Ball2(ones(3), 0.1 * i) for i in 1:3]);

julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4

julia> dim(next_set(inputset(sext), 1))
4
```

Extending a varing input set with more than one extra dimension:
1] normalize(::AffineContinuousSystem{Float64,Array{Float64,2},Array{Float64,1}}) at /home/mforets/.julia/dev/ReachabilityAnalysis/src/Continuous/normalization.jl:387
```jldoctest add_dimension_cont_sys
julia> sext = add_dimension(s, 7);

julia> statedim(sext)
10

julia> dim(next_set(inputset(sext), 1))
10
```
"""
function add_dimension(cs, m=1)
    Aext = add_dimension(cs.s.A, m)
    X0ext = add_dimension(cs.x0, m)
    if hasmethod(inputset, Tuple{typeof(cs.s)})
        Uext = map(x -> add_dimension(x, m), inputset(cs))
        s = CLCCS(Aext, Matrix(1.0I, size(Aext)), nothing, Uext)
    else
        s = LCS(Aext)
    end
    return IVP(s, X0ext)
end

# accepted types of non-deterministic inputs (non-canonical form)
const UNCF{N} = Union{<:LazySet{N}, Vector{<:LazySet{N}}, <:AbstractInput} where {N}

# accepted types for the state constraints (non-canonical form)
const XNCF{N} = Union{<:LazySet{N}, Nothing} where {N}

"""
    normalize(system::AbstractSystem)

Transform a mathematical system to a normalized (or canonical) form.

### Input

- `system` -- system; it can be discrete or continuous

### Output

Either the same system if it already conforms to a canonical form, or a new
system otherwise.

### Notes

The normalization procedure consists of transforming a given system type into a
"canonical" format that is used internally. More details are given below.

### Algorithm

The implementation of `normalize` exploits `MathematicalSystems`'s' types, which
carry information about the problem as a type parameter.

Homogeneous ODEs of the form ``x' = Ax, x ∈ \\mathcal{X}`` are canonical if the
associated problem is a `ConstrainedLinearContinuousSystem` and `A` is a matrix.
This type does not handle non-deterministic inputs.

Note that a `LinearContinuousSystem` does not consider constraints on the
state-space (such as an invariant); to specify state constraints, use a
`ConstrainedLinearContinuousSystem`. If the passed system is a `LinearContinuousSystem`
(i.e. no constraints) then the normalization fixes a universal set (`Universe`)
as the constraint set.

The generalization to canonical systems with constraints and possibly time-varying
non-deterministic inputs is considered next. These systems are of the form
``x' = Ax + u, u ∈ \\mathcal{U}, x ∈ \\mathcal{X}``. The system type is
`ConstrainedLinearControlContinuousSystem`, where `A` is a matrix, `X` is a set
and `U` is an input, that is, any concrete subtype of `AbstractInput`.

If `U` is not given as an input, normalization accepts either a `LazySet`, or
a vector of `LazySet`s. In these cases, the sets are wrapped around an appropriate
concrete input type.

If the system does not conform to a canonical form, the implementation tries
to make the transformation; otherwise an error is thrown. In particular, ODEs
of the form ``x' = Ax + Bu`` are mapped into ``x' = Ax + u, u ∈ B\\mathcal{U}``,
where now ``u`` has the same dimensions as ``x``.

The transformations described above are analogous in the discrete case, i.e.
``x_{k+1} = A x_k`` and ``x_{k+1} = Ax_{k} + u_k, u_k ∈ \\mathcal{U}, x_k ∈ \\mathcal{X}``
for the linear and affine cases respectively.
"""
function normalize(system::AbstractSystem)
    throw(ArgumentError("the system type $(typeof(system)) is currently not supported"))
end

"""
    normalize(ivp::InitialValueProblem)

Transform an initial-value problem into a normalized (or canonical) form.

### Input

- `ivp` -- initial-value problem

### Output

Either the same initial-value problem if it already conforms to a canonical form,
or a new one otherwise.

### Notes

This function extends [`normalize`](@ref normalize(::AbstractSystem)) for
initial-value problems.
"""
function normalize(ivp::InitialValueProblem)
    return IVP(normalize(system(ivp)), initial_state(ivp))
end

# "black-box" like systems are not normalized; algorithms should handle this
normalize(S::Union{BBCS, CBBCS, CBBCCS}) = S

# x' = Ax, in the continuous case
# x+ = Ax, in the discrete case
for (L_S, CL_S) in ((:LCS, :CLCS), (:LDS, :CLDS))
    @eval begin
        function normalize(system::$L_S{N, AN}) where {N, AN<:AbstractMatrix{N}}
            n = statedim(system)
            X = Universe(n)
            return $CL_S(state_matrix(system), X)
        end
    end
end

# x' = Ax, x ∈ X in the continuous case
# x+ = Ax, x ∈ X in the discrete case
for CL_S in (:CLCS, :CLDS)
    @eval begin
        function normalize(system::$CL_S{N, AN, XT}) where {N, AN<:AbstractMatrix{N}, M, XT<:XNCF{M}}
            n = statedim(system)
            X = _wrap_invariant(stateset(system), n)
            return $CL_S(state_matrix(system), X)
        end
    end
end

# x' = Ax + Bu, x ∈ X, u ∈ U in the continuous case
# x+ = Ax + Bu, x ∈ X, u ∈ U in the discrete case
for CLC_S in (:CLCCS, :CLCDS)
    @eval begin
        function normalize(system::$CLC_S{N, AN, BN, XT, UT}) where {N, AN<:AbstractMatrix{N}, NN, BN<:AbstractMatrix{NN}, M, XT<:XNCF{M}, UT<:UNCF{M}}
            n = statedim(system)
            X = _wrap_invariant(stateset(system), n)
            U = _wrap_inputs(inputset(system), input_matrix(system))
            $CLC_S(state_matrix(system), I(n, N), X, U)
        end
    end
end

# x' = Ax + b, x ∈ X, u ∈ U in the continuous case
# x+ = Ax + b, x ∈ X, u ∈ U in the discrete case
for (CA_S, CLC_S) in ((:CACS, :CLCCS), (:CADS, :CLCDS))
    @eval begin
        function normalize(system::$CA_S{N, AN, VN, XT}) where {N, AN<:AbstractMatrix{N}, VN<:AbstractVector{N}, XT<:XNCF{N}}
            n = statedim(system)
            X = _wrap_invariant(stateset(system), n)
            U = _wrap_inputs(affine_term(system))
            $CLC_S(state_matrix(system), I(n, N), X, U)
        end
    end
end

# x' = Ax + b in the continuous case
# x+ = Ax + b in the discrete case
for (A_S, CLC_S) in ((:ACS, :CLCCS), (:ADS, :CLCDS))
    @eval begin
        function normalize(system::$A_S{N, AN, VN}) where {N, AN<:AbstractMatrix{N}, VN<:AbstractVector{N}}
            n = statedim(system)
            X = Universe(n)
            U = _wrap_inputs(affine_term(system))
            $CLC_S(state_matrix(system), I(n, N), X, U)
        end
    end
end

# x' = Ax + Bu + c, x ∈ X, u ∈ U in the continuous case
# x+ = Ax + Bu + c, x ∈ X, u ∈ U in the discrete case
for (CAC_S, CLC_S) in ((:CACCS, :CLCCS), (:CACDS, :CLCDS))
    @eval begin
        function normalize(system::$CAC_S{N, AN, BN, VN, XT, UT}) where {N, AN<:AbstractMatrix{N}, BN<:AbstractMatrix{N}, VN<:AbstractVector{N}, XT<:XNCF{N}, UT<:UNCF{N}}
            n = statedim(system)
            X = _wrap_invariant(stateset(system), n)
            U = _wrap_inputs(inputset(system), input_matrix(system), affine_term(system))
            $CLC_S(state_matrix(system), I(n, N), X, U)
        end
    end
end

# fix type inference
function _normalize(ivp::IVP{LCS{N, IdentityMultiple{N}}, Interval{N, IA.Interval{N}}}) where {N}
    return IVP(CLCS(ivp.s.A, Universe(1)), ivp.x0)
end

# ===========================
# Second order systems
# ===========================

# we left the n-dimensional second order system to a 2n-dimensional first
# order system
# if derivatives_last = true (Default)
# we assum that the first n variables correspond to position and the last n
# to velocities as in x̃ = [x, x']
#
# otherwise, if derivatives_last = false
# we assum that the first n variables correspond to velocities and the last n
# to position as in x̃ = [x', x]
function _second_order_linear_matrix(M::AbstractMatrix{N}, C, K; derivatives_last=true) where {N}
    n = size(M, 1)
    M⁻¹ = inv(Matrix(M))
    Idn = Matrix(one(N)*I, n, n)
    Zn = zeros(N, n, n)

    if derivatives_last
        A = [Zn       Idn   ;
             -M⁻¹*K   -M⁻¹*C]
    else
        A = [-M⁻¹*C   -M⁻¹*K;
             Idn          Zn]
    end
    return A, M⁻¹
end

function normalize(system::SOLCS{N}; derivatives_last=true) where {N}
    n = statedim(system)
    M = mass_matrix(system)
    C = viscosity_matrix(system)
    K = stiffness_matrix(system)
    A, _ = _second_order_linear_matrix(M, C, K; derivatives_last=derivatives_last)
    return normalize(LCS(A))
end

function normalize(system::SOACS{N}; derivatives_last=true) where {N}
    n = statedim(system)
    M = mass_matrix(system)
    C = viscosity_matrix(system)
    K = stiffness_matrix(system)

    A, M⁻¹ = _second_order_linear_matrix(M, C, K; derivatives_last=derivatives_last)

    b = affine_term(system)
    if derivatives_last
        c = vcat(zeros(N, n), M⁻¹*b)
    else
        c = vcat(M⁻¹*b, zeros(N, n))
    end

    return normalize(ACS(A, c))
end

function normalize(system::SOCLCCS{N}; derivatives_last=true) where {N}
    n = statedim(system)
    M = mass_matrix(system)
    C = viscosity_matrix(system)
    K = stiffness_matrix(system)
    B = input_matrix(system)
    X = stateset(system)
    U = inputset(system)

    A, M⁻¹ = _second_order_linear_matrix(M, C, K; derivatives_last=derivatives_last)
    if derivatives_last
        B̃ = vcat(zeros(N, n), M⁻¹*B)
    else
        B̃ = vcat(M⁻¹*B, zeros(N, n))
    end
    return normalize(CLCCS(A, B̃, X, U))
end

function normalize(system::SOCACCS{N}; derivatives_last=true) where {N}
    n = statedim(system)
    M = mass_matrix(system)
    C = viscosity_matrix(system)
    K = stiffness_matrix(system)
    B = input_matrix(system)
    d = affine_term(system)
    X = stateset(system)
    U = inputset(system)

    A, M⁻¹ = _second_order_linear_matrix(M, C, K; derivatives_last=derivatives_last)

    if derivatives_last
        B̃ = vcat(zeros(N, n), M⁻¹*B)
        d̃ = vcat(zeros(N, n), M⁻¹*d)
    else
        B̃ = vcat(M⁻¹*B, zeros(N, n))
        d̃ = vcat(M⁻¹*d, zeros(N, n))
    end
    return normalize(CACCS(A, B̃, X, U, d̃))
end

# ---

@inline isidentity(B::IdentityMultiple) = B.M.λ == oneunit(B.M.λ)

_wrap_invariant(X::LazySet, n::Int) = X
_wrap_invariant(X::Nothing, n::Int) = Universe(n)

_wrap_inputs(U::AbstractInput, B::IdentityMultiple) = isidentity(B) ? U : map(u -> B*u, U)
_wrap_inputs(U::AbstractInput, B::AbstractMatrix) = map(u -> B*u, U)

_wrap_inputs(U::LazySet, B::IdentityMultiple) = isidentity(B) ? ConstantInput(U) : ConstantInput(B*U)
_wrap_inputs(U::LazySet, B::AbstractMatrix) = ConstantInput(B*U)

_wrap_inputs(U::Vector{<:LazySet}, B::IdentityMultiple) = isidentity(B) ? VaryingInput(U) : VaryingInput(map(u -> B*u, U))
_wrap_inputs(U::Vector{<:LazySet}, B::AbstractMatrix) = VaryingInput(map(u -> B*u, U))

_wrap_inputs(U::AbstractInput, B::IdentityMultiple, c::AbstractVector) = isidentity(B) ? map(u -> u ⊕ c, U) : map(u -> B*u ⊕ c, U)
_wrap_inputs(U::AbstractInput, B::AbstractMatrix, c::AbstractVector) = map(u -> B*u ⊕ c, U)

_wrap_inputs(U::LazySet, B::IdentityMultiple, c::AbstractVector) = isidentity(B) ? ConstantInput(U ⊕ c) : ConstantInput(B*U ⊕ c)
_wrap_inputs(U::LazySet, B::AbstractMatrix, c::AbstractVector) = ConstantInput(B*U ⊕ c)

_wrap_inputs(U::Vector{<:LazySet}, B::IdentityMultiple, c::AbstractVector) = isidentity(B) ? VaryingInput(map(u -> u ⊕ c, U)) : VaryingInput(map(u -> B*u ⊕ c, U))
_wrap_inputs(U::Vector{<:LazySet}, B::AbstractMatrix, c::AbstractVector) = VaryingInput(map(u -> B*u ⊕ c, U))

_wrap_inputs(c::AbstractVector) = ConstantInput(Singleton(c))

# special cases for interval systems
function MathematicalSystems.ConstrainedLinearControlContinuousSystem(A::Matrix{IT}, B::SparseMatrixCSC, X::XT, U::UT) where {XT, UT, IT<:IA.Interval}
    CLCCS(IntervalMatrix(A), IntervalMatrix(B), X, U)
end
function MathematicalSystems.ConstrainedLinearControlContinuousSystem(A::Matrix{IT}, B::Matrix{IT}, X::XT, U::UT) where {XT, UT, IT<:IA.Interval}
    CLCCS(IntervalMatrix(A), IntervalMatrix(B), X, U)
end
function MathematicalSystems.ConstrainedLinearControlDiscreteSystem(A::IntervalMatrix, B::Matrix, X::XT, U::UT) where {XT, UT}
    CLCDS(A, IntervalMatrix(B), X, U)
end

# ==========================================================
# Shared functionality for linear continuous post operators
# ==========================================================

Base.:*(im::IdentityMultiple, d::Diagonal) = im.M.λ * d  # TODO: add to MathematicalSystems

hasinput(S::AbstractSystem) = inputdim(S) > 0
isconstantinput(::ConstantInput) = true
isconstantinput(::VaryingInput) = false
isconstantinput(::LazySet) = true

# The canonical form is:
#     - If the system doesn't have input, a constrained linear continous system (CLCS)
#       x' = Ax, x ∈ X
#     - If the system has an input, a CLCCS, x' = Ax + u, x ∈ X, u ∈ U
# If the original system is unconstrained, the constraint set X is the universal set.
abstract type AbstractLinearContinuousSystem <: AbstractContinuousSystem end
abstract type AbstractNonlinearContinuousSystem <: AbstractContinuousSystem end

function _normalize(ivp::IVP{<:AbstractContinuousSystem})
    if islinear(ivp) || isaffine(ivp)
        return _normalize(ivp, AbstractLinearContinuousSystem)
    else
        return _normalize(ivp, AbstractNonlinearContinuousSystem)
    end
end

function _normalize(ivp::IVP{<:BBCS}, ::Type{AbstractNonlinearContinuousSystem})
    f = ivp.s.f # TODO add getter function in MathematicalSystems
    n = statedim(ivp)
    X = Universe(n)
    return IVP(CBBCS(f, n, X), initial_state(ivp))
end

function _normalize(ivp::IVP{<:CBBCS}, ::Type{AbstractNonlinearContinuousSystem})
    return ivp
end

function _normalize(ivp::IVP{<:AbstractContinuousSystem}, ::Type{AbstractNonlinearContinuousSystem})
    throw(ArgumentError("can't normalize this nonlinear initial-value problemof type $(typeof(ivp))"))
end

_normalize_initial_state(ivp, X0) = X0
_normalize_initial_state(ivp, X0::AbstractVector) = Singleton(X0)
_normalize_initial_state(ivp, X0::IA.Interval) = convert(Interval, X0)
_normalize_initial_state(ivp, X0::IA.IntervalBox) = convert(Hyperrectangle, X0)
_normalize_initial_state(ivp, X0::Number) = Singleton([X0])
_normalize_initial_state(ivp::IVP{<:SecondOrderSystem}, X0::Tuple{<:AbstractVector, <:AbstractVector}) = Singleton(vcat(X0[1], X0[2]))
_normalize_initial_state(ivp::IVP{<:SecondOrderSystem}, X0::Tuple{<:LazySet, <:LazySet}) = X0[1] × X0[2]

const CanonicalLinearContinuousSystem = Union{CLCS, CLCCS}

function _normalize(ivp::IVP{<:AbstractContinuousSystem}, ::Type{AbstractLinearContinuousSystem})

    # initial states normalization
    X0 = initial_state(ivp)
    X0_norm = _normalize_initial_state(ivp, X0)

    # system's normalization
    S = system(ivp)
    S_norm = normalize(S)

    if S_norm === S && X0_norm === X0
        ivp_norm = ivp
    else
        ivp_norm = IVP(S_norm, X0_norm)
    end
    return ivp_norm
end


text:
@system RespirationTracker(Q10Function) begin
    w: weather ~ ::Weather(override)
    T(w.T_air): temperature ~ track(u"°C") # should be soil temperature
    To: optimal_temperature => 20 ~ preserve(u"°C", parameter)
    Q10 => begin
        # typical Q10 value for respiration, Loomis and Amthor (1999) Crop Sci 39:1584-1596
        2
    end ~ preserve(parameter)
end

@system Carbon begin
    weather ~ hold
    pheno: phenology ~ hold

    C_to_CH2O_ratio ~ hold
    seed_mass_export_rate ~ hold
    assimilation ~ hold
    total_mass ~ hold
    green_leaf_ratio ~ hold

    dp(pheno.development_phase): development_phase ~ track::sym

    C_conc: carbon_concentration => begin
        # maize: 40% C, See Kim et al. (2007) EEB
        0.45
    end ~ preserve

    carbon_reserve_from_seed(seed_mass_export_rate, C_conc, C_to_CH2O_ratio) => begin
        seed_mass_export_rate * C_conc * C_to_CH2O_ratio
    end ~ track(u"g/d")

    #TODO: take account NSC from bulb
    #TODO: handle translocation back from carbon pool (i.e. at midnight)
    carbon_reserve(carbon_reserve_from_seed, carbon_translocation) => begin
        carbon_reserve_from_seed - carbon_translocation
    end ~ accumulate(u"g")

    carbon_translocation(carbon_reserve, carbon_translocation_rate) => begin
        carbon_reserve * carbon_translocation_rate
    end ~ track(u"g/d")

    carbon_pool(assimilation, carbon_translocation, carbon_supply) => begin
        assimilation + carbon_translocation - carbon_supply
    end ~ accumulate(u"g")

    nonstructural_carbon(carbon_pool, carbon_reserve) => begin
        carbon_pool + carbon_reserve
    end ~ track(u"g")

    nonstructural_carbon_mass(nonstructural_carbon, C_conc, CH2O_to_C_ratio) => begin
        nonstructural_carbon * CH2O_to_C_ratio / C_conc
    end ~ track(u"g")

    carbon_supply(carbon_pool, carbon_supply_rate) => begin
        carbon_pool * carbon_supply_rate
    end ~ track(u"g/d")

    # to be used by allocate_carbon()
    carbon_temperature_effect(T=nounit(weather.T_air, u"°C"), β=pheno.BF.ΔT) => begin
        #FIXME properly handle T_air
        # this needs to be f of temperature, source/sink relations, nitrogen, and probably water
        # a valve function is necessary because assimilates from CPool cannot be dumped instantaneously to parts
        # this may be used for implementing feedback inhibition due to high sugar content in the leaves
        # The following is based on Grant (1989) AJ 81:563-571

        # Normalized (0 to 1) temperature response fn parameters, Pasian and Lieth (1990)
        # Lieth and Pasian Scientifica Hortuculturae 46:109-128 1991
        # parameters were fit using rose data -
        b1 = 2.325152587
        b2 = 0.185418876 # I'm using this because it can have broad optimal region unlike beta fn or Arrhenius eqn
        b3 = 0.203535650
        Td = 48.6 #High temperature compensation point

        g1 = 1 + exp(b1 - b2 * T)
        g2 = 1 - exp(-b3 * max(0, Td - T))
        #return g2 / g1

        β
    end ~ track

    carbon_growth_factor => begin
        # translocation limitation and lag, assume it takes 1 hours to complete, 0.2=5 hrs
        # this is where source/sink (supply/demand) valve can come in to play
        # 0.2 is value for hourly interval, Grant (1989)
        1 / 5u"hr"
    end ~ preserve(u"hr^-1")

    carbon_translocating(carbon_pool) => begin
        carbon_pool < 0u"g"
    end ~ flag

    carbon_translocation_rate(carbon_temperature_effect, carbon_growth_factor) => begin
        # C_demand does not enter into equations until grain fill
        carbon_temperature_effect * carbon_growth_factor
    end ~ track(u"hr^-1", when=carbon_translocating)

    carbon_supply_rate(carbon_temperature_effect, carbon_growth_factor) => begin
        carbon_temperature_effect * carbon_growth_factor
    end ~ track(u"hr^-1")

    Rm: maintenance_respiration_coefficient => begin
        # gCH2O g-1DM day-1 at 20C for young plants, Goudriaan and van Laar (1994) Wageningen textbook p 54, 60-61
        #0.015
        #0.018 # for maize
        0.012
    end ~ preserve(u"g/g/d", parameter)

    agefn(green_leaf_ratio): carbon_age_effect => begin
        # as more leaves senesce maint cost should go down, added 1 to both denom and numer to avoid division by zero.
        #agefn = (self.p.area.green_leaf + 1) / (self.p.area.leaf + 1)
        # no maint cost for dead materials but needs to be more mechanistic, SK
        #agefn = 1.0
        # from garlic model
        #agefn = (self.p.area.green_leaf + 0.1) / (self.p.area.leaf + 0.1)
        green_leaf_ratio
    end ~ track

    # based on McCree's paradigm, See McCree(1988), Amthor (2000), Goudriaan and van Laar (1994)
    # units very important here, be explicit whether dealing with gC, gCH2O, or gCO2
    maintenance_respiration_tracker(context, weather) ~ ::RespirationTracker
    maintenance_respiration(total_mass, Rm, agefn, q=maintenance_respiration_tracker.ΔT) => begin
        total_mass * q * Rm * agefn # gCH2O dt-1, agefn effect removed. 11/17/14. SK.
    end ~ track(u"g/d")

    carbon_available(carbon_supply, maintenance_respiration) => begin
        carbon_supply - maintenance_respiration
    end ~ track(u"g/d")

    Yg: synthesis_efficiency => begin
        #1 / 1.43 # equivalent Yg, Goudriaan and van Laar (1994)
        #0.75 # synthesis efficiency, ranges between 0.7 to 0.76 for corn, see Loomis and Amthor (1999), Grant (1989), McCree (1988)
        #0.74
        0.8
    end ~ preserve(parameter)

    total_carbon(c=carbon_available, Yg) => begin
        Yg * c # gCH2O partitioned to shoot/root
    end ~ track(u"g/d")

    # pt: partitioning_table => [
    #     # root leaf sheath scape bulb
    #       0.00 0.00   0.00  0.00 0.00 ; # seed garlic before germination
    #       0.35 0.30   0.25  0.00 0.10 ; # vegetative stage between germination and scape initiation
    #       0.15 0.15   0.10  0.25 0.35 ; # period between scape initiation and scape appearance
    #       0.05 0.10   0.00  0.35 0.50 ; # period after scape appearance before removal (scape stays intact)
    #       0.05 0.00   0.00  0.00 0.95 ; # period after scape removal (scape appeared and subsequently removed)
    #       0.00 0.00   0.00  0.00 0.00 ; # dead
    # ] ~ tabulate(
    #     rows=(:seed, :vegetative, :bulb_growth_before_scape_appearance, :bulb_growth_after_scape_appearance, :bulb_growth_after_scape_removal, :dead),
    #     columns=(:root, :leaf, :sheath, :scape, :bulb),
    #     parameter
    # )

    #Jim's lettuce revision
    pt: partitioning_table => [
        # root leaf stem
          0.00 0.00 0.00 ; # seed garlic before germination
          0.35 0.30 0.25 ; # vegetative stage
          0.00 0.00 0.00 ; # dead
    ] ~ tabulate(
        rows=(:seed, :vegetative, :dead),
        columns=(:root, :leaf, :stem),
        parameter
    )
    #FIXME: stem not incorporated

    root_carbon(total_carbon, pt, dp) => begin
        total_carbon * pt[dp].root
    end ~ track(u"g/d")

    shoot_carbon(total_carbon, root_carbon) => begin
        total_carbon - root_carbon
    end ~ track(u"g/d")

    leaf_carbon(total_carbon, pt, dp) => begin
        total_carbon * pt[dp].leaf
    end ~ track(u"g/d")

    stem_carbon(total_carbon, pt, dp) => begin
        total_carbon * pt[dp].stem
    end ~ track(u"g/d")

    # sheath_carbon(total_carbon, pt, dp) => begin
    #     total_carbon * pt[dp].sheath
    # end ~ track(u"g/d")

    # scape_carbon(total_carbon, pt, dp) => begin
    #     total_carbon * pt[dp].scape
    # end ~ track(u"g/d")

    # bulb_carbon(total_carbon, pt, dp) => begin
    #     total_carbon * pt[dp].bulb
    # end ~ track(u"g/d")
end


text:
@recipe function f(band::WannierBand, data=:eigvals; ks=nothing, fermi=0, linewidth=2)
    if ks == :relative
        ks = []
        k_m = band.kpoints_cryst[div(size(band.kpoints)[1] + 1, 2)]
        for k in band.kpoints_cryst
            push!(ks, norm(k - k_m))
        end
        ks[1:div(length(ks), 2)] = -ks[1:div(length(ks), 2)]
    else
        ks = collect(1:length(band.kpoints_cryst))
    end
    eigvals = band.eigvals .- fermi
    linewidth --> linewidth
    out = []
    if data==:eigvals
        out =eigvals
        title --> "Eigenvalues"
        yguide -->"energy (eV)"
    else
        if data==:cm_x
            title --> "Center of Mass (X)"
            yguide --> "Rx (Angström)"
            for cm in band.cms
                push!(out,cm[1])
            end
        elseif data==:cm_y
            title --> "Center of Mass (Y)"
            yguide --> "Ry (Angström)"
            for cm in band.cms
                push!(out,cm[2])
            end
        elseif data==:cm_z
            title --> "Center of Mass (Z)"
            yguide --> L"R_z"*" ("*L"\AA"*")"
            for cm in band.cms
                push!(out,cm[3])
            end
        elseif data==:angmom_x
            title --> "Orbital Angular Momentum (X)"
            yguide --> L"L_x (arb. units)"
            for angmom in band.angmoms
                push!(out,angmom[1][1]+angmom[2][1])
            end
        elseif data==:angmom_y
            title --> "Orbital Angular Momentum (Y)"
            yguide --> L"L_y (arb. units)"
            for angmom in band.angmoms
                push!(out,angmom[1][2]+angmom[2][2])
            end
        elseif data==:angmom_z
            title --> "Orbital Angular Momentum (Z)"
            yguide --> "Lz (arb. units)"
            for angmom in band.angmoms
                push!(out,angmom[1][3]+angmom[2][3])
            end
        elseif data==:angmom1_x
            title --> "Orbital Angular Momentum (X)"
            yguide --> "Lx (arb. units)"
            label --> "OAM around Ge"
            for angmom in band.angmoms
                push!(out,angmom[1][1])
            end
        elseif data==:angmom1_y
            title --> "Orbital Angular Momentum (Y)"
            yguide --> "Ly (arb. units)"
            label --> "OAM around Ge"
            for angmom in band.angmoms
                push!(out,angmom[1][2])
            end
        elseif data==:angmom1_z
            title --> "Orbital Angular Momentum (Z)"
            yguide --> "Lz (arb. units)"
            label --> "OAM around Ge"
            for angmom in band.angmoms
                push!(out,angmom[1][3])
            end
        elseif data==:angmom2_x
            title --> "Orbital Angular Momentum (X)"
            yguide --> L"L_x"*" (arb. units)"
            label --> "OAM around Te"
            for angmom in band.angmoms
                push!(out,angmom[2][1])
            end
        elseif data==:angmom2_y
            title --> "Orbital Angular Momentum (Y)"
            yguide --> L"L_y"*" (arb. units)"
            label --> "OAM around Ge"
            for angmom in band.angmoms
                push!(out,angmom[2][2])
            end
        elseif data==:angmom2_z
            title --> "Orbital Angular Momentum (Z)"
            yguide --> "Lz (arb. units)"
            label --> "OAM around Ge"
            for angmom in band.angmoms
                push!(out,angmom[2][3])
            end
        elseif data==:angmom2_xy
            title --> "Orbital Angular Momentum (XY)"
            yguide --> "norm(Lx)+norm(Ly)"
            label --> "Total angmom"
            for (spin,angmom) in zip(band.spins,band.angmoms)
                push!(out,sqrt((angmom[2][1]+spin[2][1])^2+(angmom[2][2]+spin[2][2])^2))
            end
        elseif data==:spin1_x
            title --> "Spin Angular Momentum (X)"
            yguide --> "Sx (arb. units)"
            label --> "SAM around Ge"
            for spin in band.spins
                push!(out,spin[1][1])
            end
        elseif data==:spin1_y
            title --> "Spin Angular Momentum (Y)"
            yguide --> "Sy (arb. units)"
            label --> "SAM around Ge"
            for spin in band.spins
                push!(out,spin[1][2])
            end
        elseif data==:spin1_z
            title --> "Spin Angular Momentum (Z)"
            yguide --> "Sz (arb. units)"
            label --> "SAM around Ge"
            for spin in band.spins
                push!(out,spin[1][3])
            end
        elseif data==:spin2_x
            title --> "Spin Angular Momentum (X)"
            yguide --> L"S_x"*" (arb. units)"
            label --> "SAM around Te"
            for spin in band.spins
                push!(out,spin[2][1])
            end
        elseif data==:spin2_y
            title --> "Spin Angular Momentum (Y)"
            yguide --> L"S_y"*" (arb. units)"
            label --> "SAM around Te"
            for spin in band.spins
                push!(out,spin[2][2])
            end
        elseif data==:spin2_z
            title --> "Spin Angular Momentum (Z)"
            yguide --> "Sz (arb. units)"
            label --> "SAM around Te"
            for spin in band.spins
                push!(out,spin[2][3])
            end
        elseif data==:epot
            title --> "Electrostatic Potential"
            yguide --> "E (arb. units)"
            for epot in band.epots
                push!(out,epot)
            end
        end
    end
    legend --> false
    ks,out
end

@recipe function f(bands::Array{<:WannierBand,1},data::Array{Symbol,1})
    layout := length(data)
    for (i,dat) in enumerate(data)
        @series begin
            subplot := i
            bands,data[i]
        end
    end
end

@recipe function f(bands::Array{<:WannierBand,1},data::Symbol)
    for band in bands
        @series begin
            band,data
        end
    end
end

@recipe function f(bands::Array{<:WannierBand,1})
    for band in bands
        @series begin
            band
        end
    end
end

@recipe function f(dfbands::Array{DFC.Band,1},WannierBands::Array{<:WannierBand,1})
    @series begin
        label--> "DFT Calculation"
        line--> (1,1.0,:blue)
        dfbands[1]
    end
    @series begin
        label--> "Wannier Interpolation"
        line--> (2,:dot,1.0,:red)
        WannierBands[1],:eigvals
    end
    for band in dfbands[2:end]
        @series begin
            label--> ""
            line--> (1,1.0,:blue)
            band
        end
    end
    for band in WannierBands[2:end]
        @series begin
            label--> ""
            line--> (2,:dot,1.0,:red)
            band,:eigvals
        end
    end
end

@recipe function f(WannierBands::Array{<:WannierBand,1},dfbands::Array{DFC.Band,1})
    @series begin
        dfbands,WannierBands
    end
end


text:
"""
This module is for exploration of Julia, Documenter and other interesting concepts.
"""
module Learn

include("DateUtils/DayCount.jl")
include("DateUtils/Calendars.jl")
include("DateUtils/BusinessDayAdjustment.jl")
include("stuff.jl")

end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
export problem1
"""
    problem1(N=1000)

# Multiples of 3 and 5
## Problem 1

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.
"""
function problem1(N=1000)
    res = 0
    for i in 1:(N-1)
        (mod(i,3) == 0 || mod(i,5) == 0) && (res += i; true)
    end
    return res
end

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler