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

using Markdown
using InteractiveUtils

# ╔═╡ 3f5d8958-31e5-11eb-1603-55006552e57e
md" ## _Functions_
Functions are the actions done to the object. If given an input, the function does something to create the output. Functions are very important because they are the foundation of Julia.

"

# ╔═╡ c56bfffc-31e5-11eb-3dd3-0180c02087f5
md" ## Defining Functions
A function begins with `function` and ends with `end`. This is important as it clearly dictates when the function starts and finishes.
"

# ╔═╡ 4223b408-31e7-11eb-3e6b-4f63b0a05b11
function plus_two(x) #overview of the function
	x=x+2 #the action being performed
	return x #what gets returned
end

# ╔═╡ 7d46403c-31e7-11eb-0852-038f465d05e7
md" Another way of writing that same function is to assign the function directly to the action like this: `plus_two(x) = x+2`. Both yield the same result.
"

# ╔═╡ d0648672-31e5-11eb-064a-1b44500d6ced
md" ## Anonymous Functions

"

# ╔═╡ e400155c-31e5-11eb-3b15-dbaff73c7de4
md" ## Void Functions

"

# ╔═╡ eda66264-31e5-11eb-1b1a-41f28a927563
md" ## Optional Positional Arguments

"

# ╔═╡ fc7ed19a-31e5-11eb-0657-13f8cb1e6826
md" ## Key Word Arguments

"

# ╔═╡ 085c37c8-31e6-11eb-256b-79b193b33094
md" ## Function Documentation

"

# ╔═╡ Cell order:
# ╠═3f5d8958-31e5-11eb-1603-55006552e57e
# ╟─c56bfffc-31e5-11eb-3dd3-0180c02087f5
# ╠═4223b408-31e7-11eb-3e6b-4f63b0a05b11
# ╠═7d46403c-31e7-11eb-0852-038f465d05e7
# ╠═d0648672-31e5-11eb-064a-1b44500d6ced
# ╠═e400155c-31e5-11eb-3b15-dbaff73c7de4
# ╠═eda66264-31e5-11eb-1b1a-41f28a927563
# ╠═fc7ed19a-31e5-11eb-0657-13f8cb1e6826
# ╠═085c37c8-31e6-11eb-256b-79b193b33094


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

using Markdown
using InteractiveUtils

# ╔═╡ ac5c6510-3901-11eb-0c47-bb3143471434
begin
    import Pkg
	#Pkg.activate("Mage")

    using Colors
	using PlaymatSimulator
	using Serialization

	Image = ImageActor

	gd = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage"
	US = deserialize("$gd/../../PlaymatSimulator/tmp/user_selection.jls")
	GS = ingredients("$gd/Base/notebooks/game_settings.jl")

	SCREEN_WIDTH = GS.SCREEN_WIDTH
	SCREEN_HEIGHT = GS.SCREEN_HEIGHT
	SCREEN_BORDER = GS.SCREEN_BORDER

	md"""
	## GAME STAGE
	This notebook should define the `stage` of the game. The `stage` should consist of elements to be drawn below the cards and other ui elements such as background images/ANations, some ui elements, game zone/area markers, etc. Provide those objects to the game engine by defining a single dictionary object, `stage = Dict{Symbol, Any}(...)`.
	"""
end

# ╔═╡ fdf59b16-4da2-11eb-204f-0d0e1a522f6c
const SHADE_PATH = "$gd/Base/ui/zones/area_blk.png"

# ╔═╡ a2195ae4-4da2-11eb-2112-e54a53878565
const zone_shade = load(SHADE_PATH)

# ╔═╡ 235edb4c-38fe-11eb-3112-bde291f6f5b5
begin

	STAGE = OrderedDict(
		:background => Image("$(GS.BKG_NAME)",
			load(GS.BKG_PATH),
			w=SCREEN_WIDTH, h=SCREEN_HEIGHT
			),
		"Library" => Image("Library",
			load(SHADE_PATH),
			x=SCREEN_BORDER,
			y=ceil(Int32, SCREEN_HEIGHT * 0.6) + SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15) - SCREEN_BORDER,
			h=ceil(Int32, SCREEN_HEIGHT * 0.4 - 2SCREEN_BORDER),
			alpha=50,
			),
		"Battlefield" => Image("Battlefield",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.15) + SCREEN_BORDER,
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.7) - SCREEN_BORDER,
			h=ceil(Int32, SCREEN_HEIGHT - 2SCREEN_BORDER),
			alpha=50,
			),
		"Command" => Image("Command",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.85 + SCREEN_BORDER),
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - 2SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.4 - SCREEN_BORDER),
			alpha=50,
			),
		"Graveyard" => Image("Graveyard",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.85 + SCREEN_BORDER),
			y=ceil(Int32, SCREEN_HEIGHT * 0.4 + SCREEN_BORDER),
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - 2SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.6 - 2SCREEN_BORDER),
			alpha=50,
			),
		"Hand" => Image("Hand",
			load(SHADE_PATH),
			x=SCREEN_BORDER,
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.6 - SCREEN_BORDER),
			alpha=50,
			),
		)
end

# ╔═╡ Cell order:
# ╠═ac5c6510-3901-11eb-0c47-bb3143471434
# ╠═fdf59b16-4da2-11eb-204f-0d0e1a522f6c
# ╟─a2195ae4-4da2-11eb-2112-e54a53878565
# ╠═235edb4c-38fe-11eb-3112-bde291f6f5b5

ext: jl


text:
function ueberpruefeMatrixVektorb(A,m)

Akorrekt = [1 0 0 1 0 0 1 0 0;
            0 1 0 0 1 0 0 1 0;
            0 0 1 0 0 1 0 0 1;
            1 1 1 0 0 0 0 0 0;
            0 0 0 1 1 1 0 0 0;
            0 0 0 0 0 0 1 1 1;
            sqrt(2) 0 0 0 sqrt(2) 0 0 0 sqrt(2)];

mkorrekt = [13 8 3 6 7 11 3*sqrt(2)];

# homogene Lösungen berechnen
homogenkorrekt = nullspace(Akorrekt);

# inhomogene Lösung berechnen
inhomogenkorrekt = pinv(Akorrekt)*mkorrekt';

homogen = nullspace(A);
inhomogen = pinv(A)*m';

u1 = homogen[:,1];
u2 = homogen[:,2];
u3 = homogen[:,3];
p = inhomogen;

Vergleich = Akorrekt'*mkorrekt';
VergleichSuS = A'*m';

if round(norm(Vergleich);digits=2) == round(norm(VergleichSuS);digits=2)
    println("Eure Matrix und euer Vektor sind korrekt!")
    println("Die Lösung eures Gleichungssystems lautet wie folgt, wobei a,b und c beliebige relle Zahl ist.")
    q = round.(p;digits=2);
    v1 = round.(u1;digits=2);
    v2 = round.(u2;digits=2);
    v3 = round.(u3;digits=2);
    print("Für f1 lautet die Lösung: ", q[1], " + a * ", v1[1], " + b * ", v2[1], " + c * ", v3[1],"\n")
    print("Für f2 lautet die Lösung: ", q[2], " + a * ", v1[2], " + b * ", v2[2], " + c * ", v3[2],"\n")
    print("Für f3 lautet die Lösung: ", q[3], " + a * ", v1[3], " + b * ", v2[3], " + c * ", v3[3],"\n")
    print("Für f4 lautet die Lösung: ", q[4], " + a * ", v1[4], " + b * ", v2[4], " + c * ", v3[4],"\n")
    print("Für f5 lautet die Lösung: ", q[5], " + a * ", v1[5], " + b * ", v2[5], " + c * ", v3[5],"\n")
    print("Für f6 lautet die Lösung: ", q[6], " + a * ", v1[6], " + b * ", v2[6], " + c * ", v3[6],"\n")
    print("Für f7 lautet die Lösung: ", q[7], " + a * ", v1[7], " + b * ", v2[7], " + c * ", v3[7],"\n")
    print("Für f8 lautet die Lösung: ", q[8], " + a * ", v1[8], " + b * ", v2[8], " + c * ", v3[8],"\n")
    print("Für f9 lautet die Lösung: ", q[9], " + a * ", v1[9], " + b * ", v2[9], " + c * ", v3[9],"\n")
    #show(stdout,"text/plain",round.(p;digits=2))
    #show(stdout,"text/plain",round.(u1;digits=2))
    #show(stdout,"text/plain",round.(u2;digits=2))
    #show(stdout,"text/plain",round.(u3;digits=2))
    #println("p = $p")
    #println("u1 = $u1")
    #println("u2 = $u2")
    #println("u3 = $u3")
else
    println("Schaut euch noch einmal den Verlauf der Strahlen an und überdenkt eure Matrix.")
end

end


text:
using Knet: @diff, Param, value, grad
using Knet: sigm, tanh

mutable struct Layer
    wi1::Param
    wi2::Param
    wis::Param
    wk1::Param
    wk2::Param
    wks::Param
    ws1::Param
    ws2::Param
    wss::Param
end

Layer(in_size1, in_size2, layer_size) =
begin
    sq = sqrt(2/((in_size1+in_size2)/2+layer_size))
    wi1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wi2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wis = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    wk1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wk2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wks = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    ws1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    ws2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wss = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    layer = Layer(wi1, wi2, wis, wk1, wk2, wks, ws1, ws2, wss)
layer
end

(layer::Layer)(state, in_1, in_2) =
begin
    state_sq = tanh.(state)
    interm = tanh.(in_1 * layer.wi1 + in_2 * layer.wi2 + state_sq * layer.wis)
    keep   = sigm.(in_1 * layer.wk1 + in_2 * layer.wk2 + state_sq * layer.wks)
    show   = sigm.(in_1 * layer.ws1 + in_2 * layer.ws2 + state_sq * layer.wss)
    state += keep .* interm
    out    = show .* interm
[out, state]
end

mutable struct LayerS
    wi1::Param
    wi2::Param
    wis::Param
    ws1::Param
    ws2::Param
    wss::Param
end

LayerS(in_size1, in_size2, layer_size) =
begin
    sq = sqrt(2/((in_size1+in_size2)/2+layer_size))
    wi1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wi2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wis = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    ws1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    ws2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wss = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    layer = LayerS(wi1, wi2, wis, ws1, ws2, wss)
layer
end

(layer::LayerS)(state, in_1, in_2) =
begin
    show   = sigm.(in_1 * layer.ws1 + in_2 * layer.ws2 + state * layer.wss)
    out    = show .* state
    state  = tanh.(in_1 * layer.wi1 + in_2 * layer.wi2 + state * layer.wis)
[out, state]
end

mutable struct LayerX
    wf1::Param
    wf2::Param
    wfs::Param
    wi1::Param
    wi2::Param
    wd1::Param
    wd2::Param
    wds::Param
    wk1::Param
    wk2::Param
    wks::Param
    ws1::Param
    ws2::Param
    wss::Param
end

LayerX(in_size1, in_size2, layer_size) =
begin
    sq = sqrt(2/(in_size1+in_size2+layer_size))
    wf1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wf2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wfs = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    wi1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wi2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wd1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wd2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wds = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    wk1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wk2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wks = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    ws1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    ws2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wss = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    layerX = LayerX(wf1, wf2, wfs, wi1, wi2, wd1, wd2, wds, wk1, wk2, wks, ws1, ws2, wss)
layerX
end

(layer::LayerX)(state, in_1, in_2) =
begin
    st_sq   = tanh.(state)

    focus   = sigm.(in_1 * layer.wf1 + in_2 * layer.wf2 + st_sq * layer.wfs)
    interm  = tanh.(in_1 * layer.wi1 + in_2 * layer.wi2 + st_sq .* focus)

    discard = sigm.(in_1 * layer.wd1 + in_2 * layer.wd2 + st_sq * layer.wds)
    keep    = sigm.(in_1 * layer.wk1 + in_2 * layer.wk2 + st_sq * layer.wks)
    state  -= discard .* state
    state  += keep .* interm

    st_sq   = tanh.(state)
    show    = sigm.(in_1 * layer.ws1 + in_2 * layer.ws2 + st_sq * layer.wss)
    out     = show .* st_sq

[out, state]
end

struct IS
    layer_1::Layer
    layer_out::Layer
end

IS(storage_size, hm_vectors, vector_size, layer_sizes, out_size) =
begin
    is = IS(
        Layer(storage_size, hm_vectors*vector_size, layer_sizes[1]),
        Layer(storage_size, layer_sizes[end], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[is, state]
end

(is::IS)(states, in_1, in_2) =
begin
    out_1, state_1 = is.layer_1(states[1], in_1, in_2)
    out, state_out = is.layer_out(states[end], in_1, out_1)
[out, [state_1, state_out]]
end

struct GS
    layer_1::Layer
    layer_out::Layer
end

GS(in_size1, in_size2, layer_sizes, out_size) =
begin
    gs = GS(
        Layer(in_size1, in_size2, layer_sizes[1]),
        Layer(in_size1, layer_sizes[end], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[gs, state]
end

(gs::GS)(states, in_1, in_2) =
begin
    out_1, state_1 = gs.layer_1(states[1], in_1, in_2)
    out, state_out = gs.layer_out(states[end], in_1, out_1)
[out, [state_1, state_out]]
end

struct GO
    layer_1::Layer
    layer_out::Layer
end

GO(in_size1, in_size2, layer_sizes, out_size) =
begin
    go = GO(
        Layer(in_size1, in_size2, layer_sizes[1]),
        Layer(in_size1, layer_sizes[end], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[go, state]
end

(go::GO)(states, in_1, in_2) =
begin
    out_1, state_1 = go.layer_1(states[1], in_1, in_2)
    out, state_out = go.layer_out(states[end], in_1, out_1)
[out, [state_1, state_out]]
end

struct ENCODER
    is::IS
    gs::GS
    go::GO
end

ENCODER(hm_vectors, vector_size, storage_size, is_layers, gs_layers, go_layers) =
begin
    is, is_state = IS(storage_size, hm_vectors, vector_size, is_layers, storage_size)
    gs, gs_state = GS(storage_size, storage_size, gs_layers, storage_size)
    go, go_state = GO(storage_size, vector_size, go_layers, vector_size)
    encoder = ENCODER(is, gs, go)
    state   = [is_state,gs_state,[go_state for i in 1:hm_vectors]]
[encoder, state]
end

(encoder::ENCODER)(state, storage, vectors) =
begin
    is_out, is_state = encoder.is(state[1], storage, hcat(vectors...))
    gs_out, gs_state = encoder.gs(state[2], storage, is_out)
    go_out, go_state = [], []
    for (vector, state_v) in zip(vectors, state[3])
        go_out_v, go_state_v = encoder.go(state_v, is_out, vector)
        push!(go_out, go_out_v)
        push!(go_state, go_state_v)
    end
    output = [gs_out, go_out]
    state  = [is_state, gs_state, go_state]
[output, state]
end

struct DECODER
    is::IS
    gs::GS
    go::GO
end

DECODER(hm_vectors, vector_size, storage_size, is_layers, gs_layers, go_layers) =
begin
    is, is_state = IS(storage_size, hm_vectors, vector_size, is_layers, storage_size)
    gs, gs_state = GS(storage_size, storage_size, gs_layers, storage_size)
    go, go_state = GO(storage_size, vector_size, go_layers, vector_size)
    decoder = DECODER(is, gs, go)
    state   = [is_state,gs_state,[go_state for i in 1:hm_vectors]]
[decoder, state]
end

(decoder::DECODER)(state, storage, vectors, enc_storages, enc_vectors) =
begin
    attended = attend(storage, enc_storages, enc_vectors)
    is_out, is_state = decoder.is(state[1], storage, hcat(attended...))
    gs_out, gs_state = decoder.gs(state[2], storage, is_out)
    go_out, go_state = [], []
    for (vector, state_v) in zip(vectors, state[3])
        go_out_v, go_state_v = decoder.go(state_v, is_out, vector)
        push!(go_out, go_out_v)
        push!(go_state, go_state_v)
    end
    output = [gs_out, go_out]
    state  = [is_state, gs_state, go_state]
[output, state]
end

attend(dec_storage, enc_storages, enc_vectors) =
begin
    similarities = soft([sum(dec_storage .* enc_storage) for enc_storage in enc_storages])
sum([vec .* sim for (vec, sim) in zip(enc_vectors,similarities)])
end

propogate(enc, dec, enc_state, x, length_y; dec_state=nothing) =
begin
    storage_size = length(enc_state[2][end])
    enc_outs = [zeros(1, storage_size), nothing]
    enc_gs_time, enc_go_time = [], []
    for timestep in x
        enc_outs, enc_state = enc(enc_state, enc_outs[1], timestep)
        push!(enc_gs_time, enc_outs[1])
        push!(enc_go_time, enc_outs[2])
    end
    dec_outs = enc_outs
    if dec_state == nothing dec_state = enc_state end
    dec_go_time = []
    for timestep in 1:length_y
        dec_outs, dec_state = dec(dec_state, dec_outs[1], dec_outs[2], enc_gs_time, enc_go_time)
        push!(dec_go_time, dec_outs[2])
    end
dec_go_time
end

grads(result, encoder, decoder) =
begin
    grads = []
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    gradient = grad(result, getfield(layer, lfield))
                    # if gradient > 0
                    #     gradient = 1.0
                    # end
                    push!(grads, gradient)
                    # @info(mfield,nfield,lfield,norm(gradient))
                end
            end
        end
    end
grads
end

upd_rms!(encoder, decoder, grads, lr, accu_grads; alpha=0.9) =
begin
    i = 0
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    i +=1
                    accu_grads[i] = alpha .* accu_grads[i] + (1 - alpha) .* grads[i].^2
                    setfield!(layer, lfield, Param(getfield(layer, lfield) - grads[i].*lr./sqrt.(accu_grads[i].+1e-8)))
                end
            end
        end
    end
end

upd!(encoder, decoder, grads, lr) =
begin
    i = 0
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    i +=1
                    this_grad = grads[i]
                    if this_grad != nothing
                        setfield!(layer, lfield, Param(getfield(layer, lfield) - this_grad.*lr))
                    end
                end
            end
        end
    end
end

relu(x) = max(0.0,x)

soft(arr) = ((arr) -> arr/sum(arr))(exp.(arr))

sqe(out, y) = sum(out .- y) ^2

sequence_loss(out, y) =
    sum([sum(sum([sqe(out_e, y_e) for (out_e,y_e) in zip(out_t,y_t)])) for (out_t,y_t) in zip(out, y)])

drop(arr; rate=0.1) = # rand() < rate ? x : 0.0
begin
    for i in 1:length(arr)
        if rand() < rate
            array[i] *= 0
        end
    end
array
end

max_stars_repo_name: developersfox/GSTM

max_issues_repo_name: developersfox/GSTM

max_forks_repo_name: developersfox/GSTM


text:
module NVTX

using ..CUDA
using ..CUDA: libnvtx, @checked
using ..CUDA: CUstream, CUdevice, CUcontext, CUevent

using CEnum

using ExprTools

# core library
initialize_api() = return
include("libnvtx_common.jl")
include("libnvtx.jl")

include("highlevel.jl")

end


text:
function test_nt_props(nt, validator)
    propnames = propertynames(nt)
    @assert isempty(symdiff(keys(validator), propnames))
    for (pname, ptype) in validator
        @assert pname in propnames
        @assert getproperty(nt, pname) isa ptype
    end
    true
end

@testset "Config parser: $config" for config in CONFIG_FUNCTIONS
	cfg = mktemp() do path, io  # write and parse config file on-the-fly
               write(io, eval(config)())
			   flush(io)
			   parse_configuration(path)
		   end
    @test cfg isa NamedTuple

    ENVCONFIG_PROPS = Dict(:data_loader => Function,
                           :data_sampler => Function,
                           :id_key => Symbol,
                           :vectors_eltype => Type,
                           :searcher_configs => Vector,
                           :embedder_configs => Vector,
                           :config_path => String)
    @test test_nt_props(cfg, ENVCONFIG_PROPS)

    EMBEDDERSCONFIG_PROPS = Dict(:id =>String,
                                 :description => String,
                                 :language => String,
                                 :stem_words => Bool,
                                 :ngram_complexity => Int,
                                 :vectors => Symbol,
                                 :vectors_transform => Symbol,
                                 :vectors_dimension => Int,
                                 :embeddings_path => Union{Nothing, String},
                                 :embeddings_kind => Symbol,
                                 :doc2vec_method => Symbol,
                                 :glove_vocabulary => Union{Nothing, String},
                                 :oov_policy => Symbol,
                                 :embedder_kwarguments => Dict{Symbol, Any},
                                 :embeddable_fields => Union{Nothing, Vector{Symbol}},
                                 :text_strip_flags => UInt32,
                                 :sif_alpha => cfg.vectors_eltype,
                                 :borep_dimension => Int,
                                 :borep_pooling_function => Symbol,
                                 :disc_ngram => Int)
    for ec in cfg.embedder_configs
        @test test_nt_props(ec, EMBEDDERSCONFIG_PROPS)
    end

    SEARCHERCONFIG_PROPS = Dict(:id => String,
                                :id_aggregation => String,
                                :description => String,
                                :enabled => Vector{Bool},
                                :indexable_fields => Vector{Symbol},
                                :data_embedder => String,
                                :input_embedder => String,
                                :search_index => Symbol,
                                :search_index_arguments => Vector{Any},
                                :search_index_kwarguments => Dict{Symbol, Any},
                                :heuristic => Union{Nothing, Symbol},
                                :score_alpha => cfg.vectors_eltype,
                                :score_weight => cfg.vectors_eltype)
    for sc in cfg.searcher_configs
        @test test_nt_props(sc, SEARCHERCONFIG_PROPS)
    end
end


text:
# TO IMPROVE .... Perhaps via a Map?
# Do we have already this function in Gridap?
# What if there are no cells? I am assuming that there is at least one cell
# What if the number of dofs per field per cell is different among cells?
function _get_cell_dofs_field_offsets(uh::MultiFieldFEFunction)
  U = get_fe_space(uh)
  uh_dofs = get_cell_dof_values(uh)[1]
  nfields = length(U.spaces)
  dofs_field_offsets=Vector{Int}(undef,nfields+1)
  dofs_field_offsets[1]=1
  for i in 1:nfields
    dofs_field_offsets[i+1]=dofs_field_offsets[i]+length(uh_dofs.array[i])
  end
  dofs_field_offsets
end

function FESpaces._gradient(f,uh::MultiFieldFEFunction,fuh::DomainContribution)
  terms = DomainContribution()
  U = get_fe_space(uh)
  for trian in get_domains(fuh)
    g = FESpaces._change_argument(gradient,f,trian,uh)
    cell_u = lazy_map(DensifyInnerMostBlockLevelMap(),get_cell_dof_values(uh))
    cell_id = FESpaces._compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_gradient(g,cell_u,cell_id)
    monolithic_result=cell_grad
    blocks = [] # TO-DO type unstable. How can I infer the type of its entries?
    nfields = length(U.spaces)
    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
    for i in 1:nfields
      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
      block=lazy_map(x->view(x,view_range),monolithic_result)
      append!(blocks,[block])
    end
    cell_grad=lazy_map(BlockMap(nfields,collect(1:nfields)),blocks...)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

function FESpaces._jacobian(f,uh::MultiFieldFEFunction,fuh::DomainContribution)
  terms = DomainContribution()
  U = get_fe_space(uh)
  for trian in get_domains(fuh)
    g = FESpaces._change_argument(jacobian,f,trian,uh)
    cell_u = lazy_map(DensifyInnerMostBlockLevelMap(),get_cell_dof_values(uh))
    cell_id = FESpaces._compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_jacobian(g,cell_u,cell_id)
    monolithic_result=cell_grad
    blocks        = [] # TO-DO type unstable. How can I infer the type of its entries?
    blocks_coords = Tuple{Int,Int}[]
    nfields = length(U.spaces)
    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
    for j=1:nfields
      view_range_j=cell_dofs_field_offsets[j]:cell_dofs_field_offsets[j+1]-1
      for i=1:nfields
        view_range_i=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
        # TO-DO: depending on the residual being differentiated, we may end with
        #        blocks [i,j] full of zeros. I guess that it might desirable to early detect
        #        these zero blocks and use a touch[i,j]==false block in ArrayBlock.
        #        How can we detect that we have a zero block?
        block=lazy_map(x->view(x,view_range_i,view_range_j),monolithic_result)
        append!(blocks,[block])
        append!(blocks_coords,[(i,j)])
      end
    end
    cell_grad=lazy_map(BlockMap((nfields,nfields),blocks_coords),blocks...)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

function FESpaces._change_argument(
  op::typeof(jacobian),f,trian,uh::MultiFieldFEFunction)

  U = get_fe_space(uh)
  function g(cell_u)
    single_fields = GenericCellField[]
    nfields = length(U.spaces)
    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
    for i in 1:nfields
      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
      cell_values_field = lazy_map(a->view(a,view_range),cell_u)
      cf = CellField(U.spaces[i],cell_values_field)
      push!(single_fields,cf)
    end
    xh = MultiFieldCellField(single_fields)
    cell_grad = f(xh)
    cell_grad_cont_block=get_contribution(cell_grad,trian)
    bs = [cell_dofs_field_offsets[i+1]-cell_dofs_field_offsets[i] for i=1:nfields]
    lazy_map(DensifyInnerMostBlockLevelMap(),
             Fill(bs,length(cell_grad_cont_block)),
             cell_grad_cont_block)
  end
  g
end

function FESpaces._change_argument(
  op::typeof(gradient),f,trian,uh::MultiFieldFEFunction)

  U = get_fe_space(uh)
  function g(cell_u)
    single_fields = GenericCellField[]
    nfields = length(U.spaces)
    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
    for i in 1:nfields
      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
      cell_values_field = lazy_map(a->view(a,view_range),cell_u)
      cf = CellField(U.spaces[i],cell_values_field)
      push!(single_fields,cf)
    end
    xh = MultiFieldCellField(single_fields)
    cell_grad = f(xh)
    get_contribution(cell_grad,trian)
  end
  g
end

function Algebra.hessian(f::Function,uh::MultiFieldFEFunction)
  @notimplemented
end

#function FESpaces._change_argument(
#  op::typeof(hessian),f,trian,uh::MultiFieldFEFunction)
#
#  U = get_fe_space(uh)
#  function g(cell_u)
#    single_fields = GenericCellField[]
#    nfields = length(U.spaces)
#    for i in 1:nfields
#      cell_values_field = lazy_map(a->a.array[i],cell_u)
#      cf = CellField(U.spaces[i],cell_values_field)
#      cell_data = lazy_map(BlockMap((nfields,nfields),i),get_data(cf))
#      uhi = GenericCellField(cell_data,get_triangulation(cf),DomainStyle(cf))
#      push!(single_fields,uhi)
#    end
#    xh = MultiFieldCellField(single_fields)
#    cell_grad = f(xh)
#    get_contribution(cell_grad,trian)
#  end
#  g
#end


text:
# Addition
# To add, convert both decimals to the same exponent.
# (If the exponents are different, use the smaller exponent
# to make sure we're adding integers.)
function +(x::Decimal, y::Decimal)
    cx = (-1)^x.s * x.c * 10^max(x.q - y.q, 0)
    cy = (-1)^y.s * y.c * 10^max(y.q - x.q, 0)
    s = (abs(cx) > abs(cy)) ? x.s : y.s
    c = BigInt(cx) + BigInt(cy)
    norm(Decimal(s, abs(c), min(x.q, y.q)))
end
+(x::Decimal, y::Number) = +(x, decimal(y))
+(x::Number, y::Decimal) = +(decimal(x), y)
+{T<:Number}(x::Union{Number, Array{T}}, y::Array{Decimal}) = +(decimal(x), y)
+{T<:Number}(x::Array{Decimal}, y::Union{Number, Array{T}}) = +(x, decimal(y))

# Negation
-(x::Decimal) = Decimal((x.s == 1) ? 0 : 1, x.c, x.q)

# Subtraction
-(x::Decinum, y::Decinum) = +(x, -y)
-{T<:Number}(x::Union{Number, Array{T}}, y::Array{Decimal}) = broadcast(-, x, y)
-{T<:Number}(x::Array{Decimal}, y::Union{Number, Array{T}}) = broadcast(-, x, y)

# Multiplication
function *(x::Decimal, y::Decimal)
    s = (x.s == y.s) ? 0 : 1
    norm(Decimal(s, BigInt(x.c) * BigInt(y.c), x.q + y.q))
end
*(x::Decimal, y::Number) = *(x, decimal(y))
*(x::Number, y::Decimal) = *(decimal(x), y)
.*(x::Decimal, y::Array{Decimal}) = [i * x for i in y]
.*(x::Array{Decimal}, y::Decimal) = [i * y for i in x]
.*(x::Number, y::Array{Decimal}) = decimal(x) .* y
.*(x::Array{Decimal}, y::Number) = x .* decimal(y)
.*{T<:Number}(x::Union{Number, Array{T}}, y::Array{Decimal}) = broadcast(*, decimal(x), y)
.*{T<:Number}(x::Array{Decimal}, y::Union{Number, Array{T}}) = broadcast(*, x, decimal(y))
.*(x::Array{Decimal}, y::BitArray) = x .* int(y)
.*(x::Array{Decimal}, y::Union{BitArray, Array{Bool}}) = x .* int(y)
.*(x::Union{BitArray, Array{Bool}}, y::Array{Decimal}) = int(x) .* y

# Inversion
function Base.inv(x::Decimal)
    str = string(x)
    if str[1] == '-'
        str = str[2:end]
    end
    b = ('.' in str) ? length(split(str, '.')[1]) : 0
    c = round(BigInt(10)^(-x.q + DIGITS) / x.c)
    q = (x.q < 0) ? 1 - b - DIGITS : -b - DIGITS
    norm(Decimal(x.s, c, q))
end
Base.inv(x::Array{Decimal}) = map(inv, x)

# Division
/(x::Decimal, y::Decimal) = x * inv(y)
/(x::Decinum, y::Decimal) = decimal(x) * inv(y)
/(x::Decimal, y::Decinum) = x * inv(decimal(y))
./(x::Decimal, y::Array{Decimal}) = [x / i for i in y]
./(x::Array{Decimal}, y::Decimal) = [i / y for i in x]
./(x::Number, y::Array{Decimal}) = decimal(x) ./ y
./(x::Array{Decimal}, y::Number) = x ./ decimal(y)
./{T<:Number}(x::Union{Number, Array{T}}, y::Array{Decimal}) = broadcast(/, decimal(x), y)
./{T<:Number}(x::Array{Decimal}, y::Union{Number, Array{T}}) = broadcast(/, x, decimal(y))

# TODO exponentiation
# TODO max, min
# TODO inequalities


text:
using Test
using Sockets
using MbedTLS

function sslaccept(server, certfile, keyfile)
    sslconfig = MbedTLS.SSLConfig(certfile, keyfile)
    conn = accept(server)
    sslconn = MbedTLS.SSLContext()
    MbedTLS.setup!(sslconn, sslconfig)
    MbedTLS.associate!(sslconn, conn)
    MbedTLS.handshake!(sslconn)
    return sslconn
end

function sslconnect(dest, port)
    conn = connect(dest, port)
    sslconfig = MbedTLS.SSLConfig(false)
    sslconn = MbedTLS.SSLContext()
    MbedTLS.setup!(sslconn, sslconfig)
    MbedTLS.set_bio!(sslconn, conn)
    MbedTLS.handshake!(sslconn)
    return sslconn
end

function testclntsrvr(certfile, keyfile)
    outbuff = ones(UInt8, 100) * UInt8(65)
    trigger = Channel{Bool}(1)
    port = UInt16(0)

    @async begin
        (port, server) = listenany(8000)
        @info("listening on port $port")
        put!(trigger, true)
        srvrconn = sslaccept(server, certfile, keyfile)
        inbuff = read(srvrconn, 100)
        @test inbuff == outbuff
        put!(trigger, true)
    end

    take!(trigger)
    @info("connecting to port $port")
    clntconn = sslconnect("127.0.0.1", port)
    @test write(clntconn, outbuff) == 100
    @async begin
        sleep(10)
        put!(trigger, false)
    end

    @test take!(trigger)
end

testclntsrvr(joinpath(dirname(@__FILE__), "test.cert"), joinpath(dirname(@__FILE__), "test.key"))


text:
export phimv, phimv!

"""
    phimv{T}(t, A, u, vec; [tol], [m], [norm], [anorm])

Calculate the solution of a nonhomogeneous linear ODE problem with constant input
``w = e^{tA}v + tφ(tA)u`` using the Krylov subspace approximation.

# Input

- `A`   -- matrix which can be dense or sparse
- `u`   -- vector, constant input of the ODE
- `vec` -- vector on which the matrix exponential of ``A`` is applied
- `tol` -- (optional, default: 1e-7) the requested accuracy tolerance on ``w``
- `m`   -- (optional, default: min(30, n)) maximum size of the Krylov basis

---

    [w, err] = phiv( t, A, u, v, tol, m )
    PHIV computes an approximation of ``w = exp(tA)v + t φ(tA)u``
    for a general matrix A using Krylov subspace projection techniques.
    Here, ``φ(z) = (\\exp(z)-1)/z`` and ``w`` is the solution of the
    nonhomogeneous linear ODE problem ``w' = Aw + u``, ``w(0) = v``.
    It does not compute the matrix functions in isolation but instead,
    it computes directly the action of these functions on the
    operand vectors. This way of doing so allows for addressing large
    sparse problems. The matrix under consideration interacts only
    via matrix-vector products (matrix-free method).

    Roger B. Sidje (rbs@maths.uq.edu.au)
    EXPOKIT: Software Package for Computing Matrix Exponentials.
    ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
"""
function phimv( t::Number,
                   A, u::Vector{T}, vec::Vector{T};
                   tol::Real=1e-7,
                   m::Int=min(30, size(A, 1)),
                   norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

    result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
    phimv!(t, A, u, result; tol=tol, m=m, norm=norm, anorm=anorm)
    return result
end

phimv!( t::Number,
           A, u::Vector{T}, vec::Vector{T};
           tol::Real=1e-7,
           m::Int=min(30, size(A, 1)),
           norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T} = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm)

function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
                   tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

    if size(vec, 1) != size(A, 2)
        error("dimension mismatch")
    end

    # safety factors
    gamma = 0.9
    delta = 1.2

    btol = 1e-7     # tolerance for "happy-breakdown"
    maxiter = 10    # max number of time-step refinements

    rndoff = anorm*eps()

    # estimate first time-step and round to two significant digits
    #beta = norm(A*vec + u)
    avtmp = similar(u)
    mul!(avtmp, A, vec)
    avtmp .+= u
    beta = norm(avtmp)
    r = 1/m
    fact = (((m+1)/exp(1))^(m+1))*sqrt(2*pi*(m+1))
    tau = (1.0/anorm)*((fact*tol)/(4.0*beta*anorm))^r
    tau = round(tau, sigdigits=2)

    # storage for Krylov subspace vectors
    vm = Array{typeof(w)}(undef,m+1)
    for i = 1:m+1
        vm[i] = similar(w)
    end
    hm = zeros(T, m+3, m+3)

    tf = abs(t)
    tsgn = sign(t)
    tk = zero(tf)

    copyto!(w, vec)
    p = similar(w)
    k1 = 3
    mb = m

    while tk < tf
        tau = min(tf-tk, tau)

        # Arnoldi procedure
        # rmul!(copyto!(vm[1], A*w+u), 1/beta)
        rmul!(copyto!(vm[1], avtmp), 1/beta)

        for j = 1:m
            mul!(p, A, vm[j])

            for i = 1:j
                hm[i, j] = dot(vm[i], p)
                p = axpy!(-hm[i,j], vm[i], p)
            end
            s = norm(p)

            if s < btol # happy-breakdown
                tau = tf - tk
                k1 = 0
                mb = j
                break
            end
            hm[j+1,j] = s
            rmul!(copyto!(vm[j+1], p), 1/hm[j+1,j])
        end

        hm[1, mb+1] = one(T)
        if k1 != 0
            hm[m+1,m+2] = one(T); hm[m+2,m+3] = one(T)
            h = hm[m+1,m]; hm[m+1,m] = zero(T)
            # avnorm = norm(A*vm[m+1])
            avnorm = norm(mul!(p, A, vm[m+1]))
        end

        local F
        # propagate using adaptive step size
        iter = 0
        while iter <= maxiter
            mx = mb + max(1,k1)
            F = exp!(tsgn*tau*view(hm, 1:mx, 1:mx))
            if k1 == 0
                err_loc = btol
                break
            else
                F[m+1,m+1] = h*F[m,m+2]
                F[m+2,m+1] = h*F[m,m+3]
                err1 = abs(beta*F[m+1,m+1])
                err2 = abs(beta*F[m+2,m+1] * avnorm)
                if err1 > 10*err2
                    err_loc = err2
                    r = 1/m
                elseif err1 > err2
                    err_loc = (err1*err2)/(err1-err2)
                    r = 1/m
                else
                    err_loc = err1
                    r = 1/(m-1)
                end
            end

            # time-step sufficient?
            if err_loc <= delta * tau*tol
                break
            else
                tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
                tau = round(tau, sigdigits=2) # round to 2 signiﬁcant digits
                                     # to prevent numerical noise
                if iter == maxiter
                    error("Number of iterations exceeded $(maxiter). Requested tolerance might be too high.")
                end
                iter = iter + 1
            end
        end
        mx = mb + max(0, k1-2)

        for k = 1:mx
            w = axpy!(beta*F[k, mb+1], vm[k], w)
        end

        # beta = norm(A*w + u)
        mul!(avtmp, A, w)
        avtmp .+= u
        beta = norm(avtmp)
        tk = tk + tau

        tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
        tau = round(tau, sigdigits=2) # round to 2 signiﬁcant digits
                             # to prevent numerical noise

        err_loc = max(err_loc, rndoff)

        fill!(hm, zero(T))
    end
    return w
end # phimv!


text:
"""
    BinomialSample(Z, n)

An observed sample ``Z`` drawn from a Binomial distribution with `n` trials.

```math
Z \\sim \\text{Binomial}(n, p)
```

``p`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``p``.

```jldoctest
julia> BinomialSample(2, 10)          # 2 out of 10 trials successful
Z=2  | n=10
```
"""
struct BinomialSample{T,S<:Integer} <: DiscreteEBayesSample{T}
    Z::T
    n::S     #TODO: add checks that Z \in {0,...,n}
end

BinomialSample(n::Integer) = BinomialSample(missing, n)

function Base.show(io::IO, Z::BinomialSample)
    spaces_to_keep = ismissing(response(Z)) ? 1 : max(3 - ndigits(response(Z)), 1)
    spaces = repeat(" ", spaces_to_keep)
    print(io, "Z=", response(Z), spaces, "| ", "n=", ntrials(Z))
end

function Base.show(io::IO, Z::BinomialSample{<:Interval})
    Zz = response(Z)
    print(io, "Z ∈ ")
    show(IOContext(io, :compact => true), Zz)
    print(io, " | ", "n=", ntrials(Z))
end

summarize_by_default(::AbstractVector{<:BinomialSample}) = true

response(Z::BinomialSample) = Z.Z
ntrials(Z::BinomialSample) = Z.n
nuisance_parameter(Z::BinomialSample) = ntrials(Z)

likelihood_distribution(Z::BinomialSample, p) = Binomial(ntrials(Z), p)

function fill_levels(Zs::AbstractVector{<:DiscreteEBayesSample})
    skedasticity(Zs) == Homoskedastic() ||
        error("Heteroskedastic likelihood not implemented.")
    #_min, _max = extrema(response.(Zs))
    #n = ntrials(Zs[1])
    #BinomialSample.(_min:_max, n)
    #BinomialSample.(0:n, n) #TODO!  Is it responsibility of whatever passes stuff in here
    # to e.g. include 0 counts?
    sort(unique(Zs))
end

function fill_levels(Zs::MultinomialSummary{<:DiscreteEBayesSample})
    fill_levels(collect(keys(Zs)))
end

function dictfun(::Nothing, Zs_summary::VectorOrSummary{<:DiscreteEBayesSample}, f)
    skedasticity(Zs_summary) == Homoskedastic() ||
        error("Heteroskedastic likelihood not implemented.")
    Zs = fill_levels(Zs_summary)
    if !isa(f, AbstractVector)
        f = f.(Zs)
    end
    DictFunction(Zs, f)
end

#-----------------------------------------------------------------------
#---------- Beta Binomial Conjugacy-------------------------------------
#-----------------------------------------------------------------------

function default_target_computation(::BasicPosteriorTarget,::BinomialSample, ::Beta)
    Conjugate()
end

function marginalize(Z::BinomialSample, prior::Beta)
    @unpack α, β = prior
    BetaBinomial(ntrials(Z), α, β)
end

function posterior(Z::BinomialSample, prior::Beta)
    Beta(prior.α + response(Z), prior.β + ntrials(Z) - response(Z))
end

# Fit BetaBinomial
function StatsBase.fit(
    ::MethodOfMoments{<:Beta},
    Zs::VectorOrSummary{<:BinomialSample},
    ::Homoskedastic,
)
    # TODO: Let ::Homoskedastic carry type information.
    n = ntrials(Zs[1])
    μ₁ = mean(response.(Zs), weights(Zs))
    μ₂ = mean(abs2.(response.(Zs)), weights(Zs))
    denom = n * (μ₂ / μ₁ - μ₁ - 1) + μ₁
    α = (n * μ₁ - μ₂) / denom
    β = (n - μ₁) * (n - μ₂ / μ₁) / denom
    Beta(α, β)
end

function StatsBase.fit(
    method::ParametricMLE{<:Beta},
    Zs::AbstractVector{<:BinomialSample},
    skedasticity,
)
    StatsBase.fit(method, summarize(Zs), skedasticity)
end

function StatsBase.fit(
    ::ParametricMLE{<:Beta},
    Zs_summary::MultinomialSummary{<:BinomialSample},
    skedasticity,
)
    func = TwiceDifferentiable(
        params -> -loglikelihood(Zs_summary, Beta(params...)),
        [1.0; 1.0];
        autodiff = :forward,
    )
    dfc = TwiceDifferentiableConstraints([0.0; 0.0], [Inf; Inf])

    opt = optimize(func, dfc, [1.0; 1.0], IPNewton())
    Beta(Optim.minimizer(opt)...)
end

# DiscretePriorClass

function _set_defaults(
    convexclass::DiscretePriorClass,
    Zs::VectorOrSummary{<:BinomialSample};
    hints,
)
    eps = get(hints, :eps, 1e-4)
    prior_grid_length = get(hints, :prior_grid_length, 300)::Integer
    DiscretePriorClass(range(eps; stop = 1 - eps, length = prior_grid_length))
end

"""
    ChiSquaredFLocalization(α) <: FLocalization

The ``\\chi^2`` F-localization at confidence level ``1-\\alpha`` for a discrete random variable
taking values in ``\\{0,\\dotsc, N\\}``. It is equal to:
```math
f: \\sum_{x=0}^N \\frac{(n \\hat{f}_n(x) - n f(x))^2}{n f(x)} \\leq \\chi^2_{N,1-\\alpha},
```
where ``\\chi^2_{N,1-\\alpha}`` is the ``1-\\alpha`` quantile of the Chi-squared
distribution with ``N`` degrees of freedom, ``n`` is the sample size,
``\\hat{f}_n(x)`` is the proportion of samples equal to ``x`` and ``f(x)`` is then
population pmf.
"""
Base.@kwdef struct ChiSquaredFLocalization{T} <: FLocalization
    α::T = 0.05
end

vexity(::ChiSquaredFLocalization) = ConvexVexity()

struct FittedChiSquaredFLocalization{T,S,D<:AbstractDict{T,S},C} <: FittedFLocalization
    summary::D
    band::S
    chisq::C
    dof::Int
    n::Int
end

vexity(chisq::FittedChiSquaredFLocalization) = vexity(chisq.chisq)

function nominal_alpha(chisq::FittedChiSquaredFLocalization)
    nominal_alpha(chisq.chisq)
end

# TODO: Allow this to work more broadly.
function StatsBase.fit(chisq::ChiSquaredFLocalization, Zs::AbstractVector{<:BinomialSample})
    StatsBase.fit(chisq, summarize(Zs))
end

function StatsBase.fit(chisq::ChiSquaredFLocalization, Zs_summary::MultinomialSummary)
    n = nobs(Zs_summary)
    _levels = fill_levels(Zs_summary)
    _dof = length(_levels) - 1        #ntrials(_levels[1]) #again maybe Homoskedastic() should return what type of homoskedastic
    empirical_probs = Zs_summary.(_levels) ./ n
    _dict = SortedDict(keys(Zs_summary.store) .=> empirical_probs)
    α = nominal_alpha(chisq)
    band =  quantile(Chisq(_dof), 1-α)
    FittedChiSquaredFLocalization(_dict, band, chisq, _dof, n)
end

function flocalization_constraint!(
    model,
    chisq::FittedChiSquaredFLocalization,
    prior::PriorVariable,
)
    n = chisq.n

    ts = @variable(model, [1:(chisq.dof+1)])
    @constraint(model, ts .>= 0)
    band = chisq.band
    for (i, (Z, pdf_value)) in enumerate(chisq.summary)
        _pdf = pdf(prior, Z::EBayesSample)
        @constraint(model, [ts[i]; n * _pdf; n * _pdf - n * pdf_value] in RotatedSecondOrderCone())
    end
    @constraint(model, sum(ts) <= band/2)
    model
end


text:
# TODO: implement methods for n-dim arrays

"""
    Rotate90 <: Augmentor.AffineOperation

Description
--------------

Rotates the image upwards 90 degrees. This is a special case
rotation because it can be performed very efficiently by simply
rearranging the existing pixels. However, it is generally not the
case that the output image will have the same size as the input
image, which is something to be aware of.

If created using the parameter `p`, the operation will be lifted
into `Either(p=>Rotate90(), 1-p=>NoOp())`, where `p` denotes the
probability of applying `Rotate90` and `1-p` the probability for
applying [`NoOp`](@ref). See the documentation of
[`Either`](@ref) for more information.

Usage
--------------

    Rotate90()

    Rotate90(p)

Arguments
--------------

- **`p::Number`** : Optional. Probability of applying the
    operation. Must be in the interval [0,1].

See also
--------------

[`Rotate180`](@ref), [`Rotate270`](@ref), [`Rotate`](@ref),
[`Either`](@ref), [`augment`](@ref)

Examples
--------------

```jldoctest
julia> using Augmentor

julia> img = [200 150; 50 1]
2×2 Array{Int64,2}:
 200  150
  50    1

julia> img_new = augment(img, Rotate90())
2×2 Array{Int64,2}:
 150   1
 200  50
```
"""
struct Rotate90 <: AffineOperation end
Rotate90(p::Number) = Either(Rotate90(), p)

@inline supports_permute(::Type{Rotate90}) = true

toaffinemap(::Rotate90, img::AbstractMatrix) = recenter(RotMatrix(pi/2), center(img))
applyeager(::Rotate90, img::AbstractMatrix, param) = plain_array(rotl90(img))
applylazy_fallback(op::Rotate90, img::AbstractMatrix, param) = applypermute(op, img, param)

function applypermute(::Rotate90, img::AbstractMatrix{T}, param) where T
    idx = map(StepRange, axes(img))
    perm_img = PermutedDimsArray{T,2,(2,1),(2,1),typeof(img)}(img)
    view(perm_img, reverse(idx[2]), idx[1])
end

function applypermute(::Rotate90, sub::SubArray{T,2,IT,<:NTuple{2,AbstractRange}}, param) where {T,IT<:PermutedDimsArray{T,2,(2,1)}}
    img = parent(parent(sub))
    view(img, reverse(sub.indices[2]), sub.indices[1])
end

function applypermute(::Rotate90, sub::SubArray{T,2,IT,<:NTuple{2,AbstractRange}}, param) where {T,IT}
    img = parent(sub)
    perm_img = PermutedDimsArray{T,2,(2,1),(2,1),typeof(img)}(img)
    if typeof(sub.indices) <: NTuple{2, StepRange{Int64, Int64}}
        view(perm_img, reverse(sub.indices[2]), sub.indices[1])
    else
        ixs = map(StepRange, sub.indices)
        view(perm_img, reverse(ixs[2]), ixs[1])
    end
end

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

"""
    Rotate180 <: Augmentor.AffineOperation

Description
--------------

Rotates the image 180 degrees. This is a special case rotation
because it can be performed very efficiently by simply
rearranging the existing pixels. Furthermore, the output image
will have the same dimensions as the input image.

If created using the parameter `p`, the operation will be lifted
into `Either(p=>Rotate180(), 1-p=>NoOp())`, where `p` denotes the
probability of applying `Rotate180` and `1-p` the probability for
applying [`NoOp`](@ref). See the documentation of
[`Either`](@ref) for more information.

Usage
--------------

    Rotate180()

    Rotate180(p)

Arguments
--------------

- **`p::Number`** : Optional. Probability of applying the
    operation. Must be in the interval [0,1].

See also
--------------

[`Rotate90`](@ref), [`Rotate270`](@ref), [`Rotate`](@ref),
[`Either`](@ref), [`augment`](@ref)

Examples
--------------

```jldoctest
julia> using Augmentor

julia> img = [200 150; 50 1]
2×2 Array{Int64,2}:
 200  150
  50    1

julia> img_new = augment(img, Rotate180())
2×2 Array{Int64,2}:
   1   50
 150  200
```
"""
struct Rotate180 <: AffineOperation end
Rotate180(p::Number) = Either(Rotate180(), p)

@inline supports_stepview(::Type{Rotate180}) = true

toaffinemap(::Rotate180, img::AbstractMatrix) = recenter(RotMatrix(pi), center(img))
applyeager(::Rotate180, img::AbstractMatrix, param) = plain_array(rot180(img))
applylazy_fallback(op::Rotate180, img::AbstractMatrix, param) = applystepview(op, img, param)

function applystepview(::Rotate180, img::AbstractMatrix, param)
    idx = map(i->1:1:length(i), axes(img))
    indirect_view(img, (reverse(idx[1]), reverse(idx[2])))
end

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

"""
    Rotate270 <: Augmentor.AffineOperation

Description
--------------

Rotates the image upwards 270 degrees, which can also be
described as rotating the image downwards 90 degrees. This is a
special case rotation, because it can be performed very
efficiently by simply rearranging the existing pixels. However,
it is generally not the case that the output image will have the
same size as the input image, which is something to be aware of.

If created using the parameter `p`, the operation will be lifted
into `Either(p=>Rotate270(), 1-p=>NoOp())`, where `p` denotes the
probability of applying `Rotate270` and `1-p` the probability for
applying [`NoOp`](@ref). See the documentation of
[`Either`](@ref) for more information.

Usage
--------------

    Rotate270()

    Rotate270(p)

Arguments
--------------

- **`p::Number`** : Optional. Probability of applying the
    operation. Must be in the interval [0,1].

See also
--------------

[`Rotate90`](@ref), [`Rotate180`](@ref), [`Rotate`](@ref),
[`Either`](@ref), [`augment`](@ref)

Examples
--------------

```jldoctest
julia> using Augmentor

julia> img = [200 150; 50 1]
2×2 Array{Int64,2}:
 200  150
  50    1

julia> img_new = augment(img, Rotate270())
2×2 Array{Int64,2}:
 50  200
  1  150
```
"""
struct Rotate270 <: AffineOperation end
Rotate270(p::Number) = Either(Rotate270(), p)

@inline supports_permute(::Type{Rotate270}) = true

toaffinemap(::Rotate270, img::AbstractMatrix) = recenter(RotMatrix(-pi/2), center(img))
applyeager(::Rotate270, img::AbstractMatrix, param) = plain_array(rotr90(img))
applylazy_fallback(op::Rotate270, img::AbstractMatrix, param) = applypermute(op, img, param)

function applypermute(::Rotate270, img::AbstractMatrix{T}, param) where T
    idx = map(StepRange, axes(img))
    perm_img = PermutedDimsArray{T,2,(2,1),(2,1),typeof(img)}(img)
    view(perm_img, idx[2], reverse(idx[1]))
end

function applypermute(::Rotate270, sub::SubArray{T,2,IT,<:NTuple{2,AbstractRange}}, param) where {T,IT<:PermutedDimsArray{T,2,(2,1)}}
    img = parent(parent(sub))
    view(img, sub.indices[2], reverse(sub.indices[1]))
end

function applypermute(::Rotate270, sub::SubArray{T,2,IT,<:NTuple{2,AbstractRange}}, param) where {T,IT}
    img = parent(sub)
    perm_img = PermutedDimsArray{T,2,(2,1),(2,1),typeof(img)}(img)
    if typeof(sub.indices) <: NTuple{2, StepRange{Int64, Int64}}
        view(perm_img, sub.indices[2], reverse(sub.indices[1]))
    else
        ixs = map(StepRange, sub.indices)
        view(perm_img, ixs[2], reverse(ixs[1]))
    end
end

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

for deg in (90, 180, 270)
    T = Symbol(:Rotate, deg)
    @eval function showconstruction(io::IO, op::$T)
        print(io, typeof(op).name.name, "()")
    end
    @eval function Base.show(io::IO, op::$T)
        if get(io, :compact, false)
            print(io, "Rotate ", $deg, " degree")
        else
            print(io, "Augmentor.", $T, "()")
        end
    end
end

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

"""
    Rotate <: Augmentor.AffineOperation

Description
--------------

Rotate the image upwards for the given `degree`. This operation
can only be performed as an affine transformation and will in
general cause other operations of the pipeline to use their
affine formulation as well (if they have one).

In contrast to the special case rotations (e.g.
[`Rotate90`](@ref), the type `Rotate` can describe any arbitrary
number of degrees. It will always perform the rotation around the
center of the image. This can be particularly useful when
combining the operation with [`CropNative`](@ref).

Usage
--------------

    Rotate(degree)

Arguments
--------------

- **`degree`** : `Real` or `AbstractVector` of `Real` that denote
    the rotation angle(s) in degree. If a vector is provided,
    then a random element will be sampled each time the operation
    is applied.

See also
--------------

[`Rotate90`](@ref), [`Rotate180`](@ref), [`Rotate270`](@ref),
[`CropNative`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

# rotate exactly 45 degree
augment(img, Rotate(45))

# rotate between 10 and 20 degree upwards
augment(img, Rotate(10:20))

# rotate one of the five specified degrees
augment(img, Rotate([-10, -5, 0, 5, 10]))
```
"""
struct Rotate{T<:AbstractVector} <: AffineOperation
    degree::T

    function Rotate{T}(degree::T) where {T<:AbstractVector{S} where S<:Real}
        length(degree) > 0 || throw(ArgumentError("The number of different angles passed to \"Rotate(...)\" must be non-zero"))
        new{T}(degree)
    end
end
Rotate(degree::T) where {T<:AbstractVector} = Rotate{T}(degree)
Rotate(degree::Real) = Rotate(degree:degree)

@inline supports_eager(::Type{<:Rotate}) = false

randparam(op::Rotate, img) = Float64(safe_rand(op.degree))

function toaffinemap(op::Rotate, img::AbstractMatrix, angle)
    recenter(RotMatrix(deg2rad(Float64(angle))), center(img))
end

function Base.show(io::IO, op::Rotate)
    if get(io, :compact, false)
        if length(op.degree) == 1
            print(io, "Rotate ", first(op.degree), " degree")
        else
            print(io, "Rotate by θ ∈ ", op.degree, " degree")
        end
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end


text:
println("======================")
println("Tests of NormalDensity")
println("======================")

@test NormalDensity <: SymmetricDensity

###Test construction
n1 = distribution(:normal,ones(2),0.1*eye(2))
d1 = density(:normal,ones(2),0.1*eye(2))
@test isa(d1,NormalDensity)
@test issymmetric(d1)
@test mean(d1.distribution) == mean(n1) && cov(d1.distribution) == cov(n1)

###Test conditioning of the distribution
n1 = recenter(n1,2*ones(2))
condition!(d1,2*ones(2))
@test mean(d1.distribution) == mean(n1) && cov(d1.distribution) == cov(n1)

n1 = rescale(n1,1/2)
scale!(d1,1/2)
@test mean(d1.distribution) == mean(n1) && cov(d1.distribution) == cov(n1)

n1 = update(n1,3*ones(2),3*eye(2))
update!(d1,3*ones(2),3*eye(2))
@test mean(d1.distribution) == mean(n1) && cov(d1.distribution) == cov(n1)

###Test proposing and calculatinng the logprobability
for i in [1,3]
    s1 = samples(:base,2,i,Float64,Float64)
    z1 = zeros(Float64,GeneralizedMetropolisHastings._valuestuple(2,i))
    @test (srand(i+23) ; propose!(d1,s1.values)) == (srand(i+23) ; rand!(n1,z1))
    p1 = Distributions.logpdf(d1.distribution,s1.values) ; isa(p1,Number)?collect(p1):p1
    @test_approx_eq logprobability(d1,s1.values) p1
end

###Test the normal distribution with zero mean
n2 = distribution(:normal,0.1*eye(2))
d2 = density(:normal,0.1*eye(2))
@test isa(d2,NormalDensity)
@test issymmetric(d2)
@test mean(d2.distribution) == mean(n2) && cov(d2.distribution) == cov(n2)
@test_throws MethodError condition!(d2,zeros(2))
n2 = rescale(n2,2.0)
scale!(d2,2.0)
@test mean(d2.distribution) == mean(n2) && cov(d2.distribution) == cov(n2)
n2 = update(n2,eye(2))
update!(d2,eye(2))
@test mean(d2.distribution) == mean(n2) && cov(d2.distribution) == cov(n2)

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

println("=========================")
println("Tests of LogNormalDensity")
println("=========================")

@test LogNormalDensity <: ASymmetricDensity

###Test construction
l1 = distribution(:lognormal,ones(2),0.1*eye(2))
d2 = density(:lognormal,ones(2),0.1*eye(2))
@test isa(d2,LogNormalDensity)
@test !issymmetric(d2)
@test Distributions.location(d2.distribution) == Distributions.location(l1)
@test Distributions.scale(d2.distribution) == Distributions.scale(l1)

###Test conditioning of the distribution
l1 = recenter(l1,2*ones(2))
condition!(d2,2*ones(2))
@test Distributions.location(d2.distribution) == Distributions.location(l1)
@test Distributions.scale(d2.distribution) == Distributions.scale(l1)

l1 = rescale(l1,1/2)
scale!(d2,1/2)
@test Distributions.location(d2.distribution) == Distributions.location(l1)
@test Distributions.scale(d2.distribution) == Distributions.scale(l1)

l1 = update(l1,3*ones(2),3*eye(2))
update!(d2,3*ones(2),3*eye(2))
@test Distributions.location(d2.distribution) == Distributions.location(l1)
@test Distributions.scale(d2.distribution) == Distributions.scale(l1)

###Test proposing and calculatinng the logprobability
for i in [1,3]
    s2 = samples(:base,2,i,Float64,Float64)
    z2 = zeros(Float64,GeneralizedMetropolisHastings._valuestuple(2,i))
    @test (srand(i+23) ; propose!(d2,s2.values)) == (srand(i+23) ; rand!(d2.distribution,z2))
    p2 = Distributions.logpdf(d2.distribution,s2.values) ; isa(p2,Number)?collect(p2):p2
    @test_approx_eq logprobability(d2,s2.values) p2
end

println("====================================")
println("Tests of CompoundDistributionWrapper")
println("====================================")

for p in [(:uniform,()),
          (:laplace,()),
          (:triangular,()),
          (:bactrian,(:normal,0.95)),
          (:bactrian,(:laplace,0.90))]
    println(p)
    dis1 = distributions(p[1],[1.0,2.0],[1.0,0.5],p[2]...)
    den1 = density(p[1],[1.0,2.0],[1.0,0.5],p[2]...)
    @test issymmetric(den1)
    for k=1:length(dis1)
        @test mean(dis1[k]) == mean(den1.distributions[k]) && scale(dis1[k]) == scale(den1.distributions[k])
    end
    condition!(den1,2*ones(2))
    for k=1:length(dis1)
        dis1[k] = recenter(dis1[k],2.0)
        @test mean(dis1[k]) == mean(den1.distributions[k]) && scale(dis1[k]) == scale(den1.distributions[k])
    end
    scale!(den1,2.0)
    for k=1:length(dis1)
        dis1[k] = rescale(dis1[k],2.0)
        @test mean(dis1[k]) == mean(den1.distributions[k]) && scale(dis1[k]) == scale(den1.distributions[k])
    end
    update!(den1,3*ones(2),3.0*ones(2))
    for k=1:length(dis1)
        dis1[k] = update(dis1[k],3.0,3.0)
        @test mean(dis1[k]) == mean(den1.distributions[k]) && scale(dis1[k]) == scale(den1.distributions[k])
    end
    for i in [1,3]
        s3 = samples(:base,2,i,Float64,Float64)
        z3 = zeros(Float64,GeneralizedMetropolisHastings._valuestuple(2,i))
        p3 = zeros(i)
        @test (srand(i+23) ; propose!(den1,s3.values)) == (srand(i+23) ; for k=1:i for j=1:length(dis1) z3[j,k] = rand(dis1[j]) end end ; z3)
        for k=1:i for j=1:length(dis1) p3[k] += Distributions.logpdf(dis1[j],s3.values[j,k]) end end
        @test_approx_eq logprobability(den1,s3.values) p3
    end
end

###test the show functions
println()
println("====================")
println("Test show() function")
println("====================")
show(d1)
show(d2)
println("====================")
println("End  show() function")
println("====================")
println()

nothing


text:
module BeamColumn

using DiffEqOperators: CenteredDifference
using LinearAlgebra
using NLsolve

using ..Mesh

export define, solve, Model, governing_equations

Base.@kwdef mutable struct Model

   A::Union{Array{Float64}, Nothing} = nothing
   Ix::Union{Array{Float64}, Nothing} = nothing
   Iy::Union{Array{Float64}, Nothing} = nothing
   Io::Union{Array{Float64}, Nothing} = nothing
   J::Union{Array{Float64}, Nothing} = nothing
   Cw::Union{Array{Float64}, Nothing} = nothing
   E::Union{Array{Float64}, Nothing} = nothing
   ν::Union{Array{Float64}, Nothing} = nothing
   G::Union{Array{Float64}, Nothing} = nothing
   xc::Union{Array{Float64}, Nothing} = nothing
   yc::Union{Array{Float64}, Nothing} = nothing
   xs::Union{Array{Float64}, Nothing} = nothing
   ys::Union{Array{Float64}, Nothing} = nothing
   xo::Union{Array{Float64}, Nothing} = nothing  #distance from centroid to shear center
   yo::Union{Array{Float64}, Nothing} = nothing  #distance from centroid to shear center
   ax::Union{Array{Float64}, Nothing} = nothing #distance from shear center to qy applied load location
   ay::Union{Array{Float64}, Nothing} = nothing #distance from shear center to qx applied load location
   kx::Union{Array{Float64}, Nothing} = nothing
   ky::Union{Array{Float64}, Nothing} = nothing
   kϕ::Union{Array{Float64}, Nothing} = nothing
   hx::Union{Array{Float64}, Nothing} = nothing #distance from centroid to ky spring location
   hy::Union{Array{Float64}, Nothing} = nothing #distance from centroid to kx spring location

   qx::Union{Array{Float64}, Nothing} = nothing
   qy::Union{Array{Float64}, Nothing} = nothing
   P::Union{Array{Float64}, Nothing} = nothing

   end_boundary_conditions::Union{Array{Int64}, Nothing} = nothing
   supports::Union{Array{Float64}, Nothing} = nothing

   z::Union{Array{Float64}, Nothing} = nothing
   m::Union{Array{Int64}, Nothing} = nothing

   K::Union{Matrix{Float64}, Nothing} = nothing
   F::Union{Array{Float64}, Nothing} = nothing

   free_dof::Union{Array{Int64}, Nothing} = nothing

   u::Union{Array{Float64}, Nothing} = nothing
   v::Union{Array{Float64}, Nothing} = nothing
   ϕ::Union{Array{Float64}, Nothing} = nothing

end

function calculate_derivative_operators(dz)

   #Define the number of nodes.
   num_nodes = length(dz) + 1

   #Add extra dz on each end for padding nodes used in the CenteredDifference function.
   dz = [dz[1]; dz; dz[end]]

   #Calculate the 4th derivative operator.
   nth_derivative = 4
   derivative_order = 2
   Azzzz = CenteredDifference(nth_derivative, derivative_order, dz, num_nodes)

   #Convert the operator to a matrix.
   Azzzz = Array(Azzzz)

   #Trim off the ghost nodes.
   Azzzz = Azzzz[:,2:end-1]

   #Calculate the 2nd derivative operator.
   nth_derivative = 2
   derivative_order = 2
   Azz = CenteredDifference(nth_derivative, derivative_order, dz, num_nodes)
   Azz = Array(Azz)
   Azz = Azz[:,2:end-1]

   return Azzzz, Azz

end

function calculate_boundary_stencils(bc_flag, h, nth_derivative)

   #Calculate boundary conditions stencils without ghost nodes using
   #Jorge M. Souza, "Boundary Conditions in the Finite Difference Method"
   #https://cimec.org.ar/ojs/index.php/mc/article/download/2662/2607

   taylor_coeffs =  [1; h; h^2/2; h^3/6; h^4/24]

   RHS = zeros(Float64, 5)
   #row1*[A;B;C;D;E]*u2
   #row2*[A;B;C;D;E]*u2'
   #row3*[A;B;C;D;E]*u2''
   #row4*[A;B;C;D;E]*u2'''
   #row5*[A;B;C;D;E]*u2''''

   if bc_flag == 1 #simply supported end

      LHS = [1 1 1 1 0
         -1 0 1 2 0
         1 0 1 4 2
         -1 0 1 8 -6
         1 0 1 16 12]

      RHS[nth_derivative + 1] = (1/taylor_coeffs[nth_derivative + 1])
      boundary_stencil = LHS \ RHS
      boundary_stencil = ((boundary_stencil[1:4]),(zeros(4)))  #since u''=0

   elseif bc_flag == 2 #fixed end

      LHS = [1 1 1 1 0
         -1 0 1 2 1
         1 0 1 4 -2
         -1 0 1 8 3
         1 0 1 16 -4]

      RHS[nth_derivative+1] = (1/taylor_coeffs[nth_derivative + 1])
      boundary_stencil = LHS\RHS
      boundary_stencil = ((boundary_stencil[1:4]),(zeros(4)))  #since u'=0

   elseif bc_flag == 3 #free end
                #u'' u'''
      LHS = [1 1 1  0   0
           0 1 2  0   0
           0 1 4  0.5 0
           0 1 8  0   6
           0 1 16 0  0]

      RHS[nth_derivative+1] = (1/taylor_coeffs[nth_derivative+1])
      boundary_stencil1 = LHS\RHS  #at free end
      boundary_stencil1 = boundary_stencil1[1:3]   #u''=u'''=0

      # use simply supported BC to find stencil at one node in from free end
      LHS = [1 1 1 1 0
         -1 0 1 2 0
         1 0 1 4 2
         -1 0 1 8 -6
         1 0 1 16 12]

      RHS[nth_derivative+1] = (1/taylor_coeffs[nth_derivative+1])
      boundary_stencil2 = LHS\RHS #at one node over from free end
      boundary_stencil2 = boundary_stencil2[1:4]
      boundary_stencil = ((boundary_stencil1), (boundary_stencil2))  #two stencils are calculated

   end

   return boundary_stencil

end

function apply_end_boundary_conditions(A, end_boundary_conditions, nth_derivative, dz)

   #Consider the left end boundary condition.  Swap out interior stencil for boundary stencil.
   h = dz[1]
   bc_flag = end_boundary_conditions[1]
   boundary_stencil = calculate_boundary_stencils(bc_flag, h, nth_derivative)

   A[1,:] .= 0.0
   A[2,:] .= 0.0

   if (bc_flag == 1) | (bc_flag == 2)   #make this cleaner, combine
      A[2,1:length(boundary_stencil[1])] = boundary_stencil[1]
   else
      A[1,1:length(boundary_stencil[1])] = boundary_stencil[1]
      A[2,1:length(boundary_stencil[2])] = boundary_stencil[2]
   end

   #Consider the right end boundary condition.
   h = dz[end]
   bc_flag = end_boundary_conditions[2]
   boundary_stencil = calculate_boundary_stencils(bc_flag,h,nth_derivative)

   A[end,:] .= 0.0
   A[end-1,:] .= 0.0

   if (bc_flag == 1) | (bc_flag == 2)
      A[end-1,(end-(length(boundary_stencil[1])-1)):end] = reverse(boundary_stencil[1])
   else
      A[end,end-(length(boundary_stencil[1])-1):end] = reverse(boundary_stencil[1])
      A[end-1,end-(length(boundary_stencil[2])-1):end] = reverse(boundary_stencil[2])
   end

   return A

end

function discretize(member_definitions)

   #Define discretization along beam (z) and assign a member type to each node (m).
   dz, z, m = Mesh.define_line_element(member_definitions)

   return z, m

end

function define(z, m, member_definitions, section_properties, material_properties, kx, ky, kϕ, hx, hy, qx, qy, P, ax, ay, end_boundary_conditions, supports)

   #This is the user interface mapping.
   #L(1) dL(2) section_properties(3) material_properties(4)

   dz = diff(z)

   #define property vectors
   A =  Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 1)
   Ix = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 2)
   Iy = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 3)
   J = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 4)
   Cw = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 5)
   xc = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 6)
   yc = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 7)
   xs = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 8)
   ys = Mesh.create_line_element_property_array(member_definitions, m, dz, section_properties, 3, 9)

   xo = -(xc .- xs)
   yo = yc .- ys

   Io = Ix .+ Iy .+ A .* (xo.^2 + yo.^2)

   E = Mesh.create_line_element_property_array(member_definitions, m, dz, material_properties, 4, 1)
   ν = Mesh.create_line_element_property_array(member_definitions, m, dz, material_properties, 4, 2)
   G = E./(2 .*(1 .+ ν))

   model = Model(A, Ix, Iy, Io, J, Cw, E, ν, G, xc, yc, xs, ys, xo, yo, ax, ay, kx, ky, kϕ, hx, hy, qx, qy, P, end_boundary_conditions, supports, z, m, nothing, nothing, nothing, nothing, nothing, nothing)

   return model

end

function governing_equations(model)

   #Define the number of nodes.
   num_nodes = length(model.z)

   #Calculate the derivative operators.
   dz = diff(model.z)
   Azzzz,Azz = calculate_derivative_operators(dz)

   #Apply left and right end boundary condition stencils to derivative operators.
   nth_derivative = 4
   Azzzz = apply_end_boundary_conditions(Azzzz, model.end_boundary_conditions, nth_derivative, dz)

   nth_derivative = 2
   Azz = apply_end_boundary_conditions(Azz, model.end_boundary_conditions, nth_derivative, dz)

   #build identity matrix for ODE operations
   AI = Matrix(1.0I,num_nodes,num_nodes)

   #build operator matrix that doesn't update with load

   #start building operator matrix, LHS first
   A11 = zeros(num_nodes,num_nodes)
   A12 = zeros(num_nodes,num_nodes)
   A13 = zeros(num_nodes,num_nodes)
   A21 = zeros(num_nodes,num_nodes)
   A22 = zeros(num_nodes,num_nodes)
   A23 = zeros(num_nodes,num_nodes)
   A31 = zeros(num_nodes,num_nodes)
   A32 = zeros(num_nodes,num_nodes)
   A33 = zeros(num_nodes,num_nodes)

   #calculate operator quantities on LHS  AU=B
   for i = 1:num_nodes
      A11[i,:] = model.E .* model.Iy .* Azzzz[i,:] .+ model.P .* Azz[i,:] .+ model.kx .* AI[i,:]
      A13[i,:] = -model.kx .* model.hy .* AI[i,:]
      A22[i,:] = model.E .* model.Ix .* Azzzz[i,:] .+ model.P .* Azz[i,:] .+ model.ky .* AI[i,:]
      A23[i,:] = model.ky .* model.hx .* AI[i,:]
      A31[i,:] = -model.kx .* model.hy .* AI[i,:]
      A32[i,:] = model.ky .* model.hx .* AI[i,:]
      A33[i,:] = model.E .* model.Cw .* Azzzz[i,:] .-(model.G .* model.J .- (model.P .* model.Io ./ model.A)).*Azz[i,:] .+model.kx .* model.hy.^2 .*AI[i,:] .+ model.ky .* model.hx.^2 .*AI[i,:] .+ model.kϕ .* AI[i,:] .+  model.qx .* model.ax .* AI[i,:] - .+ model.qy .* model.ay .* AI[i,:]
   end

   #calculate RHS of AU=B
   B1 = model.qx
   B2 = model.qy
   B3 = model.qx .* model.ay .+ model.qy .* model.ax

   #reduce problem to free dof
   #Find fixed dof.
   fixed_dof = []
   for i = 1:length(model.supports)

      if i == 1

            fixed_dof = findall(x->abs(x-model.supports[i])<=10e-6, model.z)

      elseif i != 1

            new_fixed_dof = findall(x->abs(x-model.supports[i])<=10e-6, model.z)
            fixed_dof = [fixed_dof; new_fixed_dof]

      end

   end

   model.free_dof = setdiff(1:num_nodes,fixed_dof)

   #Define stiffness matrix.
   model.K = [A11[model.free_dof,model.free_dof] A12[model.free_dof,model.free_dof] A13[model.free_dof,model.free_dof];
      A21[model.free_dof,model.free_dof] A22[model.free_dof,model.free_dof] A23[model.free_dof,model.free_dof];
      A31[model.free_dof,model.free_dof] A32[model.free_dof,model.free_dof] A33[model.free_dof,model.free_dof]]

   #Define external force vector.
   model.F = [B1[model.free_dof]; B2[model.free_dof]; B3[model.free_dof]]

   return model

end

function residual!(R, U, K, F)

   for i=1:length(F)
    R[i] = transpose(K[i,:]) * (U) - F[i]
   end

   return R

end

function solve(model)

   #Set up solution matrices from governing equations.
   model = governing_equations(model)

   #Define the number of nodes.
   num_nodes=length(model.z)

   #Define the deformation vectors.
   u=zeros(Float64, num_nodes)
   v=zeros(Float64, num_nodes)
   ϕ=zeros(Float64, num_nodes)

   #Define the deformation initial guess for the nonlinear solver.
   deformation_guess = model.K \ model.F    #consider revising this for large systems, it might be slow...

   #Solve for the deformations.
   solution = nlsolve((R, U) ->residual!(R, U, model.K, model.F), deformation_guess)

   #Pull the displacements and twist from the solution results.
   u[model.free_dof] = solution.zero[1:length(model.free_dof)]
   v[model.free_dof] = solution.zero[length(model.free_dof)+1:2*length(model.free_dof)]
   ϕ[model.free_dof] = solution.zero[2*length(model.free_dof)+1:3*length(model.free_dof)]

   #Add deformations to data structure.
   model.u = u
   model.v = v
   model.ϕ = ϕ

   return model

end

end #module


text:
module EMLelectrical
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("electrical/wire.jl")
	include("electrical/Resistor.jl")
	include("electrical/Capacitor.jl")
	include("electrical/Indutor.jl")
	include("electrical/Supply.jl")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CSINodeDriver holds information about the specification of one CSI driver installed on a node

    IoK8sApiStorageV1CSINodeDriver(;
        allocatable=nothing,
        name=nothing,
        nodeID=nothing,
        topologyKeys=nothing,
    )

    - allocatable::IoK8sApiStorageV1VolumeNodeResources : allocatable represents the volume resources of a node that are available for scheduling. This field is beta.
    - name::String : This is the name of the CSI driver that this object refers to. This MUST be the same name returned by the CSI GetPluginName() call for that driver.
    - nodeID::String : nodeID of the node from the driver point of view. This field enables Kubernetes to communicate with storage systems that do not share the same nomenclature for nodes. For example, Kubernetes may refer to a given node as \"node1\", but the storage system may refer to the same node as \"nodeA\". When Kubernetes issues a command to the storage system to attach a volume to a specific node, it can use this field to refer to the node name using the ID that the storage system will understand, e.g. \"nodeA\" instead of \"node1\". This field is required.
    - topologyKeys::Vector{String} : topologyKeys is the list of keys supported by the driver. When a driver is initialized on a cluster, it provides a set of topology keys that it understands (e.g. \"company.com/zone\", \"company.com/region\"). When a driver is initialized on a node, it provides the same topology keys along with values. Kubelet will expose these topology keys as labels on its own node object. When Kubernetes does topology aware provisioning, it can use this list to determine which labels it should retrieve from the node object and pass back to the driver. It is possible for different nodes to use different topology keys. This can be empty if driver does not support topology.
"""
mutable struct IoK8sApiStorageV1CSINodeDriver <: SwaggerModel
    allocatable::Any # spec type: Union{ Nothing, IoK8sApiStorageV1VolumeNodeResources } # spec name: allocatable
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    nodeID::Any # spec type: Union{ Nothing, String } # spec name: nodeID
    topologyKeys::Any # spec type: Union{ Nothing, Vector{String} } # spec name: topologyKeys

    function IoK8sApiStorageV1CSINodeDriver(;allocatable=nothing, name=nothing, nodeID=nothing, topologyKeys=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("allocatable"), allocatable)
        setfield!(o, Symbol("allocatable"), allocatable)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("nodeID"), nodeID)
        setfield!(o, Symbol("nodeID"), nodeID)
        validate_property(IoK8sApiStorageV1CSINodeDriver, Symbol("topologyKeys"), topologyKeys)
        setfield!(o, Symbol("topologyKeys"), topologyKeys)
        o
    end
end # type IoK8sApiStorageV1CSINodeDriver

const _property_map_IoK8sApiStorageV1CSINodeDriver = Dict{Symbol,Symbol}(Symbol("allocatable")=>Symbol("allocatable"), Symbol("name")=>Symbol("name"), Symbol("nodeID")=>Symbol("nodeID"), Symbol("topologyKeys")=>Symbol("topologyKeys"))
const _property_types_IoK8sApiStorageV1CSINodeDriver = Dict{Symbol,String}(Symbol("allocatable")=>"IoK8sApiStorageV1VolumeNodeResources", Symbol("name")=>"String", Symbol("nodeID")=>"String", Symbol("topologyKeys")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiStorageV1CSINodeDriver }) = collect(keys(_property_map_IoK8sApiStorageV1CSINodeDriver))
Swagger.property_type(::Type{ IoK8sApiStorageV1CSINodeDriver }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1CSINodeDriver[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1CSINodeDriver }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1CSINodeDriver[property_name]

function check_required(o::IoK8sApiStorageV1CSINodeDriver)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("nodeID")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1CSINodeDriver }, name::Symbol, val)
end


text:
module TestBaseUtils

using LightGBM
using Test

@testset "deepcopy: Booster" begin

    b = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))
    bb = deepcopy(b)

    @test b != bb
    @test b !== bb
    @test LightGBM.LGBM_BoosterSaveModelToString(b) == LightGBM.LGBM_BoosterSaveModelToString(bb)

    # test for empty boosters too
    b_empty = LightGBM.Booster()
    bb_empty = deepcopy(b_empty)

    @test b_empty != bb_empty
    @test b_empty !== bb_empty

end

end # module


text:
function gadflyhistogram(bigvec, smallvec)
    using Gadfly
    plot(
        layer(
            x = bigvec,
            Geom.histogram(bincount=25),
            order=1
        ),
        layer(
            x = smallvec,
            Geom.histogram(bincount=25),
            order=2,
            Theme(default_color=color("red"))
        ),
        Guide.xlabel("Sample mean"),
        Guide.ylabel("Frequency"),
        Guide.title("Sample Mean Histogram in Dataset 1")
    );
end


text:
"""
    struct gRPCStatus
        success::Bool
        message::String
    end

`gRPCStatus` represents the status of a request. It has the following fields:

- `success`: whether the request was completed successfully.
- `message`: any error message if request was not successful
"""
struct gRPCStatus
    success::Bool
    grpc_status::Int
    message::String
    exception::Union{Nothing,Exception}
end

gRPCStatus(success::Bool, grpc_status::Int, message::AbstractString) = gRPCStatus(success, grpc_status, string(message), nothing)
function gRPCStatus(status_future)
    try
        fetch(status_future)
    catch ex
        task_exception = isa(ex, TaskFailedException) ? ex.task.exception : ex
        while isa(task_exception, TaskFailedException)
            task_exception = task_exception.task.exception
        end
        gRPCStatus(false, StatusCode.INTERNAL.code, string(task_exception), task_exception)
    end
end

"""
    struct gRPCServiceCallException
        message::String
    end

A `gRPCServiceCallException` is thrown if a gRPC request is not successful.
It has the following members:

- `message`: any error message if request was not successful
"""
struct gRPCServiceCallException <: gRPCException
    grpc_status::Int
    message::String
end

Base.show(io::IO, m::gRPCServiceCallException) = print(io, "gRPCServiceCallException: $(m.grpc_status), $(m.message)")

"""
    gRPCCheck(status; throw_error::Bool=true)

Every gRPC request returns the result and a future representing the status
of the gRPC request. Check the response of a gRPC request and raise a
`gRPCException` if it has failed. If `throw_error` is set to false, this
returns `true` or `false` indicating success instead.
"""
gRPCCheck(status_future; throw_error::Bool=true) = gRPCCheck(gRPCStatus(status_future); throw_error=throw_error)
function gRPCCheck(status::gRPCStatus; throw_error::Bool=true)
    if throw_error && !status.success
        if status.exception === nothing
            throw(gRPCServiceCallException(status.grpc_status, status.message))
        else
            throw(status.exception)
        end
    end
    status.success
end

"""
    gRPCController(;
        [ maxage::Int = 0, ]
        [ keepalive::Int64 = 60, ]
        [ negotiation::Symbol = :http2_prior_knowledge, ]
        [ revocation::Bool = true, ]
        [ request_timeout::Real = Inf, ]
        [ connect_timeout::Real = 0, ]
        [ max_message_length = DEFAULT_MAX_MESSAGE_LENGTH, ]
        [ max_recv_message_length = 0, ]
        [ max_send_message_length = 0, ]
        [ verbose::Bool = false, ]
    )

Contains settings to control the behavior of gRPC requests.
- `maxage`: maximum age (seconds) of a connection beyond which it will not
   be reused (default 180 seconds, same as setting this to 0).
- `keepalive`: interval (seconds) in which to send TCP keepalive messages on
   the connection (default 60 seconds).
- `negotiation`: how to negotiate HTTP2, can be one of `:http2_prior_knowledge`
   (no negotiation, the default), `:http2_tls` (http2 upgrade but only over
   tls), or `:http2` (http2 upgrade)
- `revocation`: whether to check for certificate recovation (default is true)
- `request_timeout`: request timeout (seconds)
- `connect_timeout`: connect timeout (seconds) (default is 300 seconds, same
   as setting this to 0)
- `max_message_length`: maximum message length (default is 4MB)
- `max_recv_message_length`: maximum message length to receive (default is
   `max_message_length`, same as setting this to 0)
- `max_send_message_length`: maximum message length to send (default is
   `max_message_length`, same as setting this to 0)
- `verbose`: whether to print out verbose communication logs (default false)
"""
struct gRPCController <: ProtoRpcController
    maxage::Clong
    keepalive::Clong
    negotiation::Symbol
    revocation::Bool
    request_timeout::Real
    connect_timeout::Real
    max_recv_message_length::Int
    max_send_message_length::Int
    verbose::Bool

    function gRPCController(;
            maxage::Integer = 0,
            keepalive::Integer = 60,
            negotiation::Symbol = :http2_prior_knowledge,
            revocation::Bool = true,
            request_timeout::Real = Inf,
            connect_timeout::Real = 0,
            max_message_length::Integer = DEFAULT_MAX_MESSAGE_LENGTH,
            max_recv_message_length::Integer = 0,
            max_send_message_length::Integer = 0,
            verbose::Bool = false
        )
        if maxage < 0 || keepalive < 0 || request_timeout < 0 || connect_timeout < 0 ||
            max_message_length < 0 || max_recv_message_length < 0 || max_send_message_length < 0
            throw(ArgumentError("Invalid gRPCController parameter"))
        end
        (max_recv_message_length == 0) && (max_recv_message_length = max_message_length)
        (max_send_message_length == 0) && (max_send_message_length = max_message_length)
        new(maxage, keepalive, negotiation, revocation, request_timeout, connect_timeout, max_recv_message_length, max_send_message_length, verbose)
    end
end

"""
    gRPCChannel(baseurl)

`gRPCChannel` represents a connection to a specific service endpoint
(service `baseurl`) of a gRPC server.

A channel also usually has a single network connection backing it and
multiple streams of requests can flow through it at any time. The number
of streams that can be multiplexed is negotiated between the client and
the server.
"""
struct gRPCChannel <: ProtoRpcChannel
    downloader::Downloader
    baseurl::String

    function gRPCChannel(baseurl::String)
        downloader = Downloader(; grace=Inf)
        Curl.init!(downloader.multi)
        Curl.setopt(downloader.multi, CURLMOPT_PIPELINING, CURLPIPE_MULTIPLEX)
        endswith(baseurl, '/') && (baseurl = baseurl[1:end-1])
        new(downloader, baseurl)
    end
end

function to_delimited_message_bytes(msg, max_message_length::Int)
    iob = IOBuffer()
    limitiob = LimitIO(iob, max_message_length)
    write(limitiob, UInt8(0))                   # compression
    write(limitiob, hton(UInt32(0)))            # message length (placeholder)
    data_len = writeproto(limitiob, msg)        # message bytes

    seek(iob, 1)                                # seek out the message length placeholder
    write(iob, hton(UInt32(data_len)))          # fill the message length
    take!(iob)
end

function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, request::T) where T <: ProtoType
    inputchannel = Channel{T}(1)
    put!(inputchannel, request)
    close(inputchannel)
    call_method(channel, service, method, controller, inputchannel)
end
call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, input::Channel{T}) where T <: ProtoType = call_method(channel, service, method, controller, input, get_response_type(method))
function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, input::Channel{T1}, ::Type{Channel{T2}}) where {T1 <: ProtoType, T2 <: ProtoType}
    call_method(channel, service, method, controller, input, Channel{T2}())
end
function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, input::Channel{T1}, ::Type{T2}) where {T1 <: ProtoType, T2 <: ProtoType}
    outchannel, status_future = call_method(channel, service, method, controller, input, Channel{T2}())
    try
        take!(outchannel), status_future
    catch ex
        gRPCCheck(status_future)    # check for core issue
        if isa(ex, InvalidStateException)
            throw(gRPCServiceCallException("Server closed connection without any response"))
        else
            rethrow()               # throw this error if there's no other issue
        end
    end
end
function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, input::Channel{T1}, outchannel::Channel{T2}) where {T1 <: ProtoType, T2 <: ProtoType}
    url = string(channel.baseurl, "/", service.name, "/", method.name)
    status_future = @async grpc_request(channel.downloader, url, input, outchannel;
        maxage = controller.maxage,
        keepalive = controller.keepalive,
        negotiation = controller.negotiation,
        revocation = controller.revocation,
        request_timeout = controller.request_timeout,
        connect_timeout = controller.connect_timeout,
        max_recv_message_length = controller.max_recv_message_length,
        max_send_message_length = controller.max_send_message_length,
        verbose = controller.verbose,
    )
    outchannel, status_future
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -25.420040452488465 <= q <= 21.65452636072683)

                     add_NL_constraint(m, :(sigmoid(-0.9974376164994152 + 0.7314791518637902*sigmoid(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.4470938857218645*sigmoid(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.10010754470216687*sigmoid(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.8594790184444285*sigmoid(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.4658863904521331*sigmoid(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + sigmoid(0.15396139316976853 + -0.4228510094486002*sigmoid(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.8035341866400141*sigmoid(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.07617073018159282*sigmoid(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.3940887684836425*sigmoid(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.6667343490308526*sigmoid(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + sigmoid(0.9665409577623398 + -0.449175249339572*sigmoid(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.2640107596849468*sigmoid(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.7284419447850534*sigmoid(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.5662248826250997*sigmoid(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.21454884110410477*sigmoid(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + sigmoid(-0.1582635851091898 + 0.5103556986884796*sigmoid(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.6017569427714062*sigmoid(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.23749573380464106*sigmoid(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.2344836196764124*sigmoid(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.07066758651932892*sigmoid(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + sigmoid(-0.47638200661148833 + 0.9126631115783219*sigmoid(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.35414400688245484*sigmoid(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.36037763135774403*sigmoid(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.18657915153735738*sigmoid(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.7872605277677782*sigmoid(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Revise
using BondGraphs
using LightGraphs
using ModelingToolkit
using Catalyst
using Test

@testset "BondGraphs.jl" begin
    @testset "Graph functions" begin include("graphfunctions_tests.jl") end
    @testset "Construction" begin include("construction_tests.jl") end
    @testset "Equations" begin include("equation_tests.jl") end
    @testset "Simulations" begin include("simulation_tests.jl") end
    @testset "Conversions" begin include("conversion_tests.jl") end
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
using Turing
using Distributions

const n_samples = 10^5

@sf test_norm() = begin
    a = normal(0, 2; sampler=langevinposterior)
    normal(a, 1; condition=3)
    a
end

as = float(sample(test_norm, n_samples))
ms = mean(as)
vs = var(as)

#Analytic Moments
va =  1/(1/4 + 1/1)
ma = va*(3/1 + 0/4)

@assert 0.95*ma < ms < 1.05*ma
@assert 0.95*va < vs < 1.05*va

#@sf test_gamma_poisson() = begin
#    rate = gamma(2, 2)
#    for i = 1:10
#        poisson(rate; condition=10)
#    end
#    rate
#end

@sf test_beta_bernoulli() = begin
    p = beta(0.5, 1.5)#, sampler=langevinposterior)
    for i = 1:10
        bernoulli(p; condition=1)
    end
    p
end

as = float(sample(test_beta_bernoulli, n_samples))
dist = Beta(9.5, 1.5)
ma = mean(dist)
va = var(dist)
ms = mean(as)
vs = var(as)

@assert 0.95*ma < ms < 1.05*ma
@assert 0.8*va < vs < 1.2*va


text:
using JLD2
using GLMakie

include("one_degree_artifacts.jl")
include("one_degree_interface_heights.jl")
z_interfaces = one_degree_interface_heights()
z_centers = (z_interfaces[1:end-1] .+ z_interfaces[2:end]) ./ 2

filename = "bathymetry-360x150-latitude-75.0.jld2"
file = jldopen(filename)
h = file["bathymetry"]
close(file)
land = h .> 0

initial_condition_path = download_initial_condition()
file = jldopen(initial_condition_path)
Tᵢ = file["T"]
Nz = size(Tᵢ, 3)

for k in 1:Nz
    z = z_centers[k]
    earth = h .> z
    view(Tᵢ, :, :, k)[earth] .= NaN
end

fig = Figure(resolution=(1800, 800))
ax = Axis(fig[2, 1])
slider = Slider(fig[3, 1:2], range=1:Nz, startvalue=Nz)
k = slider.value

depth = @lift z_centers[$k]
title = @lift string("Temperature (ᵒC) at z = ", $depth, " m")
Label(fig[1, 1:2], title)

Tᵏ = @lift view(Tᵢ, :, :, $k)

hm = heatmap!(ax, Tᵏ, colorrange=(-2, 32), colormap=:thermal)
Colorbar(fig[2, 2], hm, label="Temperature (ᵒC)")

display(fig)

record(fig, "one_degree_initial_condition.mp4", 1:Nz, framerate=12) do kk
    k[] = kk
end


text:
"""
    $(SIGNATURES)

True if the array is spin polarized, meaning it has an axis named `:spin` of dimension 2.
"""
function is_spin_polarized(array::AxisArray)
    if !is_spin_polarized(typeof(array))
        return false
    end
    @assert axisvalues(array)[findfirst(axisnames(array), :spin)] == [:↑, :↓]
    true
end
@inline is_spin_polarized{T <: AxisArray}(::Type{T}) = :spin ∈ axisnames(T)

@lintpragma("Ignore use of undeclared variable X")
@traitdef HasSpinDim{X}
@traitimpl HasSpinDim{X} <- is_spin_polarized(X)


text:
type LikelihoodModel{X, Y}
  xs :: Vector{X}
  ys :: Vector{Y}
  zs :: Vector{Y}
  datas :: Vector{Y}
  Ldx :: Matrix
  Lzx :: Matrix
end

# constructors

LikelihoodModel(xs, ys, zs, datas, measerr) =
  LikelihoodModel(xs, ys, zs, datas, likelihoodmat(datas, ys, measerr), likelihoodmat(zs, ys, measerr))

function LikelihoodModel(; xs=[], phi=Void, ys=[], datas=[], measerr=Void, prior=Void, ndata=0, zs=[], zmult=0, sfact=0, smoothkernel=Void)

  if ys == [] && phi != Void
    ys = phi.(xs)
  end

  if datas == [] && ndata > 0
    datas = phi.(rand(prior, ndata)) + rand(measerr, ndata)
  end

  if zs == [] && zmult > 0
    zs = repmat(ys, zmult) + rand(measerr, length(ys)*zmult)
  end

  if sfact > 0
    datas = smoothdata(smoothkernel, sfact, datas)
    measerr = convolute(measerr, smoothkernel)
  end

  m = LikelihoodModel(xs, ys, zs, datas, measerr)
end

# likelihood computation

likelihoodmat(xs, ys, d::Distribution) = map(a->pdf(d, a), (x-y for x in xs, y in ys))

# model wrappers

npmle(m)     = mple(m, 0)
refprior(m)  = mple(m, 1)

logl(m::LikelihoodModel, w) = logl(w, m.Ldx)
dlogl(m::LikelihoodModel, w) = dlogl(w, m.Ldx)

hz(m::LikelihoodModel, w)  =  hz(w, m.Lzx)
dhz(m::LikelihoodModel, w) = dhz(w, m.Lzx)

function mple_obj(m::LikelihoodModel, reg)
  if reg == 0
    w -> logl(m, w)
  elseif reg == 1
    w -> hz(m, w)
  else
    w -> reg*hz(m,w) + (1-reg) * logl(m,w)
  end
end

function dmple_obj(m::LikelihoodModel, reg)
  if reg == 0
    w -> dlogl(m, w)
  elseif reg == 1
    w -> dhz(m,w)
  else
    w -> reg*dhz(m, w) + (1-reg) * dlogl(m, w)
  end
end


text:
module BlockwiseInversion

  using Revise
  using Reexport

  @reexport using BandedMatrices
  @reexport using Distributed
  @reexport using LinearAlgebra
  @reexport using SharedArrays
  @reexport using SparseArrays

  import Base: length, \, size
  import Base: getindex, setindex!

  import LinearAlgebra: rank

  export FourBlockMatrix

  include("four_block_matrix.jl")
  include("array_operations.jl")

  include("getindex.jl")
  include("setindex.jl")

  include("ldiv.jl")
  include("get_sub_matrices.jl")

end


text:
"""
    $(SIGNATURES)
    quality function constructor for least square optimization
"""
function make_qfunc(model, x, y)
    function qfunc(p)
        mask = (y .!= 0)
        model_res = model(x, p)
        sum(((model_res[mask] .- y[mask]).^2) ./ y[mask])
    end
    qfunc
end

"""
    $(SIGNATURES)
    performs pre fit
    # Arguments
    - `chargedist => ChargeDist`: charge distribution on which pre fit is
                                  performed
"""
function pre_fit(chargedist::ChargeDist)
    x = chargedist.x
    y = chargedist.y
    qfunc_ped = make_qfunc(gauss, x, y)
    mxval_ped, mxidx_ped = findmax(y)
    p0_ped = [x[mxidx_ped], 0.01, mxval_ped]
    fit_ped = optimize(qfunc_ped, p0_ped, NewtonTrustRegion())
    popt_ped = Optim.minimizer(fit_ped)
    mask = (x .> (popt_ped[1] + 3 * popt_ped[2]))
    x_spe = x[mask]
    y_spe = y[mask]
    mxval_spe, mxidx_spe = findmax(y_spe)
    mask = (x_spe .< x_spe[mxidx_spe] * 2)
    x_spe = x_spe[mask]
    y_spe = y_spe[mask]
    qfunc_spe = make_qfunc(gauss, x_spe, y_spe)
    p0_spe = [x_spe[mxidx_spe], 0.1, mxval_spe]
    fit_spe = optimize(qfunc_spe, p0_spe, NewtonTrustRegion())
    popt_spe = Optim.minimizer(fit_spe)
    PreFit(popt_ped..., popt_spe...)
end

"""
    $(SIGNATURES)
    performs pmt response fit
    # Arguments
    - `charges => ChargeDist`: charge distribution on which fit is performed
    - `prefit_results => PreFitResults`: results from pre fit used as starting
                                         values for fit
"""
function pmtresp_fit(chargedist::ChargeDist, prefit::PreFit; mod=:default)
    p0 = [prefit.μₚ,
          prefit.σₚ,
          prefit.μₛ,
          prefit.σₛ,
          prefit.Aₛ / prefit.Aₚ,
          prefit.Aₚ + prefit.Aₛ]
    if mod == :uap
        push!(p0, prefit.μₛ / 5)
        push!(p0, prefit.σₛ / 5)
        push!(p0, prefit.Aₛ / 5)
        func = pmtresp_uap
        ResultStruct = PMTRespUapFit
    else
        func = pmtresp
        ResultStruct = PMTRespFit
    end
    qfunc = make_qfunc(func, chargedist.x, chargedist.y)
    fit = optimize(qfunc, p0, Newton())
    popt = Optim.minimizer(fit)

    ResultStruct(popt...)
end


text: isisogram(str::AbstractString) = !occursin(r"(?i)(\w).*\1", str)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
function defaultcb(model::IsingModel,args::Any)
  return nothing
end