text:
module ContGridWorlds

# Write your package code here.

# module FourRooms

import Random, Base.size
import Reexport: @reexport

@reexport using MinimalRLCore

using RecipesBase

"""
    Simple reward structure for goals. Can be duck typed.
"""
struct BasicRewFunc
    value::Float64
end

is_terminal(rew::BasicRewFunc) = true
get_reward(rew::BasicRewFunc) = rew.value

module FourRoomsContParams

const BASE_WALLS = [0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    1 0 1 1 1 1 0 0 0 0 0;
                    0 0 0 0 0 1 1 1 0 1 1;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;]

const GOAL_LOCS =  [0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;]

const ROOM_TOP_LEFT = 1
const ROOM_TOP_RIGHT = 2
const ROOM_BOTTOM_LEFT = 3
const ROOM_BOTTOM_RIGHT = 4

const REW_FUNCS = Dict{Int, Function}()

end

module TMazeContParams

import ..BasicRewFunc

const BASE_WALLS = [1 1 1 1 1 1 1 1 1 1 1;
                    1 0 0 1 1 1 1 1 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 1 0 0 0 1 0 0 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;]

const GOAL_LOCS = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1;
                   -1  1  1 -1 -1 -1 -1 -1  3  3 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  2  2 -1  0  0  0 -1  4  4 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;]

const REW_FUNCS = Dict([i=>BasicRewFunc(i) for i in 1:4]...)

end

include("dynamics.jl")

FourRooms(max_action_noise=0.1,
          drift_noise=0.001;
          normalized=false) =
              ContGridWorld(FourRoomsContParams.BASE_WALLS[end:-1:1, :],
                            FourRoomsContParams.GOAL_LOCS[end:-1:1, :],
                            FourRoomsContParams.REW_FUNCS,
                            max_action_noise,
                            drift_noise,
                            normalized)

TMaze(max_action_noise=0.1, drift_noise=0.001; normalized=false) =
    ContGridWorld(TMazeContParams.BASE_WALLS[end:-1:1, :],
                  TMazeContParams.GOAL_LOCS[end:-1:1, :],
                  TMazeContParams.REW_FUNCS,
                  max_action_noise,
                  drift_noise,
                  normalized)

module PlotParams

using Colors

const SIZE = 10
const BG = Colors.RGB(1.0, 1.0, 1.0)
const WALL = Colors.RGB(0.3, 0.3, 0.3)
const AC = Colors.RGB(0.69921875, 0.10546875, 0.10546875)
const GOAL = Colors.RGB(0.796875, 0.984375, 0.76953125)
const AGENT = [AC AC AC AC;
               AC AC AC AC;
               AC AC AC AC;
               AC AC AC AC;]

end

@recipe function f(env::ContGridWorld)
    ticks := nothing
    foreground_color_border := nothing
    grid := false
    legend := false
    aspect_ratio := 1
    xaxis := false
    yaxis := false

    PP = PlotParams

    s = size(env)
    screen = fill(PP.BG, (s[2] + 2)*PlotParams.SIZE, (s[1] + 2)*PlotParams.SIZE)
    screen[:, 1:PP.SIZE] .= PP.WALL
    screen[1:PP.SIZE, :] .= PP.WALL
    screen[end-(PP.SIZE-1):end, :] .= PP.WALL
    screen[:, end-(PP.SIZE-1):end] .= PP.WALL

    for i ∈ 1:s[1]
        for j ∈ 1:s[2]
            sqr_i = ((i)*PP.SIZE + 1):((i+1)*PP.SIZE)
            sqr_j = ((j)*PP.SIZE + 1):((j+1)*PP.SIZE)
            if env.walls[s[1] - i + 1, j]
                screen[sqr_i, sqr_j] .= PP.WALL
            end
            if env.goals[s[1] - i + 1, j] > 0
                screen[sqr_i, sqr_j] .= PP.GOAL
            end
        end
    end
    state = ((env.state[1] - 1, env.state[2] + 1).*PlotParams.SIZE)

    screen[collect(-1:2) .+ Int(floor(PP.SIZE*s[1] - state[1])), collect(-1:2) .+ Int(floor(state[2]))] .= PP.AGENT

    screen
end

end


text:
module Budgets

using StringBuilders, Formatting, Dates

import Base:*, +, -, print

export Percentile, Money, Event, Expense, Deduction, Payment, Budget, add!, expenses, deductions, payments, taxable, tax, nontaxable, balance, define_event

struct Percentile{T<:Real}
    p::T
end

struct Money{T<:Real}
    value::T
end

*(m::Money, x::Int) = Money(m.value*x)
*(p::Percentile, x::Money) = Money(p.p/100*x.value)
*(x::Money, p::Percentile) = p*x
+(m::Money, x::Money) = Money(m.value + x.value)

abstract type Event end

struct Expense <: Event
    date::Union{Missing, Date}
    description::String
    units::Int
    rate::Money
    sum::Money
    function Expense(date, description, units, rate::Money)
        @assert units > 0 "units must be positive"
        @assert rate.value > 0 "expense must be positive"
        new(date, description, units, rate, rate*units)
    end
end

struct Deduction <: Event
    description::String
    sum::Money
    function Deduction(description, sum)
        @assert sum.value < 0 "deduction must be negative"
        new(description, sum)
    end
end

struct Payment <: Event
    date::Date
    sum::Money
    function Payment(date, sum)
        @assert sum.value < 0 "payment must be negative"
        new(date, sum)
    end
end

+(x::Money, event::Event) = x + event.sum

Event(a, b, c, d) = Expense(a, b, c, d)
Event(a, b) = Deduction(a, b)
Event(a::Date, b) = Payment(a, b)

struct Budget
    moms::Percentile
    events::Vector{Event}
    Budget(moms::Percentile) = new(moms, Event[])
end

add!(b, x...) = push!(b.events, x...)

expenses(b) = sort(filter(x -> x isa Expense, b.events), by = x -> x.date)

deductions(b) = filter(x -> x isa Deduction, b.events)

payments(b) = sort(filter(x -> x isa Payment, b.events), by = x -> x.date)

nonexpenses(b) = filter(x -> !isa(x, Expense), b.events)

taxable(b) = reduce(+, expenses(b), init = Money(0))

tax(b) = taxable(b)*b.moms

balance(b) = taxable(b) + tax(b) + reduce(+, nonexpenses(b), init = Money(0))

# in

function define_event(l::String)
    a = match(r"^([edp]):(.*)"i, l)
    event = first(a.captures)
    fields = strip.(split(last(a.captures), ","))
    return if event == "e"
            a, b, c, d = fields
            Event(isempty(a) ? missing : Date(a), b, parse(Int, c), Money(parse(Float64, d)))
        elseif event == "d"
            a, b = fields
            Event(a, Money(parse(Float64, b)))
        elseif event == "p"
            a, b = fields
            Event(Date(a), Money(parse(Float64, b)))
        else
            error("event is neither e, d, nor p")
        end
    end
end

# out

_format(x) = format(x, commas = true, precision = 2 , stripzeros = true)

my_format(x::Number) = _format(x)
my_format(x::String) = x
my_format(x::Missing) = ""
my_format(x::Date) = Dates.format(x, "U d, Y")
function my_format(x::Percentile)
    a = my_format(x.p)
    return "$a\\%"
end
function my_format(x::Money)
    a = my_format(x.value)
    return "$a sek"
end
#=function format2md(x::Money)
    a = my_format(x.value)
    return "$a sek"
end=#
# format2md(x) = my_format(x)

function print(b::Budget)
    o = StringBuilder()
    for expense in expenses(b)
        for field in [:date, :description, :units, :rate]
            append!(o, my_format(getfield(expense, field)), " & ")
        end
        append!(o, my_format(expense.sum), " \\\\\n")
    end

    append!(o, "\\hline\n & Total & & & ", my_format(taxable(b)), "\\\\\n & VAT (", my_format(b.moms), ") & & &", my_format(tax(b)), "\\\\\n")

    for deduction in deductions(b)
        append!(o, " & ", deduction.description, " & & & ", my_format(deduction.sum), "\\\\\n")
    end

    for payment in payments(b)
        append!(o, my_format(payment.date), " & Payment & & & ", my_format(payment.sum), "\\\\\n")
    end

    append!(o, "& Balance due & & & ", my_format(balance(b)), "\\\\")
    return String(o)
end

end # module

ext: jl


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediapackage
using AWS.Compat
using AWS.UUIDs
"""
    ConfigureLogs()

Changes the Channel's properities to configure log subscription

# Required Parameters
- `id`: The ID of the channel to log subscription.

# Optional Parameters
- `egressAccessLogs`:
- `ingressAccessLogs`:
"""

configure_logs(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/configure_logs"; aws_config=aws_config)
configure_logs(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/configure_logs", args; aws_config=aws_config)

"""
    CreateChannel()

Creates a new Channel.

# Required Parameters
- `id`: The ID of the Channel. The ID must be unique within the region and it
cannot be changed after a Channel is created.

# Optional Parameters
- `description`: A short text description of the Channel.
- `tags`:
"""

create_channel(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/channels", Dict{String, Any}("id"=>id); aws_config=aws_config)
create_channel(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/channels", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    CreateHarvestJob()

Creates a new HarvestJob record.

# Required Parameters
- `endTime`: The end of the time-window which will be harvested

- `id`: The ID of the HarvestJob. The ID must be unique within the region
and it cannot be changed after the HarvestJob is submitted

- `originEndpointId`: The ID of the OriginEndpoint that the HarvestJob will harvest from.
This cannot be changed after the HarvestJob is submitted.

- `s3Destination`:
- `startTime`: The start of the time-window which will be harvested

"""

create_harvest_job(endTime, id, originEndpointId, s3Destination, startTime; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/harvest_jobs", Dict{String, Any}("endTime"=>endTime, "id"=>id, "originEndpointId"=>originEndpointId, "s3Destination"=>s3Destination, "startTime"=>startTime); aws_config=aws_config)
create_harvest_job(endTime, id, originEndpointId, s3Destination, startTime, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/harvest_jobs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "id"=>id, "originEndpointId"=>originEndpointId, "s3Destination"=>s3Destination, "startTime"=>startTime), args)); aws_config=aws_config)

"""
    CreateOriginEndpoint()

Creates a new OriginEndpoint record.

# Required Parameters
- `channelId`: The ID of the Channel that the OriginEndpoint will be associated with.
This cannot be changed after the OriginEndpoint is created.

- `id`: The ID of the OriginEndpoint.  The ID must be unique within the region
and it cannot be changed after the OriginEndpoint is created.

# Optional Parameters
- `authorization`:
- `cmafPackage`:
- `dashPackage`:
- `description`: A short text description of the OriginEndpoint.
- `hlsPackage`:
- `manifestName`: A short string that will be used as the filename of the OriginEndpoint URL (defaults to \"index\").
- `mssPackage`:
- `origination`: Control whether origination of video is allowed for this OriginEndpoint. If set to ALLOW, the OriginEndpoint
may by requested, pursuant to any other form of access control. If set to DENY, the OriginEndpoint may not be
requested. This can be helpful for Live to VOD harvesting, or for temporarily disabling origination

- `startoverWindowSeconds`: Maximum duration (seconds) of content to retain for startover playback.
If not specified, startover playback will be disabled for the OriginEndpoint.

- `tags`:
- `timeDelaySeconds`: Amount of delay (seconds) to enforce on the playback of live content.
If not specified, there will be no time delay in effect for the OriginEndpoint.

- `whitelist`: A list of source IP CIDR blocks that will be allowed to access the OriginEndpoint.
"""

create_origin_endpoint(channelId, id; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/origin_endpoints", Dict{String, Any}("channelId"=>channelId, "id"=>id); aws_config=aws_config)
create_origin_endpoint(channelId, id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/origin_endpoints", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelId"=>channelId, "id"=>id), args)); aws_config=aws_config)

"""
    DeleteChannel()

Deletes an existing Channel.

# Required Parameters
- `id`: The ID of the Channel to delete.

"""

delete_channel(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/channels/$(id)"; aws_config=aws_config)
delete_channel(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/channels/$(id)", args; aws_config=aws_config)

"""
    DeleteOriginEndpoint()

Deletes an existing OriginEndpoint.

# Required Parameters
- `id`: The ID of the OriginEndpoint to delete.

"""

delete_origin_endpoint(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/origin_endpoints/$(id)"; aws_config=aws_config)
delete_origin_endpoint(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/origin_endpoints/$(id)", args; aws_config=aws_config)

"""
    DescribeChannel()

Gets details about a Channel.

# Required Parameters
- `id`: The ID of a Channel.

"""

describe_channel(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/channels/$(id)"; aws_config=aws_config)
describe_channel(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/channels/$(id)", args; aws_config=aws_config)

"""
    DescribeHarvestJob()

Gets details about an existing HarvestJob.

# Required Parameters
- `id`: The ID of the HarvestJob.

"""

describe_harvest_job(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs/$(id)"; aws_config=aws_config)
describe_harvest_job(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs/$(id)", args; aws_config=aws_config)

"""
    DescribeOriginEndpoint()

Gets details about an existing OriginEndpoint.

# Required Parameters
- `id`: The ID of the OriginEndpoint.

"""

describe_origin_endpoint(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints/$(id)"; aws_config=aws_config)
describe_origin_endpoint(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints/$(id)", args; aws_config=aws_config)

"""
    ListChannels()

Returns a collection of Channels.

# Optional Parameters
- `maxResults`: Upper bound on number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
"""

list_channels(; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/channels"; aws_config=aws_config)
list_channels(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/channels", args; aws_config=aws_config)

"""
    ListHarvestJobs()

Returns a collection of HarvestJob records.

# Optional Parameters
- `includeChannelId`: When specified, the request will return only HarvestJobs associated with the given Channel ID.
- `includeStatus`: When specified, the request will return only HarvestJobs in the given status.
- `maxResults`: The upper bound on the number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
"""

list_harvest_jobs(; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs"; aws_config=aws_config)
list_harvest_jobs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs", args; aws_config=aws_config)

"""
    ListOriginEndpoints()

Returns a collection of OriginEndpoint records.

# Optional Parameters
- `channelId`: When specified, the request will return only OriginEndpoints associated with the given Channel ID.
- `maxResults`: The upper bound on the number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
"""

list_origin_endpoints(; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints"; aws_config=aws_config)
list_origin_endpoints(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints", args; aws_config=aws_config)

"""
    ListTagsForResource()

# Required Parameters
- `resource-arn`:

"""

list_tags_for_resource(resource_arn; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("GET", "/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    RotateChannelCredentials()

Changes the Channel's first IngestEndpoint's username and password. WARNING - This API is deprecated. Please use RotateIngestEndpointCredentials instead

# Required Parameters
- `id`: The ID of the channel to update.

"""

rotate_channel_credentials(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/credentials"; aws_config=aws_config)
rotate_channel_credentials(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/credentials", args; aws_config=aws_config)

"""
    RotateIngestEndpointCredentials()

Rotate the IngestEndpoint's username and password, as specified by the IngestEndpoint's id.

# Required Parameters
- `id`: The ID of the channel the IngestEndpoint is on.
- `ingest_endpoint_id`: The id of the IngestEndpoint whose credentials should be rotated

"""

rotate_ingest_endpoint_credentials(id, ingest_endpoint_id; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/ingest_endpoints/$(ingest_endpoint_id)/credentials"; aws_config=aws_config)
rotate_ingest_endpoint_credentials(id, ingest_endpoint_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/ingest_endpoints/$(ingest_endpoint_id)/credentials", args; aws_config=aws_config)

"""
    TagResource()

# Required Parameters
- `resource-arn`:
- `tags`:

"""

tag_resource(resource_arn, tags; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

# Required Parameters
- `resource-arn`:
- `tagKeys`: The key(s) of tag to be deleted

"""

untag_resource(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateChannel()

Updates an existing Channel.

# Required Parameters
- `id`: The ID of the Channel to update.

# Optional Parameters
- `description`: A short text description of the Channel.
"""

update_channel(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)"; aws_config=aws_config)
update_channel(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)", args; aws_config=aws_config)

"""
    UpdateOriginEndpoint()

Updates an existing OriginEndpoint.

# Required Parameters
- `id`: The ID of the OriginEndpoint to update.

# Optional Parameters
- `authorization`:
- `cmafPackage`:
- `dashPackage`:
- `description`: A short text description of the OriginEndpoint.
- `hlsPackage`:
- `manifestName`: A short string that will be appended to the end of the Endpoint URL.
- `mssPackage`:
- `origination`: Control whether origination of video is allowed for this OriginEndpoint. If set to ALLOW, the OriginEndpoint
may by requested, pursuant to any other form of access control. If set to DENY, the OriginEndpoint may not be
requested. This can be helpful for Live to VOD harvesting, or for temporarily disabling origination

- `startoverWindowSeconds`: Maximum duration (in seconds) of content to retain for startover playback.
If not specified, startover playback will be disabled for the OriginEndpoint.

- `timeDelaySeconds`: Amount of delay (in seconds) to enforce on the playback of live content.
If not specified, there will be no time delay in effect for the OriginEndpoint.

- `whitelist`: A list of source IP CIDR blocks that will be allowed to access the OriginEndpoint.
"""

update_origin_endpoint(id; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/origin_endpoints/$(id)"; aws_config=aws_config)
update_origin_endpoint(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage("PUT", "/origin_endpoints/$(id)", args; aws_config=aws_config)


text:
using .ExamplePKG
using Test

printDate(2019,12,11)

@testset "ExamplePKG.jl" begin

    @test myFun(3,2) == 3^2
    @test myFun(22,4) == 22^4

end


text:
using InteractiveUtils
include("halfedge.jl")
include("catmull-clark.jl")

code_warntype(build, (Vector{Float32}, Vector{UInt32}, Vector{UInt32}, Vector{UInt32}))

#code_warntype(get_vertex, (Vector{Float32}, UInt32))

#code_warntype(set_vertex, (Vector{Float32}, UInt32, Vector{Float32}))

#code_warntype(subdivide, (Vector{Float32}, Vector{UInt32}, Vector{UInt32}, Vector{UInt32}))

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
@testset "WeightedVariance fit!" begin
    test_fit(WeightedVariance{Float64},
             x, w,
             x -> (var(x), mean(x)),
             (x, w) -> (var(x, aweights(w), corrected = true), mean(x, weights(w))))
    test_fit(WeightedVariance{Float32},
             x, w,
             x -> (var(x), mean(x)),
             (x, w) -> (var(x, aweights(w), corrected = true), mean(x, weights(w))))

    test_fit(WeightedVariance{Float64},
             xmis, wmis,
             x -> (var(x), mean(x)),
             (x, w) -> begin
                ## doesn't work:
                # idx = findall((x) -> !ismissing(x[1]) & !ismissing(x[2]), zip(x, w))
                idx = findall(.!ismissing(x) .& .!ismissing.(w))
                xx = Float64.(x[idx])
                ww = Float64.(w[idx])
                (var(xx, aweights(ww), corrected = true), mean(xx, aweights(ww)))
             end
             )

    s, m, v = sum(x .* w), mean(x, weights(w)), var(x, aweights(w), corrected = true)
    sa, ma, va = sum(x .* w), mean(x, weights(w)), var(x, aweights(w), corrected = true)
    sf, mf, vf = sum(x .* w), mean(x, weights(w)), var(x, fweights(w), corrected = true)
    sp, mp, vp = sum(x .* w), mean(x, weights(w)), var(x, pweights(w), corrected = true)

    sval, mval, vval = fit!(WeightedVariance(), x, w)      |> x -> (sum(x), mean(x), var(x))
    szip, mzip, vzip = fit!(WeightedVariance(), zip(x, w)) |> x -> (sum(x), mean(x), var(x))

    svala, mvala, vvala = fit!(WeightedVariance(), x, w) |>
        x -> (sum(x), mean(x), var(x, corrected = true, weight_type = :analytic))
    svalf, mvalf, vvalf = fit!(WeightedVariance(), x, w) |>
        x -> (sum(x), mean(x), var(x, corrected = true, weight_type = :frequency))
    @test_throws ArgumentError fit!(WeightedVariance(), x, w) |>
        x -> (sum(x), mean(x), var(x, corrected = true, weight_type = :something))

    @test v ≈ vzip
    @test v ≈ vval
    @test m ≈ mzip
    @test m ≈ mval
    @test va ≈ vvala
    @test vf ≈ vvalf
    @test ma ≈ mvala
    @test mf ≈ mvalf

    @test s ≈ sval
    @test s ≈ szip
    @test s ≈ svala
    @test s ≈ svalf

    # After implementing :probability, these should pass/not throw any more:
    @test_throws ErrorException mvalp, vvalp = fit!(WeightedVariance(), x, w) |>
        x -> (mean(x), var(x, corrected = true, weight_type = :probability))
    @test_broken vp ≈ vvalp
    @test_broken mp ≈ mvalp
end

@testset "WeighedVariance merge!" begin
    v = var(x, aweights(w), corrected = true)

    wv = fit!(WeightedVariance(), x, w)
    ov = map(x, w) do xi, wi
        fit!(WeightedVariance(), xi, wi)
    end;
    rv = reduce(merge!, deepcopy(ov))
    rv2 = merge!(
        fit!(WeightedVariance(), x[1:end ÷ 2],       w[1:end ÷ 2]),
        fit!(WeightedVariance(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end]))

    rv_32 = reduce(merge!, deepcopy(ov), init = WeightedVariance(Float32))
    rv2_32 = merge!(
        fit!(WeightedVariance(Float32), x[1:end ÷ 2],       w[1:end ÷ 2]),
        fit!(WeightedVariance(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end]))

    @test rv.μ ≈ wv.μ
    @test rv.σ2 ≈ wv.σ2
    @test rv.W ≈ wv.W
    @test rv.W2 ≈ wv.W2

    @test rv2.μ ≈ wv.μ
    @test rv2.σ2 ≈ wv.σ2
    @test rv2.W ≈ wv.W
    @test rv2.W2 ≈ wv.W2

    @test var(rv) ≈ v
    @test var(rv2) ≈ v

    @test var(rv_32) ≈ v
    @test var(rv2_32) ≈ v

    @test mean(rv_32) ≈ wv.μ
    @test mean(rv2_32) ≈ wv.μ

    @test eltype(rv) == Float64
    @test eltype(rv2) == Float64
    @test eltype(rv_32) == Float32
    @test eltype(rv2_32) == Float32
end

@testset "WeightedVariance copy" begin
    o = fit!(WeightedVariance(), x2, w)

    o2 = copy(o)

    @test Ref(o.μ) != Ref(o.μ)
    @test Ref(o.σ2) != Ref(o.σ2)
    @test Ref(o.W) != Ref(o2.W)
    @test Ref(o.W2) != Ref(o2.W2)
    @test Ref(o.n) != Ref(o2.n)

    @test o.μ == o2.μ
    @test o.σ2 == o2.σ2
    @test o.W == o2.W
    @test o.W2 == o2.W2
    @test o.n == o2.n
end

@testset "WeightedVariance constructor" begin
    @test WeightedVariance{Float64}() == WeightedVariance()
    @test WeightedVariance{Float32}() == WeightedVariance(Float32)
    @test WeightedVariance() == WeightedVariance(0.0, 0.0, 0.0, 0.0, 0)
    @test WeightedVariance() == WeightedVariance{Float64}(0, 0, 0, 0, 0)
end


text:
using Test
using CoarseGraining3DPrisms
using LinearAlgebra

const K = 5
const x = K+1
const y = K/2

a = 0.0

#@test delta(1.0,1.0,1.0) == 1

#@show PrA

function dataPrA()
    PrA = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jd1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, ja2 in 0.:0.5:y,
                jb2 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y, jg in 0.:0.5:y
        ans = numchop( prismA(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
        PrA[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)] = ans
    end
    return PrA
end

function dataPrB2()
    PrB = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jd1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, ja2 in 0.:0.5:y,
                jb2 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y, jg in 0.:0.5:y
        ans = numchop( prismB2(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
        PrB[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)] = ans
    end
    return PrB
end

#PrA = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
function dataPrA()
    ampsInfo = []
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y, je1 in 0.:0.5:y
        if delta(jd1,jd2,je1) != 0
            for jb1 in 0.:0.5:y, jg in 0.:0.5:y
                if delta(jd1,jb1,jg) != 0
                    for ja1 in 0.:0.5:y
                        if delta(ja1,jd2,jg) != 0 && delta(ja1,jb1,je1) != 0
                            for jc in 0.:0.5:y, jb2 in 0.:0.5:y
                                if delta(jd1,jc,jb2) != 0
                                    for jf1 in 0.:0.5:y
                                        if delta(jf1,jd2,jb2) != 0 && delta(jf1,jc,je1) != 0
                                            for ja2 in 0.:0.5:y, jf2 in 0.:0.5:y
                                                if delta(jd1,ja2,jf2) != 0
                                                    for je2 in 0.:0.5:y
                                                        if delta(jc,je2,jf2) != 0 && delta(jb2,je2,ja2) != 0 #&& numchop( prismA(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg)) != 0
                                                            #if numchop( prismA(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg)) != 0
                                                            ans = numchop( prismA(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
                                                            if ans != 0
                                                                indx = [ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg]
                                                                push!(ampsInfo,(indx,ans))
                                                                #push!(lindx,)
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return ampsInfo
end

function dataB21()
    ampsInfo = []
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jd1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y, jg in 0.:0.5:y
        if numchop( prismB2(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg)) != 0
            ans = numchop( prismB2(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
            indx = [ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg]
            push!(ampsInfo,(indx,ans))
            #push!(lindx,)
        end
    end
    return ampsInfo
end

#PrB = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
#for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jd1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y, jg in 0.:0.5:y
#    ans = numchop( prismB(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
   # PrB[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)] = ans
#end

#PrB = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
#ampsB = Float64[]
#indxB = Array{Int64,1}[]
#for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jd1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y, jg in 0.:0.5:y
#    if numchop( prismB(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg)) != 0
#        ans = numchop( prismB(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
#        push!(ampsB,ans)
#        push!(indxB,[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)])
#    end
#end

#pos = findall(x-> x!=0, PrA )
#PrANZ = PrA[pos]
#posB = findall(x-> x!=0, PrB )
#PrANZ = PrB[posB]

#@show pos

#get the index for each of the shared triangular face.. return index values and amplitudes
# gives 30 elements for spins(0,0,0)

function svdRA(ampInfo::Array{Any,1},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    mat = blockReduceA(ampInfo,ja,jb,jc)[3]
    U, s, V = svd(mat)
    #V = ctranspose(V)
    return (U,s,V)
end

#test1 = svdRA(ampsA,indxA,0.,0.,0.5)[3][2]
#@time @show test1

@time datA1 = dataPrA();
#@time datPA = dataPrA();

#@time datB21 = dataB21();
#@time datPB2 = dataPrB2();

#@time for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
#    sing = svdA(datPA,ja,jb,jc)[2]
#    if sing[1] != 0
#        println(sing[1:2])
#    end
#end

@time for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
    sing = svdRA(datA1,ja,jb,jc)[2]
    if sing[1] != 0
        println(sing[1:2])
    end
end

#@time @show blockReduceA1(datA1,0.,0.,0.)

#@time println(svdA(PrA,0.0,0.5,0.5)[2][1:6])
#@time for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
#    sing = svdA(PrA,ja,jb,jc)[2]
#    if sing[1] != 0
#        println(sing[1:3])
#    end
#end

#@time fullUtensorA(PrA)
#@time fullUtensorB(PrB)
#@time fullVtensorA(PrA)
#@time prsmUAVB(PrA,PrB)
#@time println(svdR(PrA,0.,0.5,0.5)[2])
#@time println(svdRA(PrA,0.,0.5,0.5)[2])

#@time for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
#    sing = svdRA2(ampsB,indxB,ja,jb,jc)[2]
#    if sing[1] != 0
#        println(sing)
#    end
#end

#@time blockReduceA(ampsB,indxB,0.,0.5,1.)
#@time println(loppoF(PrA,0.,0.5,0.5))

CoarseGraining3DPrisms.greet()


text:
# ---
# title: 617. Merge Two Binary Trees
# id: problem617
# author: zhwang
# date: 2022-02-19
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/merge-two-binary-trees/description/>
# hidden: true
# ---
#
# Given two binary trees and imagine that when you put one of them to cover the
# other, some nodes of the two trees are overlapped while the others are not.
#
# You need to merge them into a new binary tree. The merge rule is that if two
# nodes overlap, then sum node values up as the new value of the merged node.
# Otherwise, the NOT null node will be used as the node of new tree.
#
# **Example 1:**
#
#
#
#     Input:
#     	Tree 1                     Tree 2
#               1                         2
#              / \                       / \
#             3   2                     1   3
#            /                           \   \
#           5                             4   7
#     Output:
#     Merged tree:
#     	     3
#     	    / \
#     	   4   5
#     	  / \   \
#     	 5   4   7
#
#
#
#
# **Note:** The merging process must start from the root nodes of both trees.
#
#
## @lc code=start
using LeetCode

merge_trees(::Nothing, ::Nothing) = nothing
merge_trees(::Nothing, root2::TreeNode) = root2
merge_trees(root1::TreeNode, ::Nothing) = root1
function merge_trees(root1::TreeNode, root2::TreeNode)::TreeNode
    root = TreeNode(root1.val + root2.val)
    root.left = merge_trees(root1.left, root2.left)
    root.right = merge_trees(root1.right, root2.right)
    return root
end
## @lc code=end

## @lc test=start
@testset "617.merge-two-binary-trees.jl" begin
    @test merge_trees(TreeNode{Int}([1,3,2,5]), TreeNode{Int}([2,1,3,nothing, 4, nothing, 7])) == TreeNode{Int}([3,4,5,5,4,nothing,7])
    @test merge_trees(TreeNode{Int}([1]), TreeNode{Int}([1,2])) == TreeNode{Int}([2,2])
    @test merge_trees(nothing, nothing) === nothing
end
## @lc test=end

## @lc info=start
# link: [solution to 617](https://leetcode-cn.com/problems/merge-two-binary-trees/solution/617-he-bing-er-cha-shu-by-rex-rs-tyx3/)
## @lc info=end


text:
"""
    struct DirichletFESpace <: SingleFieldFESpace
      space::SingleFieldFESpace
    end
"""
struct DirichletFESpace{B} <: SingleFieldFESpace
  space::SingleFieldFESpace
  constraint_style::Val{B}
  function DirichletFESpace(space::SingleFieldFESpace)
    cs = constraint_style(space)
    B = get_val_parameter(cs)
    new{B}(space,cs)
  end
end

constraint_style(::Type{DirichletFESpace{B}}) where B = Val{B}()

get_cell_axes(t::DirichletFESpace)= get_cell_axes(t.space)

get_cell_axes_with_constraints(t::DirichletFESpace)= get_cell_axes_with_constraints(t.space)

CellData.CellField(t::DirichletFESpace,cell_vals) = CellField(t.space,cell_vals)

get_cell_isconstrained(f::DirichletFESpace) = get_cell_isconstrained(f.space)

get_cell_constraints(f::DirichletFESpace) = get_cell_constraints(f.space)

function num_free_dofs(f::DirichletFESpace)
  num_dirichlet_dofs(f.space)
end

function zero_free_values(f::DirichletFESpace)
  zero_dirichlet_values(f.space)
end

function get_cell_dofs(f::DirichletFESpace)
  apply(elem(-),get_cell_dofs(f.space))
end

function num_dirichlet_dofs(f::DirichletFESpace)
  num_free_dofs(f.space)
end

function zero_dirichlet_values(f::DirichletFESpace)
  zero_free_values(f.space)
end

function num_dirichlet_tags(f::DirichletFESpace)
  1
end

function get_dirichlet_dof_tag(f::DirichletFESpace)
  ones(Int8,num_dirichlet_dofs(f))
end

function scatter_free_and_dirichlet_values(f::DirichletFESpace,fv,dv)
  scatter_free_and_dirichlet_values(f.space,dv,fv)
end

function gather_free_and_dirichlet_values!(fv,dv,f::DirichletFESpace,cv)
  gather_free_and_dirichlet_values!(dv,fv,f.space,cv)
  (fv,dv)
end

function TrialFESpace(f::DirichletFESpace)
  U = TrialFESpace(f.space)
  DirichletFESpace(U)
end

function get_cell_basis(f::DirichletFESpace)
  get_cell_basis(f.space)
end

function get_cell_dof_basis(f::DirichletFESpace)
  get_cell_dof_basis(f.space)
end


text:
using Revise
using AdFem
using PyPlot

n = 20
mmesh = Mesh(n, n, 1/n, degree=2)
bdedge = bcedge((x1,y1,x2,y2)->((y1<1e-5) && (y2<1e-5)), mmesh)
bdnode = bcnode((x,y)->!(y<1e-5 && x>1e-5 && x<1-1e-5), mmesh)

F1 = eval_f_on_gauss_pts((x,y)->3.0, mmesh)
F2 = eval_f_on_gauss_pts((x,y)->-1.0, mmesh)
F = compute_fem_source_term(F1, F2, mmesh)

t1 = eval_f_on_boundary_edge((x,y)->-x-y, bdedge, mmesh)
t2 = eval_f_on_boundary_edge((x,y)->2y, bdedge, mmesh)
T = compute_fem_traction_term(t1, t2, bdedge, mmesh)

D = constant(diagm(0=>[1,1,0.5]))
K = compute_fem_stiffness_matrix(D, mmesh)
rhs = T - F
bdval = [eval_f_on_boundary_node((x,y)->x^2+y^2, bdnode, mmesh);
        eval_f_on_boundary_node((x,y)->x^2-y^2, bdnode, mmesh)]
DOF = [bdnode;bdnode .+ mmesh.ndof]
K, rhs = impose_Dirichlet_boundary_conditions(K, rhs, DOF, bdval)
u = K\rhs
sess = Session(); init(sess)
S = run(sess, u)

figure(figsize=[10,4])
subplot(121)
visualize_scalar_on_fem_points(S[1:mmesh.nnode], mmesh)
subplot(122)
visualize_scalar_on_fem_points(S[mmesh.ndof+1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("numerical.png")

figure(figsize=[10,4])
x = LinRange(0, 1, 50)
X = zeros(50, 50)
Y = zeros(50, 50)
for i = 1:50
    for j = 1:50
        X[i,j] = x[i]
        Y[i,j] = x[j]
    end
end
subplot(121)
X = mmesh.nodes[:,1]
Y = mmesh.nodes[:,2]
U = (@. X^2+Y^2)
V = (@. X^2-Y^2)
subplot(121)
visualize_scalar_on_fem_points(U, mmesh)
subplot(122)
visualize_scalar_on_fem_points(V, mmesh)
savefig("exact.png")

figure(figsize=[10,4])
subplot(121)
visualize_scalar_on_fem_points(U-S[1:mmesh.nnode], mmesh)
subplot(122)
visualize_scalar_on_fem_points(V-S[mmesh.ndof+1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("difference.png")


text:
module SpikingNNFunctions

include("neurons/Neuron.jl")
include("synapses/Synapse.jl")
include("thresholds/Threshold.jl")

end


text:
using Test
using LogicCircuits

@testset "Decomposability tests" begin
    lin = LogicalΔNode[]
    ors = map(1:10) do v
        pos = LiteralNode(var2lit(Var(v)))
        push!(lin, pos)
        neg = LiteralNode(-var2lit(Var(v)))
        push!(lin, neg)
        or = ⋁Node([pos,neg])
        push!(lin, or)
        or
    end
    and = ⋀Node(ors)
    push!(lin, and)
    bias = ⋁Node([and])
    push!(lin, bias)

    @test is_decomposable(lin)

    lin2 = LogicalΔNode[]
    ors = map(1:10) do v
        pos = LiteralNode(var2lit(Var(v)))
        push!(lin2, pos)
        neg = LiteralNode(-var2lit(Var(v)))
        push!(lin2, neg)
        or = ⋁Node([pos,neg])
        push!(lin2, or)
        or
    end
    and1 = ⋀Node([ors[1], ors[2], ors[1], ors[3]])
    and2 = ⋀Node(ors[5:10])
    push!(lin2, and1)
    push!(lin2, and2)
    or1 = ⋁Node([and1, and2])
    push!(lin2, or1)
    @test !is_decomposable(lin2)

    lin3 = LogicalΔNode[]
    ors = map(1:10) do v
        pos = LiteralNode(var2lit(Var(v)))
        push!(lin3, pos)
        neg = LiteralNode(-var2lit(Var(v)))
        push!(lin3, neg)
        or = ⋁Node([pos,neg])
        push!(lin3, or)
        or
    end
    and1 = ⋀Node(ors[1:3])
    and2 = ⋀Node(ors[3:6])
    push!(lin3, and1)
    push!(lin3, and2)
    and3 = ⋀Node(ors[6:8])
    and4 = ⋀Node(ors[8:10])
    push!(lin3, and3)
    push!(lin3, and4)
    or1 = ⋁Node([and1, and2])
    or2 = ⋁Node([and3, and4])
    push!(lin3, or1)
    push!(lin3, or2)
    and5 = ⋀Node([or1, or2])
    push!(lin3, and5)

    @test !is_decomposable(lin3)

    simpletest1 = LogicalΔNode[]
    push!(simpletest1, LiteralNode(Lit(1)))
    @test is_decomposable(simpletest1)

    simpletest2 = LogicalΔNode[]
    leaf1 = LiteralNode(Lit(1))
    leaf2 = LiteralNode(Lit(-1))
    and = ⋀Node([leaf1, leaf2])
    push!(simpletest2, leaf1)
    push!(simpletest2, leaf2)
    push!(simpletest2, and)
    @test !is_decomposable(simpletest2)
end


text:
"""
Simple script to grab and install these packages into the global Julia install (easiest
way of handling this).
"""

using Pkg

Pkg.add("Distributions")
Pkg.add("Random" )
Pkg.add("Statistics" )
Pkg.add("LinearAlgebra")
Pkg.add("DataStructures")

Pkg.add("FileIO")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Plots")
Pkg.add("Measures")
Pkg.add("Pipe")
Pkg.add("PDFmerger") # *NOTE* that this will not be added on phoenix

Pkg.add("DelimitedFiles")
Pkg.add("Dates")

Pkg.add("Distributed")
Pkg.add("SharedArrays")
Pkg.add("ProgressMeter")
Pkg.add("TimerOutputs")
Pkg.add("ProgressBars")
Pkg.add("BenchmarkTools")
Pkg.add("StaticArrays")


text:
import URIParser

function getrepohttpurl(reporemoteurl::AbstractString)
    repouri = URIParser.URI(reporemoteurl)
    repopath = splitext(repouri.path)[1] # remove git suffix
    repourl = URIParser.URI("https", repouri.host, repouri.port, repopath)
    string(repourl)
end


text:
# ------- #
# Objects #
# ------- #

"""
    Object

All primitive objects must be a subtype of this. Currently there are
two primitive objects present, [`Triangle`](@ref) and [`Sphere`](@ref).
To add support for custom primitive type, two functions need to be
defined - [`intersect`](@ref) and [`get_normal`](@ref).
"""
abstract type Object end

"""
    get_color(obj::Object, pt::Vec3, sym::Val)

Computes the color at the point `pt` of the Object `obj` by dispatching
to the correct `get_color(obj.material, pt, sym, obj)` method.
"""
get_color(obj::Object, pt::Vec3, sym::Val) =
    get_color(obj.material, pt, sym, obj)

"""
    specular_exponent(obj::Object)

Returns the `specular_exponent` of the material of the Object.
"""
specular_exponent(obj::Object) = specular_exponent(obj.material)

"""
    reflection(obj::Object)

Returns the `reflection coefficient` of the material of the Object.
"""
reflection(obj::Object) = reflection(obj.material)

"""
    intersect(obj::Object, origin, direction)

Computes the intersection of the light ray with the object.
This function returns the ``t`` value where
``intersection\\_point = origin + t \\times direction``. In case
the ray does not intersect with the object `Inf` is returned.
"""
function intersect(obj::Object, origin, direction) end

"""
    get_normal(obj::Object, pt, direction)

Returns the normal at the Point `pt` of the Object `obj`.
"""
function get_normal(obj::Object, pt, direction) end

# ----------- #
# - Imports - #
# ----------- #

include("objects/sphere.jl")
include("objects/triangle.jl")
include("objects/polygon_mesh.jl")
include("objects/obj_parser.jl")


text:
@testset "Tiling" begin
    t = Tiling((1:2:5, 10:5:20))
    @test encode(t, (2, 12)) == 1
    @test encode(t, (2, 18)) == 3

    t2 = t - (1, 3)
    @test t2.ranges == (0:2:4, 7:5:17)
end


text:
module ValSplit

using ExprTools, Tricks

export @valsplit, valarg_params, valarg_has_param

include("utils.jl")

"""
    valarg_params(f, types::Type{<:Tuple}, idx::Int, ptype::Type=Any)
    valarg_params(f, idx::Int, ptypes::Type=Any)

Given a method signature `(f, types=Tuple{Vararg{Any}})`, finds all matching
methods with a concrete `Val`-typed argument in position `idx`, then returns
all parameter values for the `Val`-typed argument as a tuple. Optionally,
`ptype` can be specified to filter parameters that are instances of `ptype`.

    valarg_params(f, types::Type{<:Tuple}, idxs::Tuple, ptypes::Type)
    valarg_params(f, idxs::Tuple, ptypes::Type)

To find parameter values for multiple argument indices, `idxs` can also be
specified as a `Tuple`. To filter each argument separately, `ptypes` should be
specified as a `Tuple` type.

This function is statically compiled, and will automatically be recompiled
whenever a new method of `f` is defined.
"""
@generated function valarg_params(
    @nospecialize(f) , @nospecialize(types::Type{T}),
    @nospecialize(idx::Val{N}), @nospecialize(ptype::Type{P}=Any)
) where {T <: Tuple, N, P}
    # Check if N and P are single values or tuples
    multi_idx = N isa Tuple
    if multi_idx
        idxs = collect(Int, N)
        P = P <: Tuple ? P : Tuple{fill(P, length(idxs))...}
    else
        idxs = [N]
    end
    # Extract parameters of Val-typed argument from matching methods
    vparams = Vector{P}()
    for m in Tricks._methods(f, T)
        argtypes = fieldtypes(m.sig)[2:end]
        all(idxs .<= length(argtypes)) || continue
        types = argtypes[idxs]
        all(ty <: Val && isconcretetype(ty) for ty in types) || continue
        ps = multi_idx ? Tuple(val_param.(types)) : val_param(types[1])
        ps isa P || continue
        push!(vparams, ps)
    end
    unique!(vparams)
    vparams = Tuple(vparams)
    # Create CodeInfo and add edges so if a method is defined this recompiles
    ci = Tricks.create_codeinfo_with_returnvalue(
        [Symbol("#self#"), :f, :types, :idx, :ptype],
        [:T, :N, :P], (:T, :N, :P), :($vparams))
    ci.edges = Tricks._method_table_all_edges_all_methods(f, T)
    return ci
end
valarg_params(f, types::Type{<:Tuple}, idx::Union{Int,Tuple}, ptype::Type=Any) =
    valarg_params(f, types, Val(idx), ptype)
valarg_params(f, idx::Val{N}, ptype::Type=Any) where {N} =
    valarg_params(f, Tuple{Vararg{Any}}, idx, ptype)
valarg_params(f, idx::Union{Int,Tuple}, ptype::Type=Any) =
    valarg_params(f, Tuple{Vararg{Any}}, Val(idx), ptype)

"""
    valarg_has_param(param, f, types::Type{<:Tuple}, idx::Int, ptype::Type=Any)
    valarg_has_param(param, f, idx::Int, ptype::Type=Any)

Given a method signature `(f, types)`, returns `true` if there exists a
matching method with a `Val`-typed argument in position `idx` with parameter
`param` and parameter type `ptype`.

    valarg_has_param(params, f, types::Type{<:Tuple}, idxs::Tuple, ptypes::Type)
    valarg_has_param(params, f, idxs::Tuple, ptypes::Type)

To check parameter values for multiple argument indices, `idxs` can also be
specified as a `Tuple`. To filter each argument separately, `ptypes` should be
specified as a `Tuple` type.
"""
valarg_has_param(param::P, f, types::Type{<:Tuple}, idx, ptype::Type{P}=Any) where {P} =
    param in valarg_params(f, types, idx, ptype)
valarg_has_param(param::P, f, idx, ptype::Type{P}=Any) where {P} =
    param in valarg_params(f, Tuple{Vararg{Any}}, idx, ptype)

"""
    _valswitch(::Val{Vs}, ::Val{I}, f, default_f, args...) where {Vs, I}

Generates a switch statement with `args[N] == v` where `v` ∈ `Vs` as the
branch conditions. Each branch with value `v` calls `f(args...)`, with
`args[idx]` replaced by `Val(v)`. If all branch conditions fail, `default_f`
is called with with no arguments.
"""
@generated function _valswitch(
    vals::Val{Vs}, idx::Val{I}, f, default_f, args::Vararg{Any,N}
) where {Vs, I, N}
    vals = map(QuoteNode, Vs)
    cond_exprs = [Expr(:call, :(==), :(args[$I]), v) for v in vals]
    branch_exprs = map(vals) do v
        args = [i == I ? :(Val($v)) : :(args[$i]) for i in 1:N]
        return Expr(:call, :f, args...)
    end
    default_expr = :(default_f())
    return generate_switch_stmt(cond_exprs, branch_exprs, default_expr)
end

"""
    _valsplit(expr, idx::Int, val_idxs=[idx])

Generates function expression(s) returned by `@valsplit` macro. `idx` is the
index of the argument to split on. `val_idxs` are the indices of all arguments
that are `Val`-typed, defining the set of matching methods to switch over.
"""
function _valsplit(def::Dict{Symbol}, idx::Int, val_idxs=[idx])
    if !haskey(def, :args) || length(def[:args]) < idx
        error("Function has less than $idx arguments")
    end
    # Fill-in unnamed functions and arguments
    def[:name] = fill_unnamed(get(gensym, def, :name))
    def[:args] = map(fill_unnamed, def[:args])
    # Extract function name, signature, and argument expressions
    fname = rm_type_annotation(def[:name])
    types = args_tupletype_expr(def[:args], esc)
    argnames = collect(args_tuple_expr(def[:args]).args)
    # Extract type of argument to split on
    ptype = types.args[idx + 1]
    types.args[val_idxs .+ 1] .= :Val
    if is_vararg_expr(ptype)
        error("Cannot split Vararg arguments.")
    end
    # Escape function name and arguments
    def[:name] = esc(def[:name])
    def[:args] = map(esc, def[:args])
    if haskey(def, :whereparams)
        def[:whereparams] = map(esc, def[:whereparams])
    end
    # Generate the function body
    def[:body] = quote
        # Look up the parameters for the Val-typed argument in position idx
        vals = valarg_params($(esc(fname)), $types, $idx, $ptype)
        # Default function returns the original function body
        function default_f() $(esc(def[:body])) end
        # Generate a switch expression over the Val-type parameters
        return _valswitch(Val(vals), Val($idx), $(esc(fname)), default_f,
                          $(map(esc, argnames)...))
    end
    # Return recombined function expression
    return combinedef(def)
end
_valsplit(expr::Expr, idx::Int, val_idxs=[idx]) =
    _valsplit(splitdef(expr), idx, val_idxs)

"""
    @valsplit f
    @valsplit idx::Int f

Given a function definition `f`, the `@valsplit` macro compiles away dynamic
dispatch over methods of `f` with `Val`-typed arguments at the specified
indices (i.e., it "splits" on `Val`-typed arguments, similar to union
splitting). The resulting method is automatically recompiled whenever a new
method of `f` is defined.

Using the first form of `@valsplit`, each argument `x::T` to split on
should be annotated as `Val(x::T)`. Alternatively, an argument index `idx` can
be manually specified using the second form of the macro.

!!! note "Documentation Issues"
    Documenting a `@valsplit`-annotated function that splits on multiple
    arguments will lead to errors. Consider documenting the function separately.

# Example

Suppose we have a function `soundof` that returns how an animal sounds:

```julia
soundof(animal::Val{:dog}) = "woof"
soundof(animal::Val{:cat}) = "nyan"
```

Using `@valsplit`, we can define a new method for `soundof` that accepts a
`Symbol` argument, and branches to each instance of `soundof(::Val{T})`
(where `T` is a `Symbol`) defaulting to the function body otherwise:

```julia
@valsplit function soundof(Val(animal::Symbol))
    error("Sound not defined for animal: \$animal")
end
```

The resulting method is equivalent in behavior to the following:
```julia
function soundof(animal::Symbol)
    if animal == :dog
        return soundof(Val{:dog}())
    elseif animal == :cat
        return soundof(Val{:cat}())
    else
        error("Sound not defined for animal: \$animal")
    end
end
```
"""
macro valsplit(expr)
    def = splitdef(expr)
    # Fill-in unnamed functions and arguments
    def[:name] = fill_unnamed(get(gensym, def, :name))
    def[:args] = map(fill_unnamed, def[:args])
    # Find arguments to split by value (notated by Val(arg::T))
    split_idxs = Int[]
    unwrapped_args = []
    for (i, arg) in enumerate(get(def, :args, []))
        if Meta.isexpr(arg, :call, 2) && arg.args[1] == :Val
            push!(split_idxs, i)
            push!(unwrapped_args, arg.args[2])
        else
            push!(unwrapped_args, arg)
        end
    end
    argnames = collect(args_tuple_expr(unwrapped_args).args)
    # Generate function definition for each argument to split on
    i_def = copy(def)
    i_def[:args] = unwrapped_args
    f_exprs = Expr[]
    while !isempty(split_idxs)
        idx = first(split_idxs)
        push!(f_exprs, _valsplit(copy(i_def), idx, split_idxs))
        popfirst!(split_idxs)
        # Adjust function signature for next definition
        V_typevar = gensym(:V)
        i_def[:args][idx] = Expr(:(::), argnames[idx], V_typevar)
        push!(get!(i_def, :whereparams, []), :($V_typevar <: $(QuoteNode(Val))))
    end
    return Expr(:block, f_exprs...)
end

macro valsplit(idx::Int, expr)
    return _valsplit(expr, idx)
end

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in the file `FMI3_ext.jl` (external/additional functions)?
# - new functions, that are useful, but not part of the FMI-spec (example: `fmi3Load`)

using Libdl
using ZipFile

"""
Sets the properties of the fmu by reading the modelDescription.xml.
Retrieves all the pointers of binary functions.
Returns the instance of the FMU struct.
Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi3Load(pathToFMU::String; unpackPath=nothing)
    # Create uninitialized FMU
    fmu = FMU3()
    fmu.components = []

    pathToFMU = normpath(pathToFMU)

    # set paths for fmu handling
    (fmu.path, fmu.zipPath) = fmi3Unzip(pathToFMU; unpackPath=unpackPath) # TODO

    # set paths for modelExchangeScripting and binary
    tmpName = splitpath(fmu.path)
    pathToModelDescription = joinpath(fmu.path, "modelDescription.xml")

    # parse modelDescription.xml
    fmu.modelDescription = fmi3LoadModelDescription(pathToModelDescription) # TODO Matrix mit Dimensions
    fmu.modelName = fmu.modelDescription.modelName
    fmu.instanceName = fmu.modelDescription.modelName
    fmuName = fmi3GetModelIdentifier(fmu.modelDescription) # tmpName[length(tmpName)] TODO

    directoryBinary = ""
    pathToBinary = ""

    if Sys.iswindows()
        directories = [joinpath("binaries", "win64"), joinpath("binaries", "x86_64-windows")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$(fmuName).dll")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is Windows, but can't find valid FMU binary at `$(pathToBinary)` for path `$(fmu.path)`."
    elseif Sys.islinux()
        directories = [joinpath("binaries", "linux64"), joinpath("binaries", "x86_64-linux")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$(fmuName).so")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is Linux, but can't find valid FMU binary at `$(pathToBinary)` for path `$(fmu.path)`."
    elseif Sys.isapple()
        directories = [joinpath("binaries", "darwin64"), joinpath("binaries", "x86_64-darwin")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$(fmuName).dylib")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is macOS, but can't find valid FMU binary at `$(pathToBinary)` for path `$(fmu.path)`."
    else
        @assert false "Unsupported target platform. Supporting Windows64, Linux64 and Mac64."
    end

    lastDirectory = pwd()
    cd(directoryBinary)

    # set FMU binary handler
    fmu.libHandle = dlopen(pathToBinary)

    cd(lastDirectory)

    if fmi3IsCoSimulation(fmu)
        fmu.type = fmi3TypeCoSimulation
    elseif fmi3IsModelExchange(fmu)
        fmu.type = fmi3TypeModelExchange
    elseif fmi3IsScheduledExecution(fmu)
        fmu.type = fmi3TypeScheduledExecution
    else
        error(unknownFMUType)
    end

    if fmi3IsCoSimulation(fmu) && fmi3IsModelExchange(fmu)
        @info "fmi3Load(...): FMU supports both CS and ME, using CS as default if nothing specified." # TODO ScheduledExecution
    end

    # make URI ressource location
    tmpResourceLocation = string("file:///", fmu.path)
    tmpResourceLocation = joinpath(tmpResourceLocation, "resources")
    fmu.fmuResourceLocation = replace(tmpResourceLocation, "\\" => "/") # URIs.escapeuri(tmpResourceLocation)

    @info "fmi3Load(...): FMU resources location is `$(fmu.fmuResourceLocation)`"

    # # retrieve functions
    fmu.cInstantiateModelExchange                  = dlsym(fmu.libHandle, :fmi3InstantiateModelExchange)
    fmu.cInstantiateCoSimulation                   = dlsym(fmu.libHandle, :fmi3InstantiateCoSimulation)
    fmu.cInstantiateScheduledExecution             = dlsym(fmu.libHandle, :fmi3InstantiateScheduledExecution)
    fmu.cGetVersion                                = dlsym(fmu.libHandle, :fmi3GetVersion)
    fmu.cFreeInstance                              = dlsym(fmu.libHandle, :fmi3FreeInstance)
    fmu.cSetDebugLogging                           = dlsym(fmu.libHandle, :fmi3SetDebugLogging)
    fmu.cEnterConfigurationMode                    = dlsym(fmu.libHandle, :fmi3EnterConfigurationMode)
    fmu.cExitConfigurationMode                     = dlsym(fmu.libHandle, :fmi3ExitConfigurationMode)
    fmu.cEnterInitializationMode                   = dlsym(fmu.libHandle, :fmi3EnterInitializationMode)
    fmu.cExitInitializationMode                    = dlsym(fmu.libHandle, :fmi3ExitInitializationMode)
    fmu.cTerminate                                 = dlsym(fmu.libHandle, :fmi3Terminate)
    fmu.cReset                                     = dlsym(fmu.libHandle, :fmi3Reset)
    fmu.cEvaluateDiscreteStates                    = dlsym(fmu.libHandle, :fmi3EvaluateDiscreteStates)
    fmu.cGetNumberOfVariableDependencies           = dlsym(fmu.libHandle, :fmi3GetNumberOfVariableDependencies)
    fmu.cGetVariableDependencies                   = dlsym(fmu.libHandle, :fmi3GetVariableDependencies)

    fmu.cGetFloat32                                = dlsym(fmu.libHandle, :fmi3GetFloat32)
    fmu.cSetFloat32                                = dlsym(fmu.libHandle, :fmi3SetFloat32)
    fmu.cGetFloat64                                = dlsym(fmu.libHandle, :fmi3GetFloat64)
    fmu.cSetFloat64                                = dlsym(fmu.libHandle, :fmi3SetFloat64)
    fmu.cGetInt8                                   = dlsym(fmu.libHandle, :fmi3GetInt8)
    fmu.cSetInt8                                   = dlsym(fmu.libHandle, :fmi3SetInt8)
    fmu.cGetUInt8                                  = dlsym(fmu.libHandle, :fmi3GetUInt8)
    fmu.cSetUInt8                                  = dlsym(fmu.libHandle, :fmi3SetUInt8)
    fmu.cGetInt16                                  = dlsym(fmu.libHandle, :fmi3GetInt16)
    fmu.cSetInt16                                  = dlsym(fmu.libHandle, :fmi3SetInt16)
    fmu.cGetUInt16                                 = dlsym(fmu.libHandle, :fmi3GetUInt16)
    fmu.cSetUInt16                                 = dlsym(fmu.libHandle, :fmi3SetUInt16)
    fmu.cGetInt32                                  = dlsym(fmu.libHandle, :fmi3GetInt32)
    fmu.cSetInt32                                  = dlsym(fmu.libHandle, :fmi3SetInt32)
    fmu.cGetUInt32                                 = dlsym(fmu.libHandle, :fmi3GetUInt32)
    fmu.cSetUInt32                                 = dlsym(fmu.libHandle, :fmi3SetUInt32)
    fmu.cGetInt64                                  = dlsym(fmu.libHandle, :fmi3GetInt64)
    fmu.cSetInt64                                  = dlsym(fmu.libHandle, :fmi3SetInt64)
    fmu.cGetUInt64                                 = dlsym(fmu.libHandle, :fmi3GetUInt64)
    fmu.cSetUInt64                                 = dlsym(fmu.libHandle, :fmi3SetUInt64)
    fmu.cGetBoolean                                = dlsym(fmu.libHandle, :fmi3GetBoolean)
    fmu.cSetBoolean                                = dlsym(fmu.libHandle, :fmi3SetBoolean)

    fmu.cGetString                                 = dlsym_opt(fmu.libHandle, :fmi3GetString)
    fmu.cSetString                                 = dlsym_opt(fmu.libHandle, :fmi3SetString)
    fmu.cGetBinary                                 = dlsym_opt(fmu.libHandle, :fmi3GetBinary)
    fmu.cSetBinary                                 = dlsym_opt(fmu.libHandle, :fmi3SetBinary)

    if fmi3CanGetSetState(fmu)
        fmu.cGetFMUState                           = dlsym_opt(fmu.libHandle, :fmi3GetFMUState)
        fmu.cSetFMUState                           = dlsym_opt(fmu.libHandle, :fmi3SetFMUState)
        fmu.cFreeFMUState                          = dlsym_opt(fmu.libHandle, :fmi3FreeFMUState)
    end

    if fmi3CanSerializeFMUstate(fmu)
        fmu.cSerializedFMUStateSize                = dlsym_opt(fmu.libHandle, :fmi3SerializedFMUStateSize)
        fmu.cSerializeFMUState                     = dlsym_opt(fmu.libHandle, :fmi3SerializeFMUState)
        fmu.cDeSerializeFMUState                   = dlsym_opt(fmu.libHandle, :fmi3DeSerializeFMUState)
    end

    if fmi3ProvidesDirectionalDerivatives(fmu)
        fmu.cGetDirectionalDerivative              = dlsym_opt(fmu.libHandle, :fmi3GetDirectionalDerivative)
    end

    if fmi3ProvidesAdjointDerivatives(fmu)
        fmu.cGetAdjointDerivative              = dlsym_opt(fmu.libHandle, :fmi3GetAdjointDerivative)
    end

    # CS specific function calls
    if fmi3IsCoSimulation(fmu)
        fmu.cGetOutputDerivatives                  = dlsym(fmu.libHandle, :fmi3GetOutputDerivatives)
        fmu.cEnterStepMode                         = dlsym(fmu.libHandle, :fmi3EnterStepMode)
        fmu.cDoStep                                = dlsym(fmu.libHandle, :fmi3DoStep)
    end

    # ME specific function calls
    if fmi3IsModelExchange(fmu)
        fmu.cGetNumberOfContinuousStates           = dlsym(fmu.libHandle, :fmi3GetNumberOfContinuousStates)
        fmu.cGetNumberOfEventIndicators            = dlsym(fmu.libHandle, :fmi3GetNumberOfEventIndicators)
        fmu.cGetContinuousStates                   = dlsym(fmu.libHandle, :fmi3GetContinuousStates)
        fmu.cGetNominalsOfContinuousStates         = dlsym(fmu.libHandle, :fmi3GetNominalsOfContinuousStates)
        fmu.cEnterContinuousTimeMode               = dlsym(fmu.libHandle, :fmi3EnterContinuousTimeMode)
        fmu.cSetTime                               = dlsym(fmu.libHandle, :fmi3SetTime)
        fmu.cSetContinuousStates                   = dlsym(fmu.libHandle, :fmi3SetContinuousStates)
        fmu.cGetContinuousStateDerivatives         = dlsym(fmu.libHandle, :fmi3GetContinuousStateDerivatives)
        fmu.cGetEventIndicators                    = dlsym(fmu.libHandle, :fmi3GetEventIndicators)
        fmu.cCompletedIntegratorStep               = dlsym(fmu.libHandle, :fmi3CompletedIntegratorStep)
        fmu.cEnterEventMode                        = dlsym(fmu.libHandle, :fmi3EnterEventMode)
        fmu.cUpdateDiscreteStates                  = dlsym(fmu.libHandle, :fmi3UpdateDiscreteStates)
    end

    if fmi3IsScheduledExecution(fmu)
        fmu.cSetIntervalDecimal                    = dlsym(fmu.libHandle, :fmi3SetIntervalDecimal)
        fmu.cSetIntervalFraction                   = dlsym(fmu.libHandle, :fmi3SetIntervalFraction)
        fmu.cGetIntervalDecimal                    = dlsym(fmu.libHandle, :fmi3GetIntervalDecimal)
        fmu.cGetIntervalFraction                   = dlsym(fmu.libHandle, :fmi3GetIntervalFraction)
        fmu.cGetShiftDecimal                       = dlsym(fmu.libHandle, :fmi3GetShiftDecimal)
        fmu.cGetShiftFraction                      = dlsym(fmu.libHandle, :fmi3GetShiftFraction)
        fmu.cActivateModelPartition                = dlsym(fmu.libHandle, :fmi3ActivateModelPartition)
    end
    # initialize further variables TODO check if needed
    fmu.jac_x = zeros(Float64, fmu.modelDescription.numberOfContinuousStates)
    fmu.jac_t = -1.0
    fmu.jac_dxy_x = zeros(fmi2Real,0,0)
    fmu.jac_dxy_u = zeros(fmi2Real,0,0)

    # dependency matrix
    # fmu.dependencies

    fmu
end

"""
Create a copy of the .fmu file as a .zip folder and unzips it.
Returns the paths to the zipped and unzipped folders.
Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi3Unzip(pathToFMU::String; unpackPath=nothing)

    fileNameExt = basename(pathToFMU)
    (fileName, fileExt) = splitext(fileNameExt)

    if unpackPath === nothing
        # cleanup=true leads to issues with automatic testing on linux server.
        unpackPath = mktempdir(; prefix="fmijl_", cleanup=false)
    end

    zipPath = joinpath(unpackPath, fileName * ".zip")
    unzippedPath = joinpath(unpackPath, fileName)

    # only copy ZIP if not already there
    if !isfile(zipPath)
        cp(pathToFMU, zipPath; force=true)
    end

    @assert isfile(zipPath) ["fmi3Unzip(...): ZIP-Archive couldn't be copied to `$zipPath`."]

    zipAbsPath = isabspath(zipPath) ?  zipPath : joinpath(pwd(), zipPath)
    unzippedAbsPath = isabspath(unzippedPath) ? unzippedPath : joinpath(pwd(), unzippedPath)

    @assert isfile(zipAbsPath) ["fmi3Unzip(...): Can't deploy ZIP-Archive at `$(zipAbsPath)`."]

    numFiles = 0

    # only unzip if not already done
    if !isdir(unzippedAbsPath)
        mkpath(unzippedAbsPath)

        zarchive = ZipFile.Reader(zipAbsPath)
        for f in zarchive.files
            fileAbsPath = normpath(joinpath(unzippedAbsPath, f.name))

            if endswith(f.name,"/") || endswith(f.name,"\\")
                mkpath(fileAbsPath) # mkdir(fileAbsPath)

                @assert isdir(fileAbsPath) ["fmi3Unzip(...): Can't create directory `$(f.name)` at `$(fileAbsPath)`."]
            else
                # create directory if not forced by zip file folder
                mkpath(dirname(fileAbsPath))

                numBytes = write(fileAbsPath, read(f))

                if numBytes == 0
                    @info "fmi3Unzip(...): Written file `$(f.name)`, but file is empty."
                end

                @assert isfile(fileAbsPath) ["fmi3Unzip(...): Can't unzip file `$(f.name)` at `$(fileAbsPath)`."]
                numFiles += 1
            end
        end
        close(zarchive)
    end

    @assert isdir(unzippedAbsPath) ["fmi3Unzip(...): ZIP-Archive couldn't be unzipped at `$(unzippedPath)`."]
    @info "fmi3Unzip(...): Successfully unzipped $numFiles files at `$unzippedAbsPath`."

    (unzippedAbsPath, zipAbsPath)
end

"""
Unload a FMU.
Free the allocated memory, close the binaries and remove temporary zip and unziped FMU model description.
"""
function fmi3Unload(fmu::FMU3, cleanUp::Bool = true)

    while length(fmu.components) > 0
        fmi3FreeInstance!(fmu.components[end])
    end

    dlclose(fmu.libHandle)

    # the components are removed from the component list via call to fmi3FreeInstance!
    @assert length(fmu.components) == 0 "fmi3Unload(...): Failure during deleting components, $(length(fmu.components)) remaining in stack."

    if cleanUp
        try
            rm(fmu.path; recursive = true, force = true)
            rm(fmu.zipPath; recursive = true, force = true)
        catch e
            @warn "Cannot delete unpacked data on disc. Maybe some files are opened in another application."
        end
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable
Create a new instance of the given fmu, adds a logger if logginOn == true.
Returns the instance of a new FMU component.
For more information call ?fmi3InstantiateModelExchange
"""
function fmi3InstantiateModelExchange!(fmu::FMU3; visible::Bool = false, loggingOn::Bool = false)

    ptrLogger = @cfunction(fmi3CallbackLogger, Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Cuint, Ptr{Cchar}))

    compAddr = fmi3InstantiateModelExchange(fmu.cInstantiateModelExchange, fmu.instanceName, fmu.modelDescription.instantiationToken, fmu.fmuResourceLocation, fmi3Boolean(visible), fmi3Boolean(loggingOn), fmu.instanceEnvironment, ptrLogger)

    if compAddr == Ptr{Cvoid}(C_NULL)
        @error "fmi3InstantiateModelExchange!(...): Instantiation failed!"
        return nothing
    end
    previous_z  = zeros(fmi3Float64, fmi3GetEventIndicators(fmu.modelDescription))
    rootsFound  = zeros(fmi3Int32, fmi3GetEventIndicators(fmu.modelDescription))
    stateEvent  = fmi3False
    timeEvent   = fmi3False
    stepEvent   = fmi3False
    component = fmi3Component(compAddr, fmu, previous_z, rootsFound, stateEvent, timeEvent, stepEvent)
    push!(fmu.components, component)
    component
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable
Create a new instance of the given fmu, adds a logger if logginOn == true.
Returns the instance of a new FMU component.
For more information call ?fmi3InstantiateCoSimulation
"""
function fmi3InstantiateCoSimulation!(fmu::FMU3; visible::Bool = false, loggingOn::Bool = false, eventModeUsed::Bool = false, ptrIntermediateUpdate=nothing)

    ptrLogger = @cfunction(fmi3CallbackLogger, Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Cuint, Ptr{Cchar}))
    if ptrIntermediateUpdate === nothing
        ptrIntermediateUpdate = @cfunction(fmi3CallbackIntermediateUpdate, Cvoid, (Ptr{Cvoid}, fmi3Float64, fmi3Boolean, fmi3Boolean, fmi3Boolean, fmi3Boolean, Ptr{fmi3Boolean}, Ptr{fmi3Float64}))
    end
    if fmu.modelDescription.CShasEventMode
        mode = eventModeUsed
    else
        mode = false
    end

    compAddr = fmi3InstantiateCoSimulation(fmu.cInstantiateCoSimulation, fmu.instanceName, fmu.modelDescription.instantiationToken, fmu.fmuResourceLocation, fmi3Boolean(visible), fmi3Boolean(loggingOn),
                                            fmi3Boolean(mode), fmi3Boolean(fmu.modelDescription.CScanReturnEarlyAfterIntermediateUpdate), fmu.modelDescription.intermediateUpdateValueReferences, Csize_t(length(fmu.modelDescription.intermediateUpdateValueReferences)), fmu.instanceEnvironment, ptrLogger, ptrIntermediateUpdate)

    if compAddr == Ptr{Cvoid}(C_NULL)
        @error "fmi3InstantiateCoSimulation!(...): Instantiation failed!"
        return nothing
    end

    component = fmi3Component(compAddr, fmu)
    push!(fmu.components, component)
    component
end

# TODO not tested
"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable
Create a new instance of the given fmu, adds a logger if logginOn == true.
Returns the instance of a new FMU component.
For more information call ?fmi3InstantiateScheduledExecution
"""
function fmi3InstantiateScheduledExecution!(fmu::FMU3, ptrlockPreemption::Ptr{Cvoid}, ptrunlockPreemption::Ptr{Cvoid}; visible::Bool = false, loggingOn::Bool = false)

    ptrLogger = @cfunction(fmi3CallbackLogger, Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Cuint, Ptr{Cchar}))
    ptrClockUpdate = @cfunction(fmi3CallbackClockUpdate, Cvoid, (Ptr{Cvoid}, ))

    compAddr = fmi3InstantiateScheduledExecution(fmu.cInstantiateScheduledExecution, fmu.instanceName, fmu.modelDescription.instantiationToken, fmu.fmuResourceLocation, fmi3Boolean(visible), fmi3Boolean(loggingOn), fmu.instanceEnvironment, ptrLogger, ptrClockUpdate, ptrlockPreemption, ptrunlockPreemption)

    if compAddr == Ptr{Cvoid}(C_NULL)
        @error "fmi3InstantiateScheduledExecution!(...): Instantiation failed!"
        return nothing
    end

    component = fmi3Component(compAddr, fmu)
    push!(fmu.components, component)
    component
end


text:
using HDF5
using Test
using Pkg

@info "libhdf5 v$(HDF5.API.h5_get_libversion())"

@testset "HDF5.jl" begin

include("plain.jl")
include("compound.jl")
include("custom.jl")
include("reference.jl")
include("dataspace.jl")
include("hyperslab.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
include("external.jl")
include("swmr.jl")
include("mmap.jl")
include("properties.jl")
include("table.jl")
include("filter.jl")
include("chunkstorage.jl")
include("fileio.jl")

using MPI
if HDF5.has_parallel()
    # basic MPI tests, for actual parallel tests we need to run in MPI mode
    include("mpio.jl")
end

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: -all
#     formats: jl,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.2
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

using StatsBase
using Serialization
using Plots

# run add "https://github.com/josePereiro/Chemostat" in the julia Pkg REPL for installing the package
import Chemostat
Ch = Chemostat;

# ### Loading cache
# This notebook use the cached data from [toy_model_computation.ipynb](./toy_model_computation.ipynb)

cache_file = "toy_model_cache.jls"
bundle = deserialize(cache_file)
sort!(bundle.ξs)
sort!(bundle.βs)
println(cache_file, " loaded!!!")

println("xis: ", bundle.ξs)
println("betas: ", bundle.βs)

# ---
# ## Plots
# ---

# ### Marginals

# +
# Plot the marginals as function of βs

# An increase β implies that the solutions with greater value of the objective reaction (biom in this case)
# will be more probables (more prob density around the large values of biom). The rest of the fluxes will change
# depending in its covarianze with biom

ξ = bundle.ξs[1] # Select ξ to plot
metnet = bundle[ξ, :net]
fbaout = bundle[ξ, :fba]
ps = []
βs_ = bundle.βs[6:end] # Select βs to plot
βs_ = [βs_; reverse(βs_)]
iders_ = metnet.rxns

function plot1(iders, ξ, β)
    ps = []
    for ider in iders
        ξstr = round(ξ, digits = 1)
        βstr = round(β, digits = 1)
        p = Plots.plot(title = "$(ider)\nxi: $ξstr, beta: $βstr",
            xlabel = "flx", ylabel = "pdf",
            legend = false, yaxis = nothing)

        Ch.Plots.plot_marginal!(p, bundle, ξ, β, [:fba, :hr, :ep], ider)

        pdf_maxval = Ch.Utils.pdf_maxval(bundle, ξ, β, [:ep, :fba], ider)
        Plots.plot!(p, yaxis = [0.0, pdf_maxval * 1.2])
        push!(ps, p)
    end
    p = Plots.plot(framestyle = :none,
        legend = :topleft, legendtitlefont = 15, legendtitle = "Toy model", legendfont = 15)
    p = Plots.plot!(p, [],[], lw = 5, color = :blue, label = "FBA")
    Plots.plot!(p, [],[], lw = 5, color = :red, label = "EP")
    Plots.plot!(p, [],[], lw = 5, color = :red, label = "HR")

    push!(ps, p);

    return Plots.plot(ps..., size = [1000, 1000])
end

gif = Plots.@gif for i in eachindex(βs_)
    plot1(iders_, ξ, βs_[i])
end every 1

# saving gif
gif_file = "toy_model__all_marginals__varying_beta.gif"
cp(gif.filename, gif_file, force = true)
println(relpath(gif_file), " created!!!")
gif
# -

rm(gif.filename, force = true);

# +
# Plot the marginals as function of ξs

# Because increase xi implies a reduction in the intakes all fluxes of the network
# tends to the closest to zero limit possible

β = bundle.βs[5] # Select β to plot
ξs_ = bundle.ξs # Select ξs to plot
ξs_ = [ξs_; reverse(ξs_)]
iders_ = metnet.rxns # Select rxn iders

function plot2(iders, ξ, β, ξ0)
    ps = []

    metnet0 = bundle[ξ0, :net]

    for ider in iders
        ξstr = round(ξ, digits = 1)
        βstr = round(β, digits = 1)

        p = Plots.plot(title = "$(ider)\nxi: $ξstr, beta: $βstr",
            xlabel = "flx", ylabel = "pdf",
            legend = false,
            yaxis = nothing)

        Ch.Plots.plot_marginal!(p, bundle, ξ, β, [:fba, :hr, :ep], ider)

        lb_, ub_ = Ch.Utils.bounds(metnet0, ider)
        margin_ = abs(ub_ - lb_) * 0.1
        pdf_maxval = Ch.Utils.pdf_maxval(bundle, ξ, β, [:ep, :fba], ider)
        Plots.plot!(p,
            xaxis = [lb_ - margin_, ub_ + margin_],
            yaxis = [0.0, pdf_maxval * 1.2])

        push!(ps, p)
    end
    p = Plots.plot(framestyle = :none,
        legend = :topleft, legendtitlefont = 15, legendtitle = "Toy model", legendfont = 15)
    p = Plots.plot!(p, [],[], lw = 5, color = :blue, label = "FBA")
    Plots.plot!(p, [],[], lw = 5, color = :red, label = "EP")
    Plots.plot!(p, [],[], lw = 5, color = :black, label = "HR")

    push!(ps, p);

    return Plots.plot(ps..., size = [1000, 1000])
end

# gif
gif = Plots.@gif for i in eachindex(ξs_)
    plot2(iders_, ξs_[i], β, minimum(ξs_))
end every 1

# saving gif
gif_file = "toy_model__all_marginals__varying_xi.gif"
cp(gif.filename, gif_file, force = true)
println(relpath(gif_file), " created!!!")
gif
# -

rm(gif.filename, force = true);

# +
# Plot stoi err asfunctio of xi
β = bundle.βs[1] # Select β to plot

ξs_ = bundle.ξs#[[1, 10, 30, 40, 50]] # Select ξs to plot

p = Plots.plot(title = "Toy Model\nbeta: $β",
    xlabel = "xi", ylabel = "stoi err/ flx",
)

iders_ = metnet.mets # Select iders to plot
for ider in iders_
    errs_ = [Ch.Utils.stoi_err(bundle, ξ, β, :ep, ider) for ξ in ξs_]
    Plots.plot!(p, ξs_, errs_, label = "", color= :black)
end
Plots.plot!(p, [], [], color = :black, label = "abs stoi err", lw = 3)

flxs_ = [abs.(Ch.Utils.av(bundle, ξ, β, :ep)) for ξ in ξs_]
Plots.plot!(p, ξs_, mean.(flxs_), lw = 3, label = "", color = :white)
Plots.plot!(p, ξs_, mean.(flxs_), ls = :dash, lw = 3, label = "mean abs flx", color = :blue)
Plots.plot!(p, ξs_, minimum.(flxs_), lw = 3, label = "", color = :white)
Plots.plot!(p, ξs_, minimum.(flxs_), ls = :dash, lw = 3, label = "min abs flx", color = :red)

# gif (making a zoom)
yulims_ = Ch.Utils.logspace(-1,1, 100) |> reverse
yulims_ = [fill(maximum(yulims_), 25); yulims_]
yulims_ = [reverse(yulims_); yulims_]
yllim_ = -0.01
gif = Plots.@gif for i in eachindex(yulims_)
    iter_ = yllim_:((yulims_[i] - yllim_)/8):yulims_[i]
    Plots.plot!(p, yaxis = [yllim_, yulims_[i]], yticks = (iter_,
            string.(round.(collect(iter_), digits = 2))))
end every 1

# saving gif
gif_file = "toy_model__stoi_err_vs_xi.gif"
cp(gif.filename, gif_file, force = true)
println(relpath(gif_file), " created!!!")
gif
# -
rm(gif.filename, force = true);

# +
# Plot stoi err asfunctio of xi
ξ = bundle.ξs[10] # Select ξ to plot
ξstr = round(ξ, digits = 2)
βs_ = bundle.βs[2:end] # Select βs to plot

p = Plots.plot(title = "Toy Model\nxi: $ξstr",
    xlabel = "beta", ylabel = "stoi err",
    xaxis = :log10,
#     yaxis = [-1e-3, 1e-3]
)

iders_ = metnet.mets # Select iders to plot
for ider in iders_
    errs_ = [Ch.Utils.stoi_err(bundle, ξ, β, :ep, ider) for β in βs_]
    Plots.plot!(p, βs_, errs_, color = :black, label = "")
end

Plots.plot!(p, [], [], color = :black, label = "abs stoi err")

flxs_ = [abs.(Ch.Utils.av(bundle, ξ, β, :ep)) for β in βs_]
Plots.plot!(p, βs_, mean.(flxs_), lw = 3, label = "", color = :white)
Plots.plot!(p, βs_, mean.(flxs_), ls = :dash, lw = 3, label = "mean abs flx", color = :blue)
Plots.plot!(p, βs_, minimum.(flxs_), lw = 3, label = "", color = :white)
Plots.plot!(p, βs_, minimum.(flxs_), ls = :dash, lw = 3, label = "min abs flx", color = :red)

# gif
yulims_ = Ch.Utils.logspace(-0.5,0.9, 100) |> reverse
yulims_ = [fill(maximum(yulims_), 25); yulims_]
yulims_ = [reverse(yulims_); yulims_]
yllim_ = -0.01
gif = Plots.@gif for i in eachindex(yulims_)
    iter_ = yllim_:((yulims_[i] - yllim_)/8):yulims_[i]
    Plots.plot!(p, yaxis = [yllim_, yulims_[i]], yticks = (iter_,
            string.(round.(collect(iter_), digits = 2))))
end every 1

# saving gif
gif_file = "toy_model__stoi_err_vs_beta.gif"
cp(gif.filename, gif_file, force = true)
println(relpath(gif_file), " created!!!")
gif
# -
rm(gif.filename, force = true);

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat


text:
function (semmg)(par)
    for i in 1:length(semmg.vec)
        F += semmg.vec[i](par)
    end
end

function (semmg)(par)
    for i in 1:length(semmg.imply_vec)
        semmg.imply_vec[i](par, semmg)
    end
    for i in 1:length(semmg.loss_vec)
        F += semmg.loss_vec[i](par, semmg)
    end
end

function (semmg)(par)
    semmg.imply_column(par, semmg)
    F = semmg.loss_column(par, semmg)
end

abstract type loss_column end

struct ML <: loss_column
    vec::Vector{loss}
end

(ML)(par)

struct ML_superparallel <: loss_column

end

struct WLS_superparallel <: loss_column

end

function (WLS_superparallel)(par, model)
    for i in 1:length(imply_column)
        A[i] = model.imply_column[i].imp_cov
    end
    do A end
end

sem_fit(model)
sem_fit_superparallel_ML(model)
sem_fit_superparallel_WLS(model)
ML_superparallel(model, par)

function (semdef)(par)
    for i in 1:length(semdef.vec)
        F += semdef.vec[i](par, semdef.defvars[i])
    end
end

function (semfiml)(par)
    for i in 1:length(semfiml.vec)
        F += semfiml.vec[i](par, semfiml.filter[i])
    end
end

function (semfiml)(par, defvars)
    imply(par)
    delete(semfiml.imply)
    for loss(par, semfiml) end

    for i in 1:length(semfiml.vec)
        F += semfiml.vec[i](par, semfiml.filter[i], defvars)
    end
end

function (semel)(par)
    semel.imply(par, semel)
    F = semel.loss(par, semel)
    return F
end

function (semel)(par, filter)
    semel.imply(par, semel)
    F = semel.loss(par, semel, filter)
    return F
end

function (sematom)(par, defvars)
    semel.imply(par, semel, defvars)
    F = semel.loss(par, sematom)
    return F
end

function (semel)(par, defvars, filter)
    semel.imply(par, semel, defvars)
    F = semel.loss(par, semel, filter)
    return F
end

test

Threads.nthreads()

A = rand(100,100)

using BenchmarkTools, LinearAlgebra

@benchmark inv(A)

BLAS.get_config()

using MKL

@benchmark inv(A)

BLAS.get_config()


text:
module PencilFFTs

import AbstractFFTs
import FFTW
import MPI

using LinearAlgebra
using Reexport
using TimerOutputs

@reexport using PencilArrays

include("Transforms/Transforms.jl")
@reexport using .Transforms

import PencilArrays.Transpositions: AbstractTransposeMethod
import .Transforms: AbstractTransform, FFTReal, scale_factor

export PencilFFTPlan
export allocate_input, allocate_output, scale_factor

# Deprecated in v0.10
@deprecate get_scale_factor scale_factor

# Functions to be extended for PencilFFTs types.
import PencilArrays: get_comm, timer, topology, extra_dims

const AbstractTransformList{N} = NTuple{N, AbstractTransform} where N

include("global_params.jl")
include("plans.jl")
include("allocate.jl")
include("operations.jl")

end # module


text:
#MaxSens
using HullReachNV
using Test
using LazySets
import HullReachNV: ReLU

nnet = read_nnet("nnet/toy_nnet.nnet", last_layer_activation = ReLU())

solver = MaxSens(0.2, true)

center = [1.0, 1.0]
radius = [1.0, 1.0]

in_hyper = Hyperrectangle(center, radius)

lower = [0.0, 0.0]
upper = [1.2, 1.7]

out_hyper = Hyperrectangle(low=lower, high=upper)

problem = Problem(nnet, in_hyper, out_hyper)

print("MaxSens - test")
timed_result =@timed solve(solver, problem)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1].status)
print("\n")


text:
"""
    Farneback(Args...)

A method for dense optical flow estimation developed by Gunnar Farneback. It
computes the optical flow for all the points in the frame using the polynomial
representation of the images. The idea of polynomial expansion is to approximate
the neighbourhood of a point in a 2D function with a polynomial. Displacement
fields are estimated from the polynomial coefficients depending on how the
polynomial transforms under translation.

The different arguments are:

 -  flow_est          =  Array of SVector{2} containing estimate flow values for all points in the frame
 -  iterations        =  Number of iterations the displacement estimation algorithm is run at each
                         point
 -  window_size       =  Size of the search window at each pyramid level; the total size of the
                         window used is 2*window_size + 1
 -  σw                =  Standard deviation of the Gaussian weighting filter
 -  neighbourhood     =  size of the pixel neighbourhood used to find polynomial expansion for each pixel;
                         larger values mean that the image will be approximated with smoother surfaces,
                         yielding more robust algorithm and more blurred motion field
 -  σp                =  standard deviation of the Gaussian that is used to smooth derivatives used as a
                         basis for the polynomial expansion (Applicability)
 -  est_flag          =  true -> Use flow_est as initial estimate
                         false -> Assume zero initial flow values
 -  gauss_flag        =  false -> use box filter
                         true -> use gaussian filter instead of box filter of the same size for optical flow
                         estimation; usually, this option gives more accurate flow than with a box filter,
                         at the cost of lower speed (Default Value)

## References

Farnebäck G. (2003) Two-Frame Motion Estimation Based on Polynomial Expansion. In: Bigun J.,
Gustavsson T. (eds) Image Analysis. SCIA 2003. Lecture Notes in Computer Science, vol 2749. Springer, Berlin,
Heidelberg

Farnebäck, G.: Polynomial Expansion for Orientation and Motion Estimation. PhD thesis, Linköping University,
Sweden, SE-581 83 Linköping, Sweden (2002) Dissertation No 790, ISBN 91-7373-475-6.

"""

struct Farneback{F <: Float64, I <: Int, V <: Bool} <: OpticalFlowAlgo
    flow_est::Array{SVector{2, F}, 2}
    iterations::I
    window_size::I
    σw::F
    neighbourhood::I
    σp::F
    est_flag::V
    gauss_flag::V
end

Farneback(flow_est::Array{SVector{2, F}, 2}, iterations::I, window_size::I, σw::F, neighbourhood::I, σp::F, est_flag::V) where {F <: Float64, I <: Int, V <: Bool} = LK{F, I, V}(flow_est, iterations, window_size, σw, neighbourhood, σp, est_flag, true)

function optflow(first_img::AbstractArray{T, 2}, second_img::AbstractArray{T,2}, algo::Farneback{}) where T <: Gray
    if algo.est_flag
        @assert size(first_img) == size(algo.flow_est)
    end
    A1, B1, _ = polynomial_expansion(first_img, algo.neighbourhood, algo.σp)
    A2, B2, _ = polynomial_expansion(second_img, algo.neighbourhood, algo.σp)

    AtA = similar(A1)
    AtB = similar(B1)

    disp = Array{SVector{2, Float64}, 2}(size(first_img))
    for i = 1:size(first_img)[1]
        for j = 1:size(first_img)[2]
            disp[i,j] = SVector{2, Float64}(0.0,0.0)
        end
    end

    for k = 1:algo.iterations
        for i = 1:size(first_img)[1]
            for j = 1:size(first_img)[2]
                d = disp[i,j]

                mod_i = round(Int, i + d[2])
                mod_j = round(Int, j + d[1])

                mod_i = clamp(mod_i, 1, size(first_img)[1])
                mod_j = clamp(mod_j, 1, size(first_img)[2])

                sub_A1 = SMatrix{2,2,Float64}(@view A1[i,j,:,:])
                sub_B1 = SVector{2,Float64}(@view B1[i,j,:])
                sub_A2 = SMatrix{2,2,Float64}(@view A2[mod_i,mod_j,:,:])
                sub_B2 = SVector{2,Float64}(@view B2[mod_i,mod_j,:])

                A = (sub_A1 + sub_A2)/2
                dB = A*d - ((sub_B2 - sub_B1)/2)

                AtA[i,j,:,:] .= A'*A
                AtB[i,j,:] .= A'*dB
            end
        end

        len = floor(Int, algo.window_size/2)
        w = ImageFiltering.KernelFactors.gaussian((algo.σw),(algo.window_size))[-len:len]

        for i = 1:2
            for j = 1:2
                if algo.gauss_flag
                    AtA[:,:,i,j] = @view conv2(w, w, AtA[:,:,i,j])[1+len:size(first_img)[1]+len,1+len:size(first_img)[2]+len]
                else
                    AtA[:,:,i,j] =  @view conv2(ones(algo.window_size)./algo.window_size, ones(algo.window_size)./algo.window_size, AtA[:,:,i,j])[1+len:size(first_img)[1]+len,1+len:size(first_img)[2]+len]
                end
            end
            if algo.gauss_flag
                AtB[:,:,i] = @view conv2(w, w, AtB[:,:,i])[1+len:size(first_img)[1]+len,1+len:size(first_img)[2]+len]
            else
                AtB[:,:,i,j] = @view conv2(ones(algo.window_size)./algo.window_size, ones(algo.window_size)./algo.window_size, AtB[:,:,i,j])[1+len:size(first_img)[1]+len,1+len:size(first_img)[2]+len]
            end
        end

        for i = 1:size(first_img)[1]
            for j = 1:size(first_img)[2]
                disp[i,j] = pinv2x2(SMatrix{2,2,Float64}(AtA[i,j,:,:]))*SVector{2,Float64}(AtB[i,j,:])
            end
        end
    end

    return disp

end

"""
```
A, B, C = polynomial_expansion(img, neighbourhood, σ)
```

Returns the polynomial coefficients of the approximation of the neighbourhood of a
point in a 2D function with a polynomial.
The expansion is: f(x) = x'Ax + B'x + C

Parameters:

-  img               =  Grayscale image who polynomial expansion is to be found
-  neighbourhood     =  Size of the pixel neighbourhood used to find polynomial expansion for each pixel;
                        larger values mean that the image will be approximated with smoother surfaces,
                        yielding more robust algorithm and more blurred motion field
-  σ                 =  Standard deviation of the Gaussian that is used to smooth derivatives used as a
                        basis for the polynomial expansion (Applicability)

## References

Farnebäck, G.: Polynomial Expansion for Orientation and Motion Estimation. PhD thesis, Linköping University,
Sweden, SE-581 83 Linköping, Sweden (2002) Dissertation No 790, ISBN 91-7373-475-6.

"""

function polynomial_expansion(img::AbstractArray{T, 2}, neighbourhood::Int, σ::Real) where T <: Gray
    len = floor(Int, neighbourhood/2)
    a = ImageFiltering.Kernel.gaussian((σ,σ), (neighbourhood,neighbourhood))[-len:len, -len:len]

    c = ones(size(img))
    padded_c = parent(padarray(c, Fill(0, (len, len))))
    padded_img = parent(padarray(img, Fill(0, (len, len))))

    B0 = zeros(neighbourhood*neighbourhood, 6)
    for x = -len:len
        for y = -len:len
            B0[neighbourhood*(x + len) + (y + len) + 1, :] = [1, x, y, x*x, y*y, x*y]
        end
    end

    B = SMatrix{neighbourhood*neighbourhood, 6}(B0)
    Wa = SDiagonal{neighbourhood^2,Float64}( @view a[:])
    BtWa = B'*Wa

    # If we assume that the confidence is constant we can move several
    # matrix multiplications out of the main loop.
    Wc = SDiagonal{neighbourhood^2,Float64}(ones(1,neighbourhood^2)...)
    BtWaWc = BtWa*Wc
    G = BtWaWc*B
    G_inv = SMatrix{6,6,Float64}(pinv(convert(Array,G)))
    G_invBtWaWc = G_inv*BtWaWc

    poly_A = zeros(size(img)[1], size(img)[2], 2, 2)
    poly_B = zeros(size(img)[1], size(img)[2], 2)
    poly_C = zeros(size(img))

    for i = (len + 1):(size(img)[1] - len)
        for j = (len + 1):(size(img)[2] - len)
            sub_img = @view padded_img[(i - len):(i + len), (j - len):(j + len)]
            r = G_invBtWaWc*SVector{neighbourhood^2,Float64}(@view sub_img[:])

            poly_A[i-len, j-len, 1, 1] =  r[4]
            poly_A[i-len, j-len, 2, 1] =  r[6]/2
            poly_A[i-len, j-len, 1, 2] =  r[6]/2
            poly_A[i-len, j-len, 2, 2] =  r[5]

            poly_B[i-len, j-len, 1] = r[2]
            poly_B[i-len, j-len, 2] = r[3]

            poly_C[i-len, j-len] = r[1]
        end
    end

    return (poly_A, poly_B, poly_C)

end


text:
#=
function q(len::Int = 7, ynb::Int = 2, horiz::Int = 100, M::Int = 500)

end
=#

"""
    monte_carlo_return_period(sepp::SEPP, r::Integer = 7, horiz::Integer = 100, M::Integer = 500)::Real

Compute the retrun period of an event of *r* days by a monte calro method.

The computation is based on *M* simulation to the horizon *horiz*.
If the model is an SEMPP and magnitude is given by the user than the method only consider events above *magnitude*.
"""
function monte_carlo_return_period end

function monte_carlo_return_period(sepp::SEPP; r::Integer = 7, horiz::Integer = 100, M::Integer = 500)::Real
    sims = TS[]

    for i in 1:M
        push!(sims, discrete_simulation(sepp; start_time = 1, end_time = horiz))
    end

    warn = 0

    function count(ts::TS)::Real
        times = copy(ts.times)
        res = 0
        c = 1

        for i in 2:length(times)
            if c == 0 || times[i-1] == times[i] - 1
                c += 1
            else
                c >= r && (res += 1)
                c = 0
            end
        end

        c >= r && (res += 1)

        res == 0 && (res = 1 ; warn += 1)

        return horiz/res
    end

    ret_per = median(count.(sims))

    warn == 0 || @warn "no events within horizon for $warn simulation(s)"

    return ret_per
end

function monte_carlo_return_period(sempp::SEMPPExpKern; magnitude::Real, r::Integer = 7, horiz::Integer = 100, M::Integer = 500)::Real
    sims = MarkedTimeSeries[]

    for i in 1:M
        push!(sims, discrete_simulation(sempp; start_time = 1, end_time = horiz))
    end

    warn = 0

    function count(mts::MarkedTimeSeries)::Real
        times = copy(mts.times)
        marks = copy(mts.marks)

        res = 0
        c = 0
        marks[1] >= magnitude && (c = 1)

        for i in 2:length(times)
            if (marks[i] >= magnitude) && (c == 0 || times[i-1] == times[i] - 1)
                c+=1
            else
                c >= r && (res += 1)
                c=0
            end
        end

        c >= r && (res += 1)

        res == 0 && (res = 1 ; warn += 1)

        return horiz/res
    end

    ret_per = median(count.(sims))

    warn == 0 || @warn "no events within horizon for $warn simulation(s)"

    return ret_per
end

function rQy(sempp::SEMPPExpKern, r::Integer = 7, y::Real = 2; horiz::Union{Real, Nothing} = nothing, M::Integer = 500)
    sims = MarkedTimeSeries[]
    horiz = isnothing(horiz) ? 5*y*365 : Int(horiz*365)

    mag_min = Inf
    mag_max = -Inf

    for i in 1:M
        sim = discrete_simulation(sempp; start_time = 1, end_time = horiz)
        mag_min = minimum(vcat(sim.marks, mag_min))
        mag_max = maximum(vcat(sim.marks, mag_max))
        push!(sims, sim)
    end

    scale = mag_max - mag_min

    function count(mts::MarkedTimeSeries, magnitude::Real)::Real
        times = copy(mts.times)
        marks = copy(mts.marks)

        res = 0
        c = 0

        if !isempty(times)
            marks[1] >= magnitude && (c = 1)

            for i in 2:length(times)
                if (marks[i] >= magnitude) && (c == 0 || times[i-1] == times[i] - 1)
                    c+=1
                else
                    c >= r && (res += 1)
                    c=0
                end
            end
        end

        c >= r && (res += 1)

        res == 0 && (res = 1)

        return horiz/res
    end

    ret_per_min = median((x -> count(x, mag_min)).(sims))
    ret_per_max = median((x -> count(x, mag_max)).(sims))

    while (abs(ret_per_max - ret_per_min) >= 1) && ((mag_max-mag_min)/scale > 0.001)
        mag = 0.5 * (mag_max + mag_min)
        ret_per = median((x -> count(x, mag)).(sims))

        if ret_per > y*365
            mag_max = mag
            ret_per_max = ret_per
        else
            mag_min = mag
            ret_per_min = ret_per
        end
    end

    return 0.5*(mag_min + mag_max), ret_per_min
end

function discrete_forecast(sepp::SEPPExpKern; start_time::Integer = 1, end_time::Integer = 1000, history::Union{TimeSeries, Nothing} = nothing, M::Integer = 500)
    sims = TimeSeries[]

    if !isnothing(history)
        (first(history.times) isa TimeType) && (history = TimeSeries(Dates.value.(Date.(history.times))))
    end

    last_h = last(history.times)

    for i in 1:M
        push!(sims, discrete_simulation(sepp, start_time = start_time, end_time = end_time, history_time_series = history))
    end

    p = zeros(Float64, end_time - start_time + 1)

    real_start = last_h + start_time

    for ts in sims
        for t in ts.times
            p[t - real_start + 1] += 1
        end
    end

    p = p ./ M

    return p
end

function discrete_forecast(sempp::SEMPPExpKern; start_time::Integer = 1, end_time::Integer = 1000, history::Union{MarkedTimeSeries, Nothing} = nothing, magnitudes::AbstractVector{<:Real} = [], M::Integer = 500)
    sims = MarkedTimeSeries[]

    if !isnothing(history)
        (first(history.times) isa TimeType) && (history = MarkedTimeSeries(Dates.value.(Date.(history.times)), history.marks))
    end

    last_h = last(history.times)

    for i in 1:M
        push!(sims, discrete_simulation(sempp, start_time = start_time, end_time = end_time, history_marked_time_series = history))
    end

    p = zeros(Float64, end_time - start_time + 1, length(magnitudes) + 1 )

    real_start = last_h + start_time

    for mts in sims
        extended_marks = fill(-Inf, end_time - start_time + 1)
        extended_marks[mts.times .- real_start .+ 1] = mts.marks

        for i in 1:(end_time - start_time + 1)
            m_idx = sum(magnitudes .<= extended_marks[i])
            p[i, m_idx + 1] += 1
        end
    end

    p = p ./ M

    return p
end