text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module HumanEye

using OpticSim
using OpticSim.GlassCat
using OpticSim.Geometry
import Unitful
using Unitful.DefaultSymbols:mm,°

# Approximate average properties of the human eye

"""
# Pupil diameter as a function of scene luminance
https://jov.arvojournals.org/article.aspx?articleid=2279420
https://en.wikipedia.org/wiki/Orders_of_magnitude_(luminance)

Pupil diameter is approximately 2.8mm at 100cd/m^2. A typical overcast day is 700cd/m^2
"""

"""computes pupil diameter as a function of scene luminance `L`, in cd/m², and the angular area, `a`, over which this luminance is presented to the eye."""
𝐃sd(L,a) = 7.75 - 5.75 * ((L * a / 846)^.41) / ((L * a / 846)^.41 + 2) # the first letter of this function name is \bfD not D.

eyeradius() = 12mm
export eyeradius

"""Posterior focal length, i.e., optical distance from entrance pupil to the retina. Focal length will change depending on accomodation. This value is for focus at ∞. When the eye is focused at 25cm focal length will be ≈ 22mm. Because the index of refraction of the vitreous humor is approximately 1.33 the physical distance from the entrance pupil to the retina will be 24mm/1.33 = 18mm."""
eyefocallength() = 24mm
export eyefocallength

vc_epupil() = 3mm #distance from vertex of cornea to entrance pupil

""" distance from vertex of cornea to center of rotation"""
cornea_to_eyecenter() = 13.5mm
export cornea_to_eyecenter

entrancepupil_to_retina() = 22.9mm

"""distance from entrance pupil to center of rotation."""
entrancepupil_to_eyecenter() = entrancepupil_to_retina() - eyeradius()
export entrancepupil_to_eyecenter

"""average angle, in degrees, the eye will rotate before users will turn their head"""
comfortable_eye_rotation_angle() = 20°

"""average translation of the entrance pupil associated with comfortable eye rotation"""
comfortable_entrance_pupil_translation() = sin(comfortable_eye_rotation_angle())*entrancepupil_to_eyecenter()
export comfortable_entrance_pupil_translation

#Model eyes of varying degrees of verisimilitude

"""
    ModelEye(assembly::LensAssembly{T}, nsamples::Int = 17; pupil_radius::T = 3.0, detpixels::Int = 1000, transform::Transform{T} = identitytransform(T))

Geometrically accurate model of the human eye focused at infinity with variable `pupil_radius`.
The eye is added to the provided `assembly` to create a [`CSGOpticalSystem`](@ref) with the retina of the eye as the detector.

The eye can be positioned in the scene using the `transform` argument and the resolution of the detector specified with `detpixels`.
By default the eye is directed along the positive z-axis with the vertex of the cornea at the origin.

`nsamples` determines the resolution at which accelerated surfaces within the eye are triangulated.
"""
function ModelEye(assembly::LensAssembly{T}; nsamples::Int = 17, pupil_radius::T = 3.0, detpixels::Int = 1000, transform::Transform{T} = identitytransform(T), detdatatype::Type{D} = Float32) where {T<:Real,D<:Real}
    cornea_front = AcceleratedParametricSurface(ZernikeSurface(6.9, radius = -7.8, conic = -0.5), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.CORNEA, OpticSim.GlassCat.Air))
    cornea_rear = Plane(SVector{3,T}(0, 0, -1), SVector{3,T}(0, 0, -3.2))
    cornea = (cornea_front ∩ cornea_rear)(transform)

    anterior_chamber_front = AcceleratedParametricSurface(ZernikeSurface(6.0, radius = -6.7, conic = -0.3), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.AQUEOUS, OpticSim.GlassCat.EYE.CORNEA))
    anterior_chamber_rear = Plane(SVector{3,T}(0, 0, -1), SVector(0.0, 0.0, -3.2), vishalfsizeu = 6.0, vishalfsizev = 6.0, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.AQUEOUS, OpticSim.GlassCat.EYE.VITREOUS))
    anterior_chamber = leaf(anterior_chamber_front ∩ anterior_chamber_rear, translation(0.0, 0.0, -0.52))(transform)

    pupil = Annulus(pupil_radius, 5.9, rotate(transform, SVector{3,T}(0, 0, 1)), transform * SVector{3,T}(0.0, 0.0, -3.62))

    lens_front = AcceleratedParametricSurface(ZernikeSurface(5.1, radius = -10.0), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_rear = AcceleratedParametricSurface(ZernikeSurface(5.1, radius = 6.0, conic = -3.25), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_barrel = Cylinder(5.0, 5.0, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_csg = (lens_front - leaf(lens_rear, translation(0.0, 0.0, -3.7))) ∩ leaf(lens_barrel, translation(0.0, 0.0, -2.0))
    lens = leaf(lens_csg, translation(0.0, 0.0, -3.72))(transform)

    vitreous_chamber_csg = (
        leaf(
            Sphere(11.0, interface = FresnelInterface{T}(EYE.VITREOUS, Air, reflectance = 0.0, transmission = 0.0)),
            translation(0.0, 0.0, -13.138998863513297)
        ) ∩
        Plane(0.0, 0.0, 1.0, 0.0, 0.0, -3.72, interface = FresnelInterface{T}(EYE.VITREOUS, EYE.AQUEOUS)) ∩
        Plane(0.0, 0.0, -1.0, 0.0, 0.0, -24 + (11 - sqrt(11^2 - 10^2)), interface = NullInterface(T))
    )
    vitreous_chamber = vitreous_chamber_csg(transform)

    new_assembly = OpticSim.LensAssembly(assembly, cornea, anterior_chamber, pupil, vitreous_chamber, lens)
    retina = SphericalCap(
        11.0,
        asin(10 / 11),
        rotate(transform, SVector{3,T}(0, 0, 1)),
        transform * SVector{3,T}(0.0, 0.0, -24.0),
        interface = FresnelInterface{T}(EYE.VITREOUS, EYE.VITREOUS, reflectance = zero(T), transmission = zero(T))
    )

    return CSGOpticalSystem(new_assembly, retina, detpixels, detpixels, D)
end
export ModelEye

#! format: off
"""
    ArizonaEye(::Type{T} = Float64; accommodation::T = 0.0)

The popular Arizona eye model taken from [this definition](https://photonengr.com/wp-content/uploads/kbasefiles/ArizonaEyeModel.pdf).
The `accommodation` of the eye can be varied in this model.
Returns a `DataFrame` specifying the prescription of the eye model.
"""
function ArizonaEye(::Type{T} = Float64; accommodation::T = 0.0) where {T<:Real}
    # from https://photonengr.com/wp-content/uploads/kbasefiles/ArizonaEyeModel.pdf also in our documentation/papers directory
    return DataFrame(
        Name = ["Air", "Cornea", "Aqueous", "Lens", "Vitreous", "Retina"],
        Surface = ["Object", "Standard", "Standard", "Standard", "Standard", "Image"],
        Radius = [Inf64, 7.8, 6.5, 12.0 - 0.4*accommodation, -5.22 + 0.2*accommodation, -13.4],
        Conic = [missing, -.25, -.25, -7.52 + 1.29*accommodation, -1.35 - 0.43*accommodation, missing],
        Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.modelglass(1.377, 57.1, 0.0), OpticSim.GlassCat.modelglass(1.337, 61.3, 0.0), OpticSim.GlassCat.modelglass(1.42 + 0.0026 * accommodation - 0.00022 * accommodation^2, 51.9, 0.0), OpticSim.GlassCat.modelglass(1.336, 61.1, 0.0), missing],
        SemiDiameter = [Inf64, 5.3, 5.3, 5.3, 5.3, 5.3],
        Thickness = [Inf64, 0.55, 2.97 - 0.04*accommodation, 3.77 + 0.04accommodation, 16.713, missing],
    )
end
export ArizonaEye

end #module
export HumanEye
#! format: on


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "NativeFileDialog"
version = v"1.1.6"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mlabbe/nativefiledialog.git", "67345b80ebb429ecc2aeda94c478b3bcc5f7888e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nativefiledialog/src/
if [[ "${target}" == *-mingw* ]]; then
    c++ nfd_common.c nfd_win.cpp -lole32 -luuid -Iinclude -O2 -Wall -Wextra -fno-exceptions -shared -o "${libdir}/libnfd.${dlext}"
elif [[ "${target}" == *-apple* ]]; then
    # https://github.com/mlabbe/nativefiledialog/compare/mlabbe:67345b8...mlabbe:eb29acc
    atomic_patch -p2 ../../patches/macos-focus-bugfix.patch
    cc nfd_cocoa.m nfd_common.c -framework Foundation -framework AppKit -Iinclude -O2 -Wall -Wextra -fno-exceptions -fPIC -shared -o "${libdir}/libnfd.${dlext}"
else
    cc nfd_common.c nfd_gtk.c -Iinclude `pkg-config --cflags --libs gtk+-3.0` -O2 -Wall -Wextra -fno-exceptions -fPIC -shared -o "${libdir}/libnfd.${dlext}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libnfd", :libnfd)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GTK3_jll", uuid="77ec8976-b24b-556a-a1bf-49a033a670a6"))
    BuildDependency("Xorg_xorgproto_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

ext: jl

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
# Adapted from Transformers.jl and HTTP.jl, both released under the MIT License.
# Original files:
# https://github.com/chengchingwen/Transformers.jl/blob/27a736794487f303d7919ac03623c9d34670f744/src/datasets/download_utils.jl
# https://github.com/JuliaWeb/HTTP.jl/blob/ed02d37b373386804707262d68856c0fad329b16/src/download.jl

using HTTP
using Dates
using DataDeps
using Random:randstring

const google_base_domain = "docs.google.com"

function fetch_google_drive_id(drive_id, local_path)
    url = "https://$google_base_domain/uc?id=$drive_id&export=download"
    fetch_google_drive(url, local_path)
end

function find_google_confirmation_code(cookiejar::Dict{String,Set{HTTP.Cookies.Cookie}})
    cookies = cookiejar[google_base_domain]
    cookies = filter(
        cookie -> startswith(cookie.name, "download_warning_"),
        cookies
    )
    first(cookies).value
end

function confirm_google_download(url, cookiejar)
    HTTP.request("HEAD", url; cookies=true, cookiejar=cookiejar)
    code = find_google_confirmation_code(cookiejar)
    if isnothing(code)
        error("download failed")
    end
    "$url&confirm=$code"
end

function fetch_google_drive(url, localdir)
    cookiejar = Dict{String,Set{HTTP.Cookies.Cookie}}()
    url = confirm_google_download(url, cookiejar)

    format_progress(x) = round(x, digits=4)
    format_bytes(x) = !isfinite(x) ? "∞ B" : Base.format_bytes(x)
    format_seconds(x) = "$(round(x; digits=2)) s"
    format_bytes_per_second(x) = format_bytes(x) * "/s"

    filepaths::Array{String} = []

    HTTP.open(
        "GET",
        url,
        ["Range" => "bytes=0-"];
        cookies=true,
        cookiejar=cookiejar,
        redirect_limit=10
    ) do stream
        response = HTTP.startread(stream)
        content_disposition = HTTP.header(response, "Content-Disposition")

        filename_match = match(r"filename=\\\"(.*)\\\"", content_disposition)
        filename = filename_match !== nothing ? filename_match.captures[] : basename(tempname())
        filepath = joinpath(localdir, filename)

        parsed_total_bytes = tryparse(
            Float64,
            split(
                HTTP.header(response, "Content-Range"),
                '/'
            )[end]
        )
        total_bytes = parsed_total_bytes !== nothing ? parsed_total_bytes : NaN
        downloaded_bytes = 0
        start_time::DateTime = Dates.now()
        prev_time::DateTime = Dates.now()
        period = DataDeps.progress_update_period()

        function report_callback()
            if filename_match === nothing && isnan(total_bytes)
                # Skip reporting intermediate redirects.
                return
            end

            prev_time = Dates.now()
            taken_time = (prev_time - start_time).value / 1000
            average_speed = downloaded_bytes / taken_time
            remaining_bytes = total_bytes - downloaded_bytes
            remaining_time = remaining_bytes / average_speed
            completion_progress = downloaded_bytes / total_bytes

            @info(
                "Downloading",
                source = url,
                dest = filepath,
                progress = completion_progress |> format_progress,
                time_taken = taken_time |> format_seconds,
                time_remaining = remaining_time |> format_seconds,
                average_speed = average_speed |> format_bytes_per_second,
                downloaded = downloaded_bytes |> format_bytes,
                remaining = remaining_bytes |> format_bytes,
                total = total_bytes |> format_bytes,
            )
        end

        Base.open(filepath, "w") do fh
            while (!eof(stream))
                downloaded_bytes += write(fh, readavailable(stream))
                if !isinf(period)
                    if Dates.now() - prev_time > Millisecond(1000 * period)
                        report_callback()
                    end
                end
            end
        end
        if !isinf(period)
            report_callback()
        end

        push!(filepaths, filepath)
    end

    final_filepath = pop!(filepaths)
    for path ∈ filepaths
        # Remove intermediate redirects.
        rm(path)
    end
    final_filepath
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

find maximum volume hypercube with edges parallel to the axes inside a polyhedron
or an ellipsoid defined with l_1, l_infty, or l_2 ball constraints
=#

struct MaxVolumeJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int
    epinormeucl_constr::Bool # add an L2 ball constraint
    epinorminf_constrs::Bool # add L1 and Linfty ball constraints
end

function build(inst::MaxVolumeJuMP{T}) where {T <: Float64}
    @assert xor(inst.epinormeucl_constr, inst.epinorminf_constrs)
    n = inst.n

    A = randn(n, n)
    # ensure there will be a feasible solution
    x = randn(n)
    A = A * A' + 10I
    gamma = norm(A * x) / sqrt(n)

    model = JuMP.Model()
    JuMP.@variable(model, t)
    JuMP.@variable(model, x[1:n])
    JuMP.@objective(model, Max, t)
    JuMP.@constraint(model, vcat(t, x) in MOI.GeometricMeanCone(n + 1))

    if inst.epinormeucl_constr
        JuMP.@constraint(model, vcat(gamma, A * x) in JuMP.SecondOrderCone())
    end
    if inst.epinorminf_constrs
        JuMP.@constraint(model, vcat(gamma, A * x) in MOI.NormInfinityCone(n + 1))
        JuMP.@constraint(model, vcat(sqrt(n) * gamma, A * x) in
            MOI.NormOneCone(n + 1))
    end

    return model
end


text:
using ColorTypes
using ColorTypes.FixedPointNumbers
using ColorTypes.Random
using Test

@isdefined(CustomTypes) || include("customtypes.jl")
using .CustomTypes

@testset "comparisons" begin
    Cp3 = ColorTypes.parametric3
    for C in unique(vcat(Cp3, coloralpha.(Cp3), alphacolor.(Cp3)))
        @test C{Float64}(1,0,0) == C{Float32}(1,0,0)
        @test C{Float32}(1,0,0) != C{Float32}(1,0,0.1)
    end

    for (a, b) in ((Gray(1.0), Gray(1)),
                   (GrayA(0.8, 0.6), AGray(0.8, 0.6)),
                   (RGB(1, 0.5, 0), BGR(1, 0.5, 0)),
                   (RGBA(1, 0.5, 0, 0.8), ABGR(1, 0.5, 0, 0.8)),
                   (Gray(0.8N0f8), Gray24(0.8)))
        local a, b
        @test a !== b
        @test a == b
        @test hash(a) == hash(b)
    end
    for (a, b) in ((RGB(1, 0.5, 0), RGBA(1, 0.5, 0, 0.9)),
                   (RGB(0.5, 0.5, 0.5), Gray(0.5)),
                   (HSV(0, 0, 0.5), HSV(100, 0, 0.5)), # grays
                   (Lab(70, 0, 60), LCHab(70, 60, 90)))
        local a, b
        @test a != b
        @test hash(a) != hash(b)
    end
    # It's not obvious whether we want these to compare as equal, but
    # whatever happens, you want hashing and equality-testing to yield the
    # same result
    for (a, b) in ((RGB(1, 0.5, 0), RGBA(1, 0.5, 0, 1)),
                   (HSV(100, 0.4, 0.6), AHSV(100, 0.4, 0.6, 1)))
        local a, b
        @test (a == b) == (hash(a) == hash(b))
    end
end

@testset "BigFloat comparisons" begin
    # issue #52
    @test AGray{BigFloat}(0.5,0.25) == AGray{BigFloat}(0.5,0.25)
    @test RGBA{BigFloat}(0.5, 0.25, 0.5, 0.5) == RGBA{BigFloat}(0.5, 0.25, 0.5, 0.5)
end

@testset "isapprox" begin
    @test Gray(0.8) ≈ Gray(0.8)
    @test Gray(0.8) ≈ Gray(0.8 + eps())
    @test Gray(0.8) ≈ Gray(0.8 - eps())
    @test !(Gray(0.8) ≈ Gray(0.80000002))
    @test Gray(0.8f0) ≈ Gray(0.80000002)
    @test Gray(0.8N0f8) ≈ Gray24(0.8)

    @test Gray(0.8) ≈ 0.8 + eps()
    @test 0.8 + eps() ≈ Gray(0.8)

    @test GrayA(0.8, 0.4) ≈ GrayA(0.8 + eps(), 0.4)
    @test GrayA(0.8, 0.4) ≈ GrayA(0.8, 0.4 + eps())
    @test GrayA(0.8, 0.4) ≈ GrayA(0.8 + eps(), 0.4 + eps())

    @test RGB(0.2, 0.8, 0.4) ≈ RGB(0.2, 0.8 + eps(), 0.4)
    @test RGBA(0.2, 0.8, 0.4, 0.2) ≈ RGBA(0.2, 0.8 + eps(), 0.4, 0.2 - eps())

    @test !isapprox(Gray(1), RGB(1, 1, 1))

    c_n0f8 = RGB{N0f8}(0.2, 0.69, 0.4)
    c_f32 = RGB{Float32}(0.2, 0.69, 0.4)
    @test c_n0f8 != c_f32 && c_n0f8 ≈ c_f32

    c1 = RGB(0.2, 0.8, 0.4)
    c2 = RGB(0.2, 0.7, 0.4)
    @test c1 ≈ c2 atol=0.11
    @test !isapprox(c1, c2; atol=0.09)
    @test c1 ≈ RGBX(0.2, 0.8, 0.4)
    @test !(c1 ≈ HSV(140.0f0,0.75f0,0.8f0))  # the latter comes from convert when using Colors
    @test !(HSV(0,0,0.5) ≈ HSV(100,0,0.5))
    @test HSVA(100,0.4,0.6,0.8) ≈ HSVA{Float32}(100,0.4,0.6,0.8)
end

@testset "isless" begin
    @test Gray(0.8) < Gray(0.9)
    @test Gray(0.8) <      0.9
    @test      0.8  < Gray(0.9)
    @test Gray(0.8) <= Gray(0.9)
    @test Gray(0.9) <= Gray(0.9)
    @test Gray(0.9) > Gray(0.8)
    @test Gray(0.9) >= Gray(0.8)
    @test Gray(0.9) >= Gray(0.9)

    @test Gray(0.8f0) < Gray(0.9f0)
    @test Gray(0.8f0) <= Gray(0.9f0)
    @test Gray(0.9f0) <= Gray(0.9f0)
    @test Gray(0.9f0) > Gray(0.8f0)
    @test Gray(0.9f0) >= Gray(0.8f0)
    @test Gray(0.9f0) >= Gray(0.9f0)

    @test Gray(0.8N0f8) < Gray(0.9N0f8)
    @test Gray(0.8N0f8) <= Gray(0.9N0f8)
    @test Gray(0.9N0f8) <= Gray(0.9N0f8)
    @test Gray(0.9N0f8) > Gray(0.8N0f8)
    @test Gray(0.9N0f8) >= Gray(0.8N0f8)
    @test Gray(0.9N0f8) >= Gray(0.9N0f8)

    @test Gray(0.8) < Gray(0.9f0)
    @test Gray(0.8) <= Gray(0.9f0)
    # @test Gray(0.9) <= Gray(0.9f0) is not true due to approximation
    @test Gray(0.9) > Gray(0.8f0)
    @test Gray(0.9) >= Gray(0.8f0)
    @test Gray(0.9) >= Gray(0.9f0)

    @test Gray(0.8f0) < Gray(0.9N0f8)
    @test Gray(0.8f0) <= Gray(0.9N0f8)
    @test Gray(0.9f0) <= Gray(0.9N0f8)
    @test Gray(0.9f0) > Gray(0.8N0f8)
    @test Gray(0.9f0) >= Gray(0.8N0f8)
    # @test Gray(0.9f0) >= Gray(0.9N0f8) is not true, since 0.9N0f8 = 0.902

    @test (Gray(0.3) < Gray(NaN)) == (0.3 < NaN)
    @test (Gray(NaN) < Gray(0.3)) == (NaN < 0.3)
    @test isless(Gray(0.3), Gray(NaN)) == isless(0.3, NaN)
    @test isless(Gray(NaN), Gray(0.3)) == isless(NaN, 0.3)
    @test isless(Gray(0.3), NaN) == isless(0.3, NaN)
    @test isless(Gray(NaN), 0.3) == isless(NaN, 0.3)
    @test isless(0.3, Gray(NaN)) == isless(0.3, NaN)
    @test isless(NaN, Gray(0.3)) == isless(NaN, 0.3)

    # transparent gray doesn't support comparison
    @test_throws MethodError GrayA(0.8, 0.4) < GrayA(0.9, 0.4)
    @test_throws MethodError GrayA(0.8, 0.4) <= GrayA(0.9, 0.4)
    @test_throws MethodError GrayA(0.9, 0.4) > GrayA(0.8, 0.4)
    @test_throws MethodError GrayA(0.9, 0.4) >= GrayA(0.8, 0.4)
end

@testset "rand" begin
    function all_in_range(c::C) where {C}
        all(((f, gmin, gmax),) -> gmin <= getfield(c, f) <= gmax,
            zip(ColorTypes.colorfields(C), gamutmin(C), gamutmax(C)))
    end
    eltypes = (N0f8, N2f6, N0f16, N2f14, N0f32, N2f30, Float16, Float32, Float64)
    @testset "rand($C)" for C in (
            (Gray{T} for T in eltypes)...,
            (RGB{T} for T in eltypes)...,
            AGray{Float32}, GrayA{Float64},
            ARGB{Float32}, RGBA{N0f16}, XRGB{N0f8}, RGBX{Float64},
            BGR{Float16}, ABGR{N0f32}, BGRA{N2f14},
            HSV{Float32}, HSL{Float64}, ALab{Float32}, LCHabA{Float16},
            Gray, AGray, GrayA,
            RGB, ARGB, RGBA, BGR, ABGR, BGRA, XRGB, RGBX,
            HSV, HSL, Lab, LCHab, YIQ,
            AHSV, HSLA)
        if C <: Transparent3 && !(C <: TransparentRGB)
            @test_broken rand(C)
            continue
        end
        CC = isconcretetype(C) ? C : C{Float64}
        c = rand(C)
        @test c isa CC
        @test all_in_range(c)
        a = rand(C, (3, 5))
        @test a isa Array{CC,2}
        @test eltype(a) === CC
        @test size(a) == (3, 5)
        @test all(all_in_range, a)
        ap = a'
        @test ap[1, 1] == a[1, 1]
    end
    @testset "rand($C)" for C in (Gray24, AGray32)
        c = rand(C)
        b = c.color
        @test b&0xff == (b>>8)&0xff == (b>>16)&0xff
        @test c isa C
        a = rand(C, 3, 5)
        for el in a
            b = el.color
            @test b&0xff == (b>>8)&0xff == (b>>16)&0xff
        end
        @test eltype(a) === C
        @test size(a) == (3, 5)
    end
    @test_broken rand(RGB24)
    @test_broken rand(ARGB32)
    @test_broken rand(MersenneTwister(), RGB{N0f8})
    @test_broken rand!(Gray{Float32}[0.0, 0.1])
    @test_broken all(all_in_range, rand(ARGB{Q3f12}, 10, 10))
    @test_broken all_in_range(LCHab(50, 10, 359))
    @test_broken all_in_range(YIQ(0.5, 0.59, 0.0))
    @test_broken !all_in_range(YIQ(0.5, 0.0, -0.53))
end

@testset "mapc" begin
    @test @inferred(mapc(sqrt, Gray{N0f8}(0.04))) === Gray(sqrt(N0f8(0.04)))
    @test @inferred(mapc(sqrt, AGray{N0f8}(0.04, 0.4))) === AGray(sqrt(N0f8(0.04)), sqrt(N0f8(0.4)))
    @test @inferred(mapc(sqrt, GrayA{N0f8}(0.04, 0.4))) === GrayA(sqrt(N0f8(0.04)), sqrt(N0f8(0.4)))
    @test @inferred(mapc(x->2x, RGB{N0f8}(0.04,0.2,0.3))) === RGB(map(x->2*N0f8(x), (0.04,0.2,0.3))...)
    @test @inferred(mapc(sqrt, RGBA{N0f8}(0.04,0.2,0.3,0.7))) === RGBA(map(x->sqrt(N0f8(x)), (0.04,0.2,0.3,0.7))...)
    @test @inferred(mapc(x->1.5f0x, RGBA{N0f8}(0.04,0.2,0.3,0.4))) === RGBA(map(x->1.5f0*N0f8(x), (0.04,0.2,0.3,0.4))...)
    @test @inferred(mapc(sqrt, ACMYK{N0f8}(0.04,0.2,0.3,0.7,0.1))) === ACMYK(map(x->sqrt(N0f8(x)), (0.04,0.2,0.3,0.7,0.1))...)

    @test @inferred(mapc(max, Gray{N0f8}(0.2), Gray{N0f8}(0.3))) === Gray{N0f8}(0.3)
    @test @inferred(mapc(-, AGray{Float32}(0.3), AGray{Float32}(0.2))) === AGray{Float32}(0.3f0-0.2f0,0.0)
    @test @inferred(mapc(min, RGB{N0f8}(0.2,0.8,0.7), RGB{N0f8}(0.5,0.2,0.99))) === RGB{N0f8}(0.2,0.2,0.7)
    @test @inferred(mapc(+, RGBA{N0f8}(0.2,0.8,0.7,0.3), RGBA{Float32}(0.5,0.2,0.99,0.5))) === RGBA(0.5f0+N0f8(0.2),0.2f0+N0f8(0.8),0.99f0+N0f8(0.7),0.5f0+N0f8(0.3))
    @test @inferred(mapc(+, HSVA(0.1,0.8,0.3,0.5), HSVA(0.5,0.5,0.5,0.3))) === HSVA(0.1+0.5,0.8+0.5,0.3+0.5,0.5+0.3)
    @test @inferred(mapc(+, ACMYK(0.1,0.8,0.3,0.5,0.2), ACMYK(0.5,0.5,0.5,0.5,0.5))) === ACMYK(0.1+0.5,0.8+0.5,0.3+0.5,0.5+0.5,0.2+0.5)

    @test_throws ArgumentError mapc(min, RGB{N0f8}(0.2,0.8,0.7), BGR{N0f8}(0.5,0.2,0.99))
    @test @inferred(mapc(abs, -2)) === 2
end

@testset "reducec" begin
    @test @inferred(reducec(+, 0.0, Gray(0.3))) === 0.3
    @test @inferred(reducec(+, 1.0, Gray(0.3))) === 1.3
    @test @inferred(reducec(+, 0, Gray(0.3))) === 0.3
    @test @inferred(reducec(+, 0.0, AGray(0.3, 0.8))) === 0.3 + 0.8
    @test @inferred(reducec(+, 0.0, RGB(0.3, 0.8, 0.5))) === (0.3 + 0.8) + 0.5
    @test @inferred(reducec(+, 0.0, RGBA(0.3, 0.8, 0.5, 0.7))) === ((0.3 + 0.8) + 0.5) + 0.7
    @test @inferred(reducec(+, 0.0, ACMYK(0.3, 0.8, 0.5, 0.7, 0.2))) === (((0.3 + 0.8) + 0.5) + 0.7) + 0.2
    @test @inferred(reducec(&, true, Gray(true)))
    @test !(@inferred(reducec(&, false, Gray(true))))
    @test !(@inferred(reducec(&, true, Gray(false))))

    @test @inferred(reducec(+, 0.0, 0.3)) === 0.3
    @test @inferred(reducec(+, 0, 0.3)) === 0.3
    @test @inferred(reducec(&, true, true))
    @test !(@inferred(reducec(&, false, true)))
    @test !(@inferred(reducec(&, true, false)))
end

@testset "mapreducec" begin
    @test @inferred(mapreducec(x->x^2, +, 0.0, Gray(0.3))) === 0.3^2
    @test @inferred(mapreducec(x->x^2, +, 1.0, Gray(0.3))) === 1 + 0.3^2
    @test @inferred(mapreducec(x->x^2, +, 0, Gray(0.3))) === 0.3^2
    @test @inferred(mapreducec(x->x^2, +, 0.0, AGray(0.3, 0.8))) === 0.3^2 + 0.8^2
    @test @inferred(mapreducec(x->x^2, +, 0.0, RGB(0.3, 0.8, 0.5))) === (0.3^2 + 0.8^2) + 0.5^2
    @test @inferred(mapreducec(x->x^2, +, 0.0, RGBA(0.3, 0.8, 0.5, 0.7))) === ((0.3^2 + 0.8^2) + 0.5^2) + 0.7^2
    @test @inferred(mapreducec(x->x^2, +, 0.0, ACMYK(0.3, 0.8, 0.5, 0.7, 0.2))) === (((0.3^2 + 0.8^2) + 0.5^2) + 0.7^2) + 0.2^2

    @test !(@inferred(mapreducec(x->!x, &, true, Gray(true))))
    @test !(@inferred(mapreducec(x->!x, &, false, Gray(true))))
    @test @inferred(mapreducec(x->!x, &, true, Gray(false)))
    @test !@inferred(mapreducec(x->!x, &, false, Gray(false)))

    @test @inferred(mapreducec(x->x^2, +, 0.0, 0.3)) === 0.3^2
    @test @inferred(mapreducec(x->x^2, +, 1.0, 0.3)) === 1 + 0.3^2
    @test @inferred(mapreducec(x->x^2, +, 0, 0.3)) === 0.3^2
    @test !(@inferred(mapreducec(x->!x, &, true, true)))
    @test !(@inferred(mapreducec(x->!x, &, false, true)))
    @test @inferred(mapreducec(x->!x, &, true, false))
    @test !@inferred(mapreducec(x->!x, &, false, false))
end


text:
@testset "Bernoulli{<:LogisticLink}" begin
    test_auglik(BernoulliLikelihood(LogisticLink()))
end


text:
using ImperoHooks
using Impero
using ClimateMachine

using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Elements
import ClimateMachine.Mesh.Elements: baryweights
using GaussQuadrature
using Base.Threads
using Makie, GLMakie, AbstractPlotting
using ImageTransformations, Colors
using AbstractPlotting.MakieLayout, LinearAlgebra
using Impero: Circle
using Random

include(pwd() * "/test/gradient_test.jl")
include(pwd() * "/test/divergence_test.jl")
include(pwd() * "/test/gradient_divergence_test.jl")

include(pwd() * "/test/impero_algebra_test.jl")

include(pwd() * "/test/gridhelper_test.jl")
include(pwd() * "/test/permutation_test.jl")

include(pwd() * "/test/field_test.jl")

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
"""
inflated: Inflated component.
rest: The rest of the distribution.
prob: Probability for the inflated component.

Note that `inflated` should not be in the support of `rest`.
"""
struct HurdleModel{A<:Real, B<:UnivariateDistribution, C<:Real} <: ContinuousUnivariateDistribution
  inflated::A
  rest::B
  prob::C
  probrest::C
end
HurdleModel(inflated, rest, prob) = HurdleModel(inflated, rest, prob, 1 - prob)

function Distributions.pdf(d::HurdleModel, x::Real)
  return x == d.inflated ? d.prob : d.probrest * pdf(d.rest, x)
end

function Distributions.logpdf(d::HurdleModel, x::Real)
  if x == d.inflated
    return log(d.prob)
  else
    return log(d.probrest) + logpdf(d.rest, x)
  end
end

function Distributions.rand(rng::Random.AbstractRNG, d::HurdleModel)
  d.prob > rand(rng) ? d.inflated : rand(rng, d.rest)
end


text:
import Base: *, A_mul_B!
import Base: A_mul_Bt, A_mul_Bt!, A_mul_Bc, A_mul_Bc!
import Base: At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!
import Base: At_mul_Bt, At_mul_Bt!, Ac_mul_Bc, Ac_mul_Bc!
import Base.LinAlg.BLAS: gemm!

A_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','N',one(T),A,B,zero(T),C)
(*){T}(A::KnetMatrix{T},B::KnetMatrix{T})=A_mul_B!(similar(A,(size(A,1),size(B,2))),A,B)

A_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','T',one(T),A,B,zero(T),C)
A_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(similar(A,(size(A,1),size(B,1))),A,B)
A_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(C,A,B)
A_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt(A,B)

At_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','N',one(T),A,B,zero(T),C)
At_mul_B{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(similar(A,(size(A,2),size(B,2))),A,B)
Ac_mul_B!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(C,A,B)
Ac_mul_B{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B(A,B)

At_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','T',one(T),A,B,zero(T),C)
At_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(similar(A,(size(A,2),size(B,2))),A,B)
Ac_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(C,A,B)
Ac_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt(A,B)

function gemm!{T}(transA::Char, transB::Char, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T})
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    cublascheck(x) = (x==0 || error("CUBLAS error $x"))
    size2(x,i)=(if ndims(x)<=2; size(x,i); elseif i==1; div(length(x),size(x,ndims(x))); elseif i==2; size(x,ndims(x)); else 1; end)
    if transA == 'N'
        m=size2(A,1); k=size2(A,2)
    else
        m=size2(A,2); k=size2(A,1)
    end
    if transB == 'N'
        n=size2(B,2); k==size2(B,1) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    else
        n=size2(B,1); k==size2(B,2) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    (m == size2(C,1) && n == size(C,2)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T(alpha); beta = T(beta)
    lda = size2(A,1); ldb = size2(B,1); ldc = size2(C,1)
    if T<:Float32
        cublascheck(ccall((:cublasSgemm_v2, "libcublas"), UInt32, (Ptr{Void}, UInt32, UInt32, Cint, Cint, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Ptr{Cfloat}, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Cint), cublashandle, transa, transb, m, n, k, [alpha], A, lda, B, ldb, [beta], C, ldc))
    elseif T<:Float64
        cublascheck(ccall((:cublasDgemm_v2, "libcublas"), UInt32, (Ptr{Void}, UInt32, UInt32, Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), cublashandle, transa, transb, m, n, k, [alpha], A, lda, B, ldb, [beta], C, ldc))
    else
        error("CUBLAS does not support $T")
    end
    return C
end


text:
# Copyright 2015 Conrad Scott

function readVocab(filename :: AbstractString)
  df = readtable( filename,
                  header = false, separator = ' ',
                  skipblanks = false, ignorepadding = false,
                  names = [:Terms], eltypes = [UTF8String] )

  (lines, columns) = size(df)

  if columns != 1
    error("Vocabulary file should have a single word per line; \"$filename\" has $columns columns")
  end

  return Dictionary(df[:Terms])
end


text:
for func in [:clblasSnrm2, :clblasDnrm2, :clblasScnrm2, :clblasDznrm2]
    @eval @blasfun $func(N::Csize_t, NRM2::cl.CL_mem, offNRM2::Csize_t,
                              X::cl.CL_mem, offx::Csize_t, incx::Cint,
                              scratch_buff::cl.CL_mem,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})
    @eval @blasfun2 $func(N::Csize_t, NRM2::cl.CL_mem, offNRM2::Csize_t,
                               X::cl.CL_mem, offx::Csize_t, incx::Cint,
                               scratch_buff::cl.CL_mem)
end


text:
################################################################################
#   Utilities
################################################################################

# Masks
const DEPTH_BUFFER_BIT   = 0x00000100
const STENCIL_BUFFER_BIT = 0x00000400
const COLOR_BUFFER_BIT   = 0x00004000

function Clear(mask::Integer)
	ccall( (:glClear, lib), Void, (GLuint,), mask)
end

# Capabilities
const BLEND                    = 0x0BE2
const CULL_FACE                = 0x0B44
const DEPTH_TEST               = 0x0B71
const DITHER                   = 0x0BD0
const POLYGON_OFFSET_FILL      = 0x8037
const SAMPLE_ALPHA_TO_COVERAGE = 0x809E
const SAMPLE_COVERAGE          = 0x80A0
const SCISSOR_TEST             = 0x0C11
const STENCIL_TEST             = 0x0B90

function Disable(cap::Integer)
	ccall( (:glDisable, lib), Void, (GLenum,), cap)
end

function Enable(cap::Integer)
	ccall( (:glEnable, lib), Void, (GLenum,), cap)
end

# Error codes
const NO_ERROR          = 0
const INVALID_ENUM      = 0x0500
const INVALID_VALUE     = 0x0501
const INVALID_OPERATION = 0x0502
const OUT_OF_MEMORY     = 0x0505

function GetError()
	ret = ccall( (:glGetError, lib), GLenum, ())
	if ret != NO_ERROR
		error(ret)
	end
end

function GetString(name::Integer)
	ret = ccall( (:glGetString, lib), Ptr{GLubyte}, (GLenum,), name)
	bytestring(ret)
end

function GetStringi(name::Integer, index::Integer)
	ret = ccall( (:glGetStringi, lib), Ptr{GLubyte}, (GLenum, GLuint), name, index)
	bytestring(ret)
end

GetString(name::Integer, index::Integer) = GetStringi(name, index)

# String names
const VENDOR                   = 0x1F00
const RENDERER                 = 0x1F01
const VERSION                  = 0x1F02
const SHADING_LANGUAGE_VERSION = 0x8B8C
const EXTENSIONS               = 0x1F03

function Viewport(width::Integer, height::Integer)
	ccall( (:glViewport, lib), Void,
		(GLint, GLint, GLsizei, GLsizei),
		0, 0, width, height)
end


text:
# Helpers to gather and store some information useful to multiple diagnostics
# groups.
#

Base.@kwdef mutable struct AtmosCollectedDiagnostics
    onetime_done::Bool = false
    zvals::Union{Nothing, Array} = nothing
    repdvsr::Union{Nothing, Array} = nothing
end
const AtmosCollected = AtmosCollectedDiagnostics()

function atmos_collect_onetime(mpicomm, dg, Q)
    if !AtmosCollected.onetime_done
        FT = eltype(Q)
        grid = dg.grid
        topology = grid.topology
        N = polynomialorder(grid)
        Nq = N + 1
        Nqk = dimensionality(grid) == 2 ? 1 : Nq
        nrealelem = length(topology.realelems)
        nvertelem = topology.stacksize
        nhorzelem = div(nrealelem, nvertelem)

        localvgeo = array_device(Q) isa CPU ? grid.vgeo : Array(grid.vgeo)

        AtmosCollected.zvals = zeros(FT, Nqk * nvertelem)
        AtmosCollected.repdvsr = zeros(FT, Nqk * nvertelem)

        @visitQ nhorzelem nvertelem Nqk Nq begin
            evk = Nqk * (ev - 1) + k
            z = localvgeo[ijk, grid.x3id, e]
            MH = localvgeo[ijk, grid.MHid, e]
            AtmosCollected.zvals[evk] = z
            AtmosCollected.repdvsr[evk] += MH
        end

        # compute the full number of points on a slab
        MPI.Allreduce!(AtmosCollected.repdvsr, +, mpicomm)

        AtmosCollected.onetime_done = true
    end

    return nothing
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# SPDX-License-Identifier: X11
# 2020-07-30
# Payment (100pt)

const BILL = 1000
N = parse(Int,readline())

N % BILL == 0 ? println(0) : println(BILL - N % BILL)

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module TimeseriesSurrogates

using Distributions
using StatsBase
using InplaceOps
using AbstractFFTs
using DSP
using Interpolations
using Wavelets
using Requires

# Example systems
include("testsystems.jl")

# Periodogram interpolation
include("interpolation.jl")

include("uncertaindatasets.jl")

# The different surrogate routines
include("randomshuffle.jl")
include("randomphases.jl")
include("randomamplitudes.jl")
include("aaft.jl")
include("iaaft.jl")
include("wiaaft.jl")

export NLNS, NSAR2, AR1, randomwalk, SNLST,
        randomshuffle, randomamplitudes, randomphases, aaft, iaaft, wiaaft

# Visualization routine for time series + surrogate + periodogram/acf/histogram
using Requires
function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        # Define and export plot routines for all combinations of example processes and surrogate
        # types
        processes = (:AR1, :NSAR2, :NLNS, :randomwalk, :SNLST)
        surrogate_methods = (:randomshuffle, :randomphases, :randomamplitudes, :aaft, :iaaft)
        include("surrogate_plot.jl")
        include("plots_and_anim.jl")
    end
end

end # module


text:
module OMEinsum
export einsum, dynamic_einsum
export einsumopt

using TupleTools, LinearAlgebra
using BatchedRoutines
import LinearAlgebra: BlasFloat
const CuBlasFloat = Union{BlasFloat, Float16, ComplexF16}

include("Core.jl")
include("loop_einsum.jl")
include("utils.jl")

include("unaryrules.jl")
include("binaryrules.jl")

using Requires
function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cueinsum.jl")
end

include("interfaces.jl")
include("einsequence.jl")
include("autodiff.jl")

include("contractionorder/contractionorder.jl")

include("deprecation.jl")
end # module


text:
module ADTest

using ForwardDiff, FiniteDifferences, TensorNetworkEvolve.TensorAD, OMEinsum

function unpackargs(args, x, mask)
    start = 0
    map(args, mask) do arg, m
        if m
            @assert arg isa AbstractArray{<:TensorAD.ADTypes}
            if eltype(arg) <: Real
                res = reshape(x[start+1:start+length(arg)], size(arg)...)
                start += length(arg)
                res
            else
                rel = x[start+1:start+length(arg)]
                img = x[start+length(arg)+1:start+2*length(arg)]
                res = reshape(rel + im * img, size(arg)...)
                start += 2*length(arg)
                res
            end
        else
            arg
        end
    end
end

function packargs(args, mask)
    vecs = []
    for (arg, m) in zip(args, mask)
        if m
            @assert arg isa AbstractArray{<:TensorAD.ADTypes}
            if eltype(arg) <: Real
                push!(vecs, vec(arg))
            else
                push!(vecs, vcat(vec(real(arg)), vec(imag(arg))))
            end
        end
    end
    return vcat(vecs...)
end

function build_testfunc(f, args...; mask=map(x->x isa AbstractArray{<:TensorAD.ADTypes}, args), realpart=true, kwargs...)
    function (x)
        _args = unpackargs(args, x, mask)
        if realpart
            OMEinsum.asarray(real(f(_args...; kwargs...)), x)
        else
            OMEinsum.asarray(imag(f(_args...; kwargs...)), x)
        end
    end, packargs(args, mask)
end

function match_jacobian(f, args...; realpart=true, atol=1e-5, kwargs...)
    tf, x = build_testfunc(f, args...; realpart, kwargs...)
    j1 = ForwardDiff.jacobian(tf, x)
    j2 = TensorAD.jacobian(tf, DiffTensor(x))
    return isapprox(j1, j2.data; atol)
end

function match_random(f, args...; realpart=true, atol=1e-5, kwargs...)
    empty!(TensorAD.GLOBAL_TAPE.instructs)
    tf, x = build_testfunc(f, args...; realpart, kwargs...)
    j1 = ForwardDiff.jacobian(tf, x)
    X = DiffTensor(x)
    y = tf(X)
    gy = ndims(y) == 0 ? fill(randn(eltype(y))) : randn(eltype(y), size(y)...)
    g1 = vec(vec(gy)' * j1)

    grad_storage = Dict{UInt,Any}()
    TensorAD.accumulate_gradient!(grad_storage, TensorAD.getid(y), DiffTensor(gy, false))
    TensorAD.back!(TensorAD.GLOBAL_TAPE, grad_storage)
    g2 = vec(TensorAD.getgrad(grad_storage, X))
    @debug "g1 = $(g1), g2 = $(g2.data)"
    return isapprox(g1, g2.data; atol) && typeof(X) == typeof(g2)
end

function match_hessian(f, args...; realpart=true, atol=1e-5, kwargs...)
    tf, x = build_testfunc(f, args...; realpart, kwargs...)
    h2 = TensorAD.hessian(tf, DiffTensor(x))
    h1 = FiniteDifferences.jacobian(central_fdm(5,1), x->ForwardDiff.gradient(x->tf(x)[], x), x)[1]
    return isapprox(h1, h2.data; atol)
end

end


text:
# ---
# title: 793. Preimage Size of Factorial Zeroes Function
# id: problem793
# author: Indigo
# date: 2021-01-29
# difficulty: Hard
# categories: Binary Search
# link: <https://leetcode.com/problems/preimage-size-of-factorial-zeroes-function/description/>
# hidden: true
# ---
#
# Let `f(x)` be the number of zeroes at the end of `x!`. (Recall that `x! = 1 *
# 2 * 3 * ... * x`, and by convention, `0! = 1`.)
#
# For example, `f(3) = 0` because 3! = 6 has no zeroes at the end, while `f(11)
# = 2` because 11! = 39916800 has 2 zeroes at the end. Given `K`, find how many
# non-negative integers `x` have the property that `f(x) = K`.
#
#
#
#     **Example 1:**
#     Input: K = 0
#     Output: 5
#     Explanation: 0!, 1!, 2!, 3!, and 4! end with K = 0 zeroes.
#
#     **Example 2:**
#     Input: K = 5
#     Output: 0
#     Explanation: There is no x such that x! ends in K = 5 zeroes.
#
#
# **Note:**
#
#   * `K` will be an integer in the range `[0, 10^9]`.
#
#
## @lc code=start
using LeetCode

function preimage_size_fzf(K::Int)
    zeta(x) = x > 0 ? x ÷ 5 + zeta(x ÷ 5) : 0
    lo, hi = K, 5 * K + 1
    while lo < hi
        mid = (lo + hi) ÷ 2
        z_mid = zeta(mid)
        if z_mid == K
            return 5
        elseif z_mid < K
            lo = mid + 1
        else
            hi = mid
        end
    end
    return 0
end
## @lc code=end


text:
struct mv_group_gaussian <: distibution_sample
    μ::Vector{AbstractArray{Float64,1}}
    Σ::AbstractArray{Float64,2}
    invΣ::AbstractArray{Float64,2}
    logdetΣ::Float64
end

function log_likelihood!(r::AbstractArray,x::AbstractArray, distibution_sample::mv_group_gaussian ,group::Int64 = 1)
    z = x .- distibution_sample.μ[group]
    dcolwise_dot!(r,z, distibution_sample.invΣ * z)
    r .= -((length(distibution_sample.Σ) * Float64(log(2π)) + logdet(distibution_sample.Σ))/2) .-r
end


text:
using FFmpegPipe
using Test
using Plots
gr()

f1 = tempname()*".mp4"

s = openvideo(f1, "w", r=24)
for a in range(0, stop=pi, length=24)
    x = a .+ range(0, stop=pi, length=100)
    plt = plot(x, sin.(x))
    writeframe(s, plt)
end
close(s)

@test isfile(f1)

f2 = tempname()*".mp4"

s1 = openvideo(f1, "r")
s2 = openvideo(f2, "w", r=24)
while !eof(s1)
    img = readframe(s1)
    writeframe(s2, img[end:-1:1,:])
end
close(s2)
close(s1)

@test isfile(f2)

rm(f2)

w = 103
h = 97
s1 = openvideo(f1, "r", ss_in = 0.5, s_out = "$(w)x$h")
img = readframe(s1)
@test size(img) == (h, w)
close(s1)

rm(f1)

max_stars_repo_name: perrutquist/FfmpegPipe

max_issues_repo_name: perrutquist/FfmpegPipe

max_forks_repo_name: perrutquist/FfmpegPipe


text:
# Example for basic 3D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, opesciSLIM.TimeModeling

## Set up model structure
n = (120,100,90)	# (x,y,z) or (x,z)
d = (10.,10.,10.)
o = (0.,0.,0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.4f0
v0 = ones(Float32,n) + 0.4f0
v[:,:, Int(round(end/2)):end] = 4.0f0
rho = ones(Float32,n)

# Slowness squared [s^2/km^2]
m = (1./v).^2
m0 = (1./v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 4
model = Model(n,d,o,m)	# to include density call Model(n,d,o,m,rho)
model0 = Model(n,d,o,m0)

## Set up 3D receiver geometry by defining one receiver vector in each x and y direction
nxrec = 120
nyrec = 100
xrec = linspace(50.,1150.,nxrec)
yrec = linspace(100., 900., nyrec)
zrec = 50

# Construct 3D grid from basis vectors
(xrec, yrec, zrec) = setup_3D_grid(xrec, yrec, zrec)

# receiver sampling and recording time
timeR = 1000.	# receiver recording time [ms]
dtR = 4.	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell([250. 500. 750. 1000.])
ysrc = convertToCell([200. 400. 600. 800.])
zsrc = convertToCell([50. 60. 70. 80.])

# source sampling and number of time steps
timeS = 1000.	# source length in [ms]
dtS = 2.	# source sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01
wavelet = source(timeS,dtS,f0)

# Info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)	# no. of computational time steps
info = Info(prod(n),nsrc,ntComp)

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

# Setup operators
Pr = joProjection(info,recGeometry)
F = joModeling(info,model)
F0 = joModeling(info,model0)
Ps = joProjection(info,srcGeometry)
q = joData(srcGeometry,wavelet)

# Combined operator Pr*F*Ps'
Ffull = joModeling(info,model,srcGeometry,recGeometry)

# Nonlinear modeling
d = Pr*F*Ps'*q	# equivalent to d = Ffull*q
qad = Ps*F'*Pr'*d

# Linearized modeling
J = joJacobian(Pr*F0*Ps',q)	# equivalent to J = joSetupJacobian(Ffull,q)
dD = J*dm
rtm = J'*dD

# Plot results
figure(); imshow(d.data[1][:,1:500], vmin=-1e-1, vmax=1e-1)
rtm = reshape(rtm,model.n)
figure(); imshow(rtm[:,50,:]', ColorMap("gray"))

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
## How to reimplement event handler. In Qt we would subclass. Do the same in PySide

using PySide

## we have a hack to create a new class and add methods to an instance.
## (You would need to write python code to create a new class with methods for all instances.)

## we call enter when we enter a widget

qnew_class("Example", "QtGui.QWidget") ## QtGui -- not just Qt.
o = qnew_class_instance("Example")
qset_method(o, :enterEvent) do e
  println("Welcome stranger.")
end

qset_method(o, :leaveEvent) do e
  println("Bye bye. Y'all come back now.")
end

raise(o)


text:
import ILKBackend
import ur5
import NaiveDataset
import LinearAlgebra

mc = ur5.ModelConstants

function test_fk_dataset( ds::NaiveDataset.BinDataset )

    err_pos = 0;
    err_ori = 0;
    count   = 0;
    qd = ILKBackend.matrix(6,1)
    while ! ds.eof
        q = NaiveDataset.readVector(ds, UInt16(6));
        given_1 = NaiveDataset.readPose(ds);
        wrist_3__base,v__wrist_3__base,J_wrist_3_base = ur5.solver1(mc, q, qd);

        Rdiff = ILKBackend.orientationDistance(
                    ILKBackend.rotationView(wrist_3__base),
                    ILKBackend.rotationView(given_1))
        err_ori += Rdiff.angle;
        err_pos += LinearAlgebra.norm(ILKBackend.positionView(wrist_3__base-given_1));
        count = count + 1;
    end

    err_pos /= count
    err_ori /= count

    println("Average position error   : " * string(err_pos))
    println("Average orientation error: " * string(err_ori))
end

if abspath(PROGRAM_FILE) == @__FILE__
    if length(Base.ARGS) < 1
        print("Please provide a dataset file")
        exit(-1)
    end
    filename = Base.ARGS[1]
    ds = NaiveDataset.openBinDataset( filename )

    test_fk_dataset( ds )
end

max_stars_repo_name: ESROCOS/ilk-compiler

max_issues_repo_name: ESROCOS/ilk-compiler

max_forks_repo_name: ESROCOS/ilk-compiler


text:
struct UGrid{T <: Real, N <: Integer}
                   nxu :: N                                     # Number of frid cells in x-direction in UGrid
                   nyu :: N                                     # Number of grid cells in y-direction in UGrid
                 mask :: Array{Bool,2}                         # Mask specifying model domain wrt U grid
                    n :: N                                     # Total number of cells in model domain
                 crop :: Diagonal{T,Array{T,1}}                # Crop matrix: diagonal matrix with mask entries on diag
                 samp :: SparseMatrixCSC{T,N}                  # Sampling matrix: take full domain to model domain
               spread :: SparseMatrixCSC{T,N}                  # Spread matrix: take model domain to full domain
                 cent :: KronType{T,N}                         # Map from U grid to H grid
                   ∂x :: KronType{T,N}                         # Matrix representation of differentiation wrt x
                   ∂y :: KronType{T,N}                         # Matrix representation of differentiation wrt y
               levels :: N                                     # Number of levels in the preconditioner
                  dwt :: KronType{T,N}                         # Wavelet matrix product on u grid
                    s :: Array{T,2}                            # Ice surface elevation
                    h :: Array{T,2}                            # Ice thickness
    grounded_fraction :: Array{T,2}                            # Grid cell grounded fraction
                 βeff :: Array{T,2}                            # Effective β value on u grid(eqn 12 in Arthern 2015 JGeophysRes)
             dnegβeff :: Base.RefValue{Diagonal{T,Array{T,1}}} # Rheological operator (diagonal of βeff)
                    u :: Array{T,2}                            # Ice velocities in x direction
end

"""
    UGrid(;
            nxu,
            nyu,
            mask = trues(nxu,nyu),
            levels,
            dx,
            dy)

Construct a WAVI.jl UGrid with size (nxu,nyu)
UGrid stores fields that are defined on the problem's U grid.
(Co-ordinates of UGrid stored in a Grid under xxu and yyu fields)

Keyword arguments
=================
    - 'nxu': (required) Number of grid cells in x-direction in UGrid (should be same as grid.nx + 1)
            Note that we store the grid size here, even though it can be easily inferred from grid, to increase transparency in velocity solve.
    - 'nyu': (required) Number of grid cells in y-direction in UGrid (should be same as grid.ny)
    - 'mask': Mask specifying the model domain with respect to U grid
    - levels: (required) Number of levels in the preconditioner
    - dx: (required) Grid spacing in the x direction
    - dy: (required) Grid spacing in the y direction
"""
function UGrid(;
                nxu,
                nyu,
                mask = trues(nxu,nyu),
                levels,
                dx,
                dy)

    #check the sizes of inputs
    (size(mask) == (nxu,nyu)) || throw(DimensionMismatch("Sizes of inputs to UGrid must all be equal to nxu x nyu (i.e. $nxu x $nyu)"))

    #construct operators
    n = count(mask)
    crop = Diagonal(float(mask[:]))
    samp = crop[mask[:],:]
    spread = sparse(samp')
    cent =  spI(nyu) ⊗ c(nxu-1)
    ∂x = spI(nyu) ⊗ ∂1d(nxu-1,dx)
    ∂y =  ∂1d(nyu-1,dy) ⊗ χ(nxu-2)
    dωt = wavelet_matrix(nyu,levels,"forward" ) ⊗ wavelet_matrix(nxu,levels,"forward")

    #fields stored on UGrid
    s = zeros(nxu,nyu)
    h = zeros(nxu,nyu)
    grounded_fraction = ones(nxu,nyu)
    βeff = zeros(nxu,nyu)
    dnegβeff = Ref(crop*Diagonal(-βeff[:])*crop)
    u = zeros(nxu,nyu)

    #size assertions
    @assert size(mask)==(nxu,nyu)
    n == count(mask)
    @assert crop == Diagonal(float(mask[:]))
    @assert samp == crop[mask[:],:]
    @assert spread == sparse(samp')
    @assert size(s)==(nxu,nyu)
    @assert size(h)==(nxu,nyu)
    @assert size(grounded_fraction)==(nxu,nyu)
    @assert size(βeff)==(nxu,nyu)
    @assert size(u)==(nxu,nyu)

    #make sure boolean type rather than bitarray
    mask = convert(Array{Bool,2}, mask)

    return UGrid(
                nxu,
                nyu,
                mask,
                n,
                crop,
                samp,
                spread,
                cent,
                ∂x,
                ∂y,
                levels,
                dωt,
                s,
                h,
                grounded_fraction,
                βeff,
                dnegβeff,
                u)
end


text:
module Experimental
    using VegaLite, ColorSchemes, Setfield

    ## set branch color scheme after graph is created.  This does not currently function because the layer numbers are may
    # not correlate if a device does not exist (e.g. dclines)
    function set_branch_color_range(plot::VegaLite.VLSpec, cs::ColorSchemes.ColorScheme)
        return set_branch_color_range(plot, colorscheme2array(cs))
    end

    function set_branch_color_range(plot::VegaLite.VLSpec, cr::AbstractVector{String})
        return plot.layer[1]["encoding"]["color"]["scale"]["range"] = cr
    end

    function set_connector_color_range(plot::VegaLite.VLSpec, cs::ColorSchemes.ColorScheme)
        return set_connector_color_range(plot, colorscheme2array(cs))
    end

    function set_connector_color_range(plot::VegaLite.VLSpec, cr::AbstractVector{String})
        return plot.layer[2]["encoding"]["color"]["scale"]["range"] = cr
    end

    # function set_dcline_color_range(plot::VegaLite.VLSpec, cs::ColorSchemes.ColorScheme)
    #   return set_dcline_color_range(plot, colorscheme2array(cs))
    # end

    # function set_dcline_color_range(plot::VegaLite.VLSpec, cr::AbstractVector{String})
    #   return plot.layer[3]["encoding"]["color"]["scale"]["range"] = cr
    # end

    function set_bus_color_range(plot::VegaLite.VLSpec, cs::ColorSchemes.ColorScheme)
        return set_bus_color_range(plot, colorscheme2array(cs))
    end

    function set_bus_color_range(plot::VegaLite.VLSpec, cr::AbstractVector{String})
        return plot.layer[3]["encoding"]["color"]["scale"]["range"] = cr
    end

    function set_gen_color_range(plot::VegaLite.VLSpec, cs::ColorSchemes.ColorScheme)
        return set_gen_color_range(plot, colorscheme2array(cs))
    end

    function set_gen_color_range(plot::VegaLite.VLSpec, cr::AbstractVector{String})
        return plot.layer[4]["encoding"]["color"]["scale"]["range"] = cr
    end

    "Change the cartesian view to a geographic projection. Default is albersUsa"
    function cartesian2geo!(plot::VegaLite.VLSpec; projection_type=:albersUsa)

        # set the plot geo map projection
        @set! plot.projection=Dict{String,Any}("type"=>projection_type)

        # create lat/lon channels from x/y channels
        for i in 1:length(plot.layer)
            plot.layer[i]["encoding"]["longitude"] = plot.layer[i]["encoding"]["x"]
            delete!(plot.layer[i]["encoding"],"x")
            plot.layer[i]["encoding"]["latitude"] = plot.layer[i]["encoding"]["y"]
            delete!(plot.layer[i]["encoding"],"y")
            if haskey(plot.layer[i]["encoding"],"x2")
                plot.layer[i]["encoding"]["longitude2"] = plot.layer[i]["encoding"]["x2"]
                delete!(plot.layer[i]["encoding"],"x2")
            end
            if haskey(plot.layer[i]["encoding"],"y2")
                plot.layer[i]["encoding"]["latitude2"] = plot.layer[i]["encoding"]["y2"]
                delete!(plot.layer[i]["encoding"],"y2")
            end
        end
        return plot
    end

    "Make zoomable my modifying layer 1"
    function add_zoom!(plot::VegaLite.VLSpec)
        return plot.layer[1]["selection"]=Dict{String,Any}(
            "grid"=>Dict{String,Any}(
                "type"=>:interval,
                "resolve"=>:global,
                "bind"=>:scales,
                "translate"=>"[mousedown[!event.shiftKey], window:mouseup] > window:mousemove!",
                "zoom"=>"wheel![!event.shiftKey]"
            )
        )
    end

end


text:
# Circuitscape file types:
# - resistance/conductance map (HabitatRasterMap)
# - focal nodes (point/region)
# - short circuit region file
# - include/exclude pairs
# - current sources
# - ground points/resistances
# - raster mask file
# - source strengths

# =================================
# Habitat representations
# =================================
# TODO: support reclassification
type HabitatRasterMap
    cond::Array{Float64,2}
    hdr::FileAAGridHeader
    nodata::Real

    function HabitatRasterMap(filename::String, resistance::Bool)
        grid = FileAAGrid(filename)
        hdr = grid.hdr
        data = grid.data

        if resistance
            # convert data into conductance format
            (findfirst(data, 0) > 0) && error("Zero resistance values are not currently supported for habitat maps.  Use a short-circuit region file instead.")
            cond = 1 ./ data
        else
            cond = data
        end

        nodata = (nothing != hdr.nodata) ? -9999 : hdr.nodata
        cond[find(x->(x==nodata), data)] = 0

        new(cond, hdr, nodata)
    end
end

function show(io::IO, hab::HabitatRasterMap)
    print(io, "Habitat Raster Map as ")
    show(hab.hdr)
    println(io, "nodata: $(hab.nodata)")
    show(hab.cond)
end

# ==================================
# Habitat
# ==================================
type Habitat
    # the node map
    node_map::Array{Int,2}
    num_nodes::Int

    # conductances between nodes. store as sparse matrix?
    node1::Array{Int,1}
    node2::Array{Int,1}
    g::Array{Float64,1}

    ccs::Array{Array{Int,1},1}     # the connected component list
    component_map::Array{Int,2}    # can be Int16 probably
    num_components::Int

    # remember parameters
    connect_using_avg_resistances::Bool
    connect_four_neighbors_only::Bool
end

function Habitat(rastmap::HabitatRasterMap, connect_using_avg_resistances::Bool=false, connect_four_neighbors_only::Bool=false)
    cond = rastmap.cond
    sz = size(cond)

    # create a grid of node numbers corresponding to raster pixels with non-zero conductances
    # Note: node ids are allocated in column major order. Should not matter for the calculations though.
    node_map = zeros(Int, sz...)
    nzidxs = findn(reshape(cond, prod(sz)))[1]

    num_nodes = length(nzidxs)
    node_map[nzidxs] = 1:num_nodes

    # create the node conductance graph
    node1, node2, cond = node_conductances(cond, node_map, connect_using_avg_resistances, connect_four_neighbors_only)
    gr = simple_adjlist(num_nodes, is_directed=false)
    for idx in 1:length(node1)
        add_edge!(gr, node1[idx], node2[idx])
    end

    # find connected components and create component map
    # we can probably store only the compact form (ccs)
    ccs = connected_components(gr)
    component_map = zeros(Int, sz...)
    for compid in 1:length(ccs)
        component_map[findin(node_map, ccs[compid])] = compid
    end

    Habitat(node_map, num_nodes, node1, node2, cond, ccs, component_map, length(ccs), connect_using_avg_resistances, connect_four_neighbors_only)
end

function show(io::IO, hab::Habitat)
    println(io, "Habitat with $(hab.num_nodes) nodes, $(hab.num_components) components")
end

# Calculates conductances between adjacent nodes given a raster conductance map.
# Returns an adjacency matrix with values representing node-to-node conductance values.
function node_conductances(gmap::Array{Float64,2}, node_map::Array{Int,2}, connect_using_avg_resistances::Bool, connect_four_neighbors_only::Bool)
    s_horiz, t_horiz = _neighbors_horz(gmap)
    s_vert, t_vert = _neighbors_vert(gmap)

    s = [s_horiz, s_vert]
    t = [t_horiz, t_vert]

    g1 = gmap[s]
    g2 = gmap[t]

    if connect_using_avg_resistances
        cond = 1 ./ (((1 ./ g1) .+ (1 ./ g2)) ./ 2)
    else
        cond = (g1 .+ g2) ./ 2
    end

    if !connect_four_neighbors_only
        s_dr, t_dr = _neighbors_diag1(gmap)
        s_dl, t_dl = _neighbors_diag2(gmap)

        s_d = [s_dr, s_dl]
        t_d = [t_dr, t_dl]

        g1 = gmap[s_d]
        g2 = gmap[t_d]

        if connect_using_avg_resistances
            cond_d = 1 ./ (sqrt(2) .* ((1 ./ g1) + (1 ./ g2) ./ 2))
        else
            cond_d = (g1 .+ g2) ./ (2 * sqrt(2))
        end

        s = [s, s_d]
        t = [t, t_d]
        cond = [cond, cond_d]
    end

    node1 = node_map[s]
    node2 = node_map[t]
    (node1, node2, cond)
end

function _neighbors_horz(gmap::Array{Float64,2})
    m,n = size(gmap)

    gmap_l = gmap[:, 1:(n-1)]
    gmap_r = gmap[:, 2:n]

    gmap_lr = convert(Array{Bool,2}, gmap_l .* gmap_r)
    gmap_lr = hcat(gmap_lr, falses(m))

    s_horiz = find(gmap_lr)
    t_horiz = s_horiz .+ m

    return (s_horiz, t_horiz)
end

function _neighbors_vert(gmap::Array{Float64,2})
    m,n = size(gmap)

    gmap_u = gmap[1:(m-1), :]
    gmap_d = gmap[2:m, :]

    gmap_ud = convert(Array{Bool,2}, gmap_u .* gmap_d)
    gmap_ud = vcat(gmap_ud, falses(1,n))

    s_vert = find(gmap_ud)
    t_vert = s_vert .+ 1

    return (s_vert, t_vert)
end

function _neighbors_diag1(gmap::Array{Float64,2})
    m,n = size(gmap)

    gmap_ul = gmap[1:(m-1), 1:(n-1)]
    gmap_dr = gmap[2:m, 2:n]

    gmap_uldr = convert(Array{Bool,2}, gmap_ul .* gmap_dr)
    gmap_uldr = vcat(hcat(gmap_uldr, falses(m-1)), falses(1,n))

    s_dr = find(gmap_uldr)
    t_dr = s_dr .+ (m+1)

    return (s_dr, t_dr)
end

function _neighbors_diag2(gmap::Array{Float64,2})
    m,n = size(gmap)

    gmap_ur = gmap[1:(m-1), 2:n]
    gmap_dl = gmap[2:m, 1:(n-1)]

    gmap_urdl = convert(Array{Bool,2}, gmap_ur .* gmap_dl)
    gmap_urdl = vcat(hcat(falses(m-1), gmap_urdl), falses(1,n))

    s_dl = find(gmap_urdl)
    t_dl = s_dl .- (m-1)

    return (s_dl, t_dl)
end

# ==================================
# FocalPoints
# ==================================
type FocalPoints
    data::Array{Int,2}
    polygons::Bool

    function FocalPoints(filename::String)
        ftype = _guess_file_type(filename)

        # TODO: implement resampling to match habitat size
        if ftype == FILE_TYPE_AAGRID
            grid = FileAAGrid(filename, "int32")
            hdr = grid.hdr
            data = grid.data
            nodata = (nothing != hdr.nodata) ? -9999 : hdr.nodata
            data[find(x->(x==nodata), data)] = 0    # ignore nodata values
            data[find(x->(x < 0), data)] = 0        # ignore negative point ids

            rows,cols,vals = findnz(data)
            data = hcat(vals, rows, cols)
            data = sortrows(data)

            return new(data, polygons(data))

        elseif ftype == FILE_TYPE_TXTLIST
            error("txtlist format for focal points not supported yet")
        else
            error("Unknown file type for focal points")
        end
    end

    polygons(data::Array{Int,2}) = (length(unique(data[:,1])) < size(data)[1])
end

point_ids(fp::FocalPoints) = fp.data[:,1]
num_points(fp::FocalPoints) = size(fp.data)[1]

function show(io::IO, fp::FocalPoints)
    print(io, "FocalPoints ")
    show(fp.data)
end

# Checks to see if there are at least two focal points in a given component.
function has_focal_points(fp::FocalPoints, hab::Habitat, comp::Int)
    numfp = num_points(fp)
    node_map = hab.node_map
    fpl = fp.data
    comp_nodes = hab.ccs[comp]

    for fpid1 in 1:numfp
        # find the node id of fp 1
        nid1 = node_map[fpl[fpid1, 2:3]...]
        ((nid1 == 0) || !(nid1 in comp_nodes)) && continue
        for fpid2 in (fpid1+1):numfp
            # find the node id of fp 2
            nid2 = node_map[fpl[fpid2, 2:3]...]
            # return true if both are valid nodes, and both belong to component comp
            (nid2 > 0) && (nid2 in comp_nodes) && (return true)
        end
    end
    return false
end

function focal_points(fp::FocalPoints, hab::Habitat, comp::Int)
    numfp = num_points(fp)
    node_map = hab.node_map
    fpl = fp.data
    comp_nodes = hab.ccs[comp]

    for fpid1 in 1:numfp
        # find the node id of fp 1
        nid1 = node_map[fpl[fpid1, 2:3]...]
        ((nid1 == 0) || !(nid1 in comp_nodes)) && continue
        for fpid2 in (fpid1+1):numfp
            # find the node id of fp 2
            nid2 = node_map[fpl[fpid2, 2:3]...]
            # return true if both are valid nodes, and both belong to component comp
            (nid2 > 0) && (nid2 in comp_nodes) && produce((nid1, nid2))
        end
        produce((nid1, -1))
    end
end

type InclExclPairs
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.
#

include("FMI2.jl")

# Comfort functions for fmi2 functions using fmi2Components

"""
Set the DebugLogger for the FMU

For more information call ?fmi2SetDebugLogging
"""
function fmi2SetDebugLogging(c::fmi2Component)
    fmi2SetDebugLogging(c, fmi2False, Unsigned(0), C_NULL)
end

"""
Set the start and end time for a simulation, can only be called after fmi2Instantiate and before fmi2EnterInitializationMode

For more information call ?fmi2SetupExperiment
"""
function fmi2SetupExperiment(c::fmi2Component,
    toleranceDefined::Bool,
                tolerance::Real,
                startTime::Real,
                stopTimeDefined::Bool,
                stopTime::Real)
    c.fmu.t = startTime
    fmi2SetupExperiment(c, fmi2Boolean(toleranceDefined), fmi2Real(tolerance),
                            fmi2Real(startTime), fmi2Boolean(stopTimeDefined), fmi2Real(stopTime))
end

"""
Setup the simulation but without defining all of the parameters

For more information call ?fmi2SetupExperiment
"""
function fmi2SetupExperiment(c::fmi2Component, startTime::Real = 0.0, stopTime::Real = startTime; tolerance::Real = 0.0)

    toleranceDefined = (tolerance > 0.0)
    stopTimeDefined = (stopTime > startTime)

    fmi2SetupExperiment(c, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
end
"""
Get the values of an array of fmi2Real variables

For more information call ?fmi2GetReal
"""
function fmi2GetReal(c::fmi2Component, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Real, nvr)
    fmi2GetReal!(c, vr, nvr, values)
    values
end

"""
Get the value of a fmi2Real variable

For more information call ?fmi2GetReal
"""
function fmi2GetReal(c::fmi2Component, vr::fmi2ValueReference)
    values = fmi2GetReal(c, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Real variables by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal(c::fmi2Component, vr_string::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2GetReal(c, vr)
    end
end
"""
Get the value of a fmi2Real variable by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal(c::fmi2Component, vr_string::String)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2GetReal(c, vr)
    end
end
"""
Get the values of an array of fmi2Real variables

For more information call ?fmi2GetReal
"""
function fmi2GetReal!(c::fmi2Component, vr::Array{fmi2ValueReference}, values::Array{<:Real})
    vars = zeros(fmi2Real, length(values))
    if length(values) != length(vr)
        display("[ERROR]: Number of value references and in place array doesn't match")
    else
        fmi2GetReal!(c, vr, Csize_t(length(values)), vars)
    end
    values[:] = vars
end
"""
Get the values of an array of fmi2Real variables by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal!(c::fmi2Component, vr_string::Array{String}, values::Array{<:Real})
    vars = zeros(fmi2Real, length(values))
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    fmi2GetReal!(c, vr, values)
end
"""
Set the values of an array of fmi2Real variables

For more information call ?fmi2SetReal
"""
function fmi2SetReal(c::fmi2Component, vr::Array{fmi2ValueReference}, value::Array{<:Real})
    nvr = Csize_t(length(vr))
    fmi2SetReal(c, vr, nvr, Array{fmi2Real}(value))
end

"""
Set the value of a fmi2Real variable

For more information call ?fmi2SetReal
"""
function fmi2SetReal(c::fmi2Component, vr::fmi2ValueReference, value::Real)
    fmi2SetReal(c, [vr], [value])
end
"""
Set the values of an array of fmi2Real variables by variable name

For more information call ?fmi2SetReal
"""
function fmi2SetReal(c::fmi2Component, vr_string::Array{String}, value::Array{<:Real})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2SetReal(c, vr, value)
    end
end
"""
Set the value of a fmi2Real variable by variable name

For more information call ?fmi2SetReal
"""
function fmi2SetReal(c::fmi2Component, vr_string::String, value::Real)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2SetReal(c, vr, value)
    end
end
"""
Get the values of an array of fmi2Integer variables

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(c::fmi2Component, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Integer, nvr)

    fmi2GetInteger!(c, vr, nvr, values)
    values
end

"""
Get the value of a fmi2Integer variable

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(c::fmi2Component, vr::fmi2ValueReference)
    values = fmi2GetInteger(c, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(c::fmi2Component, vr_string::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2GetInteger(c, vr)
    end
end
"""
Get the value of a fmi2Integer variable by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(c::fmi2Component, vr_string::String)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2GetInteger(c, vr)
    end
end
"""
Get the values of an array of fmi2Integer variables

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger!(c::fmi2Component, vr::Array{fmi2ValueReference}, values::Array{<:Integer})
    vars = zeros(fmi2Integer, length(values))
    if length(values) != length(vr)
        display("[ERROR]: Number of value references and in place array doesn't match")
    else
        fmi2GetInteger!(c, vr, Csize_t(length(values)), vars)
    end
    values[:] = vars
end
"""
Get the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger!(c::fmi2Component, vr_string::Array{String}, values::Array{<:Integer})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    fmi2GetInteger!(c, vr, values)
end
"""
Set the values of an array of fmi2Integer variables

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(c::fmi2Component, vr::Array{fmi2ValueReference},value::Array{<:Integer})
    nvr = Csize_t(length(vr))
    fmi2SetInteger(c, vr, nvr, Array{fmi2Integer}(value))
end
"""
Get the value of a fmi2Integer variable

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(c::fmi2Component, vr::fmi2ValueReference, value::Integer)
    fmi2SetInteger(c,[vr], [value])
end
"""
Set the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(c::fmi2Component, vr_string::Array{String}, value::Array{<:Integer})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2SetInteger(c, vr, value)
    end
end
"""
Set the value of a fmi2Integer variable by variable name

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(c::fmi2Component, vr_string::String, value::Integer)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2SetInteger(c, vr, value)
    end
end
"""
Get the values of an array of fmi2Boolean variables

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(c::fmi2Component, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    value = Array{fmi2Boolean}(undef, nvr)
    fmi2GetBoolean!(c, vr, nvr, value)
    value
end
"""
Get the value of a fmi2Boolean variable

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(c::fmi2Component, vr::fmi2ValueReference)
    values = fmi2GetBoolean(c, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(c::fmi2Component, vr_string::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2GetBoolean(c, vr)
    end
end
"""
Get the value of a fmi2Boolean variable by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(c::fmi2Component, vr_string::String)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2GetBoolean(c, vr)
    end
end
"""
Get the values of an array of fmi2Boolean variables

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean!(c::fmi2Component, vr::Array{fmi2ValueReference}, values::Array{Bool})
    vars = Array{fmi2Boolean}(undef, length(values))

    if length(values) != length(vr)
        display("[ERROR]: Number of value references and in place array doesn't match")
    else
        fmi2GetBoolean!(c, vr, Csize_t(length(values)), vars)
    end
    values[:] = vars
end
"""
Get the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean!(c::fmi2Component, vr_string::Array{String}, values::Array{Bool})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    fmi2GetBoolean!(c, vr, values)
end
"""
Set the values of an array of fmi2Boolean variables

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(c::fmi2Component, vr::Array{fmi2ValueReference}, value::Array{Bool})
    nvr = Csize_t(length(vr))
    fmi2SetBoolean(c, vr, nvr, Array{fmi2Boolean}(value))
end
"""
Set the value of a fmi2Boolean variable

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(c::fmi2Component, vr::fmi2ValueReference, value::Bool)
    fmi2SetBoolean(c, [vr], [value])
end
"""
Set the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(c::fmi2Component, vr_string::Array{String}, value::Array{Bool})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2SetBoolean(c, vr, value)
    end
end
"""
Set the value of a fmi2Boolean variable by variable name

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(c::fmi2Component, vr_string::String, value::Bool)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2SetBoolean(c, vr, value)
    end
end
"""
Get the values of an array of fmi2String variables

For more information call ?fmi2GetString
"""
function fmi2GetString(c::fmi2Component, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    value = Vector{Ptr{Cchar}}(undef, nvr)

    fmi2GetString!(c, vr, nvr, value)
    unsafe_string.(value)
end
"""
Get the value of a fmi2String variable

For more information call ?fmi2GetString
"""
function fmi2GetString(c::fmi2Component, vr::fmi2ValueReference)
    values = fmi2GetString(c, [vr])
    values[1]
end
"""
Get the values of an array of fmi2String variables by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString(c::fmi2Component, vr_string::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2GetString(c, vr)
    end
end
"""
Get the value of a fmi2String variable by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString(c::fmi2Component, vr_string::String)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2GetString(c, vr)
    end
end
"""
Get the values of an array of fmi2String variables

For more information call ?fmi2GetString
"""
function fmi2GetString!(c::fmi2Component, vr::Array{fmi2ValueReference}, values::Array{String})
    vars = Vector{Ptr{Cchar}}(undef, length(vr))
    if length(values) != length(vr)
        display("[ERROR]: Number of value references and in place array doesn't match")
    else
        fmi2GetString!(c, vr, Csize_t(length(values)), vars)
        values[:] = unsafe_string.(vars)
    end
end
"""
Get the values of an array of fmi2String variables by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString!(c::fmi2Component, vr_string::Array{String}, values::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    fmi2GetString!(c, vr, values)
end
"""
Set the values of an array of fmi2String variables

For more information call ?fmi2SetString
"""
function fmi2SetString(c::fmi2Component, vr::Array{fmi2ValueReference}, value::Array{String})
    nvr = Csize_t(length(vr))
    ptrs = pointer.(value)
    fmi2SetString(c, vr, nvr, ptrs)
end
"""
Set the values of a fmi2String variable

For more information call ?fmi2SetString
"""
function fmi2SetString(c::fmi2Component, vr::fmi2ValueReference, value::String)
    fmi2SetString(c, [vr], [value])
end
"""
Set the values of an array of fmi2String variables by variable name

For more information call ?fmi2SetString
"""
function fmi2SetString(c::fmi2Component, vr_string::Array{String}, value::Array{String})
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if length(vr) == 0
        display("[Error]: no valueReferences could be converted")
    else
        fmi2SetString(c, vr, value)
    end
end
"""
Set the value of a fmi2String variable by variable name

For more information call ?fmi2SetString
"""
function fmi2SetString(c::fmi2Component, vr_string::String, value::String)
    vr = fmi2String2ValueReference(c.fmu, vr_string)
    if vr == nothing
        display("[ERROR]: valueReference not found")
    else
        fmi2SetString(c, vr, value)
    end
end
"""
Get the pointer to the current FMU state

For more information call ?fmi2GetFMUstate
"""
function fmi2GetFMUstate(c::fmi2Component)
    state = fmi2FMUstate()
    stateRef = Ref(state)
    fmi2GetFMUstate(c, stateRef)
    state = stateRef[]
    state
end
"""
Free the allocated memory for the FMU state

For more information call ?fmi2FreeFMUstate
"""
function fmi2FreeFMUstate(c::fmi2Component, state::fmi2FMUstate)
    stateRef = Ref(state)
    fmi2FreeFMUstate(c, stateRef)
    state = stateRef[]
end

"""
Returns the size of a byte vector the FMU can be stored in

For more information call ?fmi2SerzializedFMUstateSize
"""
function fmi2SerializedFMUstateSize(c::fmi2Component, state::fmi2FMUstate)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi2SerializedFMUstateSize(c, state, sizeRef)
    size = sizeRef[]
end

"""
Serialize the data in the FMU state pointer

For more information call ?fmi2SerzializeFMUstate
"""
function fmi2SerializeFMUstate(c::fmi2Component, state::fmi2FMUstate)
    size = fmi2SerializedFMUstateSize(c, state)
    serializedState = Array{fmi2Byte}(undef, size)
    fmi2SerializeFMUstate(c, state, serializedState, size)
    serializedState
end

"""
Deserialize the data in the serializedState fmi2Byte field

For more information call ?fmi2DeSerzializeFMUstate
"""
function fmi2DeSerializeFMUstate(c::fmi2Component, serializedState::Array{fmi2Byte})
    size = length(serializedState)
    state = fmi2FMUstate()
    stateRef = Ref(state)
    fmi2DeSerializeFMUstate(c, serializedState, Csize_t(size), stateRef)
    state = stateRef[]
end

"""
Computes directional derivatives

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(c::fmi2Component,
                                      vUnknown_ref::Array{fmi2ValueReference},
                                      vKnown_ref::Array{fmi2ValueReference},
                                      dvKnown::Array{fmi2Real} = Array{fmi2Real}([]))

    nKnown = Csize_t(length(vKnown_ref))
    nUnknown = Csize_t(length(vUnknown_ref))

    if length(dvKnown) < nKnown
        dvKnown = ones(fmi2Real, nKnown)
    end

    dvUnknown = zeros(fmi2Real, nUnknown)

    fmi2GetDirectionalDerivative!(c, vUnknown_ref, nUnknown, vKnown_ref, nKnown, dvKnown, dvUnknown)

    dvUnknown
end
"""
Computes directional derivatives

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(c::fmi2Component,
                                      vUnknown_ref::fmi2ValueReference,
                                      vKnown_ref::fmi2ValueReference,
                                      dvKnown::fmi2Real = 1.0,
                                      dvUnknown::fmi2Real = 1.0)

    fmi2GetDirectionalDerivative(c, [vUnknown_ref], [vKnown_ref], [dvKnown])[1]
end

# CoSimulation specific functions

"""
The computation of a time step is started.

For more information call ?fmi2DoStep
"""
function fmi2DoStep(c::fmi2Component, currentCommunicationPoint::Real, communicationStepSize::Real, noSetFMUStatePriorToCurrentPoint::Bool = true)
    fmi2DoStep(c, fmi2Real(currentCommunicationPoint), fmi2Real(communicationStepSize), fmi2Boolean(noSetFMUStatePriorToCurrentPoint))
end
"""
The computation of a time step is started.

For more information call ?fmi2DoStep
"""
function fmi2DoStep(c::fmi2Component, communicationStepSize::Real)
    fmi2DoStep(c, fmi2Real(c.fmu.t), fmi2Real(communicationStepSize), fmi2True)
    fmu2.t += communicationStepSize
end

# Model Exchange specific functions
"""
Set a new time instant

For more information call ?fmi2SetTime
"""
function fmi2SetTime(c::fmi2Component, time::Real)
    fmi2SetTime(c, fmi2Real(time))
end

"""
Set a new (continuous) state vector

For more information call ?fmi2SetContinuousStates
"""
function fmi2SetContinuousStates(c::fmi2Component, x::Union{Array{Float32}, Array{Float64}})
    nx = Csize_t(length(x))
    fmi2SetContinuousStates(c, Array{fmi2Real}(x), nx)
end
"""
Returns the next discrete states

For more information call ?fmi2NewDiscretestates
"""
function fmi2NewDiscreteStates(c::fmi2Component)
    eventInfo = fmi2EventInfo()
    fmi2NewDiscreteStates(c, eventInfo)
    eventInfo
end
"""
This function must be called by the environment after every completed step

For more information call ?fmi2CompletedIntegratorStep
"""
function fmi2CompletedIntegratorStep(c::fmi2Component,
                                     noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    enterEventMode = fmi2Boolean(false)
    terminateSimulation = fmi2Boolean(false)
    status = fmi2CompletedIntegratorStep!(c,
                                         noSetFMUStatePriorToCurrentPoint,
                                         enterEventMode,
                                         terminateSimulation)
    (status, enterEventMode, terminateSimulation)
end

"""
Compute state derivatives at the current time instant and for the current states.

For more information call ?fmi2GetDerivatives
"""
function  fmi2GetDerivatives(c::fmi2Component)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    derivatives = zeros(fmi2Real, nx)
    fmi2GetDerivatives(c, derivatives, nx)
    derivatives
end
"""
Returns the event indicators of the FMU

For more information call ?fmi2GetEventIndicators
"""
function fmi2GetEventIndicators(c::fmi2Component)
    ni = Csize_t(c.fmu.modelDescription.numberOfEventIndicators)
    eventIndicators = zeros(fmi2Real, ni)
    fmi2GetEventIndicators(c, eventIndicators, ni)
    eventIndicators
end
"""
Return the new (continuous) state vector x

For more information call ?fmi2GetContinuousStates
"""
function fmi2GetContinuousStates(c::fmi2Component)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    x = zeros(fmi2Real, nx)
    fmi2GetContinuousStates(c, x, nx)
    x
end

"""
Return the new (continuous) state vector x

For more information call ?fmi2GetNominalsOfContinuousStates
"""
function fmi2GetNominalsOfContinuousStates(c::fmi2Component)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    x = zeros(fmi2Real, nx)
    fmi2GetNominalsOfContinuousStates(c, x, nx)
    x
end

"""
Starts a simulation of the fmu instance for the matching fmu type, if both types are available the CoSimulation Simulation is started
"""
function fmi2Simulate(c::fmi2Component, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0, recordValues::Union{Array{fmi2ValueReference}, Array{String}} = [])

    if c.fmu.fmuType == fmi2CoSimulation::fmi2Type
        fmi2SimulateCS(c, dt, t_start, t_stop, recordValues)
    elseif c.fmu.fmuType == fmi2ModelExchange::fmi2Type
        fmi2SimulateME(c, dt, t_start, t_stop)
    else
        error(unknownFMUType)
    end
end
"""
Starts a simulation of the CoSimulation FMU instance
"""
function fmi2SimulateCS(c::fmi2Component, dt::Real, t_start::Real, t_stop::Real, recordValues::Array{fmi2ValueReference} = [])

    fmi2SetupExperiment(c, t_start, t_stop)
    fmiEnterInitializationMode(c)
    fmiExitInitializationMode(c)

    t = t_start

    sd = nothing
    record = length(recordValues) > 0

    numDigits = length(string(round(Integer, 1/dt)))

    if record
        sd = fmi2SimulationResult()
        sd.valueReferences = recordValues
        sd.dataPoints = []
        sd.fmu = c.fmu

        values = fmi2GetReal(c, sd.valueReferences)
        push!(sd.dataPoints, (t, values...))

        while t < t_stop
            fmiDoStep(c, t, dt)
            t = round(t + dt, digits=numDigits)

            values = fmi2GetReal(c, sd.valueReferences)
            push!(sd.dataPoints, (t, values...))
        end
    else
        while t < t_stop
            fmiDoStep(c, t, dt)
            t = round(t + dt, digits=numDigits)
        end
    end

    sd
end
"""
Starts a simulation of the CoSimulation FMU instance
"""
function fmi2SimulateCS(c::fmi2Component, dt::Real, t_start::Real, t_stop::Real, recordValues::Array{String})
    vr = fmi2String2ValueReference(c.fmu, recordValues)
    fmi2SimulateCS(c, dt, t_start, t_stop, vr)
end


text:
struct NNPDENS{C1,C2,O} <: NeuralNetDiffEqAlgorithm
    u0::C1
    σᵀ∇u::C2
    opt::O
end

NNPDENS(u0,σᵀ∇u;opt=Flux.ADAM(0.1)) = NNPDENS(u0,σᵀ∇u,opt)

function DiffEqBase.solve(
    prob::TerminalPDEProblem,
    pdealg::NNPDENS;
    verbose = false,
    maxiters = 300,
    trajectories = 100,
    alg,
    pabstol = 1f-6,
    save_everystep = false,
    kwargs...)

    X0 = prob.X0
    tspan = prob.tspan
    d  = length(X0)
    g,f,μ,σ,p = prob.g,prob.f,prob.μ,prob.σ,prob.p

    data = Iterators.repeated((), maxiters)

    #hidden layer
    opt = pdealg.opt
    u0 = pdealg.u0
    σᵀ∇u = pdealg.σᵀ∇u
    ps = Flux.params(u0, σᵀ∇u)

    function F(h, p, t)
        u =  h[end]
        X =  h[1:end-1].data
        _σᵀ∇u = σᵀ∇u([X;t])
        _f = -f(X, u, _σᵀ∇u, p, t)
        Flux.Tracker.collect(vcat(μ(X,p,t),[_f]))
    end

    function G(h, p, t)
        X = h[1:end-1].data
        _σᵀ∇u = σᵀ∇u([X;t])'
        Flux.Tracker.collect(vcat(σ(X,p,t),_σᵀ∇u))
    end

    function neural_sde(init_cond, F, G, tspan, args...; kwargs...)
        noise = Flux.Tracker.collect(zeros(Float32,d+1,d))
        prob = SDEProblem(F, G, init_cond, tspan, noise_rate_prototype=noise)
        map(1:trajectories) do j #TODO add Ensemble Simulation
            predict_ans = solve(prob,  args...; kwargs...)[end]
            (X,u) = (predict_ans[1:end-1], predict_ans[end])
        end
    end

    n_sde = init_cond->neural_sde(init_cond,F,G,tspan,alg;kwargs...)

    function predict_n_sde()
        _u0 = u0(X0)
        init_cond = Flux.Tracker.collect([X0;_u0])
        n_sde(init_cond)
    end

    function loss_n_sde()
        mean(sum(abs2, g(X.data) - u) for (X,u) in predict_n_sde())
    end

    iters = eltype(X0)[]

    cb = function ()
        save_everystep && push!(iters, u0(X0)[1].data)
        l = loss_n_sde()
        verbose && println("Current loss is: $l")
        l < pabstol && Flux.stop()
    end

    Flux.train!(loss_n_sde, ps, data, opt; cb = cb)

    save_everystep ? iters : u0(X0)[1].data
end #pde_solve_ns


text:
@testset "utils" begin
    # Copied from https://github.com/JuliaLang/julia/pull/39352/files
    @testset "constructing a Cholesky factor from a triangular matrix" begin
        A = [1.0 2.0; 3.0 4.0]
        let
            U = UpperTriangular(A)
            C = Cholesky(U)
            @test C isa Cholesky{Float64}
            @test C.U == U
            @test C.L == U'
        end
        let
            L = LowerTriangular(A)
            C = Cholesky(L)
            @test C isa Cholesky{Float64}
            @test C.L == L
            @test C.U == L'
        end
    end
end