text:
#= ######################################################
includes
###################################################### =#

include("./models2D.jl")

#= ######################################################
mesh2D - testing
###################################################### =#

function test_mesh2D()
    # which model
    # m = model_two_intesecting_circles(3.0, 120, 24)
    # m = model_rectangle(2.0, 4.0, 0.1)
    m = model_circle_in_circle(0.8, 1.6, 0.05)
    # m = model_circle_with_opening_line(3.0, 120, 15)
    # m = model_isosceles_triangle(2.0, 90, 0.2)
    # m = model_two_rectangles_with_holes(0.02)
    # m = model_lines()
    # m = model_cosinus(0.0, 1.0, 100)
    # m = model_trapez(2.0, 3.0, 10.0, 0.1)
    # m = model_right_triangle(3.0, 2.0, 0.05)
    # m = model_furnace1()
    # m = model_furnace2(0.195,0.5,3,0.7)
    # m = model_labyrinth([0.4, -0.4, 0.4, 0.2, 0.5, -0.3, -0.2], 0.1, 0.02)
    # m = model_square_in_square(0.8, 1.6, 0.01)
    # m = model_rect_in_rect(1.0, 3.0, 2.0, 4.0, 0.03)
    # m = model_circles_in_circle_cross(0.1, 1.8, 9, 0.01)
    # m = model_circles_in_circle_rand(0.1, 1.8, 17, 0.01)
    # m = model_circles_in_circle_rand_full(0.1, 1.8, 0.01)
    # m = model_circles_in_circle_rand_half(0.1, 1.8, 0.01)
    # m = model_circles_in_circle_rand_quarter(0.1, 1.8, 0.01)
    # element analysis
    element_analysis(m, printit = true)
    # 2D plot
    fig = Figure(resolution = (900, 900))
    ax = fig[1, 1] = Axis(fig)
    # ax.xlabelpadding = 0.0
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    plot_model(fig, ax, m, shownvec = true, shownodes = false, showcom = false)
    display(fig)
end

#= ######################################################
raycast2D - testing
###################################################### =#

function test_raycast2D()
    # 2D plot
    fig = Figure(resolution = (900, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # create tiles
    n = 10
    dx = 0.5
    dy = 0.5
    plot_empty_tiles(fig, ax, dx, dy, n)
    # create two points by random
    p1 = Point2D(rand()*n*dx,rand()*n*dy)
    p2 = Point2D(rand()*n*dx,rand()*n*dy)
    # p1 = Point2D(1*n*dx,rand()*n*dy)
    # p2 = Point2D(1*n*dx,rand()*n*dy)
    # p1 = Point2D(0.2,0.1)
    # p2 = Point2D(4.8,4.9)
    plot_points_and_connection(fig, ax, p1, p2)
    # tile walk - used for blocking calc
    max_steps = get_max_steps(n)
    tiles = Vector{Index2D{Int64}}(undef,max_steps)
    @time ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
    for i = 1:ntiles
        t = tiles[i]
        # print tile in color
        poly!(ax, Rect((t.x-1) * dx, (t.y-1) * dy, dx, dy), color = (:blue, 0.3))
    end
    display(fig)
end

function test_raycast2D_occ(;option = 3)
    # not correctly working because of missing Makie pkg in RadMod2D
    # 2D plot
    fig = Figure(resolution = (900, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # create tiles
    n = 10
    dx = 0.5
    dy = 0.5
    plot_empty_tiles(fig, ax, dx, dy, n)
    t_occ = create_randomly_occupied_tiles(n, density = 0.2)
    plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    # create two points by random
    p1 = Point2D(rand()*n*dx,rand()*n*dy)
    p2 = Point2D(rand()*n*dx,rand()*n*dy)
    # p1 = Point2D(1*n*dx,rand()*n*dy)
    # p2 = Point2D(1*n*dx,rand()*n*dy)
    # p1 = Point2D(0.2,0.1)
    # p2 = Point2D(4.8,4.9)
    plot_points_and_connection(fig, ax, p1, p2)
    # tile walk options
    if option == 1
        @time tilewalk_with_check_for_occ_tiles(fig, ax, p1, p2, dx, dy, n, t_occ)
    elseif option == 2
        @time tiles = tilewalk_with_return(fig, ax, p1, p2, dx, dy, n)
        for t in tiles
            println("tile number: ", t.x, ", ", t.y)
        end
    elseif option == 3
        # used for blocking calc
        max_steps = get_max_steps(n)
        tiles = Vector{Index2D{Int64}}(undef,max_steps)
        @time ntiles = tilewalk_with_return!(tiles, p1, p2, dx, dy, n)
        for i = 1:ntiles
            t = tiles[i]
            # print tile in color
            poly!(ax, Rect((t.x-1) * dx, (t.y-1) * dy, dx, dy), color = (:blue, 0.3))
        end
    else
        println("Option not defined...")
    end
    display(fig)
end

#= ######################################################
view2D - testing
###################################################### =#

function test_view2D_blocking_2elem()
    # test double open circle case
    m = model_two_intesecting_circles(3.0, 120, 3)
    # m = model_circle_with_opening_line(2.0, 360, 143)
    # 2D plot
    fig = Figure(resolution = (900, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # create tiles
    n = 30
    dx, dy = create_tiles(m, n)
    t_occ = check_tile_occupation(m, dx, dy, n)
    plot_empty_tiles(fig, ax, dx, dy, n)
    plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    plot_model(fig, ax, m, shownvec = false, shownodes = false, showcom = false)
    # elem pair check
    i1 = 50
    i2 = 160
    plot_model_elements(fig, ax, m, [i1,i2], shownvec = true, shownodes = false, showcom = false)
    isexisting = existing_vf_2elem(m, i1, i2)
    if isexisting
        linesegments!(ax, [Point2f0(m.elem[i1].com.x, m.elem[i1].com.y) => Point2f0(m.elem[i2].com.x, m.elem[i2].com.y)], color = :blue, linewidth = 3)
        blocking_vf_with_tiles_2elem(fig, ax, i1, i2, m, dx, dy, n, t_occ)
    end
    display(fig)
end

function test_view2D_blocking_bf_vs_tiles()
    # test blocking: brute force vs tiles
    # m = model_two_intesecting_circles(3.0, 360, 36)
    m = model_two_rectangles_with_holes(0.02)
    # blocking with brute force
    vfmat = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat)
    @time blocking_vf!(m, vfmat)
    # blocking with tilewalk
    vfmat2 = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat2)
    n = 10
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    @time blocking_vf_with_tiles!(m, vfmat2, dx, dy, n, t_occ)
    # check for difference
    diff = vfmat2 .- vfmat
    println(sum(diff))
    # findfirst(diff.==1.0)
    # 2D plot
    # fig = Figure(resolution = (1400, 900))
    # ax = fig[1, 1] = Axis(fig)
    # ax.xlabel = "X in m"
    # ax.ylabel = "Y in m"
    # ax.aspect = DataAspect()
    # plot_empty_tiles(fig, ax, dx, dy, n)
    # plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    # plot_model(fig, ax, m, shownvec = false, shownodes = false, showcom = false)
    # plot_existing(fig, ax, m, 40, vfmat2)
    # display(fig)
end

function test_view2D_vf()
    #### which model
    # m = model_two_intesecting_circles(3.0, 360, 45)
    # m = model_rectangle(2.0, 5.0, 0.1)
    # m = model_circle_in_circle(3.3, 5.2, 0.1)
    # m = model_circle_with_opening_line(2.0, 360, 143)
    # m = model_labyrinth(0.02)
    # m = model_square_in_square(0.8, 1.6, 0.01)
    # m = model_rect_in_rect(0.8, 0.8, 1.6, 1.6, 0.05)
    # m = model_circles_in_circle_rand_full(0.1, 1.8, 0.005)
    # m = model_circles_in_circle_rand_half(0.1, 1.8, 0.005)
    # m = model_circles_in_circle_rand_quarter(0.1, 1.8, 0.005)
    m = model_circles_in_circle_cross(0.1, 1.8, 9, 0.01)
    #### vfmat calculation
    vfmat = zeros(Float64, m.nelem, m.nelem)
    @time existing_vf!(m, vfmat)
    n = 20
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    @time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
    @time calculating_vf!(m, vfmat, normit = false)
    vfmatp = vfmat_to_parts(m, vfmat, normit = false)
    #### 2D plot
    fig = Figure(resolution = (1400, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # plot_empty_tiles(fig, ax, dx, dy, n)
    # plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    plot_model(fig, ax, m, shownvec = false, shownodes = false, showcom = false)
    plot_existing(fig, ax, m, 170, vfmat)
    display(fig)
    #### post processing
    # control = sum(vfmatp, dims=2)
    # save_vfmat_part_to_csv(vfmatp, "figures", "cyl_in_cyl_rand_cross_0n005", control = true)
    # vfmatp
end

function test_view2D_tiles(; n = 15)
    #### which model
    # m = model_two_intesecting_circles(3.0, 360, 45)
    # m = model_rectangle(2.0, 5.0, 0.1)
    # m = model_circle_in_circle(0.8, 1.6, 0.001)
    # m = model_circle_with_opening_line(2.0, 360, 143)
    # m = model_labyrinth(0.02)
    m = model_square_in_square(0.8, 1.6, 0.01)
    #### vfmat calculation
    vfmat = zeros(Float64, m.nelem, m.nelem)
    @time existing_vf!(m, vfmat)
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    #### 2D plot
    fig = Figure(resolution = (1400, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    plot_empty_tiles(fig, ax, dx, dy, n)
    plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    plot_model(fig, ax, m, shownvec = false, shownodes = false, showcom = false)
    display(fig)
end

function test_view2D_shadow()
    # test blocking with shadow plot
    m = model_two_rectangles_with_holes(0.02)
    # blocking with tilewalk
    vfmat = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat)
    n = 30
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    @time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
    # 2D plot
    fig = Figure(resolution = (1400, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # plot_empty_tiles(fig, ax, dx, dy, n)
    # plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
    # plot_model(fig, ax, m, shownvec = true, shownodes = false, showcom = false)
    # plot_existing(fig, ax, m, 150, vfmat)
    # plot_model_shadow_1to1(fig, ax, m, vfmat, 3, 8)
    plot_model_shadow_1toAll(fig, ax, m, vfmat, 3)
    display(fig)
end

function test_tile_occ_analysis(;n = 10)
    m = model_circle_in_circle(0.8, 1.6, 0.05)
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    # tile_occ_analysis(t_occ)
    t_max, t_min, t_mean, ratio = tile_occ_analysis(t_occ, printit = false)
end

#= ######################################################
therm2D - testing
###################################################### =#

function test_therm2D()
    # which model
    # m = model_two_intesecting_circles(3.0, 120, 24)
    m = model_rectangle(2.0, 2.5, 0.1)
    # m = model_circle_in_circle(3.3, 5.2, 0.1)
    # vfmat calculation
    vfmat = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat)
    n = 30
    dx, dy = create_tiles(m, n)
    t_occ = check_tile_occupation(m, dx, dy, n)
    blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
    calculating_vf!(m, vfmat, normit = true)
    # solve Qp
    epsilon = zeros(m.nelem,1)
    set_bc_part!(m, epsilon, 1, 0.9)
    set_bc_part!(m, epsilon, 2:4, 0.3)
    temp = zeros(m.nelem,1)
    set_bc_part!(m, temp, 1, 602)
    set_bc_part!(m, temp, 2, 600)
    set_bc_part!(m, temp, 3:4, 590)
    @time Qp, G = tempsolver(m, vfmat, temp, epsilon)
    Qp_parts = [sum(Qp[m.elem2par[i].first:m.elem2par[i].last,1]) for i = 1:m.npar]
    # 2D plot
    fig = Figure(resolution = (1400, 900))
    ax = fig[1, 1] = Axis(fig)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # plot_model(fig, ax, m, shownvec = true, shownodes = false, showcom = false)
    plot_model_with_value(fig, ax, m, Qp, "Wärmestrom in W", showcbar = false)
    display(fig)
end


text:
using MyMethods
using Test

@testset "myLS.jl" begin
    # Generate example data
    n_i = 1000
    X = 2 * randn((n_i, 2))
    y = X * [1, -2] + randn((n_i, 1))

    # Estimate the least square regression
    ls_fit = myLS(y, X)

    # Check the methods
    β_hat = coef(ls_fit)
    y_hat = predict(ls_fit)
    se_hat = inference(ls_fit, print_df = false)
    se_hc_hat = inference(ls_fit, heteroskedastic = true, print_df = false)

    # Let's check that everything is of correct type.
    @test typeof(ls_fit) == myLS
    @test typeof(β_hat) == Array{Float64,2}
    @test typeof(y_hat) == Array{Float64,2}
    @test typeof(se_hat) == NamedTuple{(:β, :se, :t, :p),
                                       Tuple{Array{Float64,2},
                                             Array{Float64,1},
                                             Array{Float64,2},
                                             Array{Float64,2}}}
    @test typeof(se_hc_hat) == NamedTuple{(:β, :se, :t, :p),
                                          Tuple{Array{Float64,2},
                                                Array{Float64,1},
                                                Array{Float64,2},
                                                Array{Float64,2}}}
end

@testset "myTSLS.jl" begin
    # Generate example data
    n_i = 1000
    ν = randn((n_i, 1))
    instrument = 2 * randn((n_i, 2))
    D = instrument * [1, -2] + ν
    y = D + 0.3 * ν + randn((n_i, 1))

    # Estimate the least square regression
    tsls_fit = myTSLS(y, D, instrument)

    # Check the methods
    β_hat = coef(tsls_fit)
    y_hat = predict(tsls_fit)
    se_hat = inference(tsls_fit, print_df = false)
    se_hc_hat = inference(tsls_fit, heteroskedastic = true, print_df = false)

    # Let's check that everything is of correct type.
    @test typeof(tsls_fit) == myTSLS
    @test typeof(β_hat) == Array{Float64,2}
    @test typeof(y_hat) == Array{Float64,2}
    @test typeof(se_hat) == NamedTuple{(:β, :se, :t, :p),
                                       Tuple{Array{Float64,2},
                                             Array{Float64,1},
                                             Array{Float64,2},
                                             Array{Float64,2}}}
    @test typeof(se_hc_hat) == NamedTuple{(:β, :se, :t, :p),
                                          Tuple{Array{Float64,2},
                                                Array{Float64,1},
                                                Array{Float64,2},
                                                Array{Float64,2}}}
end

@testset "mySieve.jl" begin
      # Generate example data
      n_i = 1000
      x = randn(n_i)
      X = get_basis(x, "Bernstein", 5, nothing)
      y = X * randn(6) + randn((n_i, 1))

      # Estimate the sieve regression
      sieve_fit = mySieve(y, x; basis="Bernstein", K=5)

      # Check the methods
      β_hat = coef(sieve_fit)
      y_hat = predict(sieve_fit)

      # Let's check that everything is of correct type.
      @test typeof(sieve_fit) == mySieve
end

ext: jl


text:
@testset "core" begin

@testset "context" begin

global global_ctx
global_ctx = GlobalContext()
local_ctx = Context()

let
    ctx = Context()
    @assert local_ctx != GlobalContext()
    dispose(ctx)
end

Context() do ctx end

end

@testset "type" begin

Context() do ctx
    typ = LLVM.Int1Type(ctx)
    @test typeof(LLVM.ref(typ)) == LLVM.API.LLVMTypeRef                 # untyped

    @test typeof(LLVM.IntegerType(LLVM.ref(typ))) == LLVM.IntegerType   # type reconstructed
    if Base.JLOptions().debug_level >= 2
        @test_throws ErrorException LLVM.FunctionType(LLVM.ref(typ))    # wrong type
    end
    @test_throws UndefRefError LLVM.FunctionType(LLVM.API.LLVMTypeRef(C_NULL))

    @test typeof(LLVM.ref(typ)) == LLVM.API.LLVMTypeRef
    @test typeof(LLVMType(LLVM.ref(typ))) == LLVM.IntegerType           # type reconstructed
    @test_throws UndefRefError LLVMType(LLVM.API.LLVMTypeRef(C_NULL))

    LLVM.IntType(8)
    @test width(LLVM.IntType(8, ctx)) == 8

    @test issized(LLVM.Int1Type(ctx))
    @test !issized(LLVM.VoidType(ctx))
end

# integer
let
    typ = LLVM.Int1Type()
    @test context(typ) == global_ctx

    show(devnull, typ)

    @test !isempty(typ)

    Context() do ctx
        typ2 = LLVM.Int1Type(ctx)
        @test context(typ2) == ctx
        @test typ2 != typ
    end
end

# floating-point

# function
Context() do ctx
    x = LLVM.Int1Type(ctx)
    y = [LLVM.Int8Type(ctx), LLVM.Int16Type(ctx)]
    ft = LLVM.FunctionType(x, y)
    @test context(ft) == ctx

    @test !isvararg(ft)
    @test return_type(ft) == x
    @test parameters(ft) == y
end

# sequential
Context() do ctx
    eltyp = LLVM.Int32Type(ctx)

    ptrtyp = LLVM.PointerType(eltyp)
    @test eltype(ptrtyp) == eltyp
    @test context(ptrtyp) == context(eltyp)

    @test addrspace(ptrtyp) == 0

    ptrtyp = LLVM.PointerType(eltyp, 1)
    @test addrspace(ptrtyp) == 1
end
Context() do ctx
    eltyp = LLVM.Int32Type(ctx)

    arrtyp = LLVM.ArrayType(eltyp, 2)
    @test eltype(arrtyp) == eltyp
    @test context(arrtyp) == context(eltyp)
    @test !isempty(arrtyp)

    @test length(arrtyp) == 2
end
Context() do ctx
    eltyp = LLVM.Int32Type(ctx)

    arrtyp = LLVM.ArrayType(eltyp, 0)
    @test isempty(arrtyp)
end
Context() do ctx
    eltyp = LLVM.Int32Type(ctx)

    vectyp = LLVM.VectorType(eltyp, 2)
    @test eltype(vectyp) == eltyp
    @test context(vectyp) == context(eltyp)

    @test size(vectyp) == 2
end

# structure
let
    st = LLVM.StructType([LLVM.VoidType()])
    @test context(st) == global_ctx
    @test !isempty(st)

    st2 = LLVM.StructType("foo")
    @test context(st2) == global_ctx
    @test isempty(st2)
end
Context() do ctx
    elem = [LLVM.Int32Type(ctx), LLVM.FloatType(ctx)]

    let st = LLVM.StructType(elem, ctx)
        @test context(st) == ctx
        @test !ispacked(st)
        @test !isopaque(st)

        let elem_it = elements(st)
            @test eltype(elem_it) == LLVMType

            @test length(elem_it) == length(elem)

            @test first(elem_it) == elem[1]
            @test last(elem_it) == elem[end]

            i = 1
            for el in elem_it
                @test el == elem[i]
                i += 1
            end

            @test collect(elem_it) == elem
        end
    end

    let st = LLVM.StructType("foo", ctx)
        @test name(st) == "foo"
        @test isopaque(st)
        elements!(st, elem)
        @test collect(elements(st)) == elem
        @test !isopaque(st)
    end
end

# other
let
    typ = LLVM.VoidType()
    @test context(typ) == global_ctx
end
Context() do ctx
    typ = LLVM.VoidType(ctx)
    @test context(typ) == ctx
end
let
    typ = LLVM.LabelType()
    @test context(typ) == global_ctx
end
Context() do ctx
    typ = LLVM.LabelType(ctx)
    @test context(typ) == ctx
end
Context() do ctx
    typ = LLVM.MetadataType(ctx)
    @test context(typ) == ctx
end
Context() do ctx
    typ = LLVM.TokenType(ctx)
    @test context(typ) == ctx
end

end

@testset "value" begin

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx)])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    entry = BasicBlock(fn, "entry")
    position!(builder, entry)
    @test name(entry) == "entry"

    typ = LLVM.Int32Type(ctx)
    val = alloca!(builder, typ, "foo")
    @test context(val) == ctx
    @test typeof(LLVM.ref(val)) == LLVM.API.LLVMValueRef                # untyped

    @test typeof(LLVM.Instruction(LLVM.ref(val))) == LLVM.AllocaInst    # type reconstructed
    if Base.JLOptions().debug_level >= 2
        @test_throws ErrorException LLVM.Function(LLVM.ref(val))        # wrong
    end
    @test_throws UndefRefError LLVM.Function(LLVM.API.LLVMValueRef(C_NULL))

    @test typeof(Value(LLVM.ref(val))) == LLVM.AllocaInst              # type reconstructed
    @test_throws UndefRefError Value(LLVM.API.LLVMValueRef(C_NULL))

    show(devnull, val)

    @test llvmtype(val) == LLVM.PointerType(typ)
    @test name(val) == "foo"
    @test !isconstant(val)
    @test !isundef(val)

    name!(val, "bar")
    @test name(val) == "bar"
end
end
end

# usage

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx)])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    entry = BasicBlock(fn, "entry")
    position!(builder, entry)

    valueinst1 = add!(builder, parameters(fn)[1],
                      ConstantInt(Int32(1), ctx))
    @test !isterminator(valueinst1)

    userinst = add!(builder, valueinst1,
                    ConstantInt(Int32(1), ctx))

    # use iteration
    let usepairs = uses(valueinst1)
        @test eltype(usepairs) == Use

        usepair = first(usepairs)
        @test value(usepair) == valueinst1
        @test user(usepair) == userinst

        for _usepair in usepairs
            @test usepair == _usepair
        end

        @test value.(collect(usepairs)) == [valueinst1]
        @test user.(collect(usepairs)) == [userinst]
    end

    valueinst2 = add!(builder, parameters(fn)[1],
                    ConstantInt(Int32(2), ctx))

    replace_uses!(valueinst1, valueinst2)
    @test user.(collect(uses(valueinst2))) == [userinst]
end
end
end

# users

Context() do ctx
    # operand iteration
    mod = parse(LLVM.Module,  """
        define void @fun(i32) {
        top:
          %1 = add i32 %0, 1
          ret void
        }""", ctx)
    fun = functions(mod)["fun"]

    for (i, instr) in enumerate(instructions(first(blocks(fun))))
        ops = operands(instr)
        @test eltype(ops) == Value
        if i == 1
            @test length(ops) == 2
            @test ops[1] == first(parameters(fun))
            @test ops[2] == ConstantInt(LLVM.Int32Type(ctx), 1)
        elseif i == 2
            @test length(ops) == 0
            @test collect(ops) == []
        end
    end
end

# constants

Context() do ctx
    @testset "constants" begin

    typ = LLVM.Int32Type(ctx)
    ptrtyp = LLVM.PointerType(typ)

    let val = null(typ)
        @test isnull(val)
    end

    let val = all_ones(typ)
        @test !isnull(val)
    end

    let val = PointerNull(ptrtyp)
        @test isnull(val)
    end

    let val = UndefValue(typ)
        @test isundef(val)
    end

    end
end

# scalar
Context() do ctx
    @testset "integer constants" begin

    # manual construction of small values
    let
        typ = LLVM.Int32Type(ctx)
        constval = ConstantInt(typ, -1)
        @test convert(Int, constval) == -1
        @test convert(UInt32, constval) == typemax(UInt32)
    end

    # manual construction of large values
    let
        typ = LLVM.Int64Type(ctx)
        constval = ConstantInt(typ, BigInt(2)^100-1)
        @test convert(Int, constval) == -1
    end

    # automatic construction
    let
        constval = ConstantInt(UInt32(1))
        @test convert(UInt, constval) == 1
    end

    # issue #81
    for T in [Int32, UInt32, Int64, UInt64]
        constval = ConstantInt(typemax(T))
        @test convert(T, constval) == typemax(T)
    end

    end

    @testset "floating point constants" begin

    t2 = LLVM.DoubleType(ctx)
    c = ConstantFP(t2, 1.1)
    @test convert(Float64, c) == 1.1

    end
end

# constant expressions
Context() do ctx
    @testset "constant expressions" begin

    # inline assembly
    let
        ft = LLVM.FunctionType(LLVM.VoidType(ctx))
        asm = InlineAsm(ft, "nop", "", false)
        @check_ir asm "void ()* asm \"nop\", \"\""
    end

    end
end

# global values
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    st = LLVM.StructType("SomeType", ctx)
    ft = LLVM.FunctionType(st, [st])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    @test isdeclaration(fn)
    @test linkage(fn) == LLVM.API.LLVMExternalLinkage
    linkage!(fn, LLVM.API.LLVMAvailableExternallyLinkage)
    @test linkage(fn) == LLVM.API.LLVMAvailableExternallyLinkage

    @test section(fn) == ""
    section!(fn, "SomeSection")
    @test section(fn) == "SomeSection"

    @test visibility(fn) == LLVM.API.LLVMDefaultVisibility
    visibility!(fn, LLVM.API.LLVMHiddenVisibility)
    @test visibility(fn) == LLVM.API.LLVMHiddenVisibility

    @test dllstorage(fn) == LLVM.API.LLVMDefaultStorageClass
    dllstorage!(fn, LLVM.API.LLVMDLLImportStorageClass)
    @test dllstorage(fn) == LLVM.API.LLVMDLLImportStorageClass

    @test !unnamed_addr(fn)
    unnamed_addr!(fn, true)
    @test unnamed_addr(fn)

    @test alignment(fn) == 0
    alignment!(fn, 4)
    @test alignment(fn) == 4
end
end

# global variables
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    @test isempty(globals(mod))
    gv = GlobalVariable(mod, LLVM.Int32Type(ctx), "SomeGlobal")
    @test !isempty(globals(mod))

    show(devnull, gv)

    @test_throws UndefRefError initializer(gv)
    init = ConstantInt(Int32(0), ctx)
    initializer!(gv, init)
    @test initializer(gv) == init

    @test !isthreadlocal(gv)
    threadlocal!(gv, true)
    @test isthreadlocal(gv)

    @test !isconstant(gv)
    constant!(gv, true)
    @test isconstant(gv)

    @test !isextinit(gv)
    extinit!(gv, true)
    @test isextinit(gv)

    @test threadlocalmode(gv) == LLVM.API.LLVMGeneralDynamicTLSModel
    threadlocalmode!(gv, LLVM.API.LLVMNotThreadLocal)
    @test threadlocalmode(gv) == LLVM.API.LLVMNotThreadLocal

    let gvars = globals(mod)
        @test gv in gvars
        unsafe_delete!(mod, gv)
        @test isempty(gvars)
    end

    # bug: PointerNull wasn't a constant, and couldn't be used as an initializer
    let ptrtyp = LLVM.PointerType(LLVM.VoidType(ctx))
        gv = GlobalVariable(mod, ptrtyp, "SomeOtherGlobal")
        init = PointerNull(ptrtyp)
        initializer!(gv, init)
    end
end
end

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    st = LLVM.StructType("SomeType", ctx)
    gv = GlobalVariable(mod, st, "SomeGlobal")

    init = null(st)
    initializer!(gv, init)
    @test initializer(gv) == init
end
end

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    gv = GlobalVariable(mod, LLVM.Int32Type(ctx), "SomeGlobal", 1)

    @test addrspace(llvmtype(gv)) == 1
end
end

end

@testset "metadata" begin

@test MDString("foo") == MDString("foo", global_ctx)

Context() do ctx
    str = MDString("foo", ctx)
    @test convert(String, str) == "foo"
end

@test MDNode([MDString("foo")]) == MDNode([MDString("foo", global_ctx)], global_ctx)

Context() do ctx
    str = MDString("foo", ctx)
    node = MDNode([str], ctx)
    ops = operands(node)
    @test length(ops) == 1
    @test ops[1] == str
end

end

@testset "module" begin

let
    mod = LLVM.Module("SomeModule")
    @test context(mod) == global_ctx

    @test name(mod) == "SomeModule"
    name!(mod, "SomeOtherName")
    @test name(mod) == "SomeOtherName"

    dispose(mod)
end

Context() do ctx
    mod = LLVM.Module("SomeModule", ctx)
    @test context(mod) == ctx
end

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    clone = LLVM.Module(mod)
    @test mod != clone
    @test context(clone) == ctx
    dispose(clone)

    show(devnull, mod)

    inline_asm!(mod, "nop")
    @test occursin("module asm", sprint(io->show(io,mod)))

    dummyTriple = "SomeTriple"
    triple!(mod, dummyTriple)
    @test triple(mod) == dummyTriple

    dummyLayout = "e-p:64:64:64"
    datalayout!(mod, dummyLayout)
    @test convert(String, datalayout(mod)) == dummyLayout

    if LLVM.version() >= v"8.0"
        md = Metadata(ConstantInt(42, ctx))

        mod_flags = flags(mod)
        mod_flags["foobar", LLVM.API.LLVMModuleFlagBehaviorError] = md

        @test occursin("!llvm.module.flags = !{!0}", sprint(io->show(io,mod)))
        @test occursin(r"!0 = !\{i\d+ 1, !\"foobar\", i\d+ 42\}", sprint(io->show(io,mod)))

        @test mod_flags["foobar"] == md
        @test_throws KeyError mod_flags["foobaz"]
    end
end
end

# type iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    st = LLVM.StructType("SomeType", ctx)
    ft = LLVM.FunctionType(st, [st])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    let ts = types(mod)
        @test keytype(ts) == String
        @test valtype(ts) == LLVMType

        @test haskey(ts, "SomeType")
        @test ts["SomeType"] == st

        @test !haskey(ts, "SomeOtherType")
        @test_throws KeyError ts["SomeOtherType"]
    end
end
end

# metadata iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    node = MDNode([MDString("SomeMDString", ctx)], ctx)

    let mds = metadata(mod)
        @test keytype(mds) == String
        @test valtype(mds) == Vector{LLVM.MetadataAsValue}

        push!(mds, "SomeMDNode", node)

        @test haskey(mds, "SomeMDNode")
        mdvals = mds["SomeMDNode"]
        @test mdvals[1] == node

        @test !haskey(mds, "SomeOtherMDNode")
        @test_throws KeyError mds["SomeOtherMDNode"]
    end
end
end

# global variable iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    dummygv = GlobalVariable(mod, LLVM.Int32Type(ctx), "SomeGlobal")

    let gvs = globals(mod)
        @test eltype(gvs) == typeof(dummygv)

        @test first(gvs) == dummygv
        @test last(gvs) == dummygv

        for gv in gvs
            @test gv == dummygv
        end

        @test collect(gvs) == [dummygv]

        @test haskey(gvs, "SomeGlobal")
        @test gvs["SomeGlobal"] == dummygv

        @test !haskey(gvs, "SomeOtherGlobal")
        @test_throws KeyError gvs["SomeOtherGlobal"]
    end
end
end

# function iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    st = LLVM.StructType("SomeType", ctx)
    ft = LLVM.FunctionType(st, [st])
    @test isempty(functions(mod))

    dummyfn = LLVM.Function(mod, "SomeFunction", ft)
    let fns = functions(mod)
        @test eltype(fns) == LLVM.Function

        @test !isempty(fns)

        @test first(fns) == dummyfn
        @test last(fns) == dummyfn

        for fn in fns
            @test fn == dummyfn
        end

        @test collect(fns) == [dummyfn]

        @test haskey(fns, "SomeFunction")
        @test fns["SomeFunction"] == dummyfn

        @test !haskey(fns, "SomeOtherFunction")
        @test_throws KeyError fns["SomeOtherFunction"]
    end
end
end

end

@testset "function" begin

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx)])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    show(devnull, fn)

    # @show personality(fn)

    @test intrinsic_id(fn) == 0

    @test callconv(fn) == LLVM.API.LLVMCCallConv
    callconv!(fn, LLVM.API.LLVMFastCallConv)
    @test callconv(fn) == LLVM.API.LLVMFastCallConv

    @test LLVM.gc(fn) == ""
    gc!(fn, "SomeGC")
    @test LLVM.gc(fn) == "SomeGC"

    let fns = functions(mod)
        @test fn in fns
        unsafe_delete!(mod, fn)
        @test isempty(fns)
    end
end
end

# attributes
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx)])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    let attrs = function_attributes(fn)
        @test eltype(attrs) == Attribute

        @test length(attrs) == 0

        let attr = EnumAttribute("sspreq", 0, ctx)
            @test kind(attr) != 0
            @test value(attr) == 0
            push!(attrs, attr)
            @test collect(attrs) == [attr]

            delete!(attrs, attr)
            @test length(attrs) == 0
        end

        let attr = EnumAttribute("sspreq")
            @test value(attr) == 0
        end

        let attr = StringAttribute("nounwind", "", ctx)
            @test kind(attr) == "nounwind"
            @test value(attr) == ""
            push!(attrs, attr)
            @test collect(attrs) == [attr]

            delete!(attrs, attr)
            @test length(attrs) == 0
        end

        let attr = StringAttribute("nounwind")
            @test value(attr) == ""
        end
    end

    for i in 1:length(parameters(fn))
        let attrs = parameter_attributes(fn, i)
            @test eltype(attrs) == Attribute
            @test length(attrs) == 0
        end
    end

    let attrs = return_attributes(fn)
        @test eltype(attrs) == Attribute
        @test length(attrs) == 0
    end
end
end

# parameter iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx)])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    let params = parameters(fn)
        @test eltype(params) == Argument

        @test length(params) == 1

        intparam = params[1]
        @test first(params) == intparam
        @test last(params) == intparam

        for param in params
            @test param == intparam
        end

        @test collect(params) == [intparam]
    end
end
end

# basic block iteration
Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)
    @test isempty(blocks(fn))

    entrybb = BasicBlock(fn, "SomeBasicBlock")
    @test entry(fn) == entrybb
    let bbs = blocks(fn)
        @test eltype(bbs) == BasicBlock

        @test !isempty(bbs)
        @test length(bbs) == 1

        @test first(bbs) == entrybb
        @test last(bbs) == entrybb

        for bb in bbs
            @test bb == entrybb
        end

        @test collect(bbs) == [entrybb]
    end
end
end

end

@testset "basic blocks" begin

Builder() do builder
LLVM.Module("SomeModule") do mod
    ft = LLVM.FunctionType(LLVM.VoidType())
    fn = LLVM.Function(mod, "SomeFunction", ft)
    bbs = blocks(fn)

    bb2 = BasicBlock(fn, "SomeOtherBasicBlock")
    @test LLVM.parent(bb2) == fn

    bb1 = BasicBlock(bb2, "SomeBasicBlock")
    @test LLVM.parent(bb2) == fn

    @test collect(bbs) == [bb1, bb2]

    move_before(bb2, bb1)
    @test collect(bbs) == [bb2, bb1]

    move_after(bb2, bb1)
    @test collect(bbs) == [bb1, bb2]

    @test bb1 in bbs
    @test bb2 in bbs
    delete!(fn, bb1)
    unsafe_delete!(fn, bb2)
    @test isempty(bbs)
end
end

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)

    bb2 = BasicBlock(fn, "SomeOtherBasicBlock", ctx)
    @test LLVM.parent(bb2) == fn
    @test isempty(instructions(bb2))

    bb1 = BasicBlock(bb2, "SomeBasicBlock", ctx)
    @test LLVM.parent(bb2) == fn
    position!(builder, bb1)
    brinst = br!(builder, bb2)
    position!(builder, bb2)
    retinst = ret!(builder)
    @test !isempty(instructions(bb2))

    @test terminator(bb1) == brinst
    @test terminator(bb2) == retinst

    # instruction iteration
    let insts = instructions(bb1)
        @test eltype(insts) == Instruction

        @test first(insts) == brinst
        @test last(insts) == brinst

        for inst in insts
            @test inst == brinst
        end

        @test collect(insts) == [brinst]
    end
end
end
end

end

@testset "instructions" begin

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)
    @test isempty(parameters(fn))

    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int1Type(ctx), LLVM.Int1Type(ctx)])
    fn = LLVM.Function(mod, "SomeOtherFunction", ft)
    @test !isempty(parameters(fn))

    bb1 = BasicBlock(fn, "entry", ctx)
    bb2 = BasicBlock(fn, "then", ctx)
    bb3 = BasicBlock(fn, "else", ctx)

    position!(builder, bb1)
    brinst = br!(builder, parameters(fn)[1], bb2, bb3)
    @test opcode(brinst) == LLVM.API.LLVMBr

    position!(builder, bb2)
    retinst = ret!(builder)

    position!(builder, bb3)
    retinst = ret!(builder)

    # terminators

    @test isterminator(brinst)
    @test isconditional(brinst)
    @test condition(brinst) == parameters(fn)[1]
    condition!(brinst, parameters(fn)[2])
    @test condition(brinst) == parameters(fn)[2]

    let succ = successors(terminator(bb1))
        @test eltype(succ) == BasicBlock

        @test length(succ) == 2

        @test succ[1] == bb2
        @test succ[2] == bb3

        @test collect(succ) == [bb2, bb3]
        @test first(succ) == bb2
        @test last(succ) == bb3

        i = 1
        for bb in succ
            @test 1 <= i <= 2
            if i == 1
                @test bb == bb2
            elseif i == 2
                @test bb == bb3
            end
            i += 1
        end

        succ[2] = bb3
        @test succ[2] == bb3
    end

    # general stuff

    @test LLVM.parent(brinst) == bb1

    # metadata
    mdval = MDNode([MDString("whatever", ctx)], ctx)
    let md = metadata(brinst)
        @test keytype(md) == LLVM.MD
        @test valtype(md) == LLVM.MetadataAsValue

        @test isempty(md)
        @test !haskey(md, LLVM.MD_dbg)

        md[LLVM.MD_dbg] = mdval
        @test md[LLVM.MD_dbg] == mdval

        @test !isempty(md)
        @test haskey(md, LLVM.MD_dbg)

        @test !haskey(md, LLVM.MD_tbaa)
        @test_throws KeyError md[LLVM.MD_tbaa]

        delete!(md, LLVM.MD_dbg)

        @test isempty(md)
        @test !haskey(md, LLVM.MD_dbg)
    end
    let val = ConstantInt(42, ctx)
        md = Metadata(val)
        @test first(operands(LLVM.Value(md, ctx))) == val
    end

    @test retinst in instructions(bb3)
    delete!(bb3, retinst)
    @test !(retinst in instructions(bb3))
    @test opcode(retinst) == LLVM.API.LLVMRet   # make sure retinst is still alive

    @test brinst in instructions(bb1)
    unsafe_delete!(bb1, brinst)
    @test !(brinst in instructions(bb1))
end
end
end

end

end


text:
module test_bukdu_system

module NAStuff

using Bukdu
struct NA <: ApplicationController; conn::Conn end
index(c::NA) = nothing
hello(c::NA) = TEST_INTERNAL_ERROR

end # module NAStuff

using Test
using Bukdu
import .Bukdu: System

function index
end

routes() do
    get("/na", NAStuff.NA, index)
    post("/hello", NAStuff.NA, NAStuff.hello)
end

result = Router.call(get, "/na")
@test result.route.action === System.not_applicable

result = Router.call(post, "/hello")
@test result.route.action === System.internal_error

Routing.empty!()

@test System._regularize_text("가1", 1) == "가"

end # module test_bukdu_system

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
function Base.nullspace{T}(A::Operator{T};tolerance=10eps(real(T)),maxlength=1_000_000)
   K=transpose_nullspace(qrfact(A'),tolerance,maxlength)
    # drop extra rows, and use QR to determine rank
    Q,R=qr(K,Val{true})
    ind=findfirst(r->abs(r)≤100tolerance,diag(R))
    Kret=ind==0?Q:Q[:,1:ind-1]
    Fun(vec(Kret'),ArraySpace(domainspace(A),1,ind-1))
end

function transpose_nullspace(QR::QROperator,tolerance,maxlength)
    T=eltype(QR)
    resizedata!(QR,:,100)

    m=size(QR.H,1)
    K=zeros(100,m-1)
    K[1:m-1,1:m-1]=eye(m-1)

    for k=1:10
        v=QR.H[:,k]
        QQ=eye(T,m)-2v*v'
        K[:,:]=K*QQ[1:end-1,2:end]
        K[k+m-1,:]=QQ[end,2:end]
    end
    k=11
    α=0.9

    while slnorm(K,floor(Int,k^α),:) >tolerance*k
        if k > maxlength
            warn("Max length of $maxlength reached.")
            break
        end

        if k ≥ QR.ncols
            resizedata!(QR,:,k+100)
        end

        v=QR.H[:,k]
        QQ=(eye(m)-2v*v')
        K[:,:]=K*QQ[1:end-1,2:end]

        if k+m-1 > size(K,1)
            K=pad(K,k+m+100,:)
        end
        K[k+m-1,:]=QQ[end,2:end]
        k+=1
    end

    K[1:floor(Int,k^α),:]
end

Base.nullspace{OO<:Operator}(A::Array{OO};kwds...) = nullspace(interlace(A);kwds...)


text:
using Test
using MultivariateBases

@testset "Orthogonal" begin
    orthogonal_test(ProbabilistsHermiteBasis, x -> [
        1,
        x,
        x^2 - 1,
        x^3 - 3x,
        x^4 - 6x^2 + 3
    ], true)
    orthogonal_test(PhysicistsHermiteBasis, x -> [
        1,
        2x,
        4x^2 - 2,
        8x^3 - 12x,
        16x^4 - 48x^2 + 12
    ], true)
end

@testset "API degree = $degree" for degree in 0:3
    api_test(ProbabilistsHermiteBasis, degree)
    api_test(PhysicistsHermiteBasis, degree)
end

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
using AdFem

function plot_cond(iter, k_chip_nodes_pred)

    k_chip_nodes = eval_f_on_fem_pts(k_exact,mesh)[chip_fem_idx_nodes]
    k_chip_nodes_out = k_chip_nodes_pred

    k_all  = k_mold * ones(nnode)
    k_all[chip_fem_idx_nodes] .= k_chip_nodes

    k_all_out = k_mold * ones(nnode)
    k_all_out[chip_fem_idx_nodes] .= k_chip_nodes_out

    figure(figsize=(15,4))
    subplot(131)
    visualize_scalar_on_fem_points(k_all, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("exact solid conductivity")

    subplot(132)
    visualize_scalar_on_fem_points(k_all_out, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("predicted solid conductivity")

    subplot(133)
    visualize_scalar_on_fem_points(k_all_out .- k_all, mesh)
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("difference in solid conductivity")

    tight_layout()
    savefig("fn$trialnum/nn_cond$iter.png")

end

function plot_velo_pres_temp_cond(iter)
    S = run(sess, S_computed)
    u_out, v_out, p_out, T_out = S[1:nnode],
                                 S[ndof+1:ndof+nnode],
                                 S[2*ndof+1:2*ndof+nelem],
                                 S[2*ndof+nelem+1:2*ndof+nelem+nnode]

    u_obs, v_obs, p_obs, T_obs = S_data[1:nnode],
                                 S_data[ndof+1:ndof+nnode],
                                 S_data[2*ndof+1:2*ndof+nelem],
                                 S_data[2*ndof+nelem+1:2*ndof+nelem+nnode]

    figure(figsize=(15,4))
    subplot(131)
    title("exact x velocity")
    visualize_scalar_on_fem_points(u_obs .* u_std, mesh)
    subplot(132)
    title("predicted x velocity")
    visualize_scalar_on_fem_points(u_out .* u_std, mesh)
    subplot(133)
    title("difference in x velocity")
    visualize_scalar_on_fem_points(u_out .* u_std .- u_obs .* u_std, mesh)
    tight_layout()
    savefig("fn$trialnum/nn_velox$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact y velocity")
    visualize_scalar_on_fem_points(v_obs .* u_std, mesh)
    subplot(132)
    title("predicted y velocity")
    visualize_scalar_on_fem_points(v_out .* u_std, mesh)
    subplot(133)
    title("difference in y velocity")
    visualize_scalar_on_fem_points(v_out .* u_std .- v_obs .* u_std, mesh)
    tight_layout()
    savefig("fn$trialnum/nn_veloy$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact pressure")
    visualize_scalar_on_fvm_points(p_obs .* p_std, mesh)
    subplot(132)
    title("predicted pressure")
    visualize_scalar_on_fvm_points(p_out .* p_std, mesh)
    subplot(133)
    title("difference in pressure")
    visualize_scalar_on_fvm_points(p_out .* p_std .- p_obs .* p_std,  mesh)
    tight_layout()
    savefig("fn$trialnum/nn_pres$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact temperature")
    visualize_scalar_on_fem_points(T_obs .* T_infty .+ T_infty, mesh)
    subplot(132)
    title("predicted temperature")
    visualize_scalar_on_fem_points(T_out .* T_infty .+ T_infty, mesh)
    subplot(133)
    title("difference in temperature")
    visualize_scalar_on_fem_points(T_out  .* T_infty .- T_obs .* T_infty, mesh)
    tight_layout()
    savefig("fn$trialnum/nn_temp$iter.png")

    #---------------------------------------------------------------------------
    k_chip_nodes = eval_f_on_fem_pts(k_exact,mesh)[chip_fem_idx_nodes]
    # k_chip_nodes_out = eval_f_on_fem_pts(k_nn,mesh)[chip_fem_idx_nodes]

    # xy = mesh.nodes
    # xy2 = zeros(mesh.nedge, 2)
    # for i = 1:mesh.nedge
    #     xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
    # end
    # xy = [xy;xy2]

    # x, y = xy[chip_fem_idx_nodes, 1], xy[chip_fem_idx_nodes, 2]
    k_chip_nodes_out = run(sess, k_chip)[1:length(chip_fem_idx_nodes)]

    k_all  = k_mold * ones(nnode)
    k_all[chip_fem_idx_nodes] .= k_chip_nodes

    k_all_out = k_mold * ones(nnode)
    k_all_out[chip_fem_idx_nodes] .= k_chip_nodes_out

    figure(figsize=(15,4))
    subplot(131)
    visualize_scalar_on_fem_points(k_all, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("exact solid conductivity")

    subplot(132)
    visualize_scalar_on_fem_points(k_all_out, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("predicted solid conductivity")

    subplot(133)
    visualize_scalar_on_fem_points(k_all_out .- k_all, mesh)
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("difference in solid conductivity")

    tight_layout()
    savefig("fn$trialnum/nn_cond$iter.png")

end

function pix_plot_velo_pres_temp_cond(iter)

    S = run(sess, S_computed)
    u_out, v_out, p_out, T_out = S[1:nnode],
                                 S[ndof+1:ndof+nnode],
                                 S[2*ndof+1:2*ndof+nelem],
                                 S[2*ndof+nelem+1:2*ndof+nelem+nnode]

    u_obs, v_obs, p_obs, T_obs = S_data[1:nnode],
                                 S_data[ndof+1:ndof+nnode],
                                 S_data[2*ndof+1:2*ndof+nelem],
                                 S_data[2*ndof+nelem+1:2*ndof+nelem+nnode]

    figure(figsize=(15,4))
    subplot(131)
    title("exact x velocity")
    visualize_scalar_on_fem_points(u_obs .* u_std, mesh)
    subplot(132)
    title("predicted x velocity")
    visualize_scalar_on_fem_points(u_out .* u_std, mesh)
    subplot(133)
    title("difference in x velocity")
    visualize_scalar_on_fem_points(u_out .* u_std .- u_obs .* u_std, mesh)
    tight_layout()
    savefig("fn$trialnum/pix_velox$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact y velocity")
    visualize_scalar_on_fem_points(v_obs .* u_std, mesh)
    subplot(132)
    title("predicted y velocity")
    visualize_scalar_on_fem_points(v_out .* u_std, mesh)
    subplot(133)
    title("difference in y velocity")
    visualize_scalar_on_fem_points(v_out .* u_std .- v_obs .* u_std, mesh)
    tight_layout()
    savefig("fn$trialnum/pix_veloy$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact pressure")
    visualize_scalar_on_fvm_points(p_obs .* p_std, mesh)
    subplot(132)
    title("predicted pressure")
    visualize_scalar_on_fvm_points(p_out .* p_std, mesh)
    subplot(133)
    title("difference in pressure")
    visualize_scalar_on_fvm_points(p_out .* p_std .- p_obs .* p_std,  mesh)
    tight_layout()
    savefig("fn$trialnum/pix_pres$iter.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact temperature")
    visualize_scalar_on_fem_points(T_obs .* T_infty .+ T_infty, mesh)
    subplot(132)
    title("predicted temperature")
    visualize_scalar_on_fem_points(T_out .* T_infty .+ T_infty, mesh)
    subplot(133)
    title("difference in temperature")
    visualize_scalar_on_fem_points(T_out  .* T_infty .- T_obs .* T_infty, mesh)
    tight_layout()
    savefig("fn$trialnum/pix_temp$iter.png")

    #---------------------------------------------------------------------------
    k_chip_nodes = eval_f_on_fem_pts(k_exact,mesh)[chip_fem_idx_nodes]

    xy = mesh.nodes
    xy2 = zeros(mesh.nedge, 2)
    for i = 1:mesh.nedge
        xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
    end
    xy = [xy;xy2]

    x, y = xy[chip_fem_idx_nodes, 1], xy[chip_fem_idx_nodes, 2]
    k_chip_nodes_out = run(sess, k_chip[1:length(chip_fem_idx_nodes)])

    k_all  = k_mold * ones(nnode)
    k_all[chip_fem_idx_nodes] .= k_chip_nodes

    k_all_out = k_mold * ones(nnode)
    k_all_out[chip_fem_idx_nodes] .= k_chip_nodes_out

    figure(figsize=(15,4))
    subplot(131)
    visualize_scalar_on_fem_points(k_all, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("exact solid conductivity")

    subplot(132)
    visualize_scalar_on_fem_points(k_all_out, mesh, vmin=minimum(k_chip_nodes), vmax=maximum(k_chip_nodes));
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("predicted solid conductivity")

    subplot(133)
    visualize_scalar_on_fem_points(k_all_out .- k_all, mesh)
    xlim(chip_left, chip_right); ylim(chip_top, chip_bottom)
    title("difference in solid conductivity")

    tight_layout()
    savefig("fn$trialnum/pix_cond$iter.png")

end


text:
import Base.length
function length(g::Union(Likelihood, Kernel))
	sum([length(g.(n)) for n in names(g)])
end


text:
export lowercase

import Base: lowercase
import Base: show

"""
    lowercase()

Creates an lowercase operator, which forces each value to be in lower case

# Producing

Stream of type `<: Subscribable{L}` where L referes to type of data of input Observable

# Examples

```jldoctest
using Rocket

source = of("Hello, world!")
subscribe!(source |> lowercase(), logger())
;

# output

[LogActor] Data: hello, world!
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
lowercase() = LowercaseOperator()

struct LowercaseOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{L}, operator::LowercaseOperator, source) where L
    return proxy(L, source, LowercaseProxy())
end

operator_right(::LowercaseOperator, ::Type{L}) where L = L

struct LowercaseProxy <: ActorProxy end

actor_proxy!(::Type{L}, proxy::LowercaseProxy, actor::A) where { L, A } = LowercaseActor{L, A}(actor)

struct LowercaseActor{L, A} <: Actor{L}
    actor :: A
end

on_next!(actor::LowercaseActor{L}, data::L) where L = next!(actor.actor, lowercase(data))
on_error!(actor::LowercaseActor, err)               = error!(actor.actor, err)
on_complete!(actor::LowercaseActor)                 = complete!(actor.actor)

Base.show(io::IO, ::LowercaseOperator)         = print(io, "LowercaseOperator()")
Base.show(io::IO, ::LowercaseProxy)            = print(io, "LowercaseProxy()")
Base.show(io::IO, ::LowercaseActor{L}) where L = print(io, "LowercaseActor($L)")


text:
#------------------------------------------------------------------------------#
#                                    Module                                    #
#------------------------------------------------------------------------------#

"""
# Description

    module HPT

High-level interface module to the Hermite Polynomial Tensor package.

# Usage

This module is automatically `include`d by the `HermitePolyTensor` package,  and
provides a common interface to the various supported backends.

```julia-repl
julia> using HermitePolyTensor

julia> typeof(HPT)
Module

```
"""
module HPT

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

import Combinatorics

#------------------------------------------------------------------------------#
#                                   Includes                                   #
#------------------------------------------------------------------------------#

include("backend/se/seHTC.jl")
import .seHTC   # Never do "using .seHTC" to avoid function name clashing

include("backend/aa/aaHTC.jl")
import .aaHTC   # Never do "using .aaHTC" to avoid function name clashing

#------------------------------------------------------------------------------#
#                                   Exports                                    #
#------------------------------------------------------------------------------#

# NOTE: This module is @reexport'ed by the main package module

export SVHP, HTC

export HT

#------------------------------------------------------------------------------#
#                                  Constants                                   #
#------------------------------------------------------------------------------#

# Valid Backends
vBE = Dict(
    "aa" => ("AA", "AbsAlg", "AbstractAlgebra"),
    "se" => ("SE", "SymEng", "SymEngine"),
)

# Default Backend
dBE = "aa"  # Since it's way faster than "se"

#------------------------------------------------------------------------------#
#                             Auxiliary Functions                              #
#------------------------------------------------------------------------------#

"""
    theBE(BE::String)::String

Returns the backend name corresponding to `BE`---the backend name or  alias,  or
the default backend `String` for invalid name/alias.
"""
function theBE(BE::String)::String
    for (k, v) in vBE
        if BE == k
            return k
        elseif BE in v
            return k
        end
    end
    return dBE
end

#------------------------------------------------------------------------------#
#                             Interface Functions                              #
#------------------------------------------------------------------------------#

"""
# Description

    SVHP(
        n::Int; ax::Int = 1, st::Bool = true, be::String = dBE
    )

Using  the  backend  specified  by  `be`,  returns  a  symbolic   `n`-th   order
Single-Variable Hermite Polynomial (SVHP) in the axis `ax`, with axis 1, 2,  and
3 corresponding to the x, y, and z directions, respectively.  A  "probabilist's"
polynomial  is  returned  if  `st`  is  `true`  (the  default),   otherwise,   a
"physicist's".

# Usage

```julia-repl
julia> using HermitePolyTensor

julia> SVHP(3, be = "SE")   # SymEngine backend
-3*x + x^3

julia> SVHP(4, ax = 2)      # Default AbstractAlgebra backend
y^4-6//1*y^2+3//1

```
"""
function SVHP(
    n::Int; ax::Int = 1, st::Bool = true, be::String = dBE
)
    # Execution
    BE = theBE(be)
    if BE == "aa"
        return aaHTC.SVHP(n, ax = ax, st = st)
    elseif BE == "se"
        return seHTC.SVHP(n, ax = ax, st = st)
    end
end

"""
# Description

    HTC(
        id::String; D::Int = 2, st::Bool = true, be::String = dBE
    )

Using the backend specified by `be`, returns a symbolic *component* of a Hermite
Tensor Polynomial.

The indices are specified by `id`, a `String` containing x's, y's,  and/or  z's;
invalid characters on `id` are ignored, while `id`  without  x's,  y's,  or  z's
represent a zero-rank  tensor  (whose  only  component  is  `1`  for  all  space
dimensions).

The Euclidean space dimensionality is specified  by  `D`.  The  tensor  rank  is
implied by the provided `id`, and is always the total  count  of  *valid*  index
characters, considering `D`.

A "probabilist's" polynomial is  returned  if  `st`  is  `true`  (the  default),
otherwise, a "physicist's".

# Usage

```julia-repl
julia> using HermitePolyTensor

julia> HTC("xzyxy", 3)      # Default AbstractAlgebra "aa" backend
x^2*y^2*z-x^2*z-y^2*z+z

julia> HTC("xzyxy", 3, be = "se")   # SymEngine "se" backend
z - x^2*z - y^2*z + x^2*y^2*z

julia> HTC("xzyxy", 3, st = false)  # Physicists's H(*x*)
32//1*x^2*y^2*z-16//1*x^2*z-16//1*y^2*z+8//1*z

```
"""
function HTC(
    id::String; D::Int = 2, st::Bool = true, be::String = dBE
)
    # Execution
    BE = theBE(be)
    if BE == "aa"
        return aaHTC.HTC(id, D = D, st = st)
    elseif BE == "se"
        return seHTC.HTC(id, D = D, st = st)
    end
end

#------------------------------------------------------------------------------#
#                          Multi-Component Functions                           #
#------------------------------------------------------------------------------#

"""
# Description

    HT(
        n::Int; D::Int = 2, st::Bool = true, be::String = dBE
    )::Dict{String,NamedTuple{(:val, :set),Tuple{Any,Set{String}}}}

Using the backend specified by the `be` kwarg, returns  a  dictionary  with  all
*unique* components of a Hermite Tensor Polynomial of order (and rank) `n`, in a
`D`-dimensional Euclidean space.

The dictionary keys are  the  "normalized"  unique  components  of  the  Hermite
Polynomial Tensor, defined by indices `"x"^ax * "y"^ay * "z"^az`,  according  to
the `HTC` function  index  conventions,  that  are  solution  to  the  following
Diophantine [1, 2] equation: `ax + ay + az = n`—there are  `((D  n))`,  or,  `D`
multichoose `n` such indices, or equivalently, `(D + n - 1) choose D` = `(D +  n
- 1)! / (D! * (n - 1)!)` indices alike [2].

Dictionary entries are `namedTuple`'s with symbols (i) `:val` and  (ii)  `:set`,
whose corresponding values are (i) the polynomial component, and (ii) the  index
`Set` to the components having the same  polynomial,  including  the  dictionary
key. Since such `Set` elements are concatenated permutations of  the  normalized
index multiset, with multiplicities `ax`, `ay`, and `az`, the number of multiset
permutations is given by the multinomial coefficient [3] and is equal to  `n!  /
(ax! * ay! * az!)`, or, conversely, by `(ax + ay + az)! / (ax! * ay! * az!)`.

"Probabilist's" polynomials are  returned  if  `st`  is  `true`  (the  default),
otherwise, "physicist's" ones.

# Usage

```julia-repl
```

# References

[1]:  Wikipedia  contributors.  "Diophantine  equation."  Wikipedia,  The   Free
Encyclopedia. Wikipedia, The Free Encyclopedia, 19 Dec. 2018. Web. 25 Jan. 2019.

[2]: Wikipedia contributors. "Combination." Wikipedia,  The  Free  Encyclopedia.
Wikipedia, The Free Encyclopedia, 22 Jan. 2019. Web. 25 Jan. 2019.

[3]: Wikipedia contributors. "Permutation." Wikipedia,  The  Free  Encyclopedia.
Wikipedia, The Free Encyclopedia, 23 Jan. 2019. Web. 25 Jan. 2019.
"""
function HT(
    n::Int; D::Int = 2, st::Bool = true, be::String = dBE
)::Dict{String,NamedTuple{(:val, :set),Tuple{Any,Set{String}}}}
    # Validations
    if n < 0
        throw(DomainError("n = $n outside the valid domain [0, ∞)"))
    elseif D <= 0 || D >= 4
        throw(DomainError("dimension D = $D is outside the valid domain [1, 3]"))
    end
    # Execution
    BE = theBE(be)
    ret = Dict{String,NamedTuple{(:val, :set),Tuple{Any,Set{String}}}}()
    pat = D == 1 ? "x"^n : D == 2 ? "x"^n * "y"^n : "x"^n * "y"^n * "z"^n
    for key in [
        join(i) for i in
        Combinatorics.multiset_combinations(pat, n)
    ]
        theVal = HTC(key, D = D, st = st, be = BE)
        theSet = Set([join(i) for i in
            Combinatorics.multiset_permutations(key, n)
        ])
        ret[key] = (val = theVal, set = theSet)
    end
    ret
end

end # module


text:
vars=25
i_train=10_000
i_test=10_000
cc1=10.0
neurons=Int64(floor(i_train/(50))) # this automatically
#generates number of neurons== number of clusters (§3.1)
_case_noise_=4

# __scale__=10
__maxima__=10
__c__=rand(MersenneTwister(0),__maxima__)
__a__=5(rand(MersenneTwister(1),__maxima__,vars).-1/2)
function ff(x)
    __ret__=0
    for i ∈ 1:length(__c__)
        __ret__+=1/(__c__[i]+sum((x.-__a__[i,:]).^2))
    end
    return __ret__
    # __scale__*length(x) + sum(x.^2 - __scale__*cos.(2.0.*π.*x)) #Rastrigin
end

xx_train=5(rand(MersenneTwister(2),i_train,vars).-1/2)
# xx_train=10rand(i_train,vars)
yy_train=zeros(i_train)
for i=1:i_train yy_train[i]=ff(xx_train[i,:]) end
scatter3d(xx_train[:,1],xx_train[:,2],yy_train,label=false,
    title="Shekel Function")
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1

if _case_noise_==1
    __noise__=rand(i_train)
elseif _case_noise_==2
    __noise__=rand(Normal(0.0, 1.0),i_train)
elseif _case_noise_==3
    __noise__=rand(GeneralizedPareto(0.0,1.0,0.25),i_train)
elseif _case_noise_==4
    __noise__=rand(LogNormal(0.25, 0.8),i_train)
elseif _case_noise_==5
    __noise__=[rand(LogNormal(0.25, 0.8),i_train);10.0.-rand(Exponential(3.0),i_train)]
    __noise__=[__noise__;15.0.+rand(Frechet(10.0,10.0),i_train)]
    __noise__=__noise__[1:3:end]
end
__noise__.-=mean(__noise__)
__noise__./=mean(abs.(__noise__))/(0.05*0.45)
histogram(__noise__,label=sum(__noise__)/sum(abs.(__noise__)),
        title=mean(abs.(__noise__)))
iso=sortperm(yy_train)
plot(yy_train[iso])
yy_train.+=__noise__
display(plot!(yy_train[iso],legend=false))

xx_test=5(rand(MersenneTwister(3),i_test,vars).-1/2)
# xx_test=10rand(i_test,vars)
yy_test=zeros(i_test)
for i=1:i_test yy_test[i]=ff(xx_test[i,:]) end
# yy_test.+=noise_fact.*rand(i_test)
yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1

## Random Forests
num_trees=100
regr_3 = RandomForestRegressor(n_trees=num_trees);
tmp_xx_train=copy(xx_train);tmp_yy_train=copy(yy_train);
tmp_xx_test=copy(xx_test);tmp_yy_test=copy(yy_test);
DecisionTree.fit!(regr_3, tmp_xx_train, tmp_yy_train[:,1]);
predl1 = DecisionTree.predict(regr_3, tmp_xx_train);
predl1_test = DecisionTree.predict(regr_3, tmp_xx_test);
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot(err[ier],label="Random Forests",color=:black,
    linestyle=:dash,legend=:topleft,linewidth=3))
__err__=string(mean(abs.(err)))

sci_model = fit!(AdaBoostRegressor(n_estimators=100),xx_train,yy_train)
predl1=ScikitLearn.predict(sci_model,xx_train)
predl1_test=ScikitLearn.predict(sci_model,xx_test)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="AdaBoost",color=:black,
    linestyle=:dashdot,legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

num_round = 10
bst = xgboost(copy(xx_train), num_round, label = copy(yy_train), eta = 1, max_depth = 7)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="XGBoost",color=:black,
    linestyle=:dot,legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

##### ANNBN RBF

inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)

a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="ANNBN",color=:black,linestyle=:dashdotdot,
    legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

println(__err__)

#


text:
module YouBotGripping

using  ModiaLang
import Modia3D
using  Unitful

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")
include("$(ModiaLang.path)/models/Translational.jl")

import Modia3D
using  Modia3D.ModiaInterface

# some constants
simplifiedContact = true  # use boxes instead of meshes for finger contact

# visual materials
vmatYellow    = VisualMaterial(color=[255,215,0])    # ground
vmatBlue      = VisualMaterial(color="DodgerBlue3")  # table
vmatGrey      = VisualMaterial(color="DarkGrey")     # work piece
vmatInvisible = VisualMaterial(transparency=1.0)     # contact boxes
tableX            = 0.3
tableY            = 0.3
tableZ            = 0.01
heigthLeg         = 0.355
widthLeg          = 0.02
diameter          = 0.05
rightFingerOffset = 0.03

# all needed paths to fileMeshes
base_frame_obj           = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/base_frame.obj")
plate_obj                = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/plate.obj")
front_right_wheel_obj    = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/front-right_wheel.obj")
front_left_wheel_obj     = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/front-left_wheel.obj")
back_right_wheel_obj     = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/back-right_wheel.obj")
back_left_wheel_obj      = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/back-left_wheel.obj")
arm_base_frame_obj       = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_base_frame.obj")
arm_joint_1_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_1.obj")
arm_joint_2_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_2.obj")
arm_joint_3_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_3.obj")
arm_joint_4_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_4.obj")
arm_joint_5_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_5.obj")
gripper_base_frame_obj   = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_base_frame.obj")
gripper_left_finger_obj  = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_left_finger.obj")
gripper_right_finger_obj = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_right_finger.obj")

# Drive train data: motor inertias and gear ratios
nullRot = nothing #SMatrix{3,3,Float64,9}(Matrix(1.0I, 3, 3))

motorInertia1 = 0.0000135 + 0.000000409
motorInertia2 = 0.0000135 + 0.000000409
motorInertia3 = 0.0000135 + 0.000000071
motorInertia4 = 0.00000925 + 0.000000071
motorInertia5 = 0.0000035 + 0.000000071
gearRatio1    = 156.0
gearRatio2    = 156.0
gearRatio3    = 100.0
gearRatio4    = 71.0
gearRatio5    = 71.0

m1=1.390
translation1 =[0.033,0,0]
rotation1 =Modia3D.rot1(180u"°")

m2=1.318
translation2=[0.155,0,0]
rotation2=Modia3D.rot123(90u"°", 0.0, -90u"°")

m3=0.821
translation3=[0,0.135,0]
rotation3=Modia3D.rot3(-90u"°")

m4=0.769
translation4=[0,0.11316,0]

m5=0.687
translation5=[0,0,0.05716]
rotation5=Modia3D.rot1(-90u"°")

### ----------------- Servo Model -----------------------
# parameters for Link
axisLink = 3
dLink  = 1.0
k1Link = 50.0
k2Link = 0.1
T2Link = 1.0
# parameters for Gripper
axisGripper = 2
k1Gripper = 50.0
k2Gripper = 0.1*100
T2Gripper = 1.0
motorInertiaGripper = 0.1
gearRatioGripper    = 1.0

#### ----------- Path Planning ------------------
referencePath1 = Modia3D.ReferencePath(
    names =    ["angle1", "angle2", "angle3", "angle4", "angle5", "gripper"],
    position = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    v_max =    [2.68512, 2.68512, 4.8879, 5.8997, 5.8997, 2.0],
    a_max =    [1.5, 1.5, 1.5, 1.5, 1.5, 0.5])

Modia3D.ptpJointSpace(referencePath = referencePath1, positions =
    [0.0  0.0    0.0       0.0   0.0  0.0;
     pi   pi/4   pi/4      0.0   0.0  diameter+0.01;
     pi   pi/4   pi/4      1.04  0.0  diameter+0.01;
     pi   pi/4   pi/4      1.04  0.0  diameter-0.002;
     pi   pi/4   pi/4      0.0   0.0  diameter-0.002;
     0.0  0.0    pi/2-0.3  0.0   0.0  diameter-0.002;
     0.0  0.3    pi/2-0.3  0.0   0.0  diameter-0.002;
     0.0  0.3    pi/2-0.3  0.0   0.0  diameter+0.01;
     0.0  0.0    pi/2      0.0   0.0  diameter+0.01;
     0.0  0.0    0.0       0.0   0.0  0.03])

getReferencePath() = referencePath1

# Controller Model
Controller = Model(
    # Interface
    refLoadAngle     = input,
    actualMotorAngle = input,
    actualMotorSpeed = input,
    refTorque        = output,
    # Components
    gainOuter     = Gain,
    angleFeedback = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadAngle, gainOuter.u)
        (gainOuter.y, angleFeedback.u1)
        (actualMotorAngle, angleFeedback.u2)
        (angleFeedback.y, gainInner.u)
        (gainInner.y, speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y, PI.u)
        (PI.y, refTorque)]
)

Drive = Model(
    # Interface
    refTorque   = input,
    actualSpeed = output,
    actualAngle = output,
    flange      = Flange,
    # Components
    torque       = UnitlessTorque,
    speedSensor  = UnitlessSpeedSensor,
    angleSensor  = UnitlessAngleSensor,
    motorInertia = Inertia,
    idealGear    = IdealGear,

    connect = :[
        (refTorque, torque.tau)
        (torque.flange, speedSensor.flange, angleSensor.flange, motorInertia.flange_a)
        (motorInertia.flange_b, idealGear.flange_a)
        (idealGear.flange_b, flange)
        (speedSensor.w, actualSpeed)
        (angleSensor.phi, actualAngle)
        ]
)

Servo = Model(
    refLoadAngle = input,
    flange       = Flange,
    ppi          = Controller,
    drive        = Drive,

    connect = :[
        # inputs of ppi
        (refLoadAngle, ppi.refLoadAngle)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorAngle, drive.actualAngle)
        # output of ppi --> input of motor
        (ppi.refTorque, drive.refTorque)
        # Flange of drive
        (drive.flange, flange)
        ]
)

# Controller Model
ControllerTrans = Model(
    # Interface
    refLoadPos     = input,
    actualMotorPos = input,
    actualMotorSpeed = input,
    refForce        = output,
    # Components
    gainOuter     = Gain,
    angleFeedback = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadPos, gainOuter.u)
        (gainOuter.y, angleFeedback.u1)
        (actualMotorPos, angleFeedback.u2)
        (angleFeedback.y, gainInner.u)
        (gainInner.y, speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y, PI.u)
        (PI.y, refForce)]
)

DriveTrans = Model(
    # Interface
    refForce    = input,
    actualSpeed = output,
    actualPos   = output,
    flange      = TranslationalFlange,
    # Components
    force       = UnitlessForce,
    speedSensor = UnitlessVelocitySensor,
    posSensor   = UnitlessPositionSensor,
    motorMass   = TranslationalMass,

    connect = :[
        (refForce, force.f)
        (force.flange, speedSensor.flange, posSensor.flange, motorMass.flange_a)
        (motorMass.flange_b, flange)
        (speedSensor.v, actualSpeed)
        (posSensor.s, actualPos)
        ]
)

ServoTrans = Model(
    refLoadPos = input,
    flange     = TranslationalFlange,
    ppi        = ControllerTrans,
    drive      = DriveTrans,

    connect = :[
        # inputs of ppi
        (refLoadPos, ppi.refLoadPos)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorPos, drive.actualPos)
        # output of ppi --> input of motor
        (ppi.refForce, drive.refForce)
        # Flange of drive
        (drive.flange, flange)
        ]
)

Ground = Model(
    ground = Object3D(parent=:world,
        translation=[0.0, 0.0, -0.005], feature = Solid(shape =Box(lengthX=2.5, lengthY=2.5, lengthZ=0.01), solidMaterial="DryWood", visualMaterial=vmatYellow))
)

Table = Model(
    plate = Object3D(parent=:world, translation=[0.0, 0.0, heigthLeg],
        feature=Solid(shape=Box(lengthX=tableX, lengthY=tableY, lengthZ=tableZ),
        collisionSmoothingRadius=0.001, solidMaterial="DryWood", collision=true, visualMaterial=vmatBlue)),
    leg1  = Object3D(parent=:plate, translation=[ tableX/2-widthLeg/2,  tableY/2-widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg),
        solidMaterial="DryWood", visualMaterial=vmatBlue)),
    leg2  = Object3D(parent=:plate, translation=[-tableX/2+widthLeg/2,  tableY/2-widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmatBlue)),
    leg3  = Object3D(parent=:plate, translation=[ tableX/2-widthLeg/2, -tableY/2+widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmatBlue)),
    leg4  = Object3D(parent=:plate, translation=[-tableX/2+widthLeg/2, -tableY/2+widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmatBlue))
)

Base = Model(
    base_frame = Object3D(parent=:world, translation=[-0.585, 0, 0.084],
        feature=Solid(shape=FileMesh(filename = base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=19.803))),
    plate = Object3D(parent=:base_frame, translation=[-0.159, 0, 0.046],
        feature=Solid(shape=FileMesh(filename = plate_obj), massProperties=MassPropertiesFromShapeAndMass(mass=2.397), contactMaterial="DryWood", collision=!simplifiedContact)),
    plate_box = Object3D(parent=:base_frame, translation=[-0.16, 0, 0.0702],# 0.06514],
        feature=Solid(shape=Box(lengthX=0.22, lengthY=0.22, lengthZ=0.01), massProperties=MassPropertiesFromShapeAndMass(mass=1.0e-9), visualMaterial=vmatInvisible, contactMaterial="DryWood", collision=simplifiedContact)),
    front_right_wheel = Object3D(parent=:base_frame, translation=[0.228, -0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = front_right_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4))),
    front_left_wheel = Object3D(parent=:base_frame, translation=[0.228, 0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = front_left_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4))),
    back_right_wheel = Object3D(parent=:base_frame, translation=[-0.228, -0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = back_right_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4))),
    back_left_wheel   = Object3D(parent=:base_frame, translation=[-0.228, 0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = back_left_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4)))
)

servoParameters1 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio1),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia1),
                          idealGear = Map(ratio=gearRatio1)
                      )
                  )

servoParameters2 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio2),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia2),
                          idealGear = Map(ratio=gearRatio2)
                      )
                  )

servoParameters3 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio3),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia3),
                          idealGear = Map(ratio=gearRatio3)
                      )
                  )

servoParameters4 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio4),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia4),
                          idealGear = Map(ratio=gearRatio4)
                      )
                  )

servoParameters5 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio5),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia5),
                          idealGear = Map(ratio=gearRatio5)
                      )
                  )

servoParameters6 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatioGripper),
                          gainInner = Map(k=k1Gripper),
                          PI        = Map(k=k2Gripper, T=T2Gripper)
                      ),
                      drive = Map(
                          motorMass = Map(m=motorInertiaGripper)
                      )
                  )

featureBody1 = Solid(shape = FileMesh(filename = arm_joint_1_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m1))
featureBody2 = Solid(shape = FileMesh(filename = arm_joint_2_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m2))
featureBody3 = Solid(shape = FileMesh(filename = arm_joint_3_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m3))
featureBody4 = Solid(shape = FileMesh(filename = arm_joint_4_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m4))
featureBody5 = Solid(shape = FileMesh(filename = arm_joint_5_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m5))

linkParameters1 = Map(
                    parent1 = Par(value = :(armBase_b)),
                    featureBody = featureBody1,
                    initRefPos = referencePath1.position[1],
                    trans = translation1,
                    rota = Par(value = :(rotation1))
)

linkParameters2 = Map(
                    parent1 = Par(value = :(link1.obj2)),
                    featureBody = featureBody2,
                    m = m2,
                    initRefPos = referencePath1.position[2],
                    trans = translation2,
                    rota = Par(value = :(rotation2))
)

linkParameters3 = Map(
                    parent1 = Par(value = :(link2.obj2)),
                    featureBody = featureBody3,
                    m = m3,
                    initRefPos = referencePath1.position[3],
                    trans = translation3,
                    rota = Par(value = :(rotation3))
)

linkParameters4 = Map(
                    parent1 = Par(value = :(link3.obj2)),
                    featureBody = featureBody4,
                    m = m4,
                    initRefPos = referencePath1.position[4],
                    trans = translation4,
                    rota = Par(value = :(nullRot))
)

linkParameters5 = Map(
                    parent1 = Par(value = :(link4.obj2)),
                    featureBody = featureBody5,
                    m = m5,
                    initRefPos = referencePath1.position[5],
                    trans = translation5,
                    rota = Par(value = :(rotation5))
)

Link = Model(
    parent1 = Par(value = :(nothing)),
    featureBody = featureBody1,
    trans = [0,0,0],
    rota = Par(value = :(nullRot)),

    obj1 = Object3D(parent=:parent1, rotation=:rota),
    body = Object3D(feature = :featureBody),
    obj2 = Object3D(parent =:body, translation = :trans),
)

Gripper = Model(
    obj1 = Object3D(parent=:(link5.obj2), rotation=Modia3D.rot3(-180u"°")),
    gripper_base_frame = Object3D(parent=:obj1, feature=Solid(shape=
        FileMesh(filename = gripper_base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.199))),
    gripper_left_finger_a = Object3D(),
    gripper_left_finger = Object3D(parent=:gripper_left_finger_a, translation=[0, 0.0082, 0],
        feature=Solid(shape= FileMesh(filename = gripper_left_finger_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.010), contactMaterial="DryWood", collision=!simplifiedContact)),
    gripper_left_finger_box = Object3D(parent=:gripper_left_finger_a, translation=[0.0, 0.0051, 0.0135],
        feature = Solid(shape=Box(lengthX=0.012, lengthY=0.01, lengthZ=0.045), massProperties=MassPropertiesFromShapeAndMass(mass=1.0e-9), visualMaterial=vmatInvisible, contactMaterial="DryWood", collision=simplifiedContact)),
    gripper_right_finger_a = Object3D(parent=:gripper_base_frame,
        translation=[0,-rightFingerOffset,0]),
    gripper_right_finger = Object3D(parent=:gripper_right_finger_a, translation=[0, -0.0082, 0],
        feature = Solid(shape = FileMesh(filename = gripper_right_finger_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.010), contactMaterial="DryWood", collision=!simplifiedContact)),
    gripper_right_finger_box = Object3D(parent=:gripper_right_finger_a, translation=[0.0, -0.0051, 0.0135],
        feature = Solid(shape=Box(lengthX=0.012, lengthY=0.01, lengthZ=0.045), massProperties=MassPropertiesFromShapeAndMass(mass=1.0e-9), visualMaterial=vmatInvisible, contactMaterial="DryWood", collision=simplifiedContact))
)

YouBot = Model(
    base = Base,
    arm_base_frame = Object3D(parent=:(base.base_frame),
        translation=[0.143, 0.0, 0.046],
        feature = Solid(shape = FileMesh(filename = arm_base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.961))),
    armBase_b = Object3D(parent=:arm_base_frame, translation=[0.024, 0, 0.115]),

    link1 = Link | linkParameters1,
    link2 = Link | linkParameters2,
    link3 = Link | linkParameters3,
    link4 = Link | linkParameters4,
    link5 = Link | linkParameters5,
    gripper = Gripper,

    rev1 = RevoluteWithFlange(obj1 = :(link1.obj1), obj2 = :(link1.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),
    rev2 = RevoluteWithFlange(obj1 = :(link2.obj1), obj2 = :(link2.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 2)), w=Var(init=0.0)),
    rev3 = RevoluteWithFlange(obj1 = :(link3.obj1), obj2 = :(link3.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 3)), w=Var(init=0.0)),
    rev4 = RevoluteWithFlange(obj1 = :(link4.obj1), obj2 = :(link4.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 4)), w=Var(init=0.0)),
    rev5 = RevoluteWithFlange(obj1 = :(link5.obj1), obj2 = :(link5.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 5)), w=Var(init=0.0)),

    prism = PrismaticWithFlange(obj1 = :(gripper.gripper_right_finger_a), obj2 = :(gripper.gripper_left_finger_a),
        axis=axisGripper, s = Var(init = getRefPathInitPosition(referencePath1, 6)) ),

    servo1 = Servo,
    servo2 = Servo,
    servo3 = Servo,
    servo4 = Servo,
    servo5 = Servo,
    servo6 = ServoTrans,

    equations=:[
        refPath = calculateRobotMovement(getReferencePath(), instantiatedModel),
        servo1.refLoadAngle = getRefPathPosition(refPath, 1),
        servo2.refLoadAngle = getRefPathPosition(refPath, 2),
        servo3.refLoadAngle = getRefPathPosition(refPath, 3),
        servo4.refLoadAngle = getRefPathPosition(refPath, 4),
        servo5.refLoadAngle = getRefPathPosition(refPath, 5),
        servo6.refLoadPos   = getRefPathPosition(refPath, 6)
    ],

    connect = :[
        (servo1.flange, rev1.flange)
        (servo2.flange, rev2.flange)
        (servo3.flange, rev3.flange)
        (servo4.flange, rev4.flange)
        (servo5.flange, rev5.flange)
        (servo6.flange, prism.flange)
        ]
)

Scenario = Model(
    gravField = UniformGravityField(g=9.81, n=[0,0,-1]),
    world = Object3D(feature=Scene(gravityField=:gravField, visualizeFrames=false, nominalLength=tableX,
    animationFile="YouBotGripping.json",
    enableContactDetection=true, elasticContactReductionFactor=1e-3)),
#   worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.2))),

    table = Table,
    ground = Ground,

    sphere = Object3D(
        feature=Solid(shape=Sphere(diameter=diameter), visualMaterial=vmatGrey, solidMaterial="DryWood", collision=true)),

    sphereJoint = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=[-0.799, 0.0, 0.1792])),

    youbot1 = YouBot,
)

modelParameters = Map(
    youbot1 = Map(
        servo1 = servoParameters1,
        servo2 = servoParameters2,
        servo3 = servoParameters3,
        servo4 = servoParameters4,
        servo5 = servoParameters5,
        servo6 = servoParameters6
    )
)

youbotModel = buildModia3D(Scenario) | modelParameters

youbot = @instantiateModel(youbotModel, unitless=true, logCode=false, log=false)

stopTime = 14.0
tolerance = 1e-6
if simplifiedContact
    requiredFinalStates = [-0.0033893963354342826, -0.005096954762027042, 0.38499191688168316, -1.445830218278958e-8, -2.7694145227833314e-8, -7.258566222266568e-13, 0.36864138584111816, -0.9931096002408847, 3.4546855331645943, -5.08129867112521e-7, 6.65198823908264e-7, -7.189333025824487e-7, -8.485351612092523e-7, 8.489451437830829e-7, -1.048206504157684e-5, 1.048341761384057e-5, -2.3863801845061742e-5, 2.3866297886359383e-5, -1.4665670255696616e-5, 1.4668022798963407e-5, 1.4943381828981332e-7, -1.4944863129434576e-7, -0.0064830265810847665, 0.16152090176322154, -0.1177121839079235, -0.05195318968695996, 0.0005198281051871344, -0.014745897597333173, 0.02969913219796793, 0.00030080034901011366]
else
    requiredFinalStates = [0.0009285649671864051, -0.00510121706010873, 0.38499183731994097, 1.3815799549216646e-8, -3.013883724849587e-8, -1.3548815880513574e-12, 0.3728261933992865, -0.996577976084585, 3.458983136911289, -5.493765692434776e-7, -2.9787026664428377e-7, -9.903133313245147e-7, -8.815312501299099e-7, 8.819099920244317e-7, -1.048429764456881e-5, 1.0485482464956499e-5, -2.38693150412509e-5, 2.3877893868562615e-5, -1.4676089351824731e-5, 1.467815727005773e-5, 1.5116246063921994e-7, -1.5117524141251032e-7, -0.00673525336167914, 0.1615038110919147, -0.1177385926739786, -0.05198945600575353, 0.0005258421886067691, -0.01468148168734971, 0.029700446414509994, 0.000299495164936174]
end
simulate!(youbot, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(youbot, ["sphereJoint.r"], figure=1)

end


text:
using Soss

m = @model σ begin
           μ ~ StudentT(3.0)
           x ~ Normal(μ,σ) ^ 10
           return x
       end;

x = rand(m(σ=3.0));

sourceSymlogdensityof(m(σ=3.0) | (;x))

s = symlogdensity(m(σ=3.0))

codegen(m(σ=3.0) | (;x))


text:
## generic operations on numbers ##

isreal(x::Real) = true
isinteger(x::Integer) = true

size(x::Number) = ()
size(x::Number,d) = convert(Int,d)<1 ? throw(BoundsError()) : 1
eltype(x::Number) = typeof(x)
eltype{T<:Number}(::Type{T}) = T
ndims(x::Number) = 0
ndims{T<:Number}(::Type{T}) = 0
length(x::Number) = 1
endof(x::Number) = 1
getindex(x::Number) = x
getindex(x::Number, i::Integer) = i == 1 ? x : throw(BoundsError())
getindex(x::Number, I::Integer...) = all([i == 1 for i in I]) ? x : throw(BoundsError())
getindex(x::Number, I::Real...) = getindex(x, to_index(i)...)
first(x::Number) = x
last(x::Number) = x

divrem(x,y) = (div(x,y),rem(x,y))
signbit(x::Real) = x < 0
sign(x::Real) = ifelse(x < 0, oftype(x,-1), ifelse(x > 0, one(x), x))
abs(x::Real) = ifelse(x < 0, -x, x)
abs2(x::Real) = x*x
copysign(x::Real, y::Real) = ifelse(signbit(x)!=signbit(y), -x, x)

conj(x::Real) = x
transpose(x::Number) = x
ctranspose(x::Number) = conj(x)
inv(x::Number) = one(x)/x
angle(z::Real) = atan2(zero(z), z)

widemul(x::Number, y::Number) = widen(x)*widen(y)

start(x::Number) = false
next(x::Number, state) = (x, true)
done(x::Number, state) = state
isempty(x::Number) = false
in(x::Number, y::Number) = x == y

reinterpret{T<:Real}(::Type{T}, x::Real) = box(T,x)

map(f::Callable, x::Number) = f(x)

zero(x::Number) = oftype(x,0)
zero{T<:Number}(::Type{T}) = oftype(T,0)
one(x::Number)  = oftype(x,1)
one{T<:Number}(::Type{T}) = oftype(T,1)

const _numeric_conversion_func_names =
    (:int,:integer,:signed,:int8,:int16,:int32,:int64,:int128,
     :uint,:unsigned,:uint8,:uint16,:uint32,:uint64,:uint128,
     :float,:float16,:float32,:float64)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
@testset "$(rpad("Gauß-Legendre",80))" begin

    for s in 1:10
        @test GaussLegendreQuadrature(s) ≈  GaussLegendreQuadrature(s; fast=true)
        @test GaussLegendreQuadrature(s) == GaussLegendreQuadrature(Float64, s)

        @test sum(weights(GaussLegendreQuadrature(s))) ≈ 1
    end

end


text:
"""
    LogBayesQuad(k::Kernel; n_candidates=100, l=1.0, σ::Real=1.0)

Tool for running Bayesian Quadrature by assuming a GP on the log integrand instead
of the integrand.

## Arguments
- `n_candidates` : will be sampled around the first samples to perform a linear approximation of the result.
- `l` : Lengthscale of the kernel
- `σ` : variance of the kernel

## Reference
- **Active Learning of Model Evidence Using Bayesian Quadrature** - Osborne et al. - NIPS 2012.
"""
struct LogBayesQuad{TK,Tl,Tσ} <: AbstractBQ{TK}
    kernel::TK
    n_candidates::Int
    l::Tl
    σ::Tσ
end

function LogBayesQuad(k::Kernel; n_candidates=100, l=1.0, σ::Real=1.0)
    k, (l, σ) = get_kernel_params(k; l, σ)
    check_kernel_parameters(l, σ)
    return LogBayesQuad(k, n_candidates, l, σ)
end

scale(bquad::LogBayesQuad) = bquad.σ

function quadrature(
    bquad::LogBayesQuad{<:SqExponentialKernel}, model::AbstractBQModel{<:MvNormal}, samples
)
    isempty(samples) && error("The collection of samples is empty")
    nsamples = length(samples)
    logf = logintegrand(model).(samples) # Evaluate integrand on samples
    normalize = false
    if normalize
        max_logf = maximum(logf) # Find maximum of all given samples
        logf .-= max_logf # Normalize log integrand
    end
    f = exp.(logf) # Compute integrand on samples

    x_c = sample_candidates(bquad, samples, bquad.n_candidates) # Sample candidates around the samples
    ke = KernelEmbedding(bquad, p_0(model))

    gp = create_gp(bquad, samples)
    f_c_0 = mean.(predict(gp, f, x_c)) # Predict integrand on x_c
    logf_c_0 = mean.(predict(gp, logf, x_c)) # Predict log-integrand on x_c
    Δ_c = exp.(logf_c_0) - f_c_0 # Compute difference of predictions

    z = kernel_mean(ke, samples) # Compute mean for the basic BQ
    K = kernelpdmat(kernel(bquad), samples) # and the kernel matrix

    z_c = kernel_mean(ke, x_c) # Compute mean for the ΔlogBQ
    K_c = kernelpdmat(kernel(bquad), x_c) # and the kernel matrix for the candidates

    m_evidence = evaluate_mean(z, K, f) # Compute m(Z|samples)
    m_correction = evaluate_mean(z_c, K_c, Δ_c) #

    C = kernel_variance(ke) # Compute the kernel variance

    var_evidence = evaluate_var(z, K, C)
    var_correction = evaluate_var(z_c, K_c, C)

    if normalize
        m = exp(log(m_evidence + m_correction) + max_logf)
        v = exp(log(var_evidence + var_correction) + 2 * max_logf)
    else
        m = m_evidence + m_correction
        v = var_evidence + var_correction
    end
    return Normal(m, sqrt(v))
end

function predict(gp::AbstractGPs.FiniteGP, f, x_c)
    f_post = posterior(gp, f) # Create a posterior given the observations f
    return mean(f_post(x_c)) # Return the predictive mean on the candidates
end

function create_gp(bquad::LogBayesQuad, samples)
    gp = GP(kernel(bquad)) # Create a GP prior
    return gp(samples) # Project it on the samples
end

"""
    sample_candidates(bquad::LogBayesQuad, samples, n_c)

Sample new candidates around the existing `samples` within an hyper-ellipse, i.e
in the space where |x_i|^n + |y_i|^n = 1.
The number of candidates is given by `bquad.n_candidates`
This is needed to estimate the correction term for using a log transformation.
See Osborne et al. 2012 and [`LogBayesQuad`](@ref)
"""
function sample_candidates(bquad::LogBayesQuad, samples, n_c)
    n_samples = length(samples)
    n_dim = length(first(samples))
    T = eltype(first(samples))
    directions = [project_on_ellipse(rand(T, n_dim) .- 0.5) for _ in 1:n_c] # Sample withing the hyper ellipse
    directions = scale_directions(bquad.l, directions) # Scale the hyperellipse
    p_indices = rand(1:n_samples, n_c) # Get n_c indices from the samples
    return directions .+= samples[p_indices] # Create candidates around the existing samples
end

function project_on_ellipse(x)
    return x / norm(x, 1)
end

scale_directions(l::Real, directions) = l .* directions
scale_directions(l::AbstractVector, directions) = hadamard(Ref(l), directions)
scale_directions(l::AbstractMatrix, directions) = Ref(l) .* directions

hadamard(x, y) = x .* y


text:
export Tria
mutable struct Tria <: AbstractElement
    id::String
    hid::Int

    node1::Node
    node2::Node
    node3::Node
    material::Material
    t::Float64
    t2::Float64
    t3::Float64

    elm_type::String
    mass_type::String

    center::Array{Float64}
    A::Float64
    T::Matrix{Float64}
end

function Tria(id,hid,node1,node2,node3,material,t,t2=0,t3=0,elm_type="TMT",mass_type="concentrate")
    #Initialize local CSys,could be optimized by using a MSE plane
    o=(node1.loc+node2.loc+node3.loc)/3
    pt1 = (node1.loc+node2.loc)/2
    pt2 = (node2.loc+node3.loc)/2
    csys = CSys(o, pt1, pt2)

    #area is considered as the average of trangles generated by splitting the quand with diagonals
    A=det([[1;1;1] node2.loc-node1.loc node3.loc-node1.loc]')/2

    E=material.E
    ν=material.ν

    #3D to local 2D
    T=zeros(18,18)
    for i in 1:6
        T[3i-2:3i,3i-2:3i]=csys.T
    end

    Tria(id,hid,node1,node2,node3,material,t,t2,t3,elm_type,mass_type,o,A,T)
end

for (root,dirs,files) in walkdir(joinpath(@__DIR__,"trias"))
    for file in files
        if file[end-2:end]==".jl"
            include(joinpath(@__DIR__,"trias",file))
        end
    end
end

function integrateK(elm::Tria)
    if elm.elm_type=="GT9"
        return K_GT9(elm)
    elseif elm.elm_type=="TMT"
        return K_TMT(elm)
    # elseif elm.elm_type=="DKT"
    #     K=K_DKT(elm)
    # elseif elm.elm_type=="DKGT"
    #     K=K_DKGT(elm)
    elseif elm.elm_type=="TMGT"
        return K_TMGT(elm)
    else
        throw("Quad element type error!")
    end
end

#WIP
function integrateKσ(elm::Tria,σ::Vector{Float64})
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V
end

#WIP
function integrateKu(elm::Tria,u)
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    n4=node4.loc+[u[7];u[8];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V
end

function integrateM!(elm::Tria)
    membrane,bending=elm.membrane,elm.bending
    ρ=elm.material.ρ
    center=elm.center
    t=elm.t
    A=elm.A
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    # M=Matrix{Float64}(ρ*t*A/18,18,18)
    M=spzeros(18,18)
    elm.Mᵉ=M
end

function integrateP(elm::Tria,elm_force)
    N=elm.N
    B=elm.B
    D=elm.D
    J=elm.J

    f,s,σ₀,ϵ₀=elmforce.f,elmforce.s,elmforce.σ₀,elmforce.ϵ₀
    a=[-1 -1]
    b=[1 1]

    function f1(x)
        N'*f*det(J)
    end

    function f2(x)
        N'*s*det(J)
    end

    function f3(x)
        x->B'*σ₀*det(J)
    end

    function f4(x)
        x->B'*D*ϵ₀*det(J)
    end

    Pᵉ_f=hquadrature(f1,a,b)
    Pᵉ_s=hquadrature(f2,[-1],[1])#ξ=-1
    Pᵉ_σ₀=hquadrature(f3,a,b)
    Pᵉ_ϵ₀=hquadrature(f4,a,b)

    Pᵉ_f,Pᵉ_s,Pᵉ_σ₀,Pᵉ_ϵ₀=calc_force(elm,elm_force) #Pᵉ=Pᵉf+Pᵉs+Pᵉσ₀+Pᵉϵ₀
    Pᵉ=Pᵉ_f+Pᵉ_s+Pᵉ_σ₀+Pᵉ_ϵ₀
end

# end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    FractionPartitioner(fraction, shuffle=true)

A method for partitioning spatial objects according to a given `fraction`.
Optionally `shuffle` elements before partitioning.
"""
struct FractionPartitioner <: AbstractPartitioner
  fraction::Float64
  shuffle::Bool

  function FractionPartitioner(fraction, shuffle)
    @assert 0 < fraction < 1 "fraction must be in interval (0,1)"
    new(fraction, shuffle)
  end
end

FractionPartitioner(fraction) = FractionPartitioner(fraction, true)

function partition(object, p::FractionPartitioner)
  n = nelms(object)
  f = round(Int, p.fraction * n)

  locs = p.shuffle ? randperm(n) : 1:n
  subsets = [locs[1:f], locs[f+1:n]]

  SpatialPartition(object, subsets)
end


text:
export fact2, member, allequal

"""
    fact2(n)

Compute the double factorial of `n`, i.e. the product of the odd (if `n` is odd) or even (if `n` is even) numbers up to `n`.
"""
function fact2(n::Int64)
    if n == 1 || n==0
        return 1;
    else
        if mod(n,2) == 0 # the number is even
            return prod(2:2:n);
        else # the number is odd
            return prod(1:2:n);
        end
    end
end

"""
    member(X, idx)

Extract the `idx` ensemble member of the ensemble matrix `X`
"""
member(X, idx) = X[:,idx]

"""
    allequal(x, r)

A function to verify that all the entries of the vector `x` are equal to `r`
"""
allequal(x, r) = all(y->y==r, x)