text:
using PkgDev
using Base.Test
using Compat
import Base.Pkg.PkgError

function temp_pkg_dir(fn::Function, remove_tmp_dir::Bool=true)
    # Used in tests below to setup and teardown a sandboxed package directory
    const tmpdir = ENV["JULIA_PKGDIR"] = joinpath(tempdir(),randstring())
    @test !isdir(Pkg.dir())
    try
        Pkg.init()
        @test isdir(Pkg.dir())
        Pkg.resolve()

        fn(Pkg.Dir.path())
    finally
        remove_tmp_dir && rm(tmpdir, recursive=true)
    end
end

temp_pkg_dir() do pkgdir

    @testset "testing a package with test dependencies causes them to be installed for the duration of the test" begin
        PkgDev.generate("PackageWithTestDependencies", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]
        @test readstring(Pkg.dir("PackageWithTestDependencies","REQUIRE")) == "julia $(PkgDev.Generate.versionfloor(VERSION))\n"

        isdir(Pkg.dir("PackageWithTestDependencies","test")) || mkdir(Pkg.dir("PackageWithTestDependencies","test"))
        open(Pkg.dir("PackageWithTestDependencies","test","REQUIRE"),"w") do f
            println(f,"Example")
        end

        open(Pkg.dir("PackageWithTestDependencies","test","runtests.jl"),"w") do f
            println(f,"using Base.Test")
            println(f,"@test haskey(Pkg.installed(), \"Example\")")
        end

        Pkg.resolve()
        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

        Pkg.test("PackageWithTestDependencies")

        @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

        # trying to pin an unregistered package errors
        try
            Pkg.pin("PackageWithTestDependencies", v"1.0.0")
            error("unexpected")
        catch err
            @test isa(err, PkgError)
            @test err.msg == "PackageWithTestDependencies cannot be pinned – not a registered package"
        end
    end

    @testset "testing a package with no runtests.jl errors" begin
        PkgDev.generate("PackageWithNoTests", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        if isfile(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
            rm(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
        end

        try
            Pkg.test("PackageWithNoTests")
            error("unexpected")
        catch err
            @test err.msg == "PackageWithNoTests did not provide a test/runtests.jl file"
        end
    end

    @testset "testing a package with failing tests errors" begin
        PkgDev.generate("PackageWithFailingTests", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        isdir(Pkg.dir("PackageWithFailingTests","test")) || mkdir(Pkg.dir("PackageWithFailingTests","test"))
        open(Pkg.dir("PackageWithFailingTests", "test", "runtests.jl"),"w") do f
            println(f,"using Base.Test")
            println(f,"@test false")
        end

        try
            Pkg.test("PackageWithFailingTests")
            error("unexpected")
        catch err
            @test err.msg == "PackageWithFailingTests had test errors"
        end
    end

    @testset "testing with code-coverage" begin
        PkgDev.generate("PackageWithCodeCoverage", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))

        src = """
module PackageWithCodeCoverage

export f1, f2, f3, untested

f1(x) = 2x
f2(x) = f1(x)
function f3(x)
    3x
end
untested(x) = 7

end"""
        linetested = [false, false, false, false, true, true, false, true, false, false]
        open(Pkg.dir("PackageWithCodeCoverage", "src", "PackageWithCodeCoverage.jl"), "w") do f
            println(f, src)
        end
        isdir(Pkg.dir("PackageWithCodeCoverage","test")) || mkdir(Pkg.dir("PackageWithCodeCoverage","test"))
        open(Pkg.dir("PackageWithCodeCoverage", "test", "runtests.jl"),"w") do f
            println(f,"using PackageWithCodeCoverage, Base.Test")
            println(f,"@test f2(2) == 4")
            println(f,"@test f3(5) == 15")
        end

        Pkg.test("PackageWithCodeCoverage")
        covdir = Pkg.dir("PackageWithCodeCoverage","src")
        covfiles = filter!(x -> contains(x, "PackageWithCodeCoverage.jl") && contains(x,".cov"), readdir(covdir))
        @test isempty(covfiles)
        Pkg.test("PackageWithCodeCoverage", coverage=true)
        covfiles = filter!(x -> contains(x, "PackageWithCodeCoverage.jl") && contains(x,".cov"), readdir(covdir))
        @test !isempty(covfiles)
        for file in covfiles
            @test isfile(joinpath(covdir,file))
            covstr = readstring(joinpath(covdir,file))
            srclines = split(src, '\n')
            covlines = split(covstr, '\n')
            for i = 1:length(linetested)
                covline = (linetested[i] ? "        1 " : "        - ")*srclines[i]
                @test covlines[i] == covline
            end
        end
    end

    @testset "testing package tags" begin
        PkgDev.generate("PackageWithTags", "MIT", config=Dict("user.name"=>"Julia Test", "user.email"=>"test@julialang.org"))
        PkgDev.tag("PackageWithTags")
        PkgDev.tag("PackageWithTags")

        # check tags
        repo = LibGit2.GitRepo(joinpath(pkgdir,"PackageWithTags"))
        try
            tags = LibGit2.tag_list(repo)
            @test ("v0.0.1" in tags) == true
            @test ("v0.0.2" in tags) == true
        finally
            finalize(repo)
        end
    end

    @testset "testing freeable" begin
        Pkg.add("Example")
        io = IOBuffer()
        f = PkgDev.freeable(io)
        @test !(any(f .== "Example") || contains(takebuf_string(io), "Example"))
        Pkg.checkout("Example")
        f = PkgDev.freeable(io)
        @test any(f .== "Example") || contains(takebuf_string(io), "Example")
    end

    @testset "testing package registration" begin
        if haskey(ENV, "CI") && lowercase(ENV["CI"]) == "true"
            info("setting git global configuration")
            run(`git config --global user.name "Julia Test"`)
            run(`git config --global user.email test@julialang.org`)
            run(`git config --global github.user JuliaTest`)
        end
        PkgDev.generate("GreatNewPackage", "MIT")
        PkgDev.register("GreatNewPackage")
        @test !isempty(readstring(joinpath(pkgdir, "METADATA", "GreatNewPackage", "url")))
    end
end

@testset "Testing package utils" begin
    @test PkgDev.getrepohttpurl("https://github.com/JuliaLang/PkgDev.jl.git") == "https://github.com/JuliaLang/PkgDev.jl"
    @test PkgDev.getrepohttpurl("git://github.com/JuliaLang/PkgDev.jl.git")  == "https://github.com/JuliaLang/PkgDev.jl"
end


text:
"""
    train!(;
        valueSelectionArray::Union{T, Array{T, 1}},
        generator::AbstractModelGenerator,
        nbEpisodes::Int64=10,
        strategy::S=DFSearch(),
        variableHeuristic::AbstractVariableSelection=MinDomainVariableSelection(),
        out_solver::Bool=false,
        verbose::Bool=true,
        evaluator::Union{Nothing, AbstractEvaluator},
        metrics::Union{Nothing,AbstractMetrics}=nothing
    ) where{ T <: ValueSelection, S <: SearchStrategy}

Training the given LearnedHeuristics and using the Basic ones to compare performances.
This function managed the training mode of the LearnedHeuristic before and after a call to `launch_experiment!`.

The function return arrays containing scores, time needed and numbers of nodes visited on each episode for every heuristic
(basic and learned heuristics)
"""
function train!(;
        valueSelectionArray::Union{T, Array{T, 1}},
        generator::AbstractModelGenerator,
        nbEpisodes::Int64=10,
        strategy::S1=DFSearch(),
        eval_strategy::S2=strategy,
        variableHeuristic::AbstractVariableSelection=MinDomainVariableSelection(),
        out_solver::Bool=false,
        verbose::Bool=true,
        evaluator::Union{Nothing, AbstractEvaluator},
        metrics::Union{Nothing,AbstractMetrics}=nothing,
        restartPerInstances = 1,
    ) where{ T <: ValueSelection, S1, S2 <: SearchStrategy}

    if isa(valueSelectionArray, T)
        valueSelectionArray = [valueSelectionArray]
    end

    for valueSelection in valueSelectionArray
        if isa(valueSelection, LearnedHeuristic)
            valueSelection.fitted_problem = typeof(generator)
            valueSelection.fitted_strategy = typeof(strategy)
            # we could add more information later ...

            # make sure it is in training mode
            testmode!(valueSelection, false)
        end
    end

    metricsArray, eval_metricsArray = launch_experiment!(
        valueSelectionArray,
        generator,
        nbEpisodes,
        strategy,
        eval_strategy,
        variableHeuristic,
        out_solver,
        verbose;
        metrics=metrics,
        evaluator=evaluator,
        restartPerInstances,
    )

    for valueSelection in valueSelectionArray
        if isa(valueSelection, LearnedHeuristic)
            # go to testing mode
            testmode!(valueSelection)

            if verbose
                print("Has been trained on : ", typeof(generator))
                print(" with strategy : ", strategy)
                print(" during ", nbEpisodes, " episodes ")
                restartPerInstances > 1 && print("with ", restartPerInstances, " restart per episode ")
                out_solver && println("out of the solver.")
                !out_solver && println("in the solver.")
                println("Training mode now desactivated !")
            end
        end
    end

    return metricsArray, eval_metricsArray
end

ext: jl


text:
using Test
using CapacityExpansion
using Clp

@testset "CEP" begin
    @testset "TX_1" begin #Compare to Merrick Testcase
        # load data
        ts_input_data = load_timeseries_data_provided("TX_1"; T=24, years=[2008])
        cep_input_data=load_cep_data_provided("TX_1")
        # run clustering
        ts_clust_res = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=1,n_clust=365) # default k-means
        # run optimization
        model = run_opt(ts_clust_res.best_results,cep_input_data,Clp.Optimizer;optimizer_config=Dict{Symbol,Any}(:LogLevel => 0))
        # compare to exact result
        exact_res=[70540.26439790576;0.0;8498.278397905757;0.0;80132.88454450261]
        @test exact_res ≈ model.variables["CAP"].data[:,1,1] atol=1
    end
    scenarios=Dict{String,OptResult}()
    #Test the workflow for single node scenarios for Germany and California
    @testset "workflow $state" for (state, years) in [["GER_1", [2016]],["CA_1", [2016]]] begin
        # laod data
        ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
        cep_data = load_cep_data_provided(state)
        ## CLUSTERING ##
        ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
        ts_full_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=30)
        ## OPTIMIZATION ##
        optimizer=Clp.Optimizer
        optimizer_config=Dict{Symbol,Any}(:LogLevel => 0)
        scenarios["$state-$years-co2"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="co2",co2_limit=1000,optimizer_config=optimizer_config)
        scenarios["$state-$years-slack"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="slack",lost_el_load_cost=1e6, lost_CO2_emission_cost=700,optimizer_config=optimizer_config)
        scenarios["$state-$years-ex"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="ex",existing_infrastructure=true,optimizer_config=optimizer_config)
        scenarios["$state-$years-simple"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="simple storage",storage="simple",optimizer_config=optimizer_config)
        scenarios["$state-$years-seasonal"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="seasonal storage",storage="seasonal",optimizer_config=optimizer_config)
        design_result=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="des&op",optimizer_config=optimizer_config)
        scenarios["$state-$years-des&op"] = run_opt(ts_full_data.best_results,cep_data,design_result.opt_config,get_cep_design_variables(design_result),optimizer;lost_el_load_cost=1e6,lost_CO2_emission_cost=700)
        end
    end
    #Test transmission for a multi-node scenario
    @testset "workflow $state" for (state, years) in [["GER_18", [2016]],["CA_14", [2016]]] begin
           # laod data
           ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
           cep_data = load_cep_data_provided(state)
           ## CLUSTERING ##
           ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
           ## OPTIMIZATION ##
           optimizer=Clp.Optimizer
           scenarios["$state-$years-trans"] = run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="trans",transmission=true,optimizer_config=Dict{Symbol,Any}(:LogLevel => 0))
       end
    end
    #Test exact values for each of the previously calculated scenarios by comparison with exact scenarios
    @load normpath(joinpath(dirname(@__FILE__),"cep_exact_data.jld2")) exact_scenarios
    @testset "exact $k" for (k,v) in scenarios begin
            #Test for each variable within this scenario
            @testset "variable $kv" for (kv,vv) in v.variables begin
                    @test sum((round.(exact_scenarios[k].variables[kv].data;sigdigits=4).==round.(vv.data;sigdigits=4)).==false)==0
                end
            end
        end
    end
end


text:
function solve_modelKARC(H, g, nlp_stop, X::PDataKARC, α::T) where {T}
    # target value should be close to satisfy αλ=||d||
    target = [(abs(α * X.shifts[i] - X.norm_dirs[i])) for i = 1:X.nshifts]
    positives = findall(X.positives)

    # pick the closest shift to the target within positive definite H+λI
    X.indmin = max(1, X.indmin)
    start = X.indmin
    bidon, indmin = findmin(target[positives[start:end]])
    X.indmin = start + indmin - 1
    p_imin = positives[X.indmin]
    X.d = X.xShift[p_imin]
    X.λ = X.shifts[p_imin]

    return X.d, X.λ
end


text:
# This file was generated, do not modify it. # hide
plot(mach_tuned_iterated_pipe, size=(600,450))

savefig(joinpath(@OUTPUT, "EX-telco-tuning.svg")) # hide

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using AstroBase
using RemoteFiles
using SPICE: furnsh, kclear, sxform, tisbod, bodfnd
using Test

pck = @RemoteFile(
    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc",
    dir=joinpath(@__DIR__(), "data"),
)
download(pck)
furnsh(path(pck))

@testset "Frames" begin
    @testset "Rotations" begin
        ep = TDBEpoch(2019, 5, 6)
        rot = Rotation(icrf, icrf, ep)
        @test origin(rot) == icrf
        @test target(rot) == icrf

        rot′ = inv(rot)
        @test rot == rot′

        r = randn(3)
        v = randn(3)

        rv = (r, v)
        r1, v1 = rot(r, v)
        rv1 = rot(rv)
        @test r1 == r
        @test v1 == v
        @test rv1 == rv

        comp = rot ∘ rot
        r1, v1 = comp(r, v)
        rv1 = comp(rv)
        @test r1 == r
        @test v1 == v
        @test rv1 == rv

        exp_path = [:IAUEarth, :ICRF, :CIRF, :TIRF, :ITRF]
        @test AstroBase.Frames.path_frames(iau_earth, itrf) == exp_path
    end
    @testset "Composition" begin
        ep = TDBEpoch(2019, 5, 6)
        et = value(seconds(j2000(ep)))
        rot = Rotation(iau_earth, iau_io, ep)
        m_act = rot.m
        m′_act = rot.m′
        s_exp = sxform("IAU_EARTH", "IAU_IO", et)
        m_exp = s_exp[1:3, 1:3]
        m′_exp = s_exp[4:6, 1:3]
        @testset for i in eachindex(m_act, m_exp)
            @test m_act[i] ≈ m_exp[i] atol = 1e-8
        end
        @testset for i in eachindex(m′_act, m′_exp)
            @test m′_act[i] ≈ m′_exp[i] atol = 1e-6
        end
    end
    @testset "IAU" begin
        ep = TDBEpoch(2019, 5, 6)
        et = value(seconds(j2000(ep)))
        @testset for body_name in AstroBase.Bodies.ALL_NAMES
            body = @eval $(Symbol(body_name))()
            id = naifid(body)
            if bodfnd(id, "PM")
                frame = @eval $(Symbol("IAU", body_name))()
                rot = Rotation(icrf, frame, ep)
                m_act = rot.m
                m′_act = rot.m′
                s_exp = tisbod("J2000", id, et)
                m_exp = s_exp[1:3, 1:3]
                m′_exp = s_exp[4:6, 1:3]
                @testset for i in eachindex(m_act, m_exp)
                    @test m_act[i] ≈ m_exp[i] atol = 1e-8
                end
                @testset for i in eachindex(m′_act, m′_exp)
                    @test m′_act[i] ≈ m′_exp[i] atol = 1e-6
                end
            end
        end
        rot1 = Rotation(iau_mercury, iau_earth, ep)
        rot2 = Rotation(iau_mercury, iau_earth, UT1Epoch(ep))
        @test rot1.m == rot2.m
        @test rot1.m′ == rot2.m′
    end
    @testset "IERS" begin
        # Reference values from Orekit (http://www.orekit.org)
        tdb = TDBEpoch(2013, 3, 18, 12, 0, 0.0)
        r = [8.59072560e+02, -4.13720368e+03, 5.29556871e+03]
        v = [7.37289205e+00, 2.08223573e+00, 4.39999794e-01]

        icrf_cirf = Rotation(icrf, cirf, tdb)
        cirf_icrf = Rotation(cirf, icrf, tdb)
        ref_icrf_cirf = [0.9999991427733881 -1.144503732190422E-8 -0.0013093710276537447;
            4.871197292135239E-8 0.9999999995949649 2.8461675919057616E-5;
            0.001309371026797659 -2.8461715302997654E-5 0.9999991423683543]
        r_cirf = [852.1380066867667, -4137.052915716736, 5296.806764985966]
        v_cirf = [7.372309565810324, 2.08224861625441, 0.4495940103957848]
        @test icrf_cirf.m ≈ ref_icrf_cirf
        @test all(icrf_cirf(r, v) .≈ (r_cirf, v_cirf))
        @test all(cirf_icrf(icrf_cirf(r, v)...) .≈ (r, v))
        @test all((cirf_icrf ∘ icrf_cirf)(r, v) .≈ (r, v))

        cirf_tirf = Rotation(cirf, tirf, tdb)
        tirf_cirf = Rotation(tirf, cirf, tdb)
        ref_cirf_tirf = [0.9972630428201841 -0.0739352650974288 0.0;
            0.0739352650974288 0.9972630428201841 0.0;
            0.0 0.0 0.9999999999999998]
        r_tirf = [1155.6798454967793, -4062.7269296118056, 5296.806764985964]
        v_tirf = [6.901921544550409, 2.537349749040834, 0.44959401039578467]
        @test cirf_tirf.m ≈ ref_cirf_tirf
        @test all(cirf_tirf(r_cirf, v_cirf) .≈ (r_tirf, v_tirf))
        @test all(tirf_cirf(cirf_tirf(r_cirf, v_cirf)...) .≈ (r_cirf, v_cirf))
        @test all((tirf_cirf ∘ cirf_tirf)(r_cirf, v_cirf) .≈ (r_cirf, v_cirf))

        tirf_itrf = Rotation(tirf, itrf, tdb)
        itrf_tirf = Rotation(itrf, tirf, tdb)
        ref_tirf_itrf = [0.9999999999999787 -3.010092334926554E-11 2.0567742763960954E-7;
            3.046034574888695E-11 0.999999999998473 -1.7475053232215554E-6;
            -2.0567742758669397E-7 1.7475053232277836E-6 0.9999999999984519]
        r_itrf = [1155.6809350526369, -4062.7361857684173, 5296.799427643569]
        v_itrf = [6.9019216369452225, 2.537348963579269, 0.4495970248578132]
        @test tirf_itrf.m ≈ ref_tirf_itrf
        @test all(tirf_itrf(r_tirf, v_tirf) .≈ (r_itrf, v_itrf))
        @test all(itrf_tirf(tirf_itrf(r_tirf, v_tirf)...) .≈ (r_tirf, v_tirf))
        @test all((itrf_tirf ∘ tirf_itrf)(r_tirf, v_tirf) .≈ (r_tirf, v_tirf))

        icrf_itrf = icrf_cirf ∘ cirf_tirf ∘ tirf_itrf
        r_exp, v_exp = tirf_itrf(cirf_tirf(icrf_cirf(r, v)))
        r_act, v_act = icrf_itrf(r, v)
        @testset for i in eachindex(r_exp, r_act)
            @test r_exp[i] ≈ r_act[i]
        end
        @testset for i in eachindex(v_exp, v_act)
            @test v_exp[i] ≈ v_act[i]
        end

        ep = TDBEpoch(2019, 5, 6)
        rot1 = Rotation(icrf, itrf, ep)
        rot2 = Rotation(icrf, itrf, UT1Epoch(ep))
        @test rot1.m == rot2.m
        @test rot1.m′ == rot2.m′
    end
end

kclear()


text:
export Variation

"""
`Variation([domainType=Float64::Type,] dim_in::Tuple)`

`Variation(dims...)`

`Variation(x::AbstractArray)`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns a matrix with its `i`th column consisting of the vectorized discretized gradient over the `i`th `direction obtained using forward finite differences.

```julia
julia> Variation(Float64,(10,2))
Ʋ  ℝ^(10, 2) -> ℝ^(20, 2)

julia> Variation(2,2,2)
Ʋ  ℝ^(2, 2, 2) -> ℝ^(8, 3)

julia> Variation(ones(2,2))*[1. 2.; 1. 2.]
4×2 Array{Float64,2}:
 0.0  1.0
 0.0  1.0
 0.0  1.0
 0.0  1.0

```

"""
struct Variation{T,N} <: LinearOperator
	dim_in::NTuple{N,Int}
end

# Constructors
#default constructor
function Variation(domainType::Type, dim_in::NTuple{N,Int}) where {N}
	N == 1 && error("use FiniteDiff instead!")
	Variation{domainType,N}(dim_in)
end

Variation(domainType::Type, dim_in::Vararg{Int}) = Variation(domainType, dim_in)
Variation(dim_in::NTuple{N,Int}) where {N} = Variation(Float64, dim_in)
Variation(dim_in::Vararg{Int}) = Variation(dim_in)
Variation(x::AbstractArray)  = Variation(eltype(x), size(x))

# Mappings

@generated function mul!(y::AbstractArray{T,2},
                         A::Variation{T,N}, b::AbstractArray{T,N}) where {T,N}

	ex = :()

	for i = 1:N
		z = zeros(Int,N)
		z[i] = 1
		z = (z...,)
		ex = :($ex; y[cnt,$i] = I[$i] == 1 ? b[I+CartesianIndex($z)]-b[I] :
		       b[I]-b[I-CartesianIndex($z)])
	end

	ex2 = quote
		cnt = 0
		for I in CartesianIndices(size(b))
			cnt += 1
			$ex
		end
		return y
	end
end

@generated function mul!(y::AbstractArray{T,N},
                              A::AdjointOperator{Variation{T,N}}, b::AbstractArray{T,2}) where {T,N}

	ex = :(y[I] = I[1] == 1  ? -(b[cnt,1] + b[cnt+1,1]) :
	              I[1] == 2  ?   b[cnt,1] + b[cnt-1,1] - b[cnt+1,1] :
	       I[1] == size(y,1) ?   b[cnt,1] : b[cnt,  1] - b[cnt+1,1]
	       )

	Nx = :(size(y,1))
	for i = 2:N
		ex = quote
			$ex
			y[I] += I[$i] == 1  ? -(b[cnt,$i] + b[cnt+$Nx,$i]) :
			        I[$i] == 2  ?   b[cnt,$i] + b[cnt-$Nx,$i] - b[cnt+$Nx,$i] :
			        I[$i] == size(y,$i) ?   b[cnt,$i] : b[cnt,  $i]   - b[cnt+$Nx,$i]
		end
		Nx = :($Nx*size(y,$i))
	end

	ex2 = quote
		cnt = 0
		for I in CartesianIndices(size(y))
			cnt += 1
			$ex
		end
		return y
	end
end

# Properties

  domainType(L::Variation{T,N}) where {T,N} = T
codomainType(L::Variation{T,N}) where {T,N} = T

size(L::Variation{T,N}) where {T,N} = ((prod(L.dim_in), N), L.dim_in)

fun_name(L::Variation)  = "Ʋ"


text:
brng_type = VSL_BRNG_MT19937
@info("Testing continuous distributions with $brng_type...")
brng = BasicRandomNumberGenerator(brng_type, 200701281)

stdout_ = stdout
pwd_ = pwd()
cd(test_folder)
mkpath("./actual")
outfile = open("./actual/continuous.out", "w")
redirect_stdout(outfile)

t_uniform = Uniform(brng, 1.0, 5.0, VSL_RNG_METHOD_STD)
@test_rand(rand(t_uniform))
@test_rand(rand(t_uniform, 2, 3))

@test_throws ErrorException rand(t_uniform, UInt64)
@test_throws ErrorException rand!(t_uniform, zeros(UInt64, 5))

t_uniform_s = Uniform(brng, 1.0f0, 5.0f0, VSL_RNG_METHOD_STD)
@test_rand(rand(t_uniform_s))
@test_rand(rand(t_uniform_s, 2, 3))

t_gaussian = Gaussian(brng, 0.0, 2.0, VSL_RNG_METHOD_BOXMULLER)
@test_rand(rand(t_gaussian))
@test_rand(rand(t_gaussian, 2, 3))

t_gaussianmv = GaussianMV(brng, 3, VSL.VSL_MATRIX_STORAGE_FULL,
                          Float64[1,2,3], Float64[1 2 3; 1 2 3; 1 2 3], VSL_RNG_METHOD_BOXMULLER)
@test_rand(rand(t_gaussianmv))

t_exponential = Exponential(brng, 1.0, 5.0, VSL_RNG_METHOD_ICDF_ACCURATE)
@test_rand(rand(t_exponential))
@test_rand(rand(t_exponential, 2, 3))

t_laplace = Laplace(brng, 1.0, 5.0, VSL_RNG_METHOD_ICDF)
@test_rand(rand(t_laplace))
@test_rand(rand(t_laplace, 2, 3))

t_weibull = Weibull(brng, 2.0, 1.0, 5.0, VSL_RNG_METHOD_ICDF)
@test_rand(rand(t_weibull))
@test_rand(rand(t_weibull, 2, 3))

t_cauchy = Cauchy(brng, 1.0, 5.0, VSL_RNG_METHOD_ICDF)
@test_rand(rand(t_cauchy))
@test_rand(rand(t_cauchy, 2, 3))

t_rayleigh = Rayleigh(brng, 1.0, 5.0, VSL_RNG_METHOD_ICDF)
@test_rand(rand(t_rayleigh))
@test_rand(rand(t_rayleigh, 2, 3))

t_lognormal = Lognormal(brng, 1.0, 5.0, 1.0, 5.0, VSL_RNG_METHOD_BOXMULLER2_ACCURATE)
@test_rand(rand(t_lognormal))
@test_rand(rand(t_lognormal, 2, 3))

t_gumbel = Gumbel(brng, 1.0, 5.0, VSL_RNG_METHOD_ICDF)
@test_rand(rand(t_gumbel))
@test_rand(rand(t_gumbel, 2, 3))

t_gamma = Gamma(brng, 2.0, 1.0, 5.0, VSL_RNG_METHOD_GNORM)
@test_rand(rand(t_gamma))
@test_rand(rand(t_gamma, 2, 3))

t_beta = Beta(brng, 3.0, 2.0, 1.0, 5.0, VSL_RNG_METHOD_CJA_ACCURATE)
@test_rand(rand(t_beta))
@test_rand(rand(t_beta, 2, 3))

close(outfile)
redirect_stdout(stdout_)

@test success(`diff ./expected/continuous.out ./actual/continuous.out`)
cd(pwd_)


text:
devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, PSI.StaticLine),
                                     :T => DeviceModel(PSY.Transformer2W, PSI.StaticTransformer),
                                     :TT => DeviceModel(PSY.TapTransformer , PSI.StaticTransformer))
services = Dict{Symbol, PSI.ServiceModel}()

@testset "Solving ED with CopperPlate" begin
    template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services);
    parameters_value = [true, false]
    systems = [c_sys5, c_sys14]
    test_results = Dict{PSY.System, Float64}(c_sys5 => 240000.0,
                                             c_sys14 => 142000.0)
    for sys in systems, p in parameters_value
        @info("Testing solve ED with CopperPlatePowerModel network")
        @testset "ED CopperPlatePowerModel model use_parameters = $(p)" begin
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = OSQP_optimizer, use_parameters = p)
        psi_checksolve_test(ED, [MOI.OPTIMAL], test_results[sys], 10000)

        end
    end
end

@testset "Solving ED with PTDF Models" begin
    template = OperationsTemplate(StandardPTDFModel, devices, branches, services);
    parameters_value = [true, false]
    systems = [c_sys5, c_sys14, c_sys14_dc]
    PTDF_ref = Dict{PSY.System, PSY.PTDF}(c_sys5 => PTDF5, c_sys14 => PTDF14, c_sys14_dc => PTDF14_dc)
    test_results = Dict{PSY.System, Float64}(c_sys5 => 340000.0,
                                             c_sys14 => 142000.0,
                                             c_sys14_dc => 142000.0)

    for sys in systems, p in parameters_value
        @info("Testing solve ED with StandardPTDFModel network")
        @testset "ED StandardPTDFModel model use_parameters = $(p)" begin
        ED = OperationsProblem(TestOpProblem, template, sys; PTDF = PTDF_ref[sys], optimizer = OSQP_optimizer, use_parameters = p)
        psi_checksolve_test(ED, [MOI.OPTIMAL], test_results[sys], 10000)
        end
    end
end

@testset "Solving ED With PowerModels with loss-less convex models" begin
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters_value = [true, false]
    networks = [DCPPowerModel,
                NFAPowerModel]
    test_results = Dict{PSY.System, Float64}(c_sys5 => 330000.0,
                                             c_sys14 => 142000.0,
                                             c_sys14_dc => 142000.0)

    for  net in networks, p in parameters_value, sys in systems
        @info("Testing solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
        template = OperationsTemplate(net, devices, branches, services);
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = ipopt_optimizer, use_parameters = p);
        #The tolerance range here is large because NFA has a much lower objective value
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], test_results[sys], 35000)
        end
    end

end

@testset "Solving ED With PowerModels with linear convex models" begin
    systems = [c_sys5, c_sys14]
    parameters_value = [true, false]
    networks = [DCPLLPowerModel,
                LPACCPowerModel]
    test_results = Dict{PSY.System, Float64}(c_sys5 => 340000.0,
                                             c_sys14 => 142000.0,
                                             c_sys14_dc => 142000.0)

    for  net in networks, p in parameters_value, sys in systems
        @info("Testing solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
        template = OperationsTemplate(net, devices, branches, services);
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = ipopt_optimizer, use_parameters = p);
        #The tolerance range here is large because NFA has a much lower objective value
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], test_results[sys], 10000)

        end
    end

end

#=
@testset "Solving ED With PowerModels with convex SOC and QC models" begin
    systems = [c_sys5, c_sys14]
    parameters_value = [true, false]
    networks = [SOCWRPowerModel,
                 QCRMPowerModel,
                 QCLSPowerModel,]
    test_results = Dict{PSY.System, Float64}(c_sys5 => 320000.0,
                                             c_sys14 => 142000.0)

    for  net in networks, p in parameters_value, sys in systems
        @info("Testing solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
        template = OperationsTemplate(net, devices, branches, services);
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = ipopt_optimizer, use_parameters = p);
        #The tolerance range here is large because Relaxations have a lower objective value
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], test_results[sys], 25000)

        end
    end

end
=#

@testset "Solving ED With PowerModels Non-Convex Networks" begin
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters_value = [true, false]
    networks = [ACPPowerModel,
                #ACRPowerModel,
                ACTPowerModel]
    test_results = Dict{PSY.System, Float64}(c_sys5 => 340000.0,
                                             c_sys14 => 142000.0,
                                             c_sys14_dc => 142000.0)

    for  net in networks, p in parameters_value, sys in systems
        @info("Testing solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
        template = OperationsTemplate(net, devices, branches, services);
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = ipopt_optimizer, use_parameters = p);
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], test_results[sys], 10000)
        end
    end

end

@testset "Solving UC Linear Networks" begin
    devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
                                        :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
    parameters_value = [true, false]
    systems = [c_sys5, c_sys5_dc]
    networks = [DCPPowerModel,
                NFAPowerModel,
                StandardPTDFModel,
                CopperPlatePowerModel]
    PTDF_ref = Dict{PSY.System, PSY.PTDF}(c_sys5 => PTDF5, c_sys5_dc => PTDF5_dc)

    for  net in networks, p in parameters_value, sys in systems
        @info("Testing solve UC with $(net) network")
        @testset "UC model $(net) and use_parameters = $(p)" begin
        template= OperationsTemplate(net, devices, branches, services);
        UC = OperationsProblem(TestOpProblem, template, sys; PTDF = PTDF_ref[sys], optimizer = GLPK_optimizer, use_parameters = p)
        psi_checksolve_test(UC, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], 340000, 100000)
        end
    end

end


text:
using BlockSystems
using Test
using ModelingToolkit

const MTK = ModelingToolkit

include("BlockSystems_test.jl")
include("transformations_test.jl")
include("function_generation_test.jl")
include("utils_test.jl")


text:
function isprime(n)
    if n == 1
        return false
    end
    if n == 2 || n == 3
        return true
    end

    for i in 2:isqrt(n)
        if n%i == 0
            return false
        else
            if i == 2
                i += 1
            else
                i += 2
            end
        end
    end
    return true
end


text:
"""
    simulate(μ, x)

Simulation of a dichotomous vector via the probit link.

# Arguments
- `μ::Vector{<:Real}`: real-valued vector of coefficients
- `x::Matrix{<:Real}`: real-valued vector or matrix of inputs

# Output
- `y::Vector{Bool}`: boolean vector of outcomes
- `p::Vector{float}`: vector of probabilities
- `z::Vector{float}`: vector of latent variables
"""

function simulate(β::Vector{T} where T<:Real, α::Real, X::Matrix{T} where T<:Real)
    N, J = size(X)
    @assert length(β) == J "μ is a $(length(β)) vector and X is a $Nx$J matrix."
    @assert α > 0 "σ must be a positive scalar."
    λ = lambda(β, α, X)
    Time = rand.(Weibull.(α, λ))
    Censor = rand.(Weibull.(α, λ))
    Indicator = Time .> Censor
    Time = ifelse.(Indicator, Censor, Time)
    return (Time, Indicator, λ, β, α)
end


text:
"""
    (oe, avg, cgi) = CpG_islands(positions::Array, refseq::NamedTuple)

Perform a sliding window scans for regions of the genome to identify areas enriched in CpG dinucleotides. The default CpG content (0.6) and average GC content (0.5) can be changed by the user.

Returns an array of the obs/exp ratios for CpG content, the average GC values, and CpG island flags by windowfor each input region in the bedfile.

kwargs: window=200, step=200, noslide=false, cpg=0.6, gcavg=0.5

"""
function CpG_islands(positions::Array, refseq::NamedTuple; window=200, step=200, noslide=false, cpg=0.6, gcavg=0.5)

    cpglim = cpg    #min CpG dinucleotide cutoff
    gclim  = gcavg  #min avg GC content

    obs = CpG_windows(positions, refseq, window=window, step=step, noslide=noslide)
    exp = CpG_windows_expected(positions, refseq, window=window, step=step, noslide=noslide)
    avg = gc_windows(positions, refseq, window=window, step=step, noslide=noslide)

    exp = map(x->round.(x, digits=4), exp)
    avg = map(x->round.(x, digits=4), avg)

    oe = Vector{Float64}[]

    for i in eachindex(obs)
        ratio = round.( (obs[i] ./ exp[i]), digits=4)
        push!(oe, ratio)
    end

    all_cgi = Vector{Any}[]

    for i in eachindex(oe)
        cgi = Int64[]
        for j in eachindex(oe[i])
            if oe[i][j] > cpglim && avg[i][j] > gclim
                push!(cgi, 1)
            else
                push!(cgi,0)
            end
        end
        push!(all_cgi, cgi)
    end

    return (oe, avg, all_cgi)

end


text:
using Test
using GRASP

@testset "GRASP.jl" begin

@eval module SymmetriesTests include("symmetries.jl") end
@eval module CSFSTests include("Configurations.jl") end
@eval module RCSFSTests include("rcsfs.jl") end

@testset "libgrasp" begin
    # Make sure that all files exist
    @test isfile(joinpath(@__DIR__, "grasp/mixing/rmix.out"))

    @test_throws ErrorException GRASP.read_rmix("invalid-file.rmix")

    @test isa(
        GRASP.read_rmix(joinpath(@__DIR__, "grasp/mixing/rmix.out")),
        GRASP.MixingFile
    )

    @test isa(
        GRASP.read_rwfn(joinpath(@__DIR__, "grasp/mixing/rwfn.out")),
        Vector{GRASP.RWFNOrbital}
    )

    # make sure that existing, but invalid files throw gracefully
    @info "Expected error output:"
    @test_throws ErrorException GRASP.read_rmix(joinpath(@__DIR__, "grasp/mixing/rwfn.out"))
    @test_throws ErrorException GRASP.read_rwfn(joinpath(@__DIR__, "grasp/mixing/rmix.out"))
    @info "End of expected error output."

    include("binaries/binaries.jl")
end

@testset "other" begin
    @test GRASP.hartree2kayser(1.0) == 219474.63137
    @test GRASP.h2k_humanize(1.0) == "219,475"
end

end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Data: 06/09/2014

module Simulator

export
    AbstractSimulator,
    simulate,
    TCASSimulator

using AbstractSimulatorImpl

using TCASSimulatorImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
abstract type IntegralOperator <: Operator end

export quadrule, elements

"""
    quadrule(operator,test_refspace,trial_refspace,p,test_element,q_trial_element, qd)

Returns an object that contains all the dynamic (runtime) information that
defines the integration strategy that will be used by `momintegrals!` to compute
the interactions between the local test/trial functions defined on the specified
geometric elements. The indices `p` and `q` refer to the position of the test
and trial elements as encountered during iteration over the output of
`geometry`.

The last argument `qd` provides access to all precomputed data required for
quadrature. For example it might be desirable to precompute all the quadrature
points for all possible numerical quadrature schemes that can potentially be
required during matrix assembly. This makes sense, since the number of point is
order N (where N is the number of faces) but these points will appear in N^2
computations. Precomputation requires some extra memory but can save a lot on
computation time.
"""
function quadrule end

"""
  elements(geo)

Create an iterable collection of the elements stored in `geo`. The order in which
this collection produces the elements determines the index used for lookup in the
data structures returned by `assemblydata` and `quaddata`.
"""
#elements(geo) = [simplex(vertices(geo,cl)) for cl in cells(geo)]
elements(geo) = [chart(geo,cl) for cl in cells(geo)]

elements(sp::Space) = elements(geometry(sp))

"""
    assemblechunk!(biop::IntegralOperator, tfs, bfs, store)

Computes the matrix of operator biop wrt the finite element spaces tfs and bfs
"""
function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store)

    test_elements, tad = assemblydata(tfs)
    bsis_elements, bad = assemblydata(bfs)

    tshapes = refspace(tfs); num_tshapes = numfunctions(tshapes)
    bshapes = refspace(bfs); num_bshapes = numfunctions(bshapes)

    qd = quaddata(biop, tshapes, bshapes, test_elements, bsis_elements)
    zlocal = zeros(scalartype(biop, tfs, bfs), num_tshapes, num_bshapes)

    assemblechunk_body!(biop,
        tshapes, test_elements, tad,
        bshapes, bsis_elements, bad,
        qd, zlocal, store)
end

function assemblerow!(biop::IntegralOperator, test_functions::Space, trial_functions::Space, store)

    test_elements = elements(geometry(test_functions))
    trial_elements, trial_assembly_data = assemblydata(trial_functions)

    test_shapes = refspace(test_functions); num_test_shapes = numfunctions(test_shapes)
    trial_shapes = refspace(trial_functions); num_trial_shapes = numfunctions(trial_shapes)

    qd = quaddata(biop, test_shapes, trial_shapes, test_elements, trial_elements)
    zlocal = zeros(scalartype(biop, test_functions, trial_functions), num_test_shapes, num_trial_shapes)

    @assert size(zlocal) == (3,3)

    @assert length(trial_elements) == numcells(geometry(trial_functions))
    @assert numfunctions(test_functions) == 1
    test_function = test_functions.fns[1]
    for shape in test_function
        p = shape.cellid
        i = shape.refid
        a = shape.coeff
        #@show shape
        tcell = test_elements[p]
        for (q,bcell) in enumerate(trial_elements)

            #@show q

            fill!(zlocal, 0)
            strat = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, qd)
            momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, strat)

            #@assert size(zlocal,2) == 3
            for j in 1:size(zlocal,2)
                for (n,b) in trial_assembly_data[q,j]
                    #n == 3 && println("bingo")
                    #@show n
                    store(a*zlocal[i,j]*b, 1, n)
                    #store(1, 1, n)
end end end end end

# function assemblerow_body!(biop::IntegralOperator,
#     tshapes, test_function,
#     bshapes, trial_elements, trial_assembly_data,
#     qd, zlocal, store)
#
#     for shape in test_function
#         p = shape.cellid
#         tcell = test_elements[p]
#         for (q,bcell) in enumerate(trial_elements)
#
#             fill!(zlocal, 0)
#             strat = quadrule(biop, tshapes, bshapes, p, tcell, q, bcell, qd)
#             momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, strat)
#
#             i = shape.refid
#             a = shape.coeff
#             for j in size(zlocal,2)
#                 for (n,b) in trial_assembly_data[q,j]
#                     store(a*zlocal[i,j]*b, m, n)
# end end end end end

function assemblechunk_body!(biop,
        test_shapes, test_elements, test_assembly_data,
        trial_shapes, trial_elements, trial_assembly_data,
        qd, zlocal, store)

    for (p,tcell) in enumerate(test_elements), (q,bcell) in enumerate(trial_elements)

        fill!(zlocal, 0)
        strat = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, qd)
        momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, strat)

        for j in 1 : size(zlocal,2), i in 1 : size(zlocal,1)
            for (n,b) in trial_assembly_data[q,j], (m,a) in test_assembly_data[p,i]
                store(a*zlocal[i,j]*b, m, n)
end end end end

immutable DoubleQuadStrategy{P,Q}
  outer_quad_points::P
  inner_quad_points::Q
end

"""
    regularcellcellinteractions!(biop, tshs, bshs, tcell, bcell, interactions, strat)

Function for the computation of moment integrals using simple double quadrature.
"""
function momintegrals!(biop, tshs, bshs, tcell, bcell, z, strat::DoubleQuadStrategy)

    # memory allocation here is a result from the type instability on strat
    # which is on purpose, i.e. the momintegrals! method is chosen based
    # on dynamic polymorphism.
    womps = strat.outer_quad_points
    wimps = strat.inner_quad_points

    M, N = size(z)

    for womp in womps
        tgeo = womp.point
        tvals = womp.value
        jx = womp.weight

        for wimp in wimps
            bgeo = wimp.point
            bvals = wimp.value
            jy = wimp.weight

            j = jx * jy
            kernel = kernelvals(biop, tgeo, bgeo)

            for m in 1 : M
                tval = tvals[m]
                for n in 1 : N
                    bval = bvals[n]

                    igd = integrand(biop, kernel, tval, tgeo, bval, bgeo)
                    z[m,n] += j * igd
                end
            end
        end
    end

    return z
end

abstract type SingularityExtractionStrategy end
regularpart_quadrule(qr::SingularityExtractionStrategy) = qr.regularpart_quadrule

function momintegrals!(op, g, f, t, s, z, strat::SingularityExtractionStrategy)

    womps = strat.outer_quad_points

    sop = singularpart(op)
    rop = regularpart(op)

    # compute the regular part
    rstrat = regularpart_quadrule(strat)
    momintegrals!(rop, g, f, t, s, z, rstrat)

    for p in 1 : length(womps)
        x = womps[p].point
        dx = womps[p].weight

        innerintegrals!(sop, x, g, f, t, s, z, strat, dx)
    end # next quadrature point

end

type QuadData{WPV1,WPV2}
  tpoints::Matrix{Vector{WPV1}}
  bpoints::Matrix{Vector{WPV2}}
end


text:
module MaxHelpingHandReskinnableFloatingDebris

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/ReskinnableFloatingDebris" ReskinnableFloatingDebris(x::Integer, y::Integer, texture::String="scenery/debris")

const placements = Ahorn.PlacementDict(
    "Floating Debris (Reskinnable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableFloatingDebris
    )
)

function Ahorn.selection(entity::ReskinnableFloatingDebris)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 6, 12, 12)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableFloatingDebris, room::Maple.Room)
    debrisSprite = Ahorn.getSprite(entity.texture, "Gameplay")
    Ahorn.drawImage(ctx, debrisSprite, -4, -4, 0, 0, 8, 8)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
type Updating_WorkingSet
  ilon :: Int64
  ilat :: Int64
  ilonlat :: Int64
  mask :: Array{Int64,2}
  ratio :: Array{Float64,2}
  plon :: Array{Float64,2}
  plat :: Array{Float64,2}
end

function update_ratio(WS :: Updating_WorkingSet,
                      wesn_range :: Array{Float64,1}; stratage :: Symbol = :default)
  (lon_westforward,lon_eastforward,lat_southforward,lat_northforward) = wesn_range ;
  plon = WS.plon;
  plat = WS.plat;
  ilon = WS.ilon
  ilat = WS.ilat

  ratio_lat = WS.ratio[ilon,:]
  ratio_lon = WS.ratio[:,ilat]

  lats = plat[ilon,:]
  lons = plon[:,ilat]

  if( :default == stratage )
    d_lon = lons - lons[ilon]
    d_lon[d_lon .> 180 ] = d_lon[d_lon .> 180] -360
    d_lon[d_lon .< -180 ] = 360 + d_lon[d_lon .< -180]
    d_lat = lats - lats[ilat]
    WS.ratio[ilon,:] = ratio_lat .* linear_update_ratio(ilat,
                                     ratio_lat, lat_northforward, lat_southforward ,d_lat)
    WS.ratio[:,ilat] = ratio_lon .* linear_update_ratio(ilon,
                                     ratio_lon, lon_westforward , lon_eastforward  ,d_lon)
    #= WS.ratio = WS.ratio - 1e-5 =#
    nothing
  else
    error("No stratage $stratage available now")
  end
end

function linear_update_ratio(i :: Int64, ratio_line :: Array{Float64},
                              up_to :: Float64, down_to :: Float64, distance :: Array{Float64})
  m_line = ones(size(ratio_line))
  index_m = falses(size(ratio_line))
  if( length(distance) != i)
    #= if( 0 != ratio_line[i+1] ) =#
    if( true)
      index_m = (distance .>= 0) & (distance .<= up_to)
    end
  elseif ((distance[1] - distance[i] > 0) & (0 != ratio_line[1] ))
    print("cat here start")
    index_m = (distance .>= 0) & (distance .<= up_to)
  end
  m_line[index_m] = abs(distance[index_m])/up_to .* m_line[index_m]

  index_m = falses(size(ratio_line))
  if (1 != i )
    #= if(0 != ratio_line[i-1]) =#
    if(true)
      index_m = (distance .<0 ) & (distance .> - down_to)
    end
  elseif ((distance[end] - distance[1] < 0 ) & (0!= ratio_line[end]) )
    #= print("cat here end") =#
    index_m = (distance .<0 ) & (distance .> - down_to)
  end
  m_line[index_m] = abs(distance[index_m]/down_to)  .* m_line[index_m]
  return m_line
end

function write_ratio_mask(parea :: AnArea_Regular, sst :: Array{Float64,2};
                          stratage :: Symbol = :default,
                          wesn_range :: Array{Float64,1} = ones(4)*10.0,
                          whole_map_cir_lon :: Bool = true)
  bask = !isnan( Selected_within(parea,sst,  map_to = :whole ))
  #= bask1 = parea.plon_selected_on_by =#
  mask = int64(bask)
  ratio = float64(mask)
  if(haskey(parea.configs, "wplon") & haskey(parea.configs,"wplat") )
    wplon = parea.configs["wplon"]
    wplat = parea.configs["wplat"]
  else
    pw = AModel_Import(filename = parea.data_based_on )
    wplon = pw.plon
    wplat = pw.plat
  end
  WS = Updating_WorkingSet(1,1,1,mask, ratio, wplon, wplat)
  bask_cw = bask[2:end,:]
  bask_ce = bask[1:end-1,:]

  bask_cs = bask[:,2:end]
  bask_cn = bask[:,1:end-1]

  bask_cs = [bask_cs falses(size(bask_cs,1),1)]
  bask_cn = [falses(size(bask_cn,1),1) bask_cn]
  if(!whole_map_cir_lon)
    bask_cw = [bask_cw ; falses(1,size(bask_cw,2)) ]
    bask_ce = [falses(1,size(bask_ce,2)) ; bask_ce ]
  else
    bask_cw = [bask_cw; bask[1,:]]
    bask_ce = [bask[end,:] ; bask_ce ]
  end
  bask_c = copy(bask)
  bask_c = bask_ce | bask_cw | bask_cn | bask_cs
  bask_c[bask] = false

  idexlon = [1:size(bask_c,1)]
  idexlat = transpose([1:size(bask_c,2)])
  idexlon = repmat(idexlon,1,size(bask_c,2))
  idexlat = repmat(idexlat,size(bask_c,1),1)

  for ijlonlat = find(bask_c[:])
    WS.ilon = idexlon[ijlonlat]
    WS.ilat = idexlat[ijlonlat]
    WS.ilonlat = ijlonlat
    update_ratio(WS, wesn_range, stratage = stratage)
  end
  return (WS.ratio, WS.mask, bask_c)
end

export write_ratio_mask


text:
struct ExtendedHamiltonian{A,T,E<:Tuple,H<:Hamiltonian{A,T},O} <: Hamiltonian{A,T}
    hamiltonian::H
    extensions::E
    terms::O
end

terms(ext::ExtendedHamiltonian) = ext.terms
basis(ext::ExtendedHamiltonian) = basis(ext.hamiltonian)
LOStructure(ext::ExtendedHamiltonian) = LOStructure(ext.terms) # TODO?
starting_address(ext::ExtendedHamiltonian) = starting_address(ext.hamiltonian)

function Base.show(io::IO, ext::ExtendedHamiltonian)
    print(io, ext.hamiltonian)
    for e in ext.extensions
        print(io, " + ", e)
    end
end

"""
    abstract type ExtensionPrototype

A prototype for a potential (or other modifier that adds terms to existing operator). When
added to an operator, it can use information from it to construct appropriate terms.

# See also

* [`HarmonicPotential`](@ref)
* [`DeltaFunctionPotential`](@ref)
* [`TranscorrelatedDeltaPotential`](@ref)
"""
abstract type ExtensionPrototype end

"""
    initialize(prototype::AbstractOperator, operator::ExtensionPrototype)

Use the information in `prototype` and `operator` to construct an appropriate term to be
added to the operator.
"""
initialize

function Base.:+(ham::Hamiltonian, proto::ExtensionPrototype)
    term = initialize(proto, ham)
    return ExtendedHamiltonian(ham, (proto,), terms(ham) + term)
end
Base.:+(proto::ExtendedHamiltonian, ham::Hamiltonian) = ham + proto

function Base.:+(ext_ham::ExtendedHamiltonian, proto::ExtensionPrototype)
    base = ext_ham.hamiltonian
    extensions = ext_ham.extensions
    term = initialize(proto, base)

    return ExtendedHamiltonian(base, (extensions..., proto), terms(ext_ham) + term)
end


text:
using Distributions, Random
Random.seed!(1)

N = 10^6
sig = 1.7

data1 = sqrt.(-(2* sig^2)*log.(rand(N)))

distG = Normal(0,sig)
data2 = sqrt.(rand(distG,N).^2 + rand(distG,N).^2)

distR = Rayleigh(sig)
data3 = rand(distR,N)

mean.([data1, data2, data3])

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#NN.jl

#Working!!!! (for odd functions with sigmoids only...)
#TODO: figure out nice way to do learning rate annealing, changing the loss function: L2/MSE, L1, softmax?, huber!?!?
#TODO: pretrain stuff, inputs that result in max activations as opt prob
#TODO: code profiling
#TODO: LSTM/RNN stuff!!!
#NOTE: can probably reduce parameter size Float64->Float32 because papers
#NOTE: might be more efficient to have explicit Integer and Float versions rather than casting everything as real

module NN

import Base.length, Base.norm, Base.sqrt
export buildNN, train, predict #, pretrain

include("NNaux.jl")
include("neurons.jl")

function buildNN(sizes,c::Float64,yl::Array{Float64,1},yu::Array{Float64,1};
                 neuron_type::String="reLu",reg_type::String="L2",verbose::Bool=True,init_type::String="zero",output_neuron_type::String="none")
  #initialize weights and build neural network type based on these inputs
  if verbose
    println("Initializing Neural Network")
  end

  if lowercase(neuron_type) == "relu"
    neuron = reLu
    gradient = GreLu
  elseif lowercase(neuron_type) == "leakyrelu"
    neuron = leakyReLu
    gradient = GleakyReLu
  elseif lowercase(neuron_type) == "tanh"
    neuron = tanh
    gradient = Gtanh
  elseif lowercase(neuron_type) == "sigm"
    neuron = sigm
    gradient = Gsigm
  elseif lowercase(neuron_type) == "softplus"
    neuron = softplus
    gradient = Gsoftplus
    error("Gradient not working")
  else
    error(string("Undefined neuron type ",neron_type))
  end

  if lowercase(output_neuron_type) == "relu"
    output_neuron = reLu
    output_neuron_gradient = GreLu
  elseif lowercase(output_neuron_type) == "leakyrelu"
    output_neuron = leakyReLu
    output_neuron_gradient = GleakyReLu
  elseif lowercase(output_neuron_type) == "tanh"
    output_neuron = tanh
    output_neuron_gradient = Gtanh
  elseif lowercase(output_neuron_type) == "sigm"
    output_neuron = sigm
    output_neuron_gradient = Gsigm
  elseif lowercase(output_neuron_type) == "softplus"
    output_neuron = softplus
    output_neuron_gradient = Gsoftplus
    error("Gradient not working")
  elseif lowercase(output_neuron_type) == "none"
    output_neuron = (x)->x
    output_neuron_gradient = (x)->ones(size(x))
  else
    error(string("Undefined neuron type ",neron_type))
  end

  assert(length(yl) == length(yu))
  assert(!( false in (yl.<yu) ) )

  #assume full connectivity
  sizes_ = [(sizes[i+1],sizes[i]) for i = 1:length(sizes)-1]

  if lowercase(init_type) == "rand"
    weights = randweights(sizes_)
  elseif lowercase(init_type) == "zero"
    weights= zero(sizes_)
  else
    error(string("no such weight initialization method: ", init_type))
  end
  biases = zero(sizes[2:end])

  nn = NeuralNetwork(weights, biases, neuron, gradient,output_neuron,output_neuron_gradient, c, reg_type, yu, yl,[0],[0], verbose)

  return nn
  #maybe do some pretraining with stacked denoising autoencoders?
end

function buildNN(sizes,c::Float64,yl::Float64,yu::Float64;
                 neuron_type::String="reLu",reg_type::String="L2",verbose::Bool=True,init_type::String="zero",output_neuron_type::String="none")
  return buildNN(sizes,c,[yl],[yu],neuron_type=neuron_type,reg_type=reg_type,verbose=verbose,init_type=init_type,output_neuron_type=output_neuron_type)
end

function pretrain(X::Array{Real,2},Y::Array{Real,2},nn::NeuralNetwork;
                  rho::Float64=0.05,alpha::Float64=0.01,beta::Float64=0.1,minibatch_size::Int=100,num_epochs::Int=10)
  #rho = sparsity parameter
  #beta = sparsity penalty

  #HOWTODO: for each layer, build a 2-layer NN, use normal feedforward, use AEBackProp
  #  BUT feedforward needs to be a sigm activation function on the end, so do pre/post scaling first
  #just going to do vanilla GD since we would have to compute rho_ for each update

  #pretraining might not be necessary if using relu/softplus units?
  if !(nn.neuron == sigm || nn.neuron == tanh)
    error("Can only pretrain with sigmoid or tanh units")
  end
  error("Not yet implemented")

  if nn.verbose
    println("Pretraining...")
  end
  #need data: X and Y
  #pretrain with AE or RBM?

  m = size(X,2)
  assert(m == size(Y,2))
  L = length(nn.biases)

  for l = 2:L-1
    #ignore input and output layers

    #something
    #pass through data to do SGD
    for epoch = 1:num_epochs
      for i = 1:m
        #minibatch GD stuff
        #for each example, need feed forward pass, and backprop pass
      end
    end #epoch

  end #l

  return nn
end

#maybe also have a train! function, probably not necessary though
function train(X::Array{Float64,2},Y::Array{Float64,2},nn::NeuralNetwork;
               alpha::Float64=0.01,mu::Float64=0.9,minibatch_size::Int=100,num_epochs::Int=10,
               update_type::String="momentum",calculate_loss::Bool=true,half_rate::Int=1,loss_fn::String="L1")
  #mu is the "mass" for nesterov/momentum, and it is also the decay rate for rmsprop
  #each column of X is a training example, same with Y--we do this to access things in memory order
  m = size(X,2)
  assert(m == size(Y,2))

  #whiten data
  #it mightttttt be necessary to have a separate preproc data pass if theres too much data for memory
  nn.xu = vec(mean(X,2))
  nn.xs = vec(std(X,2))

  losses = Float64[]
  losses_true = Float64[]

  L = length(nn.biases)

  dW = zero(nn.weights)#Dict{Int,Array{Float64,2}}([i=>zeros(size(nn.weights[i])) for i = 1:L])
  db = zero(nn.biases)#Dict{Int,Array{Float64,2}}([i=>zeros(size(nn.biases[i])) for i = 1:L])
  #momentum update stuff/cache
  vW = zero(nn.weights)
  vb = zero(nn.biases)

  if m < minibatch_size
    minibatch_size = m
    #num_epochs = 1
    #^^ probably not a good idea to use the same data multiple times if you've only got 1 batch--NVM thats just gradient descent-durr
  end

  count = 0
  batch = 1
  s = ""
  for epoch = 1:num_epochs
    println("Epoch $epoch")
    if minibatch_size != m
      p = randperm(m)
    else
      p = [1:m]
    end
    for i_ = 1:m
      i = p[i_]
      x = (X[:,i]-nn.xu)./nn.xs
      y = Y[:,i]#(Y[:,i]-((nn.yl+nn.yu)/2))./((nn.yu-nn.yl)./2)
      y_, a, z = feedforward(x,nn)
      dW_, db_ = backprop(a,z,y,nn)

      dW += dW_
      db += db_
      #add!(dW,dW_)
      #add!(db,db_)

      #increment dW, db
      if mod(count,minibatch_size) == 0
        # fold all of this into a GDupdate function?
        #calculate loss function
        if calculate_loss
          loss = 0.
          for j = 1:m
            x = (X[:,j]-nn.xu)./nn.xs
            y = (Y[:,j])#-((nn.yl+nn.yu)/2))./((nn.yu-nn.yl)./2)
            y_, a, z = feedforward(x,nn)
            loss += 0.5*norm(y-y_,1)#(y-y_)'*(y-y_)#squared error
          end

          loss /= m

          losses = [losses;loss]

          loss += 0.5*nn.c*norm(nn.weights)
          loss += 0.5*nn.c*norm(nn.biases)
          losses_true = [losses_true;loss]
        end

        if nn.verbose
          print(repeat("\b \b",length(s)))
          s = string("Minibatch: ",batch)
          print(s)
        end
        #update weights, biases (maybe include another function to make this cleaner)
        #might be a cleaner way to do this (make a function over weights than nn holds?)

        updateWeights!(nn, dW, db, vW, vb,alpha,mu,m,update_type)

        dW = zero(nn.weights)
        db = zero(nn.biases)

        #reset dW, db
        batch += 1

        if mod(batch,half_rate) == 0
          alpha /= 2
        end
      end

      count += 1
    end #i_
  end #epoch

  return nn, losses, losses_true, alpha

end

function train(X::Array{Float64,1},Y::Array{Float64,1},nn::NeuralNetwork;
               alpha::Float64=0.01,mu::Float64=0.9,minibatch_size::Int=100,num_epochs::Int=10,update_type::String="momentum",calculate_loss::Bool=true,half_rate::Int=1)
  return train(X',Y',nn,alpha=alpha,mu=mu,minibatch_size=minibatch_size,num_epochs=num_epochs,
               update_type=update_type,calculate_loss=calculate_loss,half_rate=half_rate)

end

function updateWeights!(nn::NeuralNetwork, dW::NNWeights, db::NNWeights, vW::NNWeights, vb::NNWeights,
                       alpha::Float64,mu::Float64,m::Int,update_type::String)

  if nn.reg == "L2" || nn.reg == "Ridge"
    reg = deepcopy(nn.weights)

  elseif nn.reg == "L1" || nn.reg == "Lasso"
    reg = deepcopy(nn.weights)
    for i = 1:L
      reg[i] = sign(reg[i])
      #reg[i][find(reg[i])] = 1
    end
  end
  #might be a neater way to do this
  if lowercase(update_type) == "momentum"
    #velocity update
    vW = mu*vW+(-alpha./m)*((dW) + (nn.c*reg))
    vb = mu*vb + (-alpha./m)*db
    #vW = add!(mult!(vW,mu),mult!(add!(mult!(dW,1./minibatch_size),mult!(reg,nn.c)),-alpha))
    #vb = add!(mult!(vb,mu),mult!(db,-alpha/minibatch_size))

    #"position" update
    nn.weights += vW
    nn.biases += vb
    #add!(nn.weights,vW)
    #add!(nn.biases,vb)
  elseif lowercase(update_type) == "nesterov"
    #nesterov momentum
    vW_ = deepcopy(vW)
    vb_ = deepcopy(vb)

    vW = mu*vW + (-alpha./m)*(dW + nn.c*reg)
    vb = mu*vb + (-alpha./m)*db

    #vW = add!(mult!(vW,mu),mult!(add!(mult!(dW,1./minibatch_size),mult!(reg,nn.c)),-alpha))
    #vb = add!(mult!(vb,mu),mult!(db,-alpha/minibatch_size))

    nn.weights += (-mu)*vW_+(1+mu)*vW
    nn.biases += (-mu)*vb_ + (1+mu)*vb

    #add!(nn.weights,add!(mult!(vW_,-mu),mult!(vW,1+mu)))
    #add!(nn.biases,add!(mult!(vb_,-mu),mult!(vb,1+mu)))

  elseif lowercase(update_type) == "sgd"
    #vanilla SGD
    nn.weights += (-alpha./m)*(dW+nn.c*reg)
    nn.biases += (-alpha./m)*db

    #add!(nn.weights,mult!(add!(mult!(dW,1./minibatch_size),mult!(reg,nn.c)),-alpha))
    #add!(nn.biases,mult!(db,-alpha/minibatch_size))
  elseif lowercase(update_type) == "rmsprop"
    vW = mu*vW + (1-mu)*(dW)^2
    vb = mu*vb + (1-mu)*(db)^2

    nn.weights += (-alpha/m)*dW/sqrt(vW+1e-8)
    nn.biases += (-alpha/m)*db/sqrt(vb+1e-8)

  else
    error(string("No update rule: ",update_type))
  end

end

function feedforward(x::Array{Float64,1},nn::NeuralNetwork)
  L = length(nn.biases)
  a = zero(nn.biases)#Dict{Int,Array{Float64,1}}([i=>zeros(length(nn.biases[i])) for i = 1:L])
  z = zero(nn.biases)#0 #same form as a
  a[0] = deepcopy(x)

  #i think i screwed up something here--last level should just be a linear function approximator i think
  #so for i=L, it would be teh same except a[L][:] = z[L][:]
  #maybe have something that would make it a sigmoid for classification or whatever

  for i = 1:L-1 #L
    z[i][:] = nn.weights[i]*a[i-1]+nn.biases[i]
    a[i][:] = nn.neuron(z[i])
  end
  z[L][:] = nn.weights[L]*a[L-1]+nn.biases[L]
  a[L][:] = nn.output_neuron(z[L])

  return a[L], a,z
end

function backprop(a::NNWeights,z::NNWeights,y::Array{Float64,1},nn::NeuralNetwork)
  L = length(nn.biases)
  dW = zero(nn.weights)#Dict{Int,Array{Float64,2}}([i=>zeros(size(nn.weights[i])) for i = 1:L])
  db = zero(nn.biases)#Dict{Int,Array{Float64,2}}([i=>zeros(size(nn.biases[i])) for i = 1:L])

  #indexing error with the weights--check it over
  #since im screwing up the things, I think the first gradient term should just be 1 since its just a linear classifier->f(Wx+b)=Wx+b
  #stil screwing things out, need to pen and paper fix things up

  #note: can save time if using sigm by letting the gradient be a.*(1-a) rather than sigm(z).*(1-sigm(z))
  del = zero(nn.biases)#Dict{Int,Array{Float64,2}}([i=>zeros(size(1)) for i = 1:L])
  #del[L] = -(y-a[L]).*nn.output_gradient(z[L])#nn.gradient(z[L])#ones(size(z[L]))#squared error
  del[L] = -sign(y-a[L]).*nn.output_gradient(z[L]) #L1 error
  for l = L:-1:1
    #=
    println("Layer $l")
    println(size(nn.weights[l]))
    println(size(del[l]))
    println(size(a[l]))
    println(size(z[l]))
    println(size(nn.gradient(z[l])))
    println(size(nn.weights[l]'*del[l]))
    println(size((nn.weights[l]'*del[l]).*nn.gradient(z[l])))
    println(size(del[l]*a[l]'))
    =#
    if l > 1
      #println(size(del[l-1]))
      del[l-1] = (nn.weights[l]'*del[l]).*nn.gradient(z[l-1])
    end
    dW[l] = del[l]*a[l-1]'
    db[l] = del[l]
  end

  return dW, db

end

function backprop(a::NNWeights,z::NNWeights,y::Float64,nn::NeuralNetwork)
  return backprop(a,z,[y],nn)
end

function predict(nn::NeuralNetwork, x::Float64)
  return predict(nn,[x])
end

function predict(nn::NeuralNetwork, x::Array{Float64,1})
  y, a, z = feedforward((x-nn.xu)./nn.xs,nn)
  return y#.*((nn.yu-nn.yl)/2)+((nn.yl+nn.yu)/2)
end

function predict(nn::NeuralNetwork, X::Array{Float64,2})
  L = length(nn.biases)
  m = size(X,2)
  Y = zeros(length(nn.biases[L]),m)
  for i = 1:m
    y, a, z = predict(X[:,i],nn)
    Y[:,i] = y
  end
  return Y
end

end #module


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains type definitions and constructores for those types.

# The Spec talks about Booleans, which are implemented in c as
#   `#DEFINE FALSE = 0`
#   `#DEFINE TRUE = 1`
# Which are 32 bit integers. Hence the c code expects Booleans to take up 32 bits of space
# The Julia codebase so far has implemented Bools (1 bit integers) and uses them as such.
# Swap out Bools for Cints, but make it obvious so it can be refactored later; TODO!
const CintThatIsUsedAsABoolean = Cint

"Allow users to configure their own fmi2Component2Environment"
abstract type AbstractFMI2ComponentEnvironment <: Function end

"""
Declaration of FMI2 Types
"""

# FMI2 Errors
struct FMI2Warning <: Exception
    msg::AbstractString
end
function Base.showerror(io::IO, e::FMI2Warning)
    print(io, "fmi2Warning", isempty(e.msg) ? "" : ": ", e.msg)
end

struct FMI2Discard <: Exception
    msg::AbstractString
end
function Base.showerror(io::IO, e::FMI2Discard)
    print(io, "fmi2Discard", isempty(e.msg) ? "" : ": ", e.msg)
end

struct FMI2Error <: Exception
    msg::AbstractString
end
function Base.showerror(io::IO, e::FMI2Error)
    print(io, "fmi2Error", isempty(e.msg) ? "" : ": ", e.msg)
end

struct FMI2Fatal <: Exception
    msg::AbstractString
end
function Base.showerror(io::IO, e::FMI2Fatal)
    print(io, "fmi2Fatal", isempty(e.msg) ? "" : ": ", e.msg)
end

struct FMI2Pending <: Exception
    msg::AbstractString
end
function Base.showerror(io::IO, e::FMI2Pending)
    print(io, "fmi2Pending", isempty(e.msg) ? "" : ": ", e.msg)
end

function fmiError(fmi2Status::Union{Unsigned, Integer}, message::String="")
    if fmi2Status == 1
        return FMI2Warning(message)
    elseif fmi2Status == 2
        return FMI2Discard(message)
    elseif fmi2Status == 3
        return FMI2Error(message)
    elseif fmi2Status ==  4
        return FMI2Fatal(message)
    elseif fmi2Status == 5
        return FMI2Pending(message)
    else
        return FMI2Error("Unknown error code")
    end
end

mutable struct fmi2ComponentEnvironment <: AbstractFMI2ComponentEnvironment
    logFile::String     # if not empty location of file to write logger messages
                        # defaults to stderr ??
    numWarnings::Int
    numErrors::Int
    numFatals::Int
end
fmi2ComponentEnvironment() = ("", 0, 0, 0) # not sure if this is sensible

mutable struct RealVariable
    value::Float64
    valueReference::UInt
    name::String

    # attributes
    min::Float64
    max::Float64

    # Inner constructors
    RealVariable()=new()
    function RealVariable(value, valueReference, name)
        new(value, valueReference, name, -Inf64, Inf64)
    end
end

mutable struct IntVariable
    value::Int64
    valueReference::UInt
    name::String

    # attributes
    min::Int64
    max::Int64

    # Inner constructors
    IntVariable()=new()
    function IntVariable(value, valueReference, name)
        new(value, valueReference, name, -(2^63 - 1), 2^63 - 1)
    end
end

mutable struct BoolVariable
    value::CintThatIsUsedAsABoolean
    valueReference::UInt
    name::String

    # Inner Constructors
    BoolVariable() = new()
    function BoolVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct StringVariable
    value::String
    valueReference::UInt
    name::String

    # Inner constructors
    StringVariable() = new()
    function StringVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct EnumerationVariable
    value::UInt     # Holds a number representig an enumeration
    valueReference::UInt
    name::String

    # Inner constructors
    EnumerationVariable() = new()
    function EnumerationVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct ModelVariables
    reals::Array{RealVariable,1}
    ints::Array{IntVariable,1}
    bools::Array{BoolVariable,1}
    strings::Array{StringVariable,1}
    enumerations::Array{EnumerationVariable,1}

    function ModelVariables(n_reals, n_ints, n_bools, n_strings, n_enumerations)

        reals = Array{RealVariable}(undef, n_reals)
        ints = Array{IntVariable}(undef, n_ints)
        bools = Array{BoolVariable}(undef, n_bools)
        strings = Array{StringVariable}(undef, n_strings)
        enumerations = Array{EnumerationVariable}(undef, n_enumerations)

        new(reals, ints, bools, strings, enumerations)
    end
end

mutable struct SimulationData
    time::AbstractFloat
    lastStepTime::AbstractFloat
    modelVariables::ModelVariables
    eventIndicators::Array{Float64}

    SimulationData()=new()
    function SimulationData(n_reals, n_ints, n_bools, n_strings, n_enumerations,
        numberOfEventIndicators)

        modelVariables=ModelVariables(n_reals, n_ints, n_bools, n_strings,
            n_enumerations)
        eventIndicators = Array{Float64}(undef, numberOfEventIndicators)
        new(0, 0, modelVariables, eventIndicators)
    end
end

mutable struct ModelData
    numberOfStates::Int
    numberOfDerivatives::Int    # Is always numberOfStates

    numberOfReals::Int
    numberOfInts::Int
    numberOfBools::Int
    numberOfStrings::Int
    numberOfEnumerations::Int
    numberOfExterns::Int

    numberOfEventIndicators::Int

    ModelData() = new(0,0,0,0,0,0,0,0)
end

mutable struct ExperimentData
    startTime::AbstractFloat
    stopTime::AbstractFloat
    tolerance::AbstractFloat
    stepSize::AbstractFloat

    # Inner constructors
    ExperimentData() = new(0, 1, 1e-6, 1e-6/4)
    function ExperimentData(startTime, stopTime, tolerance, stepSize)
        if startTime >= stopTime
            error("ExperimentData not valid: startTime=$startTime >= stopTime=$stopTime")
        elseif tolerance <= 0
            error("ExperimentData not valid: tolerance=$tolerance not greater than zero")
        elseif stepSize <= 0
            error("ExperimentData not valid: stepSize=$stepSize not greater than zero")
        end
        new(startTime, stopTime, tolerance, stepSize)
    end
end

struct RealAttributes
    quantity::String
    unit::String            # TODO: make types for Units and functions
    displayUnit::String     #       for unit conversion
    relativeQuantity::CintThatIsUsedAsABoolean
    min::Real
    max::Real
    nominal::Real
    unbound::CintThatIsUsedAsABoolean

    # Inner constructor
    RealAttributes() = new()

    function RealAttributes(quantity, unit, displayUnit, relativeQuantity,
        min, max, nominal, unbound)

        if nominal <= 0
            error("Nominal > 0.0 required")
        elseif min > max
            error("Minimum is greater than maximum")
        end

        new(quantity, unit, displayUnit, relativeQuantity, min, max, nominal,
            unbound)
    end

    function RealAttributes(quantity, unit, displayUnit, min, max, nominal)

        if nominal <= 0
            error("Nominal > 0.0 required")
        elseif min > max
            error("Minimum is greater than maximum")
        end

        new(quantity, unit, displayUnit, false, min, max, nominal, false)
    end
end

struct IntegerAttributes
    quantity::String
    min::Int
    max::Int

    # Inner constructors
    IntegerAttributes()=new()
    IntegerAttributes(quantity, min, max)=new(quantity, min, max)
end

"""
Item of an enumeartion attribute
"""
struct EnumerationItemAttribute
    name::String
    "Unique number in the same enumeration"
    value::Int
    description::String

    # Constructor
    EnumerationItemAttribute() = new()
    EnumerationItemAttribute(name, value, description) = new(name, value, description)
end

"""
Attributes for simpleType enumeration
"""
struct EnumerationAttributes
    "Physical quantity of the variable"
    quantity::String
    "Items of an enumeration"
    items::Array{EnumerationItemAttribute}

    # Constructor
    EnumerationAttributes() = new()
    EnumerationAttributes(quantity, items) = new(quantity, items)
end

struct RealProperties
    declaredType::String
    variableAttributes::RealAttributes
    start::Float64
    derivative::UInt
    reinit::CintThatIsUsedAsABoolean

    # Inner constructors
    RealProperties() = new()
    RealProperties(declaredType, variableAttributes, start, derivative, reinit) = new(declaredType, variableAttributes, start, derivative, reinit)
    function RealProperties(declaredType, variableAttributes, start, derivative, reinit_in::String)
        if reinit_in == "true"
            reinit = true
        elseif reinit_in == "false"
            reinit = false
        else
            error("Could not parste input reinit=$reinit_in to Bool.")
        end

        return new(declaredType, variableAttributes, start, derivative, reinit)
    end
end

struct IntegerProperties
    declaredType::String
    variableAttributes::IntegerAttributes
    start::Int

    # Inner constructors
    IntegerProperies() = new()
    IntegerProperties(declaredType, variableAttributes, start) = new(declaredType, variableAttributes, start)
end

struct BooleanProperties
    declaredType::String
    start::CintThatIsUsedAsABoolean

    BooleanProperties() = new()
    BooleanProperties(declaredType, start) = new(declaredType, start)
end

struct StringProperties
    declaredType::String
    start::String

    StringProperties() = new()
    StringProperties(declaredType, start) = new(declaredType, start)
end

struct EnumerationProperties
    declaredType::String
    quantity::String
    min::Int
    max::Int
    start::Int

    EnumerationProperties() = new()
    EnumerationProperties(declaredType, quantity, min, max, start) = new(declaredType, quantity, min, max, start)
end

struct ScalarVariable
    name::String
    valueReference::Unsigned

    # Optional
    description::String
    causality::String           # TODO: Change to enumeration??
    variability::String         # TODO: Change to enumeration??
    initial::String             # TODO: Change to enumeration??
    canHandleMultipleSetPerTimelnstant::CintThatIsUsedAsABoolean

    # Type specific properties of ScalarVariable
    typeSpecificProperties::Union{RealProperties, IntegerProperties, BooleanProperties, StringProperties, EnumerationProperties}

    # Inner constructors
    function ScalarVariable(name, valueReference, description, typeSpecificProperties)
        if isempty(strip(name))
            error("ScalarVariable $name not valid: name can't be empty or only whitespace")
        elseif valueReference < 0
            error("ScalarVariable $name not valid: valueReference=$valueReference not unsigned")
        end
        new(name, Unsigned(valueReference), description, "local", "continous", "", false, typeSpecificProperties)
    end

    function ScalarVariable(name, valueReference, description, causality,
        variability, initial, canHandleMultipleSetPerTimelnstant, typeSpecificProperties)

        # Check name
        if isempty(strip(name))
            error("ScalarVariable $name not valid: name can't be empty or only whitespace")
        end

        # Check valueReference
        if valueReference < 0
            error("ScalarVariable $name not valid: valueReference=$valueReference not unsigned")
        end

        #check if causality and variability are correct
        if isempty(variability)
            variability = "continuous"
        elseif !in(variability, ["constant", "fixed", "tunable", "discrete", "continuous"])
            error("ScalarVariable $name not valid: variability has to be one of \"constant\", \"fixed\", \"tunable\", \"discrete\" or \"continuous\" but is \"$variability\"")
        end

        if isempty(causality)
            causality = "local"
        elseif !in(causality,["parameter", "calculatedParameter","input", "output", "local", "independent"])
          error("ScalarVariable $name not valid: causality has to be one of \"parameter\", \"calculatedParameter\", \"input\", \"output\", \"local\", \"independent\" but is \"$causality\"")
        elseif causality == "parameter"
            if (variability!="fixed" && variability!="tunable" && variability!="constant")
                error("ScalarVariable $name not valid: causality is \"parameter\", so variability has to be \"fixed\" or \"tunable\" but is \"$variability\"")
            end
            if isempty(initial)
                initial = "exact"
            elseif initial!="exact"
                error("ScalarVariable $name not valid: causality is \"parameter\", so initial has to be \"exact\" or empty but is \"$initial\"")
            end
        elseif causality=="calculatedParameter"
            if (variability!="fixed" && variability!="tunable")
                error("ScalarVariable $name not valid: causality is \"calculatedParameter\", so variability has to be \"fixed\" or \"tunable\" but is \"$variability\"")
            end
            if isempty(initial)
                initial = "calculated"
            elseif !in(initial, ["approx", "calculated"])
                error("ScalarVariable $name not valid: causality is \"calculatedParameter\", so initial has to be \"approx\", \"calculated\" or empty but is \"$initial\"")
            end
        elseif causality=="input"
            if !isempty(initial)
                error("ScalarVariable $name not valid: causality is \"input\", so initial has to be empty but is \"$initial\"")
            end
        elseif causality=="independent"
            if variability!="continuous"
                error("ScalarVariable $name not valid: causality is \"independent\", so variability has to be \"continuous\" but is \"$causality\"")
            end
        end

        new(name, Unsigned(valueReference), description, causality, variability, initial, canHandleMultipleSetPerTimelnstant, typeSpecificProperties)
    end
end

"""
Definition of a simple type according to fmi2SimpleType
"""
struct SimpleType
    "Unique name with respect to all other element of TypeDefinitions in ModelDescription"
    name::String
    description::String

    # ToDo: Add Bool and String to this, maybe as empty Attributes?
    attributes::Union{RealAttributes, IntegerAttributes, EnumerationAttributes}

    # Constructor
    SimpleType() = new()
    SimpleType(name, description, attributes)=new(name, description, attributes)
end

struct LogCategory
    name::String
    description::String

    LogCategory(name) = new(name, "")
    LogCategory(name, description) = new(name, description)
end

"""
Containing all informations from modelDescription.xml
"""
mutable struct ModelDescription
    # FMI model description
    fmiVersion::String
    modelName::String
    guid::String
    description::String
    author::String
    version::String
    copyright::String
    license::String
    generationTool::String
    generationDateAndTime::String
    variableNamingConvention::NamingConvention
    numberOfEventIndicators::Int

    # Model exchange
    isModelExchange::CintThatIsUsedAsABoolean
    # Co-Simulation
    isCoSimulation::CintThatIsUsedAsABoolean
    modelIdentifier::String

    # Unit definitions
    # Type definitions
    typeDefinitions::Union{Array{SimpleType}, Nothing}

    logCategories::Union{Array{LogCategory}, Nothing}

    # Default experiment
    defaultExperiment::ExperimentData

    # Vendor annotations

    # Model variables
    modelVariables::Array{ScalarVariable,1}

    # Model structure
    modelStructure

    # Constructor for uninitialized struct
    function ModelDescription()
        md = new()
        md.isModelExchange = false
        md.isCoSimulation = false
        return md
    end
end

"""
Functional Mockupt Unit (FMU) struct.
"""
mutable struct FMU
    modelName::String
    instanceName::String
    FMUPath::String                     # TODO: find better type for paths
    fmuResourceLocation::String         # is URI
    fmuGUID::String

    modelDescription::ModelDescription

    modelData::ModelData

    simulationData::SimulationData

    experimentData::ExperimentData

    eventInfo::EventInfo

    modelState::ModelState

    # CSV and log file
    csvFile::IOStream
    logFile::IOStream

    # Other stuff
    libHandle::Ptr{Nothing}
    libLoggerHandle::Ptr{Nothing}
    tmpFolder::String
    fmuType::FMUType
    fmi2CallbackFunctions::CallbackFunctions
    fmi2Component::Ptr{Nothing}

    # Constructor
    FMU() = new()
end

max_stars_repo_name: AnHeuermann/JuliaFMU

max_issues_repo_name: AnHeuermann/JuliaFMU

max_forks_repo_name: jwscook/JuliaFMI


text:
module TestExpansionErrors

using ContextualMacros
using Test

@contextualmacro @f

macro expansion_error(ex)
    quote
        err = try
            $Base.@eval $__module__ $ex
            nothing
        catch _err
            _err
        end
        $Test.@test err !== nothing
        err
    end
end

msgof(err) = sprint(showerror, err)

function with_fake_context(f)
    orig = copy(ContextualMacros.DEFINITION_CONSUMERS)
    try
        ContextualMacros.register("@fake_context", :f)
        f()
    finally
        copyto!(empty!(ContextualMacros.DEFINITION_CONSUMERS), orig)
    end
end

function with_no_contexts(f)
    with_fake_context() do
        ContextualMacros._cleanup()
        f()
    end
end

with_fake_context() do
    @testset "NotInContextError (with `@fake_context`)" begin
        err = @expansion_error @f
        @test occursin("Macro `@f` used outside context.", msgof(err))
        @test occursin("@fake_context", msgof(err))
    end
end

with_no_contexts() do
    @testset "NotInContextError (no contexts)" begin
        err = @expansion_error @f
        @test occursin("Macro `@f` used outside context.", msgof(err))
        @test occursin("No package implements `@f`.", msgof(err))
    end
end

end  # module


text:
const NULLMOVE_REDUCTION=2.0
const MAXTHINKINGTIME=(1000000000*14) # n sec

const TT_EXACT = 1
const TT_ALPHA = 2
const TT_BETA  = 3
const TT_INVALID = -1

struct TTKEY
    square::Array{Int,1} # incrementally updated, this array is usefull if we want to
    WhitePiecesInHands::Array{Int,1}
    BlackPiecesInHands::Array{Int,1}
    nextMove::Int # SENTE or GOTE
end

function tt_probe(depth::Float64, alpha::Int, beta::Int, gs::GameStatus)
    val::Int = TT_INVALID
    best::Move = Move(0,0,0,0,0,0)
    key = TTKEY(gs.board.square, gs.board.WhitePiecesInHands, gs.board.BlackPiecesInHands, gs.board.nextMove)
    contents = get( gs.tt, hash(key), -1)
    fl::Int = TT_INVALID

    if contents != -1
        # match!
        best = contents.best

        if contents.depth >= depth
            if contents.flags == TT_EXACT
                val = contents.val
                fl = TT_EXACT
            end
            if (contents.flags == TT_ALPHA)&&(contents.val <= alpha)
                val = alpha
                fl = TT_ALPHA
            end
            if (contents.flags == TT_BETA)&&(contents.val >= beta)
                val = beta
                fl = TT_BETA
            end
        end
    end
    return val, best, fl
end

function tt_save(depth::Float64, val::Int, flags::Int, best::Move, gs::GameStatus)
    key = TTKEY(gs.board.square, gs.board.WhitePiecesInHands, gs.board.BlackPiecesInHands, gs.board.nextMove)
    contents = get( gs.tt, hash(key), -1)
    if contents != -1
        if contents.depth > depth
            return
        end
    end
    newcont = TransP(hash(key),best,depth,flags,val)
    gs.tt[hash(key)] = newcont
end

#@iprofile begin
function qsearch( gs::GameStatus, ply::Int, alpha::Int, beta::Int)
    # quiescence search
    val::Int = 0
    val2::Int = 0
    teban::Int = ((ply & 0x1) == 0) ? gs.side : (gs.side$1)
    ev = (teban == SENTE) ? 1 : -1
    bestValue = -Infinity

    if gs.timedout
        return 0
    end
    gs.triangularLength[ply+1] = ply
    if in_check( teban, gs.board)
        return PVS( gs, ply, 1.0, alpha, beta,true)
    end
    val = ev * EvalBonanza( SENTE, gs.board, gs)

    if val >= beta
        return val
    end
    if val > alpha
        alpha = val
        bestValue = val
    end
    if (gs.inodes) & 1023 == 0
        now = time_ns()
        period::Int = now - gs.nsStart
        if period > gs.maxThinkingTime #MAXTHINKINGTIME
            gs.timedout = true
            return val
        end
    end
    if float(ply) > (gs.depth + 4.0) # changed!
        return alpha
    end

    gs.MoveBeginIndex = gs.moveBufLen[ply+1]
    gs.moveBufLen[ply+1+1] = generateQBB(gs.board, gs.moveBuf, teban, gs.moveBufLen[ply+1], gs)
    for i = gs.moveBufLen[ply+1]+1:gs.moveBufLen[ply+2]
        makeMove( gs.board, i, gs.moveBuf, teban)
        oote = in_check( teban, gs.board)
        if oote||(!oote && seeMoveFlag(gs.moveBuf[i]) == 0)
            #println("check!")
            takeBack( gs.board, i, gs.moveBuf, teban)
        else

	    gs.inodes += 1
	    val = -qsearch( gs, ply+1, -beta, -alpha)
            takeBack( gs.board, i, gs.moveBuf, teban)
            if bestValue < val
                # both sides want to maximize from *their* perspective
                gs.pvmovesfound += 1

                gs.triangularArray[ply+1,ply+1] = gs.moveBuf[i]
		# save this move
                for j = (ply +1+1):gs.triangularLength[ply+1+1]
                    gs.triangularArray[ply+1,j] = gs.triangularArray[ply+1+1,j]
                    # and append the latest best PV from deeper plies
                end
                gs.triangularLength[ply+1] = gs.triangularLength[ply + 1+1]
                rememberPV(gs)
            end
            bestValue=max(val,bestValue)
            alpha = max( alpha, val)
            if alpha >= beta
                return beta
            end
        end
    end
    bestValue
end

#end # @iprofile

function selectHash( gs::GameStatus, ply::Int, i::Int, depth::Float64, best::Move)
    # re-orders the move list so that the best move is selected as the next move to try.
    j::Int = 0
    k::Int = 0
    temp::Move = Move(0,0,0,0,0,0)

    for j = i:gs.moveBufLen[ply+1+1]
        if best.move == gs.moveBuf[j].move
            temp = gs.moveBuf[j]
	    gs.moveBuf[j] = gs.moveBuf[i]
	    gs.moveBuf[i] = temp
	    return
        end
    end
end

function selectmove( gs::GameStatus, ply::Int, i::Int, depth::Float64)
    # re-orders the move list so that the best move is selected as the next move to try.
    best::Int = 0
    j::Int = 0
    k::Int = 0
    temp::Move = Move(0,0,0,0,0,0)

    if gs.followpv && depth > 1.0
        for j = i:gs.moveBufLen[ply+1+1]
            if gs.moveBuf[j].move == gs.lastPV[ply+1].move
                temp = gs.moveBuf[j]
		gs.moveBuf[j] = gs.moveBuf[i]
		gs.moveBuf[i] = temp
		return
            end
        end
    end

    if gs.board.nextMove == GOTE
	best = gs.blackHeuristics[seeMoveFrom(gs.moveBuf[i])+1,seeMoveTo(gs.moveBuf[i])+1]
	j = i
	for k = (i + 1):gs.moveBufLen[ply+1+1]
            if gs.blackHeuristics[seeMoveFrom(gs.moveBuf[k])+1,seeMoveTo(gs.moveBuf[k])+1] > best
	        best = gs.blackHeuristics[seeMoveFrom(gs.moveBuf[k])+1,seeMoveTo(gs.moveBuf[k])+1]
                j = k
            end
        end
	if j > i
            temp = gs.moveBuf[j]
	    gs.moveBuf[j] = gs.moveBuf[i]
	    gs.moveBuf[i] = temp
        end
    else
	best = gs.whiteHeuristics[seeMoveFrom(gs.moveBuf[i])+1,seeMoveTo(gs.moveBuf[i])+1]
	j = i
	for k = (i + 1):gs.moveBufLen[ply+1+1]
            if gs.whiteHeuristics[seeMoveFrom(gs.moveBuf[k])+1,seeMoveTo(gs.moveBuf[k])+1] > best
	        best = gs.whiteHeuristics[seeMoveFrom(gs.moveBuf[k])+1,seeMoveTo(gs.moveBuf[k])+1]
                j = k
            end
        end
	if j > i
            temp = gs.moveBuf[j]
	    gs.moveBuf[j] = gs.moveBuf[i]
	    gs.moveBuf[i] = temp
        end
    end
end

#@iprofile begin
function PVS( gs::GameStatus, ply::Int, depth::Float64, alpha::Int, beta::Int, isID::Bool)
    movesfound::Int = 0
    gs.pvmovesfound = 0
    bestValue = -Infinity
    teban::Int = ((ply & 0x1) == 0) ? gs.side : gs.side $ 1
    val::Int = 0
    tt_flag::Int = TT_ALPHA
    tt_val::Int = 0
    tt_bestMove = Move(0,0,0,0,0,0)
    gs.triangularLength[ply+1] = ply
    if (depth <= 0.0)
	gs.followpv = false
        va::Int = qsearch(gs, ply, alpha, beta)
	return va
    end

    if (gs.inodes) & 1023 == 0
        now = time_ns()
        period::Int = now - gs.nsStart
        if period > gs.maxThinkingTime #MAXTHINKINGTIME
            gs.timedout = true
            return 0
        end
    end

    val, best, fl = tt_probe( depth, alpha, beta, gs)
    #if val != TT_INVALID
    #    return val
    #end
    # if fl == TT_EXACT
    #     return val
    # elseif fl == TT_ALPHA
    #     alpha = max(val,alpha)
    #     #return alpha
    # elseif fl == TT_BETA
    #     beta = min(beta,val)
    #     return beta
    # end

    if (!gs.followpv) && gs.allownull
        if true
            if in_check( teban, gs.board)
                gs.allownull = false
                gs.inodes += 1
                gs.board.nextMove $= 1
                val = -PVS( gs, ply, depth - NULLMOVE_REDUCTION, -beta, -beta+1,isID)
                gs.board.nextMove $= 1
                if gs.timedout
                    return 0
                end
                gs.allownull = true
                if val >= beta
                    return val
                end
            else
            end
        end
    end

    # IID
    if (!gs.followpv) && gs.allownull
        gs.allownull = false
        gs.inodes += 1
        gs.board.nextMove $= 1
        val = -PVS( gs, ply, depth - NULLMOVE_REDUCTION, -beta, -alpha,isID)
        gs.board.nextMove $= 1
        if gs.timedout
            return 0
        end
        gs.allownull = true
        if val >= beta
            return val
        end
    end
    ev = (teban == SENTE) ? 1 : -1
    tmpeval = ev * EvalBonanza( SENTE, gs.board, gs)

    if (!gs.followpv) && gs.allownull && depth < 3.0 && ((tmpeval - 100.0 * depth) >= beta)
        return tmpeval - 100.0*depth
    end

    gs.allownull = true
    # 0.5手延長うまくいかない。qsearchがPVSをも一回呼ぶため。
    movesfound = 0
    gs.pvmovesfound = 0
    gs.MoveBeginIndex = gs.moveBufLen[ply+1]
    gs.moveBufLen[ply+1+1] = generateBB(gs.board, gs.moveBuf, teban, gs.moveBufLen[ply+1], gs)

    beg = gs.moveBufLen[ply+1]+1
    for i = beg:gs.moveBufLen[ply+2]
        if (i == beg)&&(best.move != 0)
            selectHash( gs, ply, i, depth, best)
        else
            selectmove( gs, ply, i, depth)
        end
        makeMove( gs.board, i, gs.moveBuf, teban)
        ext::Float64 = 0.0
        #if i == beg
        #    ext = 0.5
        #end
        if in_check( teban, gs.board)
            #println("check!")
            takeBack( gs.board, i, gs.moveBuf, teban)
        else
	    gs.inodes += 1
	    movesfound += 1

	    if i > beg
                val = -PVS( gs, ply+1, depth-1.0+ext, -alpha-1, -alpha,isID)
                if (val > alpha) && (val < beta)
                    val = -PVS( gs, ply+1, depth-1.0+ext, -beta, -alpha,isID)
                end
            else
                val = -PVS( gs, ply+1, depth-1.0+ext, -beta, -alpha,isID)
            end
            takeBack( gs.board, i, gs.moveBuf, teban)
            if gs.timedout
                return 0
            end
            if bestValue < val
                # both sides want to maximize from *their* perspective
                gs.pvmovesfound += 1
                gs.triangularArray[ply+1,ply+1] = gs.moveBuf[i]
		# save this move
                for j = (ply +1+1):gs.triangularLength[ply+1+1]
                    gs.triangularArray[ply+1,j] = gs.triangularArray[ply+1+1,j]
                    # and append the latest best PV from deeper plies
                end
                gs.triangularLength[ply+1] = gs.triangularLength[ply + 1+1]
                rememberPV(gs)

                tt_flag = TT_ALPHA
                tt_val = alpha
                tt_bestMove = gs.moveBuf[i]
                tt_save( depth, tt_val, tt_flag, tt_bestMove, gs)
            end
            bestValue=max(val,bestValue)
            alpha = max( alpha, val)
            if alpha >= beta
                if gs.board.nextMove == GOTE
                    gs.blackHeuristics[seeMoveFrom(gs.moveBuf[i])+1,seeMoveTo(gs.moveBuf[i])+1] += Int(depth*depth)
                else
                    gs.whiteHeuristics[seeMoveFrom(gs.moveBuf[i])+1,seeMoveTo(gs.moveBuf[i])+1] += Int(depth*depth)
                end

                tt_flag = TT_BETA
                tt_val = beta
                tt_bestMove = gs.moveBuf[i]
                tt_save( depth, tt_val, tt_flag, gs.moveBuf[i], gs)
                return beta
            end
        end
    end

    if gs.pvmovesfound > 0
        if gs.board.nextMove == GOTE
            gs.blackHeuristics[seeMoveFrom(gs.triangularArray[ply+1,ply+1])+1,seeMoveTo(gs.triangularArray[ply+1,ply+1])+1] += Int(depth*depth)
        else
            gs.whiteHeuristics[seeMoveFrom(gs.triangularArray[ply+1,ply+1])+1,seeMoveTo(gs.triangularArray[ply+1,ply+1])+1] += Int(depth*depth)
        end
        tt_flag = TT_EXACT
        tt_save( depth, alpha, tt_flag, tt_bestMove, gs)
    end
    if movesfound == 0
	if in_check(teban, gs.board)
            #tt_save( depth, -Infinity+ply-1, tt_flag, tt_bestMove, gs)
            return -Infinity+ply-1
        end
    end
    bestValue
end

#end # @iprofile

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
# # Recursive net on IMDB sentiment treebank

# In this example, we create a recursive neural network to perform sentiment analysis using
# IMDB data.
# This type of model can be used
# for learning tree-like structures (directed acyclic graphs).
# It computes compositional vector representations for prhases of variable length
# which are used as features for performing classification.

# ![treebank](../treebank/docs/treebank.png)

# [Source](https://nlp.stanford.edu/~socherr/EMNLP2013_RNTN.pdf)

# This example uses the [Standford Sentiment Treebank dataset
# (SST)](https://nlp.stanford.edu/sentiment/index.html) which is often used
# as one of the benchmark datasets to test new language models.
# It has five different classes (very negative to very positive) and the
# goal is to perform sentiment analysis.

# To run this example, we need the following packages:

using Flux
using Flux: logitcrossentropy, throttle
using Flux.Data: Tree, children, isleaf
using Parameters: @with_kw

# The script `data.jl` contains the function `getdata` that obtains
# and process the SST dataset.

include("data.jl")

# We set default values for the hyperparameters:

@with_kw mutable struct Args
    lr::Float64 = 1e-3    ## Learning rate
    N::Int = 300
    throttle::Int = 10    ## Throttle timeout
end

# ## Build the model

# The function `train` loads the data, builds and trains the model.
# For more information on how the recursive neural network works, see
# section 4 of [Recursive Deep Models for Semantic Compositionality
# Over a Sentiment Treebank](https://nlp.stanford.edu/~socherr/EMNLP2013_RNTN.pdf).

function train(; kws...)
    ## Initialize HyperParameters
    args = Args(; kws...)

    ## Load data
    @info("Loading Data...")
    train_data, alphabet = getdata()

    @info("Constructing model....")
    embedding = randn(Float32, args.N, length(alphabet))

    @info "Size of the embedding" size(embedding)

    W = Dense(2*args.N, args.N, tanh)
    combine(a, b) = W([a; b])

    sentiment = Chain(Dense(args.N, 5))

    function forward(tree)
      if isleaf(tree)
        token, sent = tree.value
        phrase = embedding * token
        phrase, logitcrossentropy(sentiment(phrase), sent)
      else
        _, sent = tree.value
        c1, l1 = forward(tree[1])
        c2, l2 = forward(tree[2])
        phrase = combine(c1, c2)
        phrase, l1 + l2 + logitcrossentropy(sentiment(phrase), sent)
      end
    end

    loss(tree) = forward(tree)[2]

    opt = ADAM(args.lr)
    ps = params(embedding, W, sentiment)
    evalcb = () -> @show loss(train_data[1])
    @info("Training Model...")
    Flux.train!(loss, ps, zip(train_data), opt,cb = throttle(evalcb, args.throttle))
end

# ## Train the model

cd(@__DIR__)
train()


text:
mutable struct LatentFactorModelPrecisionOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    lfm::LatentFactorModelXMLElement
    weight::Real

    function LatentFactorModelPrecisionOperatorXMLElement(
                        lfm::LatentFactorModelXMLElement)
        return new(nothing, lfm, 1.0)
    end
end

function make_xml(op::LatentFactorModelPrecisionOperatorXMLElement)
    el = new_element(bn.LAT_FAC_PREC_OP)
    set_attribute(el, bn.WEIGHT, op.weight)
    set_attribute(el, bn.RANDOM_SCAN, bn.FALSE)

    make_xml(op.lfm)
    add_ref_el(el, op.lfm.el)
    add_ref_el(el, op.lfm.precision_prior_el)
    op.el = el
    return el
end


text:
disabled() = get(ENV, "R123_DISABLE_AESNI", "") != ""

const __m128i = NTuple{2, VecElement{UInt64}}
has_aesni() = try
    ccall(
        "llvm.x86.aesni.aeskeygenassist", llvmcall, __m128i, (__m128i, UInt8),
        __m128i((0x0123456789123450, 0x9876543210987654)), 0x1
    ) ≡ __m128i((0x857c266f7c266e85, 0x2346382146382023))
catch e
    false
end

const filename = joinpath(dirname(@__FILE__), "aes-ni")
isfile(filename) && rm(filename)

if disabled()
    @info "AES-NI is disabled."
elseif has_aesni()
    @info "AES-NI is enabled."
    touch(filename)
else
    @info "AES-NI is not supported."
end