text:
const depsfile = normpath(@__DIR__, "..", "deps", "deps.jl")

if isfile(depsfile)
    include(depsfile)
    gccworks = try
        success(`$gcc -v`)
    catch
        false
    end
    if !gccworks
        error("GCC wasn't found. Please make sure that gcc is on the path and run Pkg.build(\"PackageCompiler\")")
    end
else
    error("Package wasn't built correctly. Please run Pkg.build(\"PackageCompiler\")")
end

system_compiler = gcc
executable_ext = Sys.iswindows() ? ".exe" : ""

if Sys.iswindows()
    function run_PATH(cmd)
        bindir = dirname(cmd.exec[1])
        run(setenv(cmd, ["PATH=" * bindir * ";" * ENV["PATH"]]))
    end
else
    const run_PATH = run
end

"""
    static_julia(juliaprog::String; kw_args...)

compiles the Julia file at path `juliaprog` with keyword arguments:

    cprog                            C program to compile (required only when building an executable, if not provided a minimal driver program is used)
    verbose                          increase verbosity
    quiet                            suppress non-error messages
    builddir                         build directory
    outname                          output files basename
    snoopfile                        specify script calling functions to precompile
    clean                            remove build directory
    autodeps                         automatically build required dependencies
    object                           build object file
    shared                           build shared library
    init_shared                      add `init_jl_runtime` and `exit_jl_runtime` to shared library for runtime initialization
    executable                       build executable file
    rmtemp                           remove temporary build files
    copy_julialibs                   copy Julia libraries to build directory
    copy_files                       copy user-specified files to build directory (either `nothing` or a string array)
    release                          build in release mode, implies `-O3 -g0` unless otherwise specified
    Release                          perform a fully automated release build, equivalent to `-atjr`
    sysimage <file>                  start up with the given system image file
    home <dir>                       set location of `julia` executable
    startup_file {yes|no}            load `~/.julia/config/startup.jl`
    handle_signals {yes|no}          enable or disable Julia's default signal handlers
    sysimage_native_code {yes|no}    use native code from system image if available
    compiled_modules {yes|no}        enable or disable incremental precompilation of modules
    depwarn {yes|no|error}           enable or disable syntax and method deprecation warnings
    warn_overwrite {yes|no}          enable or disable method overwrite warnings
    compile {yes|no|all|min}         enable or disable JIT compiler, or request exhaustive compilation
    cpu_target <target>              limit usage of CPU features up to <target> (implies default `--sysimage_native_code=no`)
    optimize {0,1,2,3}               set the optimization level
    debug <level>                    enable / set the level of debug info generation
    inline {yes|no}                  control whether inlining is permitted
    check_bounds {yes|no}            emit bounds checks always or never
    math_mode {ieee,fast}            disallow or enable unsafe floating point optimizations
    cc                               system C compiler
    cc_flags <flags>                 pass custom flags to the system C compiler when building a shared library or executable (either `nothing` or a string array)
"""
function static_julia(
        juliaprog;
        cprog = nothing, verbose = false, quiet = false, builddir = nothing, outname = nothing, snoopfile = nothing,
        clean = false, autodeps = false, object = false, shared = false, init_shared = false, executable = false, rmtemp = false,
        copy_julialibs = false, copy_files = nothing, release = false, Release = false,
        sysimage = nothing, home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )

    cprog == nothing && (cprog = normpath(@__DIR__, "..", "examples", "program.c"))
    builddir == nothing && (builddir = "builddir")
    outname == nothing && (outname = splitext(basename(juliaprog))[1])
    cc == nothing && (cc = system_compiler)

    verbose && quiet && (quiet = false)

    if Release
        autodeps = rmtemp = copy_julialibs = release = true
    end

    if autodeps
        executable && (shared = true)
        shared && (object = true)
    end

    if release
        optimize == nothing && (optimize = "3")
        debug == nothing && (debug = "0")
    end

    if juliaprog != nothing
        juliaprog = abspath(juliaprog)
        isfile(juliaprog) || error("Cannot find file: \"$juliaprog\"")
        quiet || println("Julia program file:\n  \"$juliaprog\"")
    end

    if executable
        cprog = abspath(cprog)
        isfile(cprog) || error("Cannot find file: \"$cprog\"")
        quiet || println("C program file:\n  \"$cprog\"")
    end

    quiet || println("Build directory:\n  \"$builddir\"")

    if [clean, object, shared, executable, rmtemp, copy_julialibs, copy_files] == [false, false, false, false, false, false, nothing]
        quiet || println("Nothing to do")
        return
    end

    if clean
        if isdir(builddir)
            verbose && println("Remove build directory")
            rm(builddir, recursive = true)
        else
            verbose && println("Build directory does not exist")
        end
    end

    if [object, shared, executable, rmtemp, copy_julialibs, copy_files] == [false, false, false, false, false, nothing]
        quiet || println("Clean completed")
        return
    end

    if !isdir(builddir)
        verbose && println("Make build directory")
        mkpath(builddir)
    end

    o_file = outname * ".a"
    s_file = outname * ".$(Libdl.dlext)"
    e_file = outname * executable_ext

    if object
        if snoopfile != nothing
            snoopfile = abspath(snoopfile)
            verbose && println("Executing snoopfile: \"$snoopfile\"")
            precompfile = joinpath(builddir, "precompiled.jl")
            snoop(nothing, nothing, snoopfile, precompfile)
            jlmain = joinpath(builddir, "julia_main.jl")
            open(jlmain, "w") do io
                # this file will get included via @eval Module() include(...)
                # but the juliaprog should be included in the main module
                juliaprog != nothing && println(io, "Base.include(Main, $(repr(relpath(juliaprog, builddir))))")
                println(io, "Base.include(@__MODULE__, $(repr(relpath(precompfile, builddir))))")
            end
            juliaprog = jlmain
        end
        build_object(
            juliaprog, o_file, builddir, verbose,
            sysimage, home, startup_file, handle_signals, sysimage_native_code, compiled_modules,
            depwarn, warn_overwrite, compile, cpu_target, optimize, debug, inline, check_bounds, math_mode
        )
    end

    shared && build_shared(s_file, o_file, init_shared, builddir, verbose, optimize, debug, cc, cc_flags)

    executable && build_exec(e_file, cprog, s_file, builddir, verbose, optimize, debug, cc, cc_flags)

    rmtemp && remove_temp_files(builddir, verbose)

    copy_julialibs && copy_julia_libs(builddir, verbose)

    copy_files != nothing && copy_files_array(copy_files, builddir, verbose, "Copy user-specified files to build directory:")

    quiet || println("All done")
end

function julia_flags(optimize, debug, cc_flags)
    allflags = Base.shell_split(PackageCompiler.allflags())
    bitness_flag = Sys.ARCH == :arm ? "-mbe32" : Sys.ARCH == :aarch64 ? `` : Int == Int32 ? "-m32" : "-m64"
    allflags = `$allflags $bitness_flag`
    optimize == nothing || (allflags = `$allflags -O$optimize`)
    debug == 2 && (allflags = `$allflags -g`)
    cc_flags == nothing || isempty(cc_flags) || (allflags = `$allflags $cc_flags`)
    allflags
end

function build_julia_cmd(
        sysimage, home, startup_file, handle_signals, sysimage_native_code, compiled_modules,
        depwarn, warn_overwrite, compile, cpu_target, optimize, debug, inline, check_bounds, math_mode
    )
    # TODO: `sysimage_native_code` and `compiled_modules` may be removed in future, see: https://github.com/JuliaLang/PackageCompiler.jl/issues/47
    sysimage_native_code == nothing && cpu_target != nothing && (sysimage_native_code = "no")
    compiled_modules == nothing && (compiled_modules = "no")
    # TODO: `startup_file` may be removed in future with `julia-compile`, see: https://github.com/JuliaLang/julia/issues/15864
    startup_file == nothing && (startup_file = "no")
    julia_cmd = `$(Base.julia_cmd())`
    function get_flag_idx(julia_cmd, flag_str)
        findfirst(f->(length(f) > length(flag_str) &&
                      f[1:length(flag_str)]==flag_str), julia_cmd.exec)
    end
    function set_flag(julia_cmd, flag_str, val)
        flag_idx = get_flag_idx(julia_cmd, flag_str)
        if flag_idx != nothing
            julia_cmd.exec[flag_idx] = "$flag_str$val"
        else
            push!(julia_cmd.exec, "$flag_str$val")
        end
    end
    sysimage == nothing || set_flag(julia_cmd, "-J", sysimage)
    home == nothing || set_flag(julia_cmd, "-H=", home)
    startup_file == nothing || set_flag(julia_cmd, "--startup-file=", startup_file)
    handle_signals == nothing || set_flag(julia_cmd, "--handle-signals=", handle_signals)
    sysimage_native_code == nothing || set_flag(julia_cmd, "--sysimage-native-code=", sysimage_native_code)
    compiled_modules == nothing || set_flag(julia_cmd, "--compiled-modules=", compiled_modules)
    depwarn == nothing || set_flag(julia_cmd, "--depwarn=", depwarn)
    warn_overwrite == nothing || set_flag(julia_cmd, "--warn-overwrite=", warn_overwrite)
    compile == nothing || set_flag(julia_cmd, "--compile=", compile)
    cpu_target == nothing || set_flag(julia_cmd, "-C", cpu_target)
    optimize == nothing || set_flag(julia_cmd, "-O", optimize)
    debug == nothing || set_flag(julia_cmd, "-g", debug)
    inline == nothing || set_flag(julia_cmd, "--inline=", inline)
    check_bounds == nothing || set_flag(julia_cmd, "--check-bounds=", check_bounds)
    math_mode == nothing || set_flag(julia_cmd, "--math-mode=", math_mode)
    # Disable incompatible flags
    set_flag(julia_cmd, "--code-coverage=", "none")

    julia_cmd
end

function build_object(
        juliaprog, o_file, builddir, verbose,
        sysimage, home, startup_file, handle_signals, sysimage_native_code, compiled_modules,
        depwarn, warn_overwrite, compile, cpu_target, optimize, debug, inline, check_bounds, math_mode
    )
    # TODO really refactor this :D
    build_object(
        juliaprog, o_file, builddir, verbose;
        sysimage = sysimage, startup_file = startup_file,
        handle_signals = handle_signals, sysimage_native_code = sysimage_native_code,
        compiled_modules = compiled_modules,
        depwarn = depwarn, warn_overwrite = warn_overwrite,
        compile = compile, cpu_target = cpu_target, optimize = optimize,
        debug_level = debug, inline = inline, check_bounds = check_bounds, math_mode = math_mode
    )
end

function build_object(
        juliaprog, o_file, builddir, verbose; julia_flags...
    )
    Sys.iswindows() && (juliaprog != nothing) && (juliaprog = replace(juliaprog, "\\" => "\\\\"))
    command = ExitHooksStart() * InitBase() * InitREPL()
    if juliaprog != nothing
        command *= Include(juliaprog)
    end
    command *= ExitHooksEnd()
    verbose && println("Build static library $(repr(o_file)):\n  $command")
    cd(builddir) do
        run_julia(
            command; julia_flags...,
            output_o = o_file, track_allocation = "none", code_coverage = "none", startup_file = "no",
        )
    end
end

function build_shared(s_file, o_file, init_shared, builddir, verbose, optimize, debug, cc, cc_flags)
    if init_shared
        i_file = joinpath(abspath(builddir), "lib_init.c")
        open(i_file, "w") do io
            print(io, """
        // Julia headers (for initialization and gc commands)
        #include "uv.h"
        #include "julia.h"
        void init_jl_runtime() // alternate name for jl_init_with_image, with hardcoded library name
        {
            // JULIAC_PROGRAM_LIBNAME defined on command-line for compilation
            const char rel_libname[] = JULIAC_PROGRAM_LIBNAME;
            jl_init_with_image(NULL, rel_libname);
        }
        void exit_jl_runtime(int retcode) // alternate name for jl_atexit_hook
        {
            jl_atexit_hook(retcode);
        }
        """
            )
        end
        i_file = `$i_file`
    else
        i_file = ``
    end
    # Prevent compiler from stripping all symbols from the shared lib.
    if Sys.isapple()
        o_file = `-Wl,-all_load $o_file`
    else
        o_file = `-Wl,--whole-archive $o_file -Wl,--no-whole-archive`
    end
    command = `$cc -shared -DJULIAC_PROGRAM_LIBNAME=$(repr(escape_string(joinpath(builddir, s_file)))) -o $s_file $o_file $i_file $(julia_flags(optimize, debug, cc_flags))`
    if Sys.isapple()
        command = `$command -Wl,-install_name,@rpath/$(joinpath(builddir, s_file))`
    elseif Sys.iswindows()
        command = `$command -Wl,--export-all-symbols`
    end
    verbose && println("Build shared library $(repr(s_file)):\n  $command")
    cd(builddir) do
        run_PATH(command)
    end
end

function build_exec(e_file, cprog, s_file, builddir, verbose, optimize, debug, cc, cc_flags)
    command = `$cc -DJULIAC_PROGRAM_LIBNAME=\"$s_file\" -o $e_file $cprog $s_file $(julia_flags(optimize, debug, cc_flags))`
    if Sys.iswindows()
        # functionality doesn't readily exist on this platform
    elseif Sys.isapple()
        command = `$command -Wl,-rpath,@executable_path`
    else
        command = `$command -Wl,-rpath,\$ORIGIN`
    end
    if Int == Int32
        # TODO: this was added because of an error with julia on win32 that suggested this line, it seems to work but I'm not sure if it's correct
        command = `$command -march=pentium4`
    end
    verbose && println("Build executable $(repr(e_file)):\n  $command")
    cd(builddir) do
        run_PATH(command)
    end
end

function remove_temp_files(builddir, verbose)
    verbose && println("Remove temporary files:")
    remove = false
    for tmp in filter(x -> endswith(x, ".a") || startswith(x, "cache_ji_v"), readdir(builddir))
        verbose && println("  $tmp")
        rm(joinpath(builddir, tmp), recursive = true)
        remove = true
    end
    verbose && !remove && println("  none")
end

function copy_files_array(files_array, builddir, verbose, message)
    verbose && println(message)
    copy = false
    for src in files_array
        isfile(src) || error("Cannot find file: \"$src\"")
        dst = joinpath(builddir, basename(src))
        if filesize(src) != filesize(dst) || ctime(src) > ctime(dst) || mtime(src) > mtime(dst)
            verbose && println("  $(basename(src))")
            cp(src, dst, force = true, follow_symlinks = false)
            copy = true
        end
    end
    verbose && !copy && println("  none")
end

function copy_julia_libs(builddir, verbose)
    # TODO: these flags should probably be emitted also by `julia-config.jl` and `compiler_flags.jl`
    shlibdir = Sys.iswindows() ? Sys.BINDIR : joinpath(Sys.BINDIR, Base.LIBDIR)
    private_shlibdir = joinpath(Sys.BINDIR, Base.PRIVATE_LIBDIR)
    libfiles = String[]
    dlext = "." * Libdl.dlext
    for dir in (shlibdir, private_shlibdir)
        if Sys.iswindows() || Sys.isapple()
            append!(libfiles, joinpath.(dir, filter(x -> endswith(x, dlext) && !startswith(x, "sys"), readdir(dir))))
        else
            append!(libfiles, joinpath.(dir, filter(x -> occursin(r"^lib.+\.so(?:\.\d+)*$", x), readdir(dir))))
        end
    end
    filter!(v -> !occursin(r"debug", v), libfiles)
    copy_files_array(libfiles, builddir, verbose, "Copy Julia libraries to build directory:")
end


text:
"Random Variable of Unknown distribution `ω -> f(ω[id], args...)`"
struct URandVar{F, TPL <: Tuple} <: RandVar
  f::F
  args::TPL
  id::ID
  URandVar(f::F, args::TPL = (), id = uid()) where {F, TPL} =
    new{F, TPL}(f, args, id)
end

id(rv::URandVar) = rv.id

"Parameters of the Random Variable"
params(x::URandVar) = x.args

"Function of the Random Variable"
func(x::URandVar) = x.f

"Name of a random variable"
name(rv::URandVar) = string(rv.f)

"`constant(c)` Constant random variable which always outputs `c`"
constant(c) = URandVar(ω -> c)

@inline ppapl(rv::URandVar, ωπ) =  rv.f(ωπ, rv.args...)

@inline (rv::URandVar)(ω::Ω) = apl(rv, ω)

ext: jl


text:
# Inspired from SOSTools
export getmonomialsforcertificate, randpsd, randsos

cfld(x::NTuple{2,Int}, n) = (cld(x[1], n), fld(x[2], n))

# TODO sparse with Newton polytope (Polyhedra.jl for convex hull)
function _getmonomialsforcertificate(X::AbstractVector{<:AbstractMonomial}, sparse=:No)
    if sparse == :No
        # Cheap approximation of the convex hull as the approximation of:
        #
        # z such that mindegree < sum(z) < maxdegree
        # |\
        # |#\ <-------- sum(z) = maxdegree
        # |##\
        # |\<-\-------- sum(z) = mindegree
        # | \##\
        # +---------
        #
        # and:
        #
        # z such that minmultideg < z < maxmultideg
        # | +----+ <--- maxmultidegree
        # | |####|
        # | |####|
        # | +----+
        # | ^---------- minmultidegree
        # +---------
        mindeg, maxdeg = cfld(extdegree(X), 2)
        n = nvariables(X)
        minmultideg, maxmultideg = Vector{Int}(undef, n), Vector{Int}(undef, n)
        for i in 1:n
            exponent_i(m) = exponents(m)[i]
            minmultideg[i] = cld(mapreduce(exponent_i, min, X), 2)
            maxmultideg[i] = fld(mapreduce(exponent_i, max, X), 2)
        end
        monomials(variables(X), mindeg:maxdeg,
                  m -> all(minmultideg .<= exponents(m) .<= maxmultideg))
    else
        error("Not supported yet :(")
    end
end
getmonomialsforcertificate(X::AbstractVector, sparse=:No) = _getmonomialsforcertificate(monovec(X), sparse)

function randpsd(n; r=n, eps=0.1)
    Q = randn(n,n)
    d = zeros(Float64, n)
    d[1:r] = eps .+ abs.(randn(r))
    Q' * Diagonal(d) * Q
end

function _randsos(X::AbstractVector{<:AbstractMonomial}; r=-1, monotype=:Classic, eps=0.1)
    if monotype == :Classic
        x = getmonomialsforcertificate(X)
    elseif monotype == :Gram
        x = X
    else
        throw(ArgumentError("Monotype $monotype not known"))
    end
    n = length(x)
    if r < 0
        r = n
    end
    MatPolynomial(randpsd(n, r=r, eps=eps), x)
end

randsos(X::AbstractVector; kws...) = _randsos(monovec(X); kws...)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis: calculate_interpolation_polynomials,
                calculate_interpolation_polynomial_derivatives

p = :(1 + u + v)
A = [1.0 0.0 0.0; 1.0 1.0 0.0; 1.0 0.0 1.0]
basis = calculate_interpolation_polynomials(p, A)
@test basis[1] == :(1.0 + -u + -v)
@test basis[2] == :(+u)
@test basis[3] == :(+v)

dbasis = calculate_interpolation_polynomial_derivatives(basis, 2)
@test isapprox(dbasis[1,1], -1.0)
@test isapprox(dbasis[2,1], -1.0)
@test isapprox(dbasis[1,2], 1.0)
@test isapprox(dbasis[2,2], 0.0)
@test isapprox(dbasis[1,3], 0.0)
@test isapprox(dbasis[2,3], 1.0)


text: setalarm(employed, vacation) = employed && !vacation

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@testset "flatland" begin
    import DSJulia.Flatland

    @testset "rectangle" begin
        rect = Flatland.Rectangle((1, 1), l=2, w=2)
        square = Flatland.Square((2.0, 1), l=2)

        @test Flatland.ncorners(rect) == Flatland.ncorners(square) == 4
        @test Flatland.center(rect) == (1,1) !=  Flatland.center(square)

        @test Flatland.xlim(rect) == (0, 2)
        @test Flatland.ylim(rect) == Flatland.ylim(square) == (0, 2)

        @test Flatland.area(rect) ≈ Flatland.area(square) ≈ 4

        Flatland.move!(square, (-1, 0))
        @test Flatland.xlim(rect) == Flatland.xlim(square)

        @test (0.1, 0.3) ∈ rect
        @test (7, 0.1) ∉ square

        @test rect ∩ Flatland.Rectangle((2, 2), l=2, w=2)
        @test !(rect ∩ Flatland.Rectangle((4, 2), l=2, w=2))
    end

    @testset "circle" begin
        circle = Flatland.Circle((1, 1), R=1)

        @test Flatland.ncorners(circle) == 0
        @test Flatland.center(circle) == (1, 1)

        @test Flatland.xlim(circle) == (0, 2)
        @test Flatland.ylim(circle) == (0, 2)

        @test Flatland.area(circle) ≈ π

        @test (0.5, 0.5) ∈ circle
        @test (7, 0.1) ∉ circle

        @test circle ∩ Flatland.Circle((2, 2), R=0.5)
        @test !(circle ∩ Flatland.Circle((5, 3), R=0.5))
    end

    @testset "pregular polygons" begin
        pent = Flatland.RegularPolygon((1, 1), 5, R=1)

        @test Flatland.ncorners(pent) == 5
        @test Flatland.center(pent) == (1, 1)

        @test Flatland.xlim(pent)[1] ≥ 0
        @test Flatland.ylim(pent)[1] ≥ 0

        @test Flatland.xlim(pent)[2] ≤ 2
        @test Flatland.ylim(pent)[2] ≤ 2

        @test Flatland.area(Flatland.RegularPolygon((1, 1), 1000_000, R=1)) ≈ π

        @test (0.5, 0.5) ∈ pent
        @test (7, 0.1) ∉ pent

        @test pent ∩ Flatland.RegularPolygon((2, 2), 5, R=1)
        @test !(pent ∩ Flatland.RegularPolygon((8, 1), 5, R=1/2))
end

    @testset "triangle" begin
        triangle = Flatland.Triangle((0, 0), (1, 0), (0.5, 0.5))

        @test Flatland.ncorners(triangle) == 3
        @test all(Flatland.center(triangle) .≈ (0.5, 0.5/3))

        @test all(Flatland.xlim(triangle) .≈ (0, 1))
        @test all(Flatland.ylim(triangle) .≈ (0, 0.5))

        @test Flatland.area(triangle) ≈ 0.25

        @test (0.5, 0.3) ∈ triangle
        @test (0.1, 0.2) ∉ triangle

        @test triangle ∩ Flatland.Triangle((0.25, 0.25), (0.5, 1.5), (1.5, 1.5))
        @test !(triangle ∩ Flatland.Triangle((3, 4), (2, 2), (1.5, 1.5)))
    end
end


text:
function solve_vf_terminal!(ddm::DDM_AbstractVF, dograd)
    evs = value_function(ddm)
    ev = EV(evs)
    dev = dEV(evs)

    exp_trm = exploratory_terminal(statespace(ddm))

     ev[:, :, end    ] .= 0
     ev[:, :, exp_trm] .= 0

     if dograd
         dev[:, :, :, end    ] .= 0
         dev[:, :, :, exp_trm] .= 0
    end
    return nothing
end

# ---------------------------------------------

function solve_vf_infill!(t::DCDPTmpVars, ddm::DDM_AbstractVF, θ, ichar, dograd; kwargs...)

    evs = value_function(ddm)
    wp = statespace(ddm)

    for i in ind_inf(wp)
        idxd, idxs, horzn = dp1space(wp,i), collect(sprimes(wp,i)), _horizon(wp,i)

        tvw = view(t, idxd)
        @views EV0  = view( EV(evs), :, :,    i)
        @views dEV0 = view(dEV(evs), :, :, :, i)

        update_static_payoffs!(tvw, ddm, θ, i, ichar, dograd)
        _ubv = ubV(tvw)
        _ubv .+= discount(ddm) .* view(EV(evs), :, :, idxs)

        if dograd
            fill!(dEV0, 0)
            _dubv = dubVperm(tvw)
            _dubv .+= discount(ddm) .* view(dEV(evs), :, :, :, idxs)
        end

        if horzn == :Finite
            if dograd
                vfit!(EV0, dEV0, tvw, ddm)
            else
                vfit!(EV0, tvw, ddm)
            end

        elseif horzn == :Infinite
            converged, iter, bnds =  solve_inf_vfit_pfit!(EV0, tvw, ddm; kwargs...)
            converged || @warn "Did not converge at state $i after $iter pfit. McQueen-Porteus bnds: $bnds. θ = $θ"

            if dograd
                # TODO: only allows 0-payoff if no action
                _ubv[:,:,1] .= discount(ddm) .* EV0
                gradinf!(dEV0, tvw, ddm)   # note: destroys ubV & dubV
            end
        else
            throw(error("i = $i, horzn = $horzn but must be :Finite or :Infinite"))
        end
    end
    return nothing
end

# ---------------------------------------------

function learningUpdate!(t::DCDPTmpVars, ddm::DDM_AbstractVF, θ, ichar, dograd)
    evs = value_function(ddm)
    wp = statespace(ddm)
    lrn2inf = inf_fm_lrn(wp)
    exp2lrn = exploratory_learning(wp)

    EV1   = view( EV(evs), :, :,    lrn2inf)
    EV2   = view( EV(evs), :, :,    exp2lrn)
    dEV1  = view(dEV(evs), :, :, :, lrn2inf)
    dEV2  = view(dEV(evs), :, :, :, exp2lrn)

    idx_σ = idx_ρ(reward(ddm))
    dEVσ  = view(dEV(evs), :, :, idx_σ, exp2lrn)

    _βΠψ!(t, ddm, θ)
    A_mul_B_md!(EV2, Πψtmp(t), EV1, 2)

    if dograd
        # dubVtilde/dθ = du/dθ[:,:,:,2:dmaxp1] + β * Πψ ⊗ I * dEV/dθ[:,:,:,2:dmaxp1]
        A_mul_B_md!(dEV2, Πψtmp(t), dEV1, 2)

        # ∂EVtilde/∂σ[:,:,2:dmaxp1] = ∂u/∂σ[:,:,2:dmaxp1] + β * dΠψ/dσ ⊗ I * EV[:,:,2:dmaxp1]
        _βΠψdθρ!(t, ddm, θ)
        A_mul_B_md!(dEVσ, Πψtmp(t), EV1, 2)
    end
    return nothing
end

# ---------------------------------------------

function solve_vf_explore!(t::DCDPTmpVars, ddm::DDM_AbstractVF, θ, ichar, dograd; kwargs...)

    evs = value_function(ddm)
    wp = statespace(ddm)
    dmaxp1 = _dmax(wp)+1
    exp2lrn = exploratory_learning(wp)

    # Views of ubV so we can efficiently access them
    tvw = view(t, OneTo(dmaxp1))

     ubV0 = view(ubV(tvw),      :, :,    1)
     ubV1 = view(ubV(tvw),      :, :,    2:dmaxp1)
    dubV0 = view(dubVperm(tvw), :, :, :, 1)
    dubV1 = view(dubVperm(tvw), :, :, :, 2:dmaxp1)
    βEV1  = view( EV(evs),      :, :,    exp2lrn)
    βdEV1 = view(dEV(evs),      :, :, :, exp2lrn)

    # do VFI
    for i in ind_exp(wp)
        ip0 = sprime(wp,i,0)    # s_{t+1} if d_t = 0
        horzn = _horizon(wp,i)

         EV0_ip = view( EV(evs), :, :,    ip0)
        dEV0_ip = view(dEV(evs), :, :, :, ip0)

        # compute u + βEV(d) ∀ d ∈ actionspace(wp,i)
        update_static_payoffs!(tvw, ddm, θ, i, ichar, dograd)
        ubV0 .+= discount(ddm) .* EV0_ip
        ubV1 .+= βEV1 # β already baked in

        if dograd
            dubV0 .+= discount(ddm) .* dEV0_ip
            dubV1 .+= βdEV1  # β already baked in
        end

        if horzn == :Finite
            EV0_i = view(EV(evs), :, :, i)
            if dograd
                dEV0_i = view(dEV(evs), :, :, :, i)
                @views vfit!(EV0_i, dEV0_i, tvw, ddm)
            else
                @views vfit!(EV0_i, tvw, ddm)
            end

        elseif horzn == :Infinite
            converged, iter, bnds =  solve_inf_vfit_pfit!(EV0_ip, tvw, ddm; kwargs...)
            converged || @warn "Did not converge at state $i after $iter pfit. McQueen-Porteus bnds: $bnds. θt = $θt, σ = $σ"

            if dograd
                # TODO: only allows 0-payoff if no action
                ubV(tvw)[:,:,1] .= discount(ddm) .* EV0_ip
                gradinf!(dEV0_ip, tvw, ddm)
            end
        else
            throw(error("i = $i, horzn = $horzn but must be :Finite or :Infinite"))
        end
    end
    return nothing
end

# ---------------------------------------------

function solve_vf_all!(t, ddm::DDM_AbstractVF, θ, ichar, dograd; kwargs...)
    solve_vf_terminal!(   ddm,           dograd)
    solve_vf_infill!(  t, ddm, θ, ichar, dograd; kwargs...)
    learningUpdate!(   t, ddm, θ, ichar, dograd)
    solve_vf_explore!( t, ddm, θ, ichar, dograd; kwargs...)
    return nothing
end

@deprecate solve_vf_terminal!(evs,    ddm,           dograd)             solve_vf_terminal!(  ddm,           dograd)
@deprecate solve_vf_infill!(  evs, t, ddm, θ, ichar, dograd; kwargs...)  solve_vf_infill!( t, ddm, θ, ichar, dograd; kwargs...)
@deprecate learningUpdate!(   evs, t, ddm, θ, ichar, dograd)             learningUpdate!(  t, ddm, θ, ichar, dograd)
@deprecate solve_vf_explore!( evs, t, ddm, θ, ichar, dograd; kwargs...)  solve_vf_explore!(t, ddm, θ, ichar, dograd; kwargs...)
@deprecate solve_vf_all!(     evs, t, ddm, θ, ichar, dograd; kwargs...)  solve_vf_all!(    t, ddm, θ, ichar, dograd; kwargs...)

# need to reset VF to 0 for benchmarking
function solve_vf_all_timing!(t, ddm, θ, ichar, dograd; kwargs...)
    fill!(value_function(ddm), 0)
    solve_vf_all!(t, ddm, θ, ichar, dograd; kwargs...)
    return nothing
end

function solve_vf_infill_timing!(t, ddm, θ, ichar, dograd; kwargs...)
    fill!(value_function(ddm), 0)
    solve_vf_infill!(t, ddm, θ, ichar, dograd; kwargs...)
    return nothing
end

function solve_vf_and_update_itp!(ddm::DDM_AbstractVF, θ, ichar, dograd; kwargs...)
    t = DCDPTmpVars(ddm)
    solve_vf_all!(t, ddm, θ, ichar, dograd; kwargs...)
    vf = value_function(ddm)
    update_interpolation!(vf, dograd)
    return nothing
end

solve_vf_and_update_itp!(args...; kwargs...) = nothing


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

mutable struct HloInstructionProto_SliceDimensions <: ProtoType
    start::Int64
    limit::Int64
    stride::Int64
    HloInstructionProto_SliceDimensions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloInstructionProto_SliceDimensions

mutable struct HloInstructionProto <: ProtoType
    name::AbstractString
    opcode::AbstractString
    shape::ShapeProto
    metadata::OpMetadata
    literal::LiteralProto
    parameter_number::Int64
    fusion_kind::AbstractString
    tuple_index::Int64
    dimensions::Base.Vector{Int64}
    window::Window
    convolution_dimension_numbers::ConvolutionDimensionNumbers
    feature_group_count::Int64
    slice_dimensions::Base.Vector{HloInstructionProto_SliceDimensions}
    exponent_bits::Int32
    mantissa_bits::Int32
    dynamic_slice_sizes::Base.Vector{Int64}
    padding_config::PaddingConfig
    outfeed_config::Array{UInt8,1}
    distribution::Int32
    epsilon::Float32
    feature_index::Int64
    channel_id::Int64
    infeed_config::Array{UInt8,1}
    custom_call_target::AbstractString
    custom_call_opaque::AbstractString
    outfeed_shape::ShapeProto
    dot_dimension_numbers::DotDimensionNumbers
    fft_type::Int32
    fft_length::Base.Vector{Int64}
    gather_dimension_numbers::GatherDimensionNumbers
    gather_slice_sizes::Base.Vector{Int64}
    channel_name::AbstractString
    cost_estimate_ns::Int64
    id::Int64
    operand_ids::Base.Vector{Int64}
    control_predecessor_ids::Base.Vector{Int64}
    called_computation_ids::Base.Vector{Int64}
    sharding::OpSharding
    backend_config::AbstractString
    replica_groups::Base.Vector{ReplicaGroup}
    all_reduce_id::Int64
    cross_replica_sum_barrier::AbstractString
    is_host_transfer::Bool
    scatter_dimension_numbers::ScatterDimensionNumbers
    precision_config::PrecisionConfig
    source_target_pairs::Base.Vector{SourceTarget}
    domain_entry_sharding::OpSharding
    domain_exit_sharding::OpSharding
    constrain_layout::Bool
    operand_shapes_with_layout::Base.Vector{ShapeProto}
    HloInstructionProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloInstructionProto
const __fnum_HloInstructionProto = Int[1,2,3,7,8,9,11,13,14,15,16,50,17,18,19,20,21,22,23,24,25,26,27,28,53,29,30,31,32,33,34,41,42,35,36,37,38,40,43,49,45,46,47,48,51,52,54,55,56,57]
const __pack_HloInstructionProto = Symbol[:dimensions,:dynamic_slice_sizes,:fft_length,:gather_slice_sizes,:operand_ids,:control_predecessor_ids,:called_computation_ids]
meta(t::Type{HloInstructionProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_HloInstructionProto, ProtoBuf.DEF_VAL, true, __pack_HloInstructionProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct HloComputationProto <: ProtoType
    name::AbstractString
    instructions::Base.Vector{HloInstructionProto}
    program_shape::ProgramShapeProto
    id::Int64
    root_id::Int64
    HloComputationProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloComputationProto
const __fnum_HloComputationProto = Int[1,2,4,5,6]
meta(t::Type{HloComputationProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_HloComputationProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct HloScheduleProto_InstructionSequence <: ProtoType
    instruction_ids::Base.Vector{Int64}
    HloScheduleProto_InstructionSequence(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloScheduleProto_InstructionSequence
const __pack_HloScheduleProto_InstructionSequence = Symbol[:instruction_ids]
meta(t::Type{HloScheduleProto_InstructionSequence}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_HloScheduleProto_InstructionSequence, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct HloScheduleProto_SequencesEntry <: ProtoType
    key::Int64
    value::HloScheduleProto_InstructionSequence
    HloScheduleProto_SequencesEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloScheduleProto_SequencesEntry (mapentry)

mutable struct HloScheduleProto <: ProtoType
    sequences::Base.Dict{Int64,HloScheduleProto_InstructionSequence} # map entry
    HloScheduleProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloScheduleProto

mutable struct HloInputOutputAliasProto_AliasEntryProto <: ProtoType
    output_shape_index::Base.Vector{Int64}
    parameter_number::Int64
    parameter_shape_index::Base.Vector{Int64}
    HloInputOutputAliasProto_AliasEntryProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloInputOutputAliasProto_AliasEntryProto
const __pack_HloInputOutputAliasProto_AliasEntryProto = Symbol[:output_shape_index,:parameter_shape_index]
meta(t::Type{HloInputOutputAliasProto_AliasEntryProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_HloInputOutputAliasProto_AliasEntryProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

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

mutable struct DynamicParameterBindingProto_Binding <: ProtoType
    dynamic_param_num::Int64
    dynamic_param_index::Base.Vector{Int64}
    target_param_num::Int64
    target_param_index::Base.Vector{Int64}
    target_param_dim_num::Int64
    DynamicParameterBindingProto_Binding(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DynamicParameterBindingProto_Binding
const __pack_DynamicParameterBindingProto_Binding = Symbol[:dynamic_param_index,:target_param_index]
meta(t::Type{DynamicParameterBindingProto_Binding}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_DynamicParameterBindingProto_Binding, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

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

mutable struct HloModuleProto <: ProtoType
    name::AbstractString
    entry_computation_name::AbstractString
    entry_computation_id::Int64
    computations::Base.Vector{HloComputationProto}
    host_program_shape::ProgramShapeProto
    id::Int64
    schedule::HloScheduleProto
    input_output_alias::HloInputOutputAliasProto
    dynamic_parameter_binding::DynamicParameterBindingProto
    HloModuleProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloModuleProto
const __fnum_HloModuleProto = Int[1,2,6,3,4,5,7,8,9]
meta(t::Type{HloModuleProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_HloModuleProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct LogicalBufferProto_Location <: ProtoType
    computation_name::AbstractString
    instruction_name::AbstractString
    shape_index::Base.Vector{Int64}
    LogicalBufferProto_Location(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LogicalBufferProto_Location
const __pack_LogicalBufferProto_Location = Symbol[:shape_index]
meta(t::Type{LogicalBufferProto_Location}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_LogicalBufferProto_Location, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct LogicalBufferProto <: ProtoType
    id::Int64
    size::Int64
    defined_at::LogicalBufferProto_Location
    color::Int64
    LogicalBufferProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LogicalBufferProto

mutable struct BufferAllocationProto_Assigned <: ProtoType
    logical_buffer_id::Int64
    offset::Int64
    size::Int64
    BufferAllocationProto_Assigned(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BufferAllocationProto_Assigned

mutable struct BufferAllocationProto <: ProtoType
    index::Int64
    size::Int64
    is_thread_local::Bool
    is_tuple::Bool
    is_entry_computation_parameter::Bool
    is_constant::Bool
    parameter_number::Int64
    parameter_shape_index::Base.Vector{Int64}
    maybe_live_out::Bool
    color::Int64
    assigned::Base.Vector{BufferAllocationProto_Assigned}
    BufferAllocationProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BufferAllocationProto
const __fnum_BufferAllocationProto = Int[1,2,3,11,5,12,6,10,7,8,9]
const __pack_BufferAllocationProto = Symbol[:parameter_shape_index]
meta(t::Type{BufferAllocationProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_BufferAllocationProto, ProtoBuf.DEF_VAL, true, __pack_BufferAllocationProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_HeapSimulatorTrace_Event_Kind <: ProtoEnum
    ALLOC::Int32
    FREE::Int32
    SHARE_WITH::Int32
    __enum_HeapSimulatorTrace_Event_Kind() = new(0,1,2)
end #struct __enum_HeapSimulatorTrace_Event_Kind
const HeapSimulatorTrace_Event_Kind = __enum_HeapSimulatorTrace_Event_Kind()

mutable struct HeapSimulatorTrace_Event <: ProtoType
    kind::Int32
    buffer_id::Int64
    computation_name::AbstractString
    instruction_name::AbstractString
    share_with_canonical_id::Int64
    HeapSimulatorTrace_Event(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HeapSimulatorTrace_Event

mutable struct HeapSimulatorTrace <: ProtoType
    events::Base.Vector{HeapSimulatorTrace_Event}
    whole_module_simulation::Bool
    HeapSimulatorTrace(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HeapSimulatorTrace

mutable struct HloModuleGroupProto <: ProtoType
    name::AbstractString
    hlo_modules::Base.Vector{HloModuleProto}
    HloModuleGroupProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloModuleGroupProto

mutable struct BufferAssignmentProto_BufferAlias <: ProtoType
    source_buffer_id::Int64
    location::LogicalBufferProto_Location
    BufferAssignmentProto_BufferAlias(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BufferAssignmentProto_BufferAlias

mutable struct BufferAssignmentProto <: ProtoType
    logical_buffers::Base.Vector{LogicalBufferProto}
    buffer_aliases::Base.Vector{BufferAssignmentProto_BufferAlias}
    buffer_allocations::Base.Vector{BufferAllocationProto}
    heap_simulator_traces::Base.Vector{HeapSimulatorTrace}
    BufferAssignmentProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BufferAssignmentProto

mutable struct HloProto <: ProtoType
    hlo_module::HloModuleProto
    buffer_assignment::BufferAssignmentProto
    HloProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloProto
const __fnum_HloProto = Int[1,3]
meta(t::Type{HloProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_HloProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct HloSnapshot <: ProtoType
    hlo::HloProto
    arguments::Base.Vector{LiteralProto}
    result::LiteralProto
    execution_platform::AbstractString
    HloSnapshot(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HloSnapshot

export HloInstructionProto_SliceDimensions, HloInstructionProto, HloComputationProto, HloScheduleProto_InstructionSequence, HloScheduleProto_SequencesEntry, HloScheduleProto, HloInputOutputAliasProto_AliasEntryProto, HloInputOutputAliasProto, DynamicParameterBindingProto_Binding, DynamicParameterBindingProto, HloModuleProto, LogicalBufferProto_Location, LogicalBufferProto, BufferAllocationProto_Assigned, BufferAllocationProto, HeapSimulatorTrace_Event_Kind, HeapSimulatorTrace_Event, HeapSimulatorTrace, HloModuleGroupProto, BufferAssignmentProto_BufferAlias, BufferAssignmentProto, HloProto, HloSnapshot
# mapentries: "HloScheduleProto_SequencesEntry" => ("Int64", "HloScheduleProto_InstructionSequence")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

typealias MortarElements3D Union{Tri3, Tri6, Quad4}

function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
    return p - dot(p-x0, n0)*n0
end

function inv3(P::Matrix)
    n, m = size(P)
    @assert n == m == 3
    a, b, c, d, e, f, g, h, i = P
    A = e*i - f*h
    B = -d*i + f*g
    C = d*h - e*g
    D = -b*i + c*h
    E = a*i - c*g
    F = -a*h + b*g
    G = b*f - c*e
    H = -a*f + c*d
    I = a*e - b*d
    return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I]
end

function vertex_inside_polygon(q, P; atol=1.0e-6)
    N = length(P)
    angle = 0.0
    for i=1:N
        A = P[i] - q
        B = P[mod(i,N)+1] - q
        c = norm(A)*norm(B)
        isapprox(c, 0.0; atol=atol) && return true
        cosa = dot(A,B)/c
        isapprox(cosa, 1.0; atol=atol) && return false
        isapprox(cosa, -1.0; atol=atol) && return true
        try
            angle += acos(cosa)
        catch
            info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
            info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
            info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
            info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
            rethrow()
        end
    end
    return isapprox(angle, 2*pi; atol=atol)
end

function calculate_centroid(P)
    N = length(P)
    P0 = P[1]
    areas = [norm(1/2*cross(P[i]-P0, P[mod(i,N)+1]-P0)) for i=2:N]
    centroids = [1/3*(P0+P[i]+P[mod(i,N)+1]) for i=2:N]
    C = 1/sum(areas)*sum(areas.*centroids)
    return C
end

function get_cells(P, C)
    N = length(P)
    cells = Vector[]
    # shared edge etc.
    N < 3 && return cells
    # trivial case, polygon already triangle / quadrangle
    #N == 3 && return Vector[P]
    #N == 4 && return Vector[P]
    #V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N])
    #A = 1/2*abs(dot(n, V))
    #info("A = $A")
    cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N]
    return cells

    maxa = 0.0
    maxj = 0
    for i=1:N
        A = P[i] - C
        B = P[mod(i,N)+1] - C
        theta = acos(dot(A,B)/(norm(A)*norm(B)))
        if theta > maxa
            maxa = theta
            maxj = i
        end
    end
    info("max angle $(maxa/pi*180) at index $maxj, N=$N")
    indices = mod(collect(maxj:maxj+N), N)
    info("indices = $indices")
end

""" Test does P contain q. """
function contains{T}(P::Vector{T}, q::T; check_is_close=true, rtol=1.0e-5)
    if q in P
        return true
    end
    if check_is_close
        for p in P
            if isapprox(p, q; rtol=rtol)
                return true
            end
        end
    end
    return false
end

function get_polygon_clip(xs, xm, n; debug=false)
    # objective: search does line xm1 - xm2 clip xs
    nm = length(xm)
    ns = length(xs)
    P = Vector{Float64}[]

    # 1. test is master point inside slave, if yes, add to clip
    for i=1:nm
        if vertex_inside_polygon(xm[i], xs)
            debug && info("1. $(xm[i]) inside S -> push")
            push!(P, xm[i])
        end
    end

    # 2. test is slave point inside master, if yes, add to clip
    for i=1:ns
        if vertex_inside_polygon(xs[i], xm)
            contains(P, xs[i]) && continue
            debug && info("2. $(xs[i]) inside M -> push")
            push!(P, xs[i])
        end
    end

    for i=1:nm
        # 2. find possible intersection
        xm1 = xm[i]
        xm2 = xm[mod(i,nm)+1]
        #info("intersecting line $xm1 -> $xm2")
        for j=1:ns
            xs1 = xs[j]
            xs2 = xs[mod(j,ns)+1]
            #info("clipping polygon edge $xs1 -> $xs2")
            tnom = dot(cross(xm1-xs1, xm2-xm1), n)
            tdenom = dot(cross(xs2-xs1, xm2-xm1), n)
            isapprox(tdenom, 0) && continue
            t = tnom/tdenom
            (0 <= t <= 1) || continue
            q = xs1 + t*(xs2 - xs1)
            #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
            if vertex_inside_polygon(q, xm)
                contains(P, q) && continue
                debug && info("3. $q inside M -> push")
                push!(P, q)
            end
        end
    end

    return P
end

function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
    element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
    basis(xi) = get_basis(element, xi, time)
    dbasis(xi) = get_dbasis(element, xi, time)
    f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p
    L(theta) = inv3([dbasis(theta[1:2])*x -n0])
#   L2(theta) = inv(ForwardDiff.get_value([dbasis(theta[2:3])*x -n0]))
    # FIXME: for some reason forwarddiff gives NaN's here.
    theta = zeros(3)
    dtheta = zeros(3)
    for i=1:max_iterations
        dtheta = L(theta) * f(theta)
        theta -= dtheta
        if norm(dtheta) < iter_tol
            return theta[1:2], theta[3]
        end
    end

    info("failed to project vertex from auxiliary plane back to surface")
    info("element type: $E")
    info("element connectivity: $(get_connectivity(element))")
    info("auxiliary plane: x0 = $x0, n0 = $n0")
    info("element geometry: $(x.data)")
    info("vertex to project: $p")
    info("parameter vector before giving up: $theta")
    info("increment in parameter vector before giving up: $dtheta")
    info("norm(dtheta) before giving up: $(norm(dtheta))")
    info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
    info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")

    info("iterations:")
    theta = zeros(3)
    dtheta = zeros(3)
    for i=1:max_iterations
        info("iter $i, theta = $theta")
        info("f = $(f(theta))")
        info("L = $(L(theta))")
        dtheta = L(theta) * f(theta)
        info("dtheta = $(dtheta)")
        theta -= dtheta
    end

    error("project_point_to_surface: did not converge in $max_iterations iterations!")
end

function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false)
    normals = Dict{Int64, Vector{Float64}}()
    for element in elements
        conn = get_connectivity(element)
        J = transpose(element([0.0, 0.0], time, Val{:Jacobian}))
        normal = cross(J[:,1], J[:,2])
        for nid in conn
            if haskey(normals, nid)
                normals[nid] += normal
            else
                normals[nid] = normal
            end
        end
    end
    # normalize to unit normal
    S = collect(keys(normals))
    for j in S
        normals[j] /= norm(normals[j])
    end
    if rotate_normals
        for j in S
            normals[j] = -normals[j]
        end
    end
    return normals
end

function check_orientation!(P, n; debug=false)
    C = mean(P)
    np = length(P)
    s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np]
    all(s .< 0) && return
    debug && info("polygon not in ccw order, fixing")
    # project points to new orthogonal basis Q and sort there
    t1 = (P[1]-C)/norm(P[1]-C)
    t2 = cross(n, t1)
    Q = [n t1 t2]
    sort!(P, lt=(A, B) -> begin
        A_proj = Q'*(A-C)
        B_proj = Q'*(B-C)
        a = atan2(A_proj[3], A_proj[2])
        b = atan2(B_proj[3], B_proj[2])
        return a > b
    end)
end

function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}; debug=true)

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)
    area = 0.0

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals = calculate_normals(slave_elements, time, Val{2};
                                rotate_normals=props.rotate_normals)
    update!(slave_elements, "normal", normals)

    # 2. loop all slave elements
    for slave_element in slave_elements

        slave_element_nodes = get_connectivity(slave_element)
        nsl = length(slave_element)
        X1 = slave_element("geometry", time)
        n1 = Field([normals[j] for j in slave_element_nodes])

        # project slave nodes to auxiliary plane (x0, Q)
        #xi = get_reference_element_midpoint(slave_element)
        xi = [1/3, 1/3]
        N = vec(get_basis(slave_element, xi, time))
        x0 = N*X1
        n0 = N*n1
        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]

        # 3. loop all master elements
        for master_element in slave_element("master elements", time)

            master_element_nodes = get_connectivity(master_element)
            nm = length(master_element)
            X2 = master_element("geometry", time)

            # 3.1 project master nodes to auxiliary plane and create polygon clipping
            M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
            P = get_polygon_clip(S, M, n0)
            length(P) < 3 && continue # no clipping or shared edge (no volume)
            check_orientation!(P, n0)
            C0 = calculate_centroid(P)

            De = zeros(nsl, nsl)
            Me = zeros(nsl, nm)
            ge = zeros(field_dim*nsl)

            # 4. loop integration cells
            for cell in get_cells(P, C0)
                virtual_element = Element(Tri3)
                update!(virtual_element, "geometry", cell)
                #x_cell = Field(cell)

                # 5. loop integration point of integration cell
                for ip in get_integration_points(virtual_element, 3)
                    N = vec(get_basis(virtual_element, ip, time))
                    #dN = vec(get_dbasis(virtual_element, ip, time))
                    #JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
                    #wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
                    detJ = virtual_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ

                    # project gauss point from auxiliary plane to master and slave element
                    #x_gauss = N*x_cell
                    x_gauss = virtual_element("geometry", ip, time)
                    if isnan(x_gauss[1])
                        info("is nan")
                        info("x_gauss = $x_gauss")
                        info("cell = $cell")
                        info("C0 = $C0")
                        info("P = $P")
                        info("S = $S")
                        info("M = $M")
                        info("n0 = $n0")
                        error("nan, unable to continue")
                    end
                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                    xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)

                    # add contributions
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    N2 = vec(get_basis(master_element, xi_m, time))
                    De += w*N1*N1'
                    Me += w*N1*N2'
                    if props.adjust
                        u1 = slave_element("displacement", time)
                        u2 = master_element("displacement", time)
                        x_s = N1*(X1+u1)
                        x_m = N2*(X2+u2)
                        ge += w*vec((x_m-x_s)*N1')
                    end
                    area += w
                end # integration points done

            end # integration cells done

            # 6. add contribution to contact virtual work
            sdofs = get_gdofs(problem, slave_element)
            mdofs = get_gdofs(problem, master_element)

            for i=1:field_dim
                lsdofs = sdofs[i:field_dim:end]
                lmdofs = mdofs[i:field_dim:end]
                add!(problem.assembly.C1, lsdofs, lsdofs, De)
                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
                add!(problem.assembly.C2, lsdofs, lsdofs, De)
                add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
            end
            add!(problem.assembly.g, sdofs, ge)

        end # master elements done

    end # slave elements done, contact virtual work ready

    debug && info("area of interface: $area")

end


text:
module TimeAxisArrays

using Reexport, Iterators
@reexport using AxisArrays

export TimeAxisArray, RegularTimeAxisArray,
timestamps, interval,
collapse, downsample, moving,
lead, lag, percentchange,
dropif, dropnan,
readtimeaxisarray, date, datetime

symbolize(x) = x |> string |> Symbol
defaultaxisnames = map(symbolize, ["Timestamp"; "Column"; map(i -> "Axis$i", 3:10)])
defaultaxisvalues = map(symbolize, 'A':'Z')

# Typealiases

typealias Timestamp Union{Number,Dates.AbstractTime}
typealias AbstractTimeVector{T<:Timestamp} AbstractVector{T}
typealias CategoryVector Vector{Symbol}

typealias CategoricalAxis{T} Axis{T, CategoryVector}
typealias TimeAxis{T<:AbstractTimeVector} Axis{defaultaxisnames[1], T}
typealias RegularTimeAxis{T<:Range} Axis{defaultaxisnames[1], T}

typealias TimeAxisArray{T,N,D,Ax<:Tuple{TimeAxis,Vararg{CategoricalAxis}}} AxisArray{T,N,D,Ax}
typealias RegularTimeAxisArray{T,N,D,Ax<:Tuple{RegularTimeAxis,Vararg{CategoricalAxis}}} AxisArray{T,N,D,Ax}

# Helper functions

maketimeaxis(timestamps::AbstractTimeVector) =
    Axis{defaultaxisnames[1]}(timestamps)

makecategoricalaxes(axesvals) = #TODO: Precise typing
    map((n,vs) -> Axis{n}(vs), defaultaxisnames[2:length(axesvals)+1], axesvals)

# Constructors

## 1-D TimeAxisArray

(::Type{TimeAxisArray})(data::AbstractArray, timestamps::AbstractTimeVector) =
    AxisArray(data, maketimeaxis(timestamps))

(::Type{RegularTimeAxisArray})(data::AbstractArray, timestamps::Range) =
    TimeAxisArray(data, timestamps)

## Custom higher-axis values with automatic names

(::Type{TimeAxisArray})(data::AbstractArray, timestamps::AbstractTimeVector, axesvals::CategoryVector...) =
    AxisArray(data, maketimeaxis(timestamps), makecategoricalaxes(axesvals)...)

(::Type{RegularTimeAxisArray})(data::AbstractArray, timestamps::Range, axesvals::CategoryVector...) =
    TimeAxisArray(data, timestamps, axesvals...)

## Custom higher-axis values with custom higher-axis names

(::Type{TimeAxisArray})(data::AbstractArray, timestamps::AbstractTimeVector, axes::CategoricalAxis...) =
    AxisArray(data, maketimeaxis(timestamps), axes...)

(::Type{RegularTimeAxisArray})(data::AbstractArray, timestamps::Range, axes::CategoricalAxis...) =
    TimeAxisArray(data, timestamps, axes...)

include("operations.jl")
include("io.jl")
include("utilities.jl")

end # module

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS118
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This problem has a quadratic objective with linear constraints.
#############################################################################

using JuMP
using Base.Test

let

m = Model()

L = zeros(15)
L[1] =  8.0
L[2] = 43.0
L[3] =  3.0

U = zeros(15)
U[1] = 21.0
U[2] = 57.0
U[3] = 16.0
for k in 1:4
    U[3*k+1] =  90.0
    U[3*k+2] = 120.0
    U[3*k+3] =  60.0
end

@variable(m, L[i] <= x[i=1:15] <= U[i])

@NLobjective(m, Min,
    sum{2.3     * x[3*k+1]   +
        0.0001  * x[3*k+1]^2 +
        1.7     * x[3*k+2]   +
        0.0001  * x[3*k+2]^2 +
        2.2     * x[3*k+3] +
        0.00015 * x[3*k+3]^2, k=0:4})

# setobjective(m, :Min,
#   sum([(2.3     * x[3*k+1]   +
#         0.0001  * x[3*k+1]^2 +
#         1.7     * x[3*k+2]   +
#         0.0001  * x[3*k+2]^2 +
#         2.2     * x[3*k+3]^2 +
#         0.00015 * x[3*k+3]^2) for k in 0:4 ]))

# constr1
for j in 1:4
    @constraint(m, x[3*j+1] - x[3*j-2] + 7 <= 13)
    @constraint(m, x[3*j+1] - x[3*j-2] + 7 >=  0)
end

# constr2
for j in 1:4
    @constraint(m, x[3*j+2] - x[3*j-1] + 7 <= 14)
    @constraint(m, x[3*j+2] - x[3*j-1] + 7 >=  0)
end

# constr3
for j in 1:4
    @constraint(m, x[3*j+3] - x[3*j  ] + 7 <= 13)
    @constraint(m, x[3*j+3] - x[3*j  ] + 7 >=  0)
end

@constraint(m, x[1] + x[2] + x[3]    >= 60)
@constraint(m, x[4] + x[5] + x[6]    >= 50)
@constraint(m, x[7] + x[8] + x[9]    >= 70)
@constraint(m, x[10] + x[11] + x[12] >= 85)
@constraint(m, x[13] + x[14] + x[15] >= 100)

# Initial solution (could also use 'start' keyword in @variable)
setvalue(x[1], 20.0)
setvalue(x[2], 55.0)
setvalue(x[3], 15.0)
setvalue(x[4], 20.0)
setvalue(x[5], 60.0)
setvalue(x[6], 20.0)
setvalue(x[7], 20.0)
setvalue(x[8], 60.0)
setvalue(x[9], 20.0)
setvalue(x[10], 20.0)
setvalue(x[11], 60.0)
setvalue(x[12], 20.0)
setvalue(x[13], 20.0)
setvalue(x[14], 60.0)
setvalue(x[15], 20.0)

solve(m)

#println(getvalue(x))

@test_approx_eq_eps getvalue(x[1]) 8.0  1e-5
@test_approx_eq_eps getvalue(x[2]) 49.0 1e-5
@test_approx_eq_eps getvalue(x[3]) 3.0  1e-5
@test_approx_eq_eps getvalue(x[4]) 1.0  1e-5
@test_approx_eq_eps getobjectivevalue(m) 664.82045 1e-5

end


text:
module TestQuasinormalModes

include("../src/CustomTypes.jl")
using CustomTypes

const tolerance = 1e-4 ## tolerance we compare to

##============================================================
import QuasinormalModes as QNM

nr = 48
nl = 16

avals = [0.0, 0.354, 0.7, 0.99] ## same as in generate.py

for s=[-2,-1,0]
   for m=[-2,0,3]
      lmin = max(abs(s),abs(m))
      for l=[lmin, lmin+1, lmin+4]
         for n=[0,1,2]
            for a in avals
               om, la qnm(n,s,m,l,a)
               println("s=$s, (n=$n,l=$l,m=$m), a=$a, ω=$om, Λ=$la")

               f_om, f_la, f_vs, f_vr, rvals = QNM.compute_om(
                                         myI(nr),
                                         myI(nl),
                                         myI(s),
                                         myI(l),
                                         myI(m),
                                         myF(a),
                                         myC(om+0.1*rand(myC)),
                                         verbose=false
                                        )
               @test abs(f_om-om)/max(1,abs(om)) < tolerance
            end
         end
      end
   end
end

end


text:
@testset "Find" begin
    seq = dna"ACGNA"
    @test findnext(DNA_A, seq, 1) == 1
    @test findnext(DNA_C, seq, 1) == 2
    @test findnext(DNA_G, seq, 1) == 3
    @test findnext(DNA_N, seq, 1) == 4
    @test findnext(DNA_T, seq, 1) == nothing
    @test findnext(DNA_A, seq, 2) == 5

    @test_throws BoundsError findnext(DNA_A, seq, 0)
    @test_throws BoundsError findnext(DNA_A, seq, 6)

    @test findprev(DNA_A, seq, 4) == 1
    @test findprev(DNA_C, seq, 4) == 2
    @test findprev(DNA_G, seq, 4) == 3
    @test findprev(DNA_N, seq, 4) == 4
    @test findprev(DNA_T, seq, 4) == nothing
    @test findprev(DNA_G, seq, 2) == nothing

    @test_throws BoundsError findprev(DNA_A, seq, 0)
    @test_throws BoundsError findprev(DNA_A, seq, 6)

    seq = dna"ACGNAN"
    @test findfirst(DNA_A, seq) == 1
    @test findfirst(DNA_N, seq) == 4
    @test findfirst(DNA_T, seq) == nothing

    @test findlast(DNA_A, seq) == 5
    @test findlast(DNA_N, seq) == 6
    @test findlast(DNA_T, seq) == nothing
end


text:
struct LatentVariable
    user::Array{Float64,2}
    item::Array{Float64,2}
end

function svd(dataset::Persa.CFDatasetAbstract, features::Int)::LatentVariable
    matrix = Persa.getmatrix(dataset)
    (result, -) = svds(matrix, nsv = features)
    U = result.U
    V = result.Vt
    return LatentVariable(U, V')
end

Base.length(lv::LatentVariable) = maximum(size(lv))

user(lv::LatentVariable) = lv.user
item(lv::LatentVariable) = lv.item

user(lv::LatentVariable, id::Int) = user(lv)[id,:]
item(lv::LatentVariable, id::Int) = item(lv)[id,:]

Base.size(lv::LatentVariable) = (size(lv.user)[2], size(lv.item)[2])


text:
export Epoch, start!, step!, release!

type Epoch{A <: OptAlgo}
    algo::A
    itertaionNo::Int
    started::Bool
end

Epoch{A}(algo::A) = Epoch{A}(algo, 0, false)

Epoch{A}(::Type{A}, pars...) = Epoch{A}(A(pars...), 0, false)

@generated function release!{A}(epoch::Epoch{A})
    if length(methods(release!, (A, ))) > 0
        :( release!(epoch.algo) )
    else
        :( )
    end
end

function start!(epoch::Epoch)
    start!(epoch.algo)
    epoch.itertaionNo = 0
    epoch.started = true
end

function step!(epoch::Epoch)
    !epoch.started && error("Epoch is not started.")
    step!(epoch.algo)
    epoch.itertaionNo += 1
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
## Types

# Projection context.  TODO: Will this be exposed?
#type Context
#    rep::Ptr{Void} # Pointer to internal projCtx struct
#end

"""
Cartographic projection type
"""
type Projection
    #ctx::Context   # Projection context object
    rep::Ptr{Void} # Pointer to internal projPJ struct

    # [geod]: a structure containing the parameters of the spheroid
    # some of the fields in [geod] are mildly duplicative of the information
    # available in [rep], which can be exposed only through _get_spheroid_defn
    geod::_geodesic
end

function Projection(proj_ptr::Ptr{Void})
    proj = Projection(proj_ptr, null_geodesic())
    finalizer(proj, freeProjection)
    proj
end

"""
Construct a projection from a string in proj.4 "plus format"

The projection string `proj_string` is defined in the proj.4 format,
with each part of the projection specification prefixed with '+' character.
For example:

    `wgs84 = Projection("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")`
"""
Projection(proj_string::String) = Projection(_init_plus(proj_string))

function freeProjection(proj::Projection)
    _free(proj.rep)
    proj.rep = C_NULL
end

# Pretty printing
Base.print(io::IO, proj::Projection) = print(io, strip(_get_def(proj.rep)))
Base.show(io::IO, proj::Projection) = print(io, "Projection(\"$proj\")")


text:
using Flux, HDF5, Random

include("BayesianSZNet.jl")
using .BayesianSZNet

PS = ["05", "01", "005", "001"]
N_SAMPLES = 256
N_VAL = 50000
N_MODELS = 5

Random.seed!(72)
h5open("data/dr12q_superset.hdf5", "r+") do datafile
    X_va = gpu(read(datafile, "X_va"))
    ẑs = zeros(N_SAMPLES, N_VAL, N_MODELS)

    fcnn_group = create_group(datafile, "fcnn")
    sznet_group = create_group(datafile, "sznet")
    for p in PS
        for i in 1:N_MODELS
            ẑs[:, :, i] = sample(FCNN("models/fcnn_$(i)_0.$(p)_1.0e-5.bson"), X_va)
        end
        write_dataset(fcnn_group, "zs_pred_$(p)", ẑs)

        for i in 1:N_MODELS
            ẑs[:, :, i] = sample(SZNet("models/sznet_$(i)_0.$(p)_1.0e-9.bson"), X_va)
        end
        write_dataset(sznet_group, "zs_pred_$(p)", ẑs)
    end
end

max_stars_repo_name: podondra/bdl-redshift

max_issues_repo_name: podondra/bdl-redshift

max_forks_repo_name: podondra/bdl-redshift


text:
# This file includes code that was formerly a part of Julia.
# License is MIT: LICENSE.md

using MurmurHash3
using MurmurHash3: mmhash128_4, mmhash128_8_c, mmhash128_8_a, mmhash128_8_u

using Test

p1 = SubString("--hello--",3,7)
p2 = "hello"

const MaybeSub = Union{String,SubString{String}}
const mhseed32 = 0x56419c81
const mhseed64 = 0x71e729fd56419c81

# 32-bit MurmurHash (calling C version in Base)
mh32c(siz::Int, ptr::Ptr{UInt8}, seed=0%UInt32) =
    ccall(:memhash32_seed, UInt32, (Ptr{UInt8}, Csize_t, UInt32), ptr, siz, seed)
mh32c(str::MaybeSub, seed=0%UInt32) = mh32c(sizeof(str), pointer(str), seed)

# Calling Julia version in MurmurHash.jl
mh32j(str::MaybeSub, seed=0%UInt32) = mmhash32(sizeof(str), pointer(str), seed)
mh32j(len::Integer, val::Unsigned, seed=0%UInt32) = mmhash32(len, val, seed)

# 128-bit MurmurHash (either 32-bit or 64-bit implementation, C version in Base), lower 64-bits
mh128c(siz::Int, ptr::Ptr{UInt8}, seed=0) =
    ccall(:memhash_seed, UInt64, (Ptr{UInt8}, Csize_t, UInt32), ptr, siz, seed)
mh128c(s::MaybeSub, h=0) = mh128c(sizeof(s), pointer(s), h%UInt32)

mh128j4(s::MaybeSub, h=0) = mmhash128_4(sizeof(s), pointer(s), h%UInt32)
mh128j8a(s::MaybeSub, h=0) = mmhash128_8_a(sizeof(s), pointer(s), h%UInt32)
mh128j8u(s::MaybeSub, h=0) = mmhash128_8_u(sizeof(s), pointer(s), h%UInt32)
mh128j8c(s::MaybeSub, h=0) = mmhash128_8_c(s, h%UInt32)

# 128-bit MurmurHash (Julia version)
mh128j(str::MaybeSub, seed=0) = mmhash128_a(sizeof(str), pointer(str), seed%UInt32)
mh128j(len::Int, val::Unsigned, seed=0) = mmhash128_a(len, val, seed%UInt32)
mh128j_c(str::AbstractString, seed=0) = mmhash128_c(str, seed%UInt32)

load_u64(p) = unsafe_load(reinterpret(Ptr{UInt64}, pointer(p1)))

const sizp2 = sizeof(p2)
const unsp2 = load_u64(p2)

pd(msg, v) = (print(stderr, msg); dump(stderr, v))

pd("32c:   ", mh32c(p2))
pd("32j:   ", mh32j(p2))
pd("32u:   ", mh32j(sizp2, unsp2))
pd("128c:  ", mh128c(p2))
pd("128j:  ", mh128j(p2))
pd("128u:  ", mh128j(sizp2, unsp2))
pd("128j4: ", mh128j4(p2))
pd("128ja: ", mh128j8a(p2))
pd("128ju: ", mh128j8u(p2))
pd("128jc: ", mh128j8c(p2))

@testset "MurmurHash3" begin
    @testset "32-bit MurmurHash" begin
        @test mh32j(p1) == mh32c(p1)
        @test mh32j(p2) == mh32c(p2)
        @test mh32j(sizp2, unsp2) == mh32c(p2)
    end
    @static if UInt === UInt64
    @testset "Aligned vs unaligned" begin
        @test mh128j_c(p1) == mh128j(p2)
    end
    @testset "128-bit MurmurHash" begin
        @test last(mh128j(sizp2, unsp2)) == mh128c(p1)
        @static if UInt === UInt64
            @test last(mh128j(p1)) == mh128c(p1)
            @test last(mh128j_c(p2)) == mh128c(p2)
            @test last(mh128j_c(p2)) == mh128c(p1)
        else
            @test_broken last(mh128j(p1)) == mh128c(p1)
            @test_broken last(mh128j_c(p2)) == mh128c(p2)
            @test_broken last(mh128j_c(p2)) == mh128c(p1)
        end
    end
    end
end


text:
"""
    cheb(xmin, xmax, N)

Return a Chebyshev-Lobatto grid, together with its first and second
differentiation matrices.

# Arguments
* `xmin::Real`: rightmost grid point
* `xmax::Real`: leftmost grid point
* `N::Integer`: total number of grid points
"""
function cheb(xmin::T, xmax::T, N::Integer) where {T<:Real}
    x, D, D2 = cheb(N)
    x = 0.5 * (xmax + xmin .+ (xmax - xmin) * x)
    D  ./= 0.5 * (xmax - xmin)
    D2 ./= 0.25 * (xmax - xmin)^2
    x, D, D2
end

# taken from Trefethen (2000), "Spectral Methods in MatLab"

"When omitting grid limits, default to [-1,1] interval"
function cheb(N::Integer)
    @assert(N > 1, "number of points should be greater than 1...")
    M = N - 1
    x = -cos.(pi*(0:M)/M)

    c  = [2; ones(M-1, 1); 2] .* (-1).^(0:M)
    X  = repeat(x, 1, N)
    dX = X - X'

    D = (c * (1 ./ c)') ./ (dX + Matrix(I, N, N))  # off-diagonal entries
    D = D - diagm(0 => sum(D', dims=1)[:])         # diagonal entries

    x, D, D*D
end

"""
    fourier(xsize, N)

Return the first and second derivative matrices on a Fourier grid

# Arguments
* `xsize::Real`: size of (periodic) domain
* `N::Integer`: total number of grid points
"""
function fourier(xsize::T, N::Integer) where {T<:Real}
    D, D2 = fourier(N)
    pi_ = T(pi)
    D  ./= xsize / (2*pi_)
    D2 ./= xsize^2 / (4*pi_^2)
    D, D2
end

# taken from Trefethen (2000), "Spectral Methods in MatLab"

"When omitting `xsize`, default to 2π"
function fourier(N::Integer)
    @assert(mod(N,2)==0, "number of points needs to be even")
    h = 2*pi/N

    column = [0; 0.5*(-1).^(1:N-1) .* cot.((1:N-1)*h/2)]
    tmp = [circshift(column, i) for i in 0:N-1]
    D   = hcat(tmp...)

    column2 = [-pi^2/(3*h^2) - 1/6; -0.5*(-1).^(1:N-1) ./ sin.((1:N-1)*h/2).^2]
    tmp = [circshift(column2, i) for i in 0:N-1]
    D2  = hcat(tmp...)

    D, D2
end

struct ChebInterpolator{T,A,FT<:FFTW.r2rFFTWPlan}
    xmin     :: T
    xmax     :: T
    c        :: A
    fft_plan :: FT
end
"""
    ChebInterpolator(xmin, xmax, N)

Build an interpolator to act on functions defined on a Chebyshev-Lobatto grid.
Uses `FFTW`.

# Arguments
* `xmin::Real`: rightmost grid point
* `xmax::Real`: leftmost grid point
* `N::Integer`: total number of grid points
"""
function ChebInterpolator(xmin::T, xmax::T, N::Int) where {T<:Real}
    M  = N - 1
    x  = Vector{T}(undef, N)

    # use only one thread for the FFTW! by default it allocates a bunch of them,
    # and it hurts performance (by a huge amount!) when using this inside loops.
    # it's much better to thread the loops themselves, which is what we do.
    FFTW.set_num_threads(1)

    # Create the FFT plan for the DCT-I
    fft_plan = FFTW.plan_r2r(x, FFTW.REDFT00)

    c  = M * [2; ones(Int,M-1); 2] .* (-1).^(0:M)

    ChebInterpolator{T,typeof(c),typeof(fft_plan)}(xmin, xmax, c, fft_plan)
end
"""
    ChebInterpolator(xp::Vector)

Build it directly from a `Vector` with the grid points
"""
ChebInterpolator(xp::Vector) = ChebInterpolator(xp[1], xp[end], length(xp))

"""
# Examples

```
julia> xx, = Jecco.cheb(0.0, 2.0, 16);

julia> f = xx.^2;

julia> interp = Jecco.ChebInterpolator(xx);

julia> f_interp = interp(f);

julia> f_interp(0.2)
0.03999999999999987
```
"""
function (interp::ChebInterpolator)(fp)
    # compute the DCT-I of the coefficients
    fft_fp = interp.fft_plan * fp

    # compute the spectral coefficients
    spec_coeff = fft_fp ./ interp.c

    function (x0::T) where {T<:Real}
        @assert interp.xmin <= x0 <= interp.xmax
        X = (2 * x0 - (interp.xmin + interp.xmax)) / (interp.xmax - interp.xmin)
        sum_l = zero(T)
        @fastmath @inbounds for i in LinearIndices(spec_coeff)
            sum_l += spec_coeff[i] * cos( (i-1)*acos(X) )
        end
        sum_l
    end
end


text:
# https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, -5 ≤ x[1:2] ≤ 10)
@disjunction(
    m,
    begin
        con1[i=1:2], 0 ≤ x[i] ≤ [3,4][i]
    end,
    begin
        con2[i=1:2], [5,4][i] ≤ x[i] ≤ [9,6][i]
    end,
    reformulation = :big_m,
    name = :y
)
print(m)

# ┌ Warning: [con1[1] : x[1] in [0.0, 3.0], con1[2] : x[2] in [0.0, 4.0]] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
# ┌ Warning: [con2[1] : x[1] in [5.0, 9.0], con2[2] : x[2] in [4.0, 6.0]] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
# Feasibility
# Subject to
#  XOR(disj_y) : y[1] + y[2] == 1.0         <- XOR constraint
#  con1[1,lb] : -x[1] + 5 y[1] <= 5.0      <- left-side of con1[1]
#  con1[1,ub] : x[1] + 7 y[1] <= 10.0      <- right-side of con1[1]
#  con1[2,lb] : -x[2] + 5 y[1] <= 5.0      <- left-side of con1[2]
#  con1[2,ub] : x[2] + 6 y[1] <= 10.0      <- right-side of con1[2]
#  con2[1,lb] : -x[1] + 10 y[2] <= 5.0     <- left-side of con2[1]
#  con2[1,ub] : x[1] + y[2] <= 10.0        <- right-side of con2[1]
#  con2[2,lb] : -x[2] + 9 y[2] <= 5.0      <- left-side of con2[2]
#  con2[2,ub] : x[2] + 4 y[2] <= 10.0      <- right-side of con2[2]
#  x[1] >= -5.0                             <- varaible bounds
#  x[2] >= -5.0                             <- variable bounds
#  x[1] <= 10.0                             <- variable bounds
#  x[2] <= 10.0                             <- variable bounds
#  y[1] >= 0.0                              <- lower bound on binary
#  y[2] >= 0.0                              <- lower bound on binary
#  y[1] <= 1.0                              <- upper bound on binary
#  y[2] <= 1.0                              <- upper bound on binary
#  y[1] binary                              <- indicator variable (1st disjunct) is binary
#  y[2] binary                              <- indicator variable (2nd disjunct) is binary


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# benchmark/benchmarks.jl
# Runs benchmarks on all supported McCormick operators.
# Allocations are expected for the following operations:
# - MC(Float64, Interval, 1): Nonstatic definition of static type.
# - exp2, exp10, tanh, atanh, acosh: Underlying interval operation may allocate.
#############################################################################

using BenchmarkTools, IntervalArithmetic

# don't really want to measure interval performance, so we're using the
# fastest interval rounding mode
setrounding(Interval, :none)
using McCormick

const SUITE = BenchmarkGroup()

function test_comp1(f, x, n)
    z = f(x)
    for i = 1:n
        z = f(z)
    end
    z
end

function test_comp2l(f, x, y, n)
    z = f(x, y)
    for i = 1:n
        z = f(z, y)
    end
    z
end

function test_comp2r(f, x, y, n)
    z = f(x, y)
    for i = 1:n
        z = f(x, z)
    end
    z
end

for T in (NS, ) #(NS, Diff, MV)
        begin
            S = SUITE["Constructors $(T)"] = BenchmarkGroup()
            S["MC(Float64)"] = @benchmarkable MC{5,$T}(2.1)
            S["MC(Float64,Float64)"] = @benchmarkable MC{5,$T}(0.1,2.0)
            S["MC(Interval{Float64})"] = @benchmarkable MC{5,$T}($(Interval{Float64}(0.1,2.0)))
            S["MC(Float64, Interval{Float64}, Int)"] = @benchmarkable MC{5,$T}(1.0,$(Interval{Float64}(0.1,2.0)),2)
            S = SUITE["OOP $(T) McCormick"] = BenchmarkGroup()
        end
    begin
        for op in (+, -, exp, expm1, log, log2, log10, log1p,
                   abs, one, zero, real, max, min, inv, cosh, sqrt,
                   isone, isnan)
            S[string(op)*"(x)"] = @benchmarkable $(test_comp1)($op, a, 1000) setup = (a = MC{5,$T}(1.0,$(Interval{Float64}(0.1,2.0)),2))
        end
        # domain violations asin, asind, acos, acosd, atan, atand, tan, tand, sec, secd,
        for op in (sin, sind, cos, cosd, step, sign, sinh)
            S[string(op)*"(x), x > 0"] = @benchmarkable $(test_comp1)($op, a, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
            S[string(op)*"(x), 0 ∈ x"] = @benchmarkable $(test_comp1)($op, a, 1000) setup = (a = MC{5,$T}(0.1,$(Interval{Float64}(-0.4,0.3)),2))
            S[string(op)*"(x), x < 0"] = @benchmarkable $(test_comp1)($op, a, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        end
        for op in (min, max, *, -, +, /)
            S[string(op)*"(x, Float64)"] =  @benchmarkable $(test_comp2l)($op, x, q, 10000) setup = (x = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2); q = 1.34534)
            S[string(op)*"(Float64, x)"] =  @benchmarkable $(test_comp2r)($op, q, x, 10000) setup = (x = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2); q = 1.34534)
            S[string(op)*"(x, y)"] =   @benchmarkable $(test_comp2l)($op, x, y, 1000) setup = (x = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2);
                                                                            y = MC{5,$T}(0.5,$(Interval{Float64}(0.3,0.9)),2);)
        end
        S["x^2, x > 0"] = @benchmarkable test_comp2l($^, a, 2, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^2, 0 ∈ x"] = @benchmarkable test_comp2l($^, a, 2, 1000) setup = (a = MC{5,$T}(0.1,$(Interval{Float64}(-0.4,0.3)),2))
        S["x^2, x < 0"] = @benchmarkable test_comp2l($^, a, 2, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        S["x^3, x > 0"] = @benchmarkable test_comp2l($^, a, 3, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^3, 0 ∈ x"] = @benchmarkable test_comp2l($^, a, 3, 1000) setup = (a = MC{5,$T}(0.1,$(Interval{Float64}(-0.4,0.3)),2))
        S["x^3, x < 0"] = @benchmarkable test_comp2l($^, a, 3, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        S["x^4, x > 0"] = @benchmarkable test_comp2l($^, a, 4, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^4, 0 ∈ x"] = @benchmarkable test_comp2l($^, a, 4, 1000) setup = (a = MC{5,$T}(0.1,$(Interval{Float64}(-0.4,0.3)),2))
        S["x^4, x < 0"] = @benchmarkable test_comp2l($^, a, 4, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        S["x^-1, x > 0"] = @benchmarkable test_comp2l($^, a, -1, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^-1, x < 0"] = @benchmarkable test_comp2l($^, a, -1, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        S["x^-2, x > 0"] = @benchmarkable test_comp2l($^, a, -2, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^-2, x < 0"] = @benchmarkable test_comp2l($^, a, -2, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
        S["x^-3, x > 0"] = @benchmarkable test_comp2l($^, a, -3, 1000) setup = (a = MC{5,$T}(0.4,$(Interval{Float64}(0.1,0.9)),2))
        S["x^-3, x < 0"] = @benchmarkable test_comp2l($^, a, -3, 1000) setup = (a = MC{5,$T}(-0.5,$(Interval{Float64}(-0.9,-0.1)),2))
    end
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


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

JLLWrappers.@generate_main_file_header("SISL")
JLLWrappers.@generate_main_file("SISL", UUID("d647f7b9-f016-52c0-8515-7d1c032fadec"))
end  # module SISL_jll


text:
"""
Make a symbol with extra bits easily:

   S"dsaf asdfl" => Symbol("dsaf asdfl")
"""
macro S_str(string) Base.Meta.quot(Symbol(string)) end

"""
    @get use_forward_model = opts | 123

is sugar for

    use_forward_model = get(opts, :use_forward_model, 123)
"""
macro get(ex)
    @assert ex.head==:(=)
    @assert ex.args[2].head==:call
    @assert ex.args[2].args[1]==:|
    out = ex.args[1]
    var, def = ex.args[2].args[2:3]
    esc(:($out = get($var, $(QuoteNode(out)), $def)))
end

"""
    @get! use_forward_model = opts | false

is sugar for

    use_forward_model = get!(opts, :use_forward_model, false)
"""
macro get!(ex)
    @assert ex.head==:(=)
    @assert ex.args[2].head==:call
    @assert ex.args[2].args[1]==:|
    out = ex.args[1]
    var, def = ex.args[2].args[2:3]
    esc(:($out = get!($var, $(QuoteNode(out)), $def)))
end

"""
    volume_area_mean_h(gl::Glacier)
    volume_area_mean_h(area::Number; isicecap::Bool=false)
    volume_area_mean_h(area::Number,c::Number,gamma::Number)

Volume Area scaling formula.  Returns mean thickness in meters.
"""
function volume_area_mean_h(area::Number; isicecap::Bool=false)
    if isicecap
        c = 0.034
        gamma = 1.36
    else
        c = 0.054
        gamma = 1.25
    end
    volume_area_mean_h(area,c,gamma)
end
volume_area_mean_h(area::Number,c::Number,gamma::Number) = c*(area/1e6)^(gamma-1) * 1e3


text:
include(joinpath(dirname(@__FILE__), "MyPetStore", "src", "MyPetStore.jl"))

include("test_UserApi.jl")
include("test_StoreApi.jl")
include("test_PetApi.jl")

const server = "http://petstore.swagger.io/v2"
TestUserApi.test(server)
TestStoreApi.test(server)
TestPetApi.test(server)


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ f497faa1-3b8e-4891-b0e5-679fd29500d7
using Pkg

# ╔═╡ 91ab0997-63de-4ed0-af8d-2221568bef7c
Pkg.activate("..")

# ╔═╡ 0f496f40-14bc-45f7-a890-dbe531e2500d
using NLS_Fit, Plots, DelimitedFiles

# ╔═╡ c703884c-de01-4539-b8a2-9fb133d0cad6
md"# Environment & packages"

# ╔═╡ c798b6b4-dd3d-4c33-ab73-8e95eac97613
md"# Raw data plot"

# ╔═╡ 2d2d35af-2e14-495f-a6c6-dc348a60fc2b
begin
	XY = readdlm("../data/recalibration.txt")
	X = XY[:,1]
	Y = XY[:,2]
	plot(X,Y, seriestype = :scatter, label = "raw data", title = "Recalibration")
end

# ╔═╡ 4dcbd1a8-b79a-43a1-9a8a-6fc5c373dc73
md"# Create model"

# ╔═╡ 29ccb74f-8a62-478d-a849-5610cb2e62c5
md"## Uncalibrable model"

# ╔═╡ 44a9850b-4ecc-42c6-8270-2ffad6903b09
begin
	model = Gaussian_Peak() + Gaussian_Peak() + Gaussian_Peak()

	θ1 = Float64[1,5,1]
	θ2 = Float64[2,10,1]
	θ3 = Float64[1,20,2]
	θ_init_model = vcat(θ1,θ2,θ3)
end

# ╔═╡ 21353220-0e6b-40e4-8930-4889f9bb757a
md"## Calibrable model

Complete the previous model with a parameterized transformation used as calibration map
"

# ╔═╡ 5ceea856-87f0-4874-860d-57cf1f31c5e8
begin
	recal_map = Map_Affine(X[1]=>X[1],X[end]=>X[end])
	recal_model = Model2Fit_Recalibration(model,recal_map)

	θc = Float64[1,1]
	θ_init_recal_model = vcat(θ_init_model, θc)

	Y_init = eval_y(recal_model,X,θ_init_recal_model)

	plot!(X,Y_init, label = "model θ_init")
end

# ╔═╡ fed24429-8e20-4601-8bf8-5a4b401b2534
 NLS_Fit.eval_calibrated_x(recal_model,X,θ_init_recal_model)

# ╔═╡ f331ea35-4d56-4ced-a76e-de621de25593
md"Define parameter bounds"

# ╔═╡ fcc0e63e-d873-4af8-bfe2-55b41e2b9f1a
begin
	ε = eps(Float64)
	lower_bound = Float64[0,5,ε,0,10,ε,0,20,ε,0.5,0.0]
	upper_bound = Float64[+Inf,5,2.5,+Inf,10,2.5,+Inf,20,2.5,1.5,2.0]
	bc = NLS_Solver.BoundConstraints(lower_bound,upper_bound)
end

# ╔═╡ e61704f1-5eae-4a35-9d5f-0f734e46ef03
md"# Solve the nonlinear least squares problem"

# ╔═╡ e33fe5dd-db88-4de9-8816-b4c36d6765d1
begin
	nls = NLS_ForwardDiff_From_Model2Fit(recal_model,X,Y)
	conf = NLS_Solver.LevenbergMarquardt_BC_Conf()
	result = NLS_Solver.solve(nls,θ_init_recal_model,bc,conf)
end

# ╔═╡ b2e98c3f-5de5-4c23-bbb7-80d39cf219fd
md"Converged?"

# ╔═╡ fb6f9b74-1687-4846-b070-58581a1b234b
NLS_Solver.converged(result)

# ╔═╡ 4bb856a1-5696-42c8-8e63-fd7abe0eb276
md"Objective function value"

# ╔═╡ 9fa86072-9483-4b63-8dc9-e94469c1a8e9
NLS_Solver.objective_value(result)

# ╔═╡ 3367e5bd-3159-4327-90a4-01b7799e7064
solution(result)

# ╔═╡ c30df4ad-de11-446c-84ed-33755a78e5aa
md"# Plot calibrated model"

# ╔═╡ e7284fe8-79b8-4684-922c-f004cb98333e
begin
	θ_fit_recal_model = NLS_Solver.solution(result)
	Y_fit_recal_model = eval_y(recal_model,X,θ_fit_recal_model)
	plot!(X,Y_fit_recal_model, label = "model θ_fit")
end

# ╔═╡ Cell order:
# ╟─c703884c-de01-4539-b8a2-9fb133d0cad6
# ╠═f497faa1-3b8e-4891-b0e5-679fd29500d7
# ╠═91ab0997-63de-4ed0-af8d-2221568bef7c
# ╠═0f496f40-14bc-45f7-a890-dbe531e2500d
# ╟─c798b6b4-dd3d-4c33-ab73-8e95eac97613
# ╠═2d2d35af-2e14-495f-a6c6-dc348a60fc2b
# ╠═4dcbd1a8-b79a-43a1-9a8a-6fc5c373dc73
# ╠═29ccb74f-8a62-478d-a849-5610cb2e62c5
# ╠═44a9850b-4ecc-42c6-8270-2ffad6903b09
# ╟─21353220-0e6b-40e4-8930-4889f9bb757a
# ╠═5ceea856-87f0-4874-860d-57cf1f31c5e8
# ╠═fed24429-8e20-4601-8bf8-5a4b401b2534
# ╠═f331ea35-4d56-4ced-a76e-de621de25593
# ╠═fcc0e63e-d873-4af8-bfe2-55b41e2b9f1a
# ╠═e61704f1-5eae-4a35-9d5f-0f734e46ef03
# ╠═e33fe5dd-db88-4de9-8816-b4c36d6765d1
# ╟─b2e98c3f-5de5-4c23-bbb7-80d39cf219fd
# ╠═fb6f9b74-1687-4846-b070-58581a1b234b
# ╟─4bb856a1-5696-42c8-8e63-fd7abe0eb276
# ╠═9fa86072-9483-4b63-8dc9-e94469c1a8e9
# ╠═3367e5bd-3159-4327-90a4-01b7799e7064
# ╠═c30df4ad-de11-446c-84ed-33755a78e5aa
# ╠═e7284fe8-79b8-4684-922c-f004cb98333e


text:
using CollegeStratBase, Test

function present_value_test()
    x = 3.9;
    T = 7;
    R = 1.2;
    pvFactor = pv_factor(R, T);
    @test pvFactor > 1.0
    @test pvFactor < T

    xV = fill(x, T);
    pValue = present_value(xV, R);
    @test isapprox(x * pvFactor, pValue)

    xV = LinRange(-0.5, 3.2, T);
    pValue = present_value(xV, R);
    pv = sum(xV .* ((1/R) .^ (0 : (T-1))));
    @test isapprox(pv, pValue)
end

@testset "Helpers" begin
    present_value_test();
end

# -----------


text:
module TetMesh

using StaticArrays: MVector, @MVector, SVector
using WriteVTK: vtk_grid, MeshCell, VTKCellTypes
using Lazy: @as
using Base: setindex

struct NodeInfo
    coords::MVector{2, Float64}
    adjacent_nodes::Vector{Int}
end

struct ElementInfo{NodesPerFace, InternalNodes}
    control_nodes::MVector{3, Int}
    faces::MVector{3, Int}
    face_nodes::MVector{3, MVector{NodesPerFace, Int}}
    internal_nodes::MVector{InternalNodes, Int}
    adjacent_elements::Vector{Int}
end

struct FaceDetails{NodesPerFace}
    number::Int
    element::Int
    nodes::MVector{NodesPerFace, Int}
end

struct BoundaryRepresentation{NodesPerFace}
    nodes::Vector{Int}
    faces::Vector{FaceDetails{NodesPerFace}}
end

struct Mesh{NodesPerFace, InternalNodes, BoundaryNodesPerFace}
    nodes::Vector{NodeInfo}
    elements::Vector{ElementInfo{NodesPerFace, InternalNodes}}
    boundaries::Vector{BoundaryRepresentation{BoundaryNodesPerFace}}
    boundary_tags::Vector{Vector{String}}
end

function deserialize(stream::IO)
    max_element_adjacencies = read(stream, Int)
    max_node_adjacencies = read(stream, Int)
    NodesPerFace = read(stream, Int)
    InternalNodes = read(stream, Int)

    # Read nodes
    count = read(stream, Int)
    nodes = [deserialize(stream, NodeInfo) for i in 1:count]

    # Read elements
    count = read(stream, Int)
    elements = [
        deserialize(stream, ElementInfo{NodesPerFace, InternalNodes})
        for i in 1:count
    ]

    # Read boundaries
    count = read(stream, Int)
    boundaries = [
        deserialize(stream, BoundaryRepresentation{NodesPerFace+2})
        for i in 1:count
    ]

    # Read tags
    count = read(stream, Int)
    tags = [deserialize(stream, Vector{String}) for i in 1:count]

    Mesh{NodesPerFace, InternalNodes, NodesPerFace+2}(nodes, elements, boundaries, tags)
end

"""
Returns three items:
  - points: An array of the coordinates of vertices of triangles in the mesh;
  - elements: An array of WriteVTK.MeshCell with indices of the vertices of
    triangles in `points`;
  - node_map: `node_map[i]` maps a node `i` in the triangle mesh to a node index
    in the full (possibly higher-order) mesh.
"""
function points_and_cells(msh::Mesh)
    node_indices = zeros(Int, length(msh.nodes))
    coords = Vector{SVector{2, Float64}}(undef, length(msh.nodes))
    vertices = NTuple{3, Int}[]
    node_map = Int[]

    node_index = 1
    for element in msh.elements
        verts = (0, 0, 0)
        i = 1
        for node in (element.control_nodes .+ 1)
            if node_indices[node] == 0
                node_indices[node] = node_index
                coords[node_index] = msh.nodes[node].coords
                node_index += 1
                push!(node_map, node)
            end
            verts = setindex(verts, node_indices[node], i)
            i += 1
        end
        push!(vertices, verts)
    end

    nnodes = node_index - 1
    resize!(coords, nnodes)

    (@as v coords reinterpret(Float64, v) reshape(v, 2, :)),
    [MeshCell(VTKCellTypes.VTK_TRIANGLE, verts) for verts in vertices],
    node_map
end

function initialize_vtk(target, msh)
    points, cells, node_map = points_and_cells(msh)
    vtk_grid(target, points, cells), node_map
end

function add_scalar_data(target, msh::Mesh, node_map::Vector{Int}, data::AbstractVector, name)
    @assert length(data) == length(msh.nodes)
    target[name] = @view data[node_map]
    target
end

function add_vector_data(target, msh::Mesh, node_map::Vector{Int}, data::AbstractMatrix, name)
    @assert size(data, 2) == length(msh.nodes)
    target[name] = @view data[:, node_map]
    target
end

function deserialize(stream::IO, ::Type{NodeInfo})
    coords = deserialize(stream, MVector{2, Float64})
    adjacent = deserialize(stream, Vector{Int})
    NodeInfo(coords, adjacent)
end

function deserialize(stream::IO, ::Type{MVector{N, T}}) where {N, T}
    x = @MVector zeros(T, N)
    read!(stream, x)
end

function deserialize(stream::IO, ::Type{Vector{T}}) where T
    sz = read(stream, Int)
    x = Vector{T}(undef, sz)
    read!(stream, x)
end

function deserialize(stream::IO, ::Type{ElementInfo{NodesPerFace, InternalNodes}}) where {NodesPerFace, InternalNodes}
    control_nodes = deserialize(stream, MVector{3, Int})
    faces = deserialize(stream, MVector{3, Int})
    face_nodes = @MVector [deserialize(stream, MVector{NodesPerFace, Int}) for i in 1:3]
    internal_nodes = deserialize(stream, MVector{InternalNodes, Int})
    adjacent_elements = deserialize(stream, Vector{Int})
    ElementInfo(control_nodes, faces, face_nodes, internal_nodes, adjacent_elements)
end

function deserialize(stream::IO, ::Type{BoundaryRepresentation{N}}) where N
    nodes = deserialize(stream, Vector{Int})
    faces = Vector{FaceDetails{N}}(undef, read(stream, Int))
    for i in eachindex(faces)
        faces[i] = deserialize(stream, FaceDetails{N})
    end
    BoundaryRepresentation{N}(nodes, faces)
end

function deserialize(stream::IO, ::Type{FaceDetails{N}}) where N
    number = read(stream, Int)
    element = read(stream, Int)
    nodes = deserialize(stream, MVector{N, Int})
    FaceDetails{N}(number, element, nodes)
end

function deserialize(stream::IO, ::Type{Vector{String}})
    count = read(stream, Int)
    [deserialize(stream, String) for i in 1:count]
end

function deserialize(stream::IO, ::Type{String})
    count = read(stream, Int)
    String(read!(stream, Vector{UInt8}(undef, count)))
end

end # module TetMesh

max_stars_repo_name: slmcbane/Le-TO

max_issues_repo_name: slmcbane/Le-TO

max_forks_repo_name: slmcbane/Le-TO


text:
include("differentiation.jl")
using Plots

h = 2.0 .^ (-2:-1:-52)
ffderr = []
bfderr = []
cfderr = []
f(x) = exp(x)
x = 1.0
dfexact = exp(1.0)
for h0 in h
    append!(ffderr, abs(ffd(f, x, h0 )-dfexact))
    append!(bfderr, abs(bfd(f, x, h0 )-dfexact))
    append!(cfderr, abs(cfd(f, x, h0 )-dfexact))
end

plot(h, ffderr, xaxis=:log, yaxis=:log)
plot!(h, bfderr, xaxis=:log, yaxis=:log)
plot!(h, cfderr, xaxis=:log, yaxis=:log)


text:
module BasicBlockRewriter

using Core: GotoNode, GotoIfNot, ReturnNode, SlotNumber, SSAValue, NewvarNode, SimpleVector
using Core: svec
using Base: mapany

export fragments

function relocatable_fragment(stmts, ssa1::Integer)
    toslot = Dict{Union{SlotNumber,SSAValue},SlotNumber}()
    stmtso = Any[]
    for (i, stmt) in enumerate(stmts)
        if isa(stmt, GotoNode) || isa(stmt, GotoIfNot) || isa(stmt, ReturnNode)
            break
        end
        push!(stmtso, rewrite_statement(stmt, ssa1:ssa1+length(stmts)-1, i+ssa1-1, toslot))
    end
    return stmtso
end

function rewrite_statement(stmt, rng, ssa, toslot)
    # Pass over some concretely-typed values that appear in code
    # These are not exhaustive, but for the time being I'd rather hear about unhandled things
    # than to have surprises.
    stmt === nothing && return stmt
    isa(stmt, Symbol) && return stmt
    isa(stmt, String) && return stmt
    isa(stmt, Char) && return stmt
    isa(stmt, LineNumberNode) && return nothing   # FIXME
    isa(stmt, Module) && return stmt
    isa(stmt, UnionAll) && return stmt
    isa(stmt, VersionNumber) && return stmt
    isa(stmt, Regex) && return stmt
    isa(stmt, Core.Box) && return stmt
    # Handle some "nontrival" code elements
    isa(stmt, NewvarNode) && return NewvarNode(rewrite_statement(stmt.slot, rng, ssa, toslot))
    isa(stmt, SlotNumber) && return get!(() -> SlotNumber(length(toslot)+1), toslot, stmt)
    if isa(stmt, SSAValue)
        id = stmt.id
        return id ∈ rng ? SSAValue(id - first(rng) + 1) : get!(() -> SlotNumber(length(toslot)+1), toslot, stmt)
    end
    isa(stmt, QuoteNode) && return QuoteNode(rewrite_statement(stmt.value, rng, ssa, toslot))
    if isa(stmt, Expr)
        out = Expr(stmt.head)
        out.args = mapany(arg -> rewrite_statement(arg, rng, ssa, toslot), stmt.args)
        return out
    end
    if isa(stmt, GlobalRef)
        return GlobalRef(rewrite_statement(stmt.mod, rng, ssa, toslot), rewrite_statement(stmt.name, rng, ssa, toslot))
    end
    isa(stmt, SimpleVector) && return svec([rewrite_statement(stmt[i], rng, ssa, toslot) for i = 1:length(stmt)]...)
    # The remainder are not concretely typed.
    # Put these last to avoid slow subtyping as often as possible.
    isa(stmt, Tuple) && return map(item -> rewrite_statement(item, rng, ssa, toslot), stmt)
    # More passing-over (not concretely typed)
    isa(stmt, Number) && return stmt
    (isa(stmt, AbstractString) || isa(stmt, Base.CodeUnits)) && return stmt
    (isa(stmt, DataType) || stmt === Union{}) && return stmt
    (isa(stmt, Val) || isa(stmt, Enum)) && return stmt
    isa(stmt, Function) && return stmt
    error("unhandled statement ", stmt, " of type ", typeof(stmt))
end

function fragments(m::Method)
    if isdefined(m, :generator) && m.generator isa Core.GeneratedFunctionStub
        @warn "skipping generated function $m"
        return Any[]
    end
    src = Base.uncompressed_ast(m)
    cfg = Core.Compiler.compute_basic_blocks(src.code)
    return Any[relocatable_fragment(src.code[rng.start:rng.stop], rng.start) for rng in map(bb->bb.stmts, cfg.blocks)]
end

end