text:
"""
    @unroll1 for-loop

Unroll the first iteration of a `for`-loop.
Set `$first` to true in the first iteration.

Example:
```
    @unroll1 for i in 1:10
        if $first
            a, state = iterate('A':'Z')
        else
            a, state = iterate('A':'Z', state)
        end
        println(i => a)
    end
```
"""
macro unroll1(expr)
    @assert expr isa Expr
    @assert expr.head == :for
    iterspec = expr.args[1]

    @assert iterspec isa Expr
    @assert  iterspec.head == :(=)

    i_ = iterspec.args[1]
    i = esc(i_)
    iter = esc(iterspec.args[2])
    body = esc(expr.args[2])

    body_1 = eval(Expr(:let, :(first = true), Expr(:quote, body)))
    body_i = eval(Expr(:let, :(first = false), Expr(:quote, body)))

    if i_ isa Expr && i_.head == :tuple
        decl = esc(quote local $(i_.args...) end)
    else
        decl = quote local $i end
    end

    quote
        local st
        $decl
        @goto enter
        while true
            @goto exit
            @label enter
                ϕ = iterate($iter)
                ϕ === nothing && break
                $i, st = ϕ
                $(body_1)
            @label exit
            while true
                ϕ = iterate($iter, st)
                ϕ === nothing && break
                $i, st = ϕ
                $(body_i)
            end
            break
        end
    end
end


text:
using Multilane
using POMDPs
using POMDPToolbox
using BenchmarkTools

nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=200.) #2.=>col_length=8
_discount = 1.
nb_cars=10

rmodel = NoCrashRewardModel()

dmodel = NoCrashIDMMOBILModel(nb_cars, pp)

mdp = NoCrashMDP{typeof(rmodel), typeof(dmodel.behaviors)}(dmodel, rmodel, _discount, true);

rng = MersenneTwister(5)

s = initial_state(mdp::NoCrashMDP, rng)

policy = Multilane.BehaviorPolicy(mdp, Multilane.NORMAL, false, rng)

hr = HistoryRecorder(rng=rng, max_steps=10000)

hist = simulate(hr, mdp, policy, s)
@time hist = simulate(hr, mdp, policy, s)
@show n_steps(hist)
hrbm = @benchmark simulate(hr, mdp, policy, s)
@show hrbm

ro = RolloutSimulator(rng=rng, max_steps=10000)
simulate(ro, mdp, policy, s)
@time simulate(ro, mdp, policy, s)
robm = @benchmark simulate(ro, mdp, policy, s)
@show robm

ext: jl


text:
using Test
using AbstractTrees

using FlatZincParser
using FlatZincParser: TokenStream, AbstractNode, DataNode, Node, next!, reset!, context
using FlatZincParser: match!, match_many!, match_token, match_many

@testset "lexer" begin
    @testset "inany" begin
        using FlatZincParser: inany
        @test !inany('a', 'b':'d', 'A':'Z')
        @test inany('a', 'a':'d', 'A':'Z')
        @test inany('B', 'a':'d', 'A':'Z')
    end

    @testset "read while" begin
        using FlatZincParser: read_while, innone
        iob = IOBuffer("123abc#12")
        @test read_while(iob, ('1':'3', 'a':'z')) == "123abc"
        @test read(iob, Char) == '#'

        iob = IOBuffer("123abc#12")
        @test read_while(iob, 'a':'z') == ""
        @test read_while(iob, ('1':'3', 'a':'z', '#')) == "123abc#12"
        @test eof(iob)
        iob = IOBuffer("123abc#12!")
        @test read_while(iob, ('1':'3', 'a':'z', '#')) == "123abc#12"

        iob = IOBuffer("123a\"bc#12")
        @test read_while(iob, '"'; condition=innone) == "123a"
        @test read(iob, Char) == '\"'
    end

    @testset "peek2" begin
        using FlatZincParser: peekn
        io = IOBuffer("12345")
        pos = position(io)
        @test peekn(1, io) == '1'
        @test peekn(2, io) == '2'
        @test position(io) == pos
    end

    @testset "integer parsing" begin
        tok = tokenize("0x23fA -0o123456 9 -15;")
        @test tok[1].lexme == "0x23fA"
        @test tok[2].lexme == "-0o123456"
        @test tok[3].lexme == "9"
        @test tok[4].lexme == "-15"
    end

    @testset "float parsing" begin
        tok = tokenize("1.0 -2.3 -1e-15 8e12 -1.23E+4;")
        @test tok[1].lexme == "1.0"
        @test tok[2].lexme == "-2.3"
        @test tok[3].lexme == "-1e-15"
        @test tok[4].lexme == "8e12"
        @test tok[5].lexme == "-1.23E+4"
    end

    @testset "string literal parsing" begin
        tok = tokenize("1.0 \"\"\"this is my string\"\"\" -1e-15 8e12 -1.23E+4;")
        @test tok[2].type === :string_literal
        @test tok[2].lexme == "this is my string"
    end

    @testset "eol parsing" begin
        ts = tokenize("1.2 2.6 0.2 ;")
        @test length(ts) == 4
    end

    @testset "tokenize files" begin
        small_file = joinpath(@__DIR__,"files", "queens3_4.fzn")
        big_file = joinpath(@__DIR__,"files", "2018_test-scheduling_t100m10r3-2.fzn")
        # huge_file = joinpath(@__DIR__,"files", "oocsp_racks_030_f7.fzn")
        for file in [small_file, big_file]
            @time tokens = open(tokenize, file)
            i, b, f = 0, 0, 0
            for t in tokens
                if t.type === :int_literal
                    parse(Int, t.lexme)
                    i += 1
                end
                if t.type === :bool_literal
                    parse(Bool, t.lexme)
                    b += 1
                end
                if t.type === :float_literal
                    parse(Float64, t.lexme)
                    f += 1
                end
            end
            println("Parsed $i integers, $b bools and $f floats")
        end
    end

    # using BenchmarkTools
    # big_file = joinpath(@__DIR__,"files", "2018_test-scheduling_t100m10r3-2.fzn")
    # string = open(big_file) do io
    #     read(io, String)
    # end;
    # @btime tokenize($string)
    # 23.758 ms (239903 allocations: 17.54 MiB)
end


text:
using VLKeggSDK
using DataFrames

# lookup the following compounds -
test_compund_array = [
    "cpd:C20820"
    "cpd:C19376"
    "cpd:C01799"
    "cpd:C00092"
    "cpd:C00085"
    "cpd:C00031"
    "cpd:C00002"
    "cpd:C00001"
    "cpd:C00003"
    "cpd:C00004"
    "cpd:C00092"
    "cpd:C00244"
];

# get the data for these compounds -
record_components = (get_compound_records(test_compund_array) |> check) |> DataFrame


text:
#using SatelliteToolbox # main.jlにてincludeしているため不要

#事前準備関数定義
function set_SPCal(Name, t0,dt)
    tles = read_tle(Name) #TLE読み込み
    #println("readTLE_Complete")
    orbp = init_orbit_propagator(Val{:sgp4}, tles[1]) #SGP4での軌道モデルの読み込み
    #println("OrbitModelDownload_Complete")
    eop_IAU2000A = get_iers_eop(:IAU2000A)  #IAU-2000Aでの地球モデルの読み込み
    #println("EarthModelDownload_Complete")
    o, r, v = propagate_to_epoch!(orbp, t0 - dt/(24*60*60))
    #println("1stStepOrbitCalcurate_Complete")
    return orbp, eop_IAU2000A, o, r, v, tles
end

#ループ毎処理関数定義
function SPCal(n, t0, dt, orbp, eop_IAU2000A)
    #=  衛星位置の赤道面座標系での計算　=#
    o, r, v = step!(orbp, dt)
    #=  衛星位置の地心直交座標系への変換  =#
    xg = rECItoECEF(GCRF(), ITRF(), t0 + n*dt/(24*60*60), eop_IAU2000A) * r
    vg = rECItoECEF(GCRF(), ITRF(), t0 + n*dt/(24*60*60), eop_IAU2000A) * v
    return xg, vg
end

function orbit_cal(DataNum,dt,start_time,TLEFileName,with_print=false)
	t0 = DatetoJD(start_time)
	JD_log = zeros(DataNum)
	x_ecef_log = zeros(DataNum, 3)    #衛星位置＠地心直交座標系
	x_geod_log = zeros(DataNum, 3)    #衛星位置＠測地座標系
	v_ecef_log = zeros(DataNum, 3)    #衛星速度＠地心直交座標系

	#シミュレーションループ
	orbp, eop_IAU2000A, o, r, v, tles = set_SPCal(TLEFileName, t0,dt)

	for i=1:DataNum
		JD_log[i] = DatetoJD(start_time + Second(i * dt))

		x_sat, v_sat = SPCal(i, t0, dt, orbp, eop_IAU2000A)
		x_ecef_log[i,:] = x_sat
		v_ecef_log[i,:] = v_sat

		if with_print == true
			#println(x_ecef_log[i,:])
			print((i-1)*dt)
			print(", ")
			print(x_ecef_log[i,1])
			print(", ")
			print(x_ecef_log[i,2])
			print(", ")
			println(x_ecef_log[i,3])
		end
		x_geo = ECEFtoGeodetic([x_sat[1]; x_sat[2]; x_sat[3]])
		x_geod_log[i,1]=x_geo[1]
		if x_geo[2]>0
			x_geod_log[i,2]=x_geo[2]
		else
			x_geod_log[i,2]=x_geo[2]+2π
		end
		x_geod_log[i,2]=x_geo[2]
		x_geod_log[i,3]=x_geo[3]
		#println(x_geod_log[i,:])
	end

	return JD_log,x_ecef_log, v_ecef_log, x_geod_log, tles, eop_IAU2000A
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
"""
    This provides units and creates a non-dimensionalization object
"""
module Units
using Unitful
import Unitful: superscript
using Parameters
using Setfield # allows modifying fields in immutable struct

import Base: show, isapprox, isequal, convert, length, size, getindex, setindex!, getproperty
import Base.Broadcast: broadcasted
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, PerpleX_LaMEM_Diagram

# Define additional units that are useful in geodynamics
@unit    Myrs  "Myrs"   MillionYears    1000000u"yr"    false

function __init__()
    Unitful.register(Units)
end

Unitful.register(Units)

# Define a number of useful units
const km    = u"km"
const m     = u"m"
const cm    = u"cm"
const mm    = u"mm"
const Myrs  = u"Myrs"
const yr    = u"yr"
const s     = u"s"
const kg    = u"kg"
const g     = u"g"
const Pa    = u"Pa"
const MPa   = u"MPa"
const kbar  = u"kbar"
const Pas   = u"Pa*s"
const K     = u"K"
const C     = u"°C"
const mol   = u"mol"
const kJ    = u"kJ"
const J     = u"J"
const Watt  = u"W"
const μW    = u"μW"
const μm    = u"μm"

export
    km, m, cm, mm, μm, Myrs, yr, s, MPa, Pa, kbar, Pas, K, C, g, kg, mol, J, kJ, Watt, μW,
    GeoUnit, GeoUnits, GEO_units, SI_units, NO_units, AbstractGeoUnit,
    nondimensionalize, dimensionalize,
    superscript, upreferred, GEO, SI, NONE, isDimensional, Value, NumValue, Unit, UnitValue,
    isdimensional,
    compute_units

include("unpack.jl")    # adds macros for unpacking GeoUnit variables with or w/out units

"""
AbstractGeoUnit

Abstract supertype for geo units.
"""
abstract type AbstractGeoUnit{TYPE, DIMENSIONAL} <: Number end

abstract type AbstractUnitType end

struct GEO <: AbstractUnitType end
struct SI  <: AbstractUnitType end
struct NONE<: AbstractUnitType end

# The GeoUnit struct encodes dimensional info in the type info
"""
    Structure that holds a GeoUnit parameter, and encodes the units and whether it is dimensional or not in the name.

    Having that is useful, as non-dimensionalization removes the units from a number and we thus no longer know how to transfer it back to the correct units.
    With the GeoUnit struct, this information is retained, and we can thus redimensionalize it at a later StepRange

"""
struct GeoUnit{T,U}
    val     ::  T
    unit    ::  U
    isdimensional::Bool
end

# Different ways of specifying the GeoUnit:
GeoUnit(val) = GeoUnit{typeof(ustrip.(val)), typeof(unit(val[1]))}(  ustrip.(val), unit(val[1]),
                           isa(val[1], Union{Unitful.FreeUnits, Unitful.Quantity} ) )

GeoUnit{T}(val) where {T,U}   = GeoUnit{T,typeof(unit(val[1]))}( T.(ustrip.(val)), unit(val[1]),
                           isa(val[1], Union{Unitful.FreeUnits, Unitful.Quantity} ) )

GeoUnit{T,U}(val) where {T,U}   = GeoUnit{T,U}( T.(ustrip.(val)), unit(val[1]),
                                            isa(val[1], Union{Unitful.FreeUnits, Unitful.Quantity} ) )

Base.convert(t::Type{GeoUnit{T,U}}, x::Quantity{T,U,X})  where {T,U,X} = GeoUnit(x)
Base.convert(t::Type{GeoUnit{T,U}}, x::T)  where {T,U} = GeoUnit(x)

# Define methods to deal with cases when the input has integers
GeoUnit(val::Union{Int64, AbstractArray{Int64}}) = GeoUnit(Float64.(val))
GeoUnit(val::Union{Int32, AbstractArray{Int32}}) = GeoUnit(Float32.(val))
GeoUnit(val::Union{Quantity{Int64}, AbstractArray{<:Quantity{<:Int64}}}) = GeoUnit(Float64.(val))
GeoUnit(val::Union{Quantity{Int32}, AbstractArray{<:Quantity{<:Int32}}}) = GeoUnit(Float32.(val))

# helper functions
Unit(v::GeoUnit{T,U}) where {T,U} =  Unitful.unit(v.unit*1)
isdimensional(v::GeoUnit{T,U})  where {T,U}     =   v.isdimensional                 # is it a nondimensional number or not?
isdimensional(v::Number)                        =   false                           # nope
NumValue(v::GeoUnit)                            =   v.val                           # numeric value, with no units
NumValue(v::Number)                             =   v                               # numeric value
NumValue(v::AbstractArray)                      =   v                               # numeric value
Value(v::GeoUnit)                               =   Unitful.Quantity.(v.val,v.unit)  # value, with units

function UnitValue(v::GeoUnit{T,U}) where {T,U}
    if v.isdimensional
        return Value(v)             # returns value with units
    else
        return NumValue(v.val)
    end
end

Base.isapprox(x::GeoUnit, y::Number; kwargs...) = Base.isapprox(x.val, y; kwargs...)
Base.isequal(x::GeoUnit,  y::Number) = Base.isequal(x.val, y)
Base.isequal(x::GeoUnit,  y::AbstractArray) = Base.isequal(x.val, y)
Base.isequal(x::GeoUnit,  y::GeoUnit) = Base.isequal(x.val, y.val)

Base.convert(::Type{<:AbstractArray},   v::GeoUnit)          =   v.val
Base.convert(::Type{<:Number},          v::GeoUnit)          =   v.val
Base.convert(::Type{GeoUnit},  v::Number)           =   GeoUnit(v)
Base.convert(::Type{GeoUnit},  v::Int32)            =   GeoUnit(Float32(v))
Base.convert(::Type{GeoUnit},  v::Int64)            =   GeoUnit(Float64(v))
Base.convert(::Type{GeoUnit},  v::Quantity)         =   GeoUnit(v)
Base.convert(::Type{GeoUnit},  v::AbstractArray)    =   GeoUnit(v)
Base.convert(::Type{GeoUnit},  v::AbstractArray{Int32}) =   GeoUnit(ustrip.(Float32.(v))*unit(v[1]))
Base.convert(::Type{GeoUnit},  v::AbstractArray{Int64}) =   GeoUnit(ustrip.(Float64.(v))*unit(v[1]))

Base.convert(::Type{GeoUnit{T}},  v::Quantity)      where T     =   GeoUnit(T.(v))
Base.convert(::Type{GeoUnit{T}},  v::Number)        where T     =   GeoUnit(T(v))
Base.convert(::Type{GeoUnit{T}},  v::AbstractArray) where T     =   GeoUnit(T.(v))

#Base.convert(::Type{GeoUnit{T,U}},  v::T)    where {T,U}    =   GeoUnit{T,typeof(unit(v[1]))}(v)

Base.promote_rule(::Type{GeoUnit}, ::Type{Quantity}) = GeoUnit

function Base.show(io::IO, x::GeoUnit{T,U})  where {T,U} # output
    val = x.val
    if x.isdimensional == true
        println("GeoUnit{dimensional, $(x.unit)}, ")
    else
        println("GeoUnit{nondimensional, $(x.unit)}, ")
    end
    show(io, MIME("text/plain"), val)
end

# define a few basic routines so we can easily operate with GeoUnits
Base.length(v::GeoUnit)         =   length(v.val)
Base.size(v::GeoUnit)           =   size(v.val)
Base.getindex(A::GeoUnit{T,U}, inds::Vararg{Int,N}) where {T,U,N} = A.val[inds...]

# Multiply with numnber
Base.:*(x::GeoUnit, y::Number)  =   x.val*y
Base.:+(x::GeoUnit, y::Number)  =   x.val+y
Base.:/(x::GeoUnit, y::Number)  =   x.val/y
Base.:-(x::GeoUnit, y::Number)  =   x.val-y

Base.:*(x::Number, y::GeoUnit)  = y.val*x
Base.:+(x::Number, y::GeoUnit)  = y.val+x
Base.:/(x::Number, y::GeoUnit)  = x/y.val
Base.:-(x::Number, y::GeoUnit)  = x-y.val

# Multiplying a GeoUnit with another one, returns a GeoUnit
Base.:+(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2} = GeoUnit(Value(x) + Value(y))
Base.:/(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2} = GeoUnit(Value(x) / Value(y))
Base.:*(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2} = GeoUnit(Value(x) * Value(y))
Base.:-(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2} = GeoUnit(Value(x) - Value(y))
Base.:^(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2} = GeoUnit(Value(x) ^ Value(y))

Base.:*(x::GeoUnit, y::Quantity)   = UnitValue(x).*y
Base.:/(x::GeoUnit, y::Quantity)   = UnitValue(x)./y
Base.:+(x::GeoUnit, y::Quantity)   = UnitValue(x) .+ y
Base.:-(x::GeoUnit, y::Quantity)   = UnitValue(x) .- y

Base.:*(x::Quantity, y::GeoUnit)   = x.*UnitValue(y)
Base.:/(x::Quantity, y::GeoUnit)   = x./UnitValue(y)
Base.:+(x::Quantity, y::GeoUnit)   = x .+ UnitValue(y)
Base.:-(x::Quantity, y::GeoUnit)   = x .- UnitValue(y)

# If we multiply a GeoUnit with an abstract array, we only return values, not units (use GeoUnits for that)
Base.:*(x::GeoUnit, y::AbstractArray)   = NumValue(x) .* y
Base.:/(x::GeoUnit, y::AbstractArray)   = broadcast(/,NumValue(x), y)
Base.:+(x::GeoUnit, y::AbstractArray)   = broadcast(+,NumValue(x), y)
Base.:-(x::GeoUnit, y::AbstractArray)   = broadcast(-,NumValue(x), y)

Base.:*(y::AbstractArray, x::GeoUnit)   = broadcast(*,y, NumValue(x))
Base.:/(y::AbstractArray, x::GeoUnit)   = broadcast(/,y, NumValue(x))
Base.:+(y::AbstractArray, x::GeoUnit)   = broadcast(+,y, NumValue(x))
Base.:-(y::AbstractArray, x::GeoUnit)   = broadcast(-,y, NumValue(x))

Base.:*(x::GeoUnit, y::AbstractArray{<:Quantity})   = broadcast(*,Value(x), y)
Base.:/(x::GeoUnit, y::AbstractArray{<:Quantity})   = broadcast(/,Value(x), y)
Base.:+(x::GeoUnit, y::AbstractArray{<:Quantity})   = broadcast(+,Value(x), y)
Base.:-(x::GeoUnit, y::AbstractArray{<:Quantity})   = broadcast(-,Value(x), y)

Base.:*(y::AbstractArray{<:Quantity}, x::GeoUnit)   = broadcast(*,y, Value(x))
Base.:/(y::AbstractArray{<:Quantity}, x::GeoUnit)   = broadcast(/,y, Value(x))
Base.:+(y::AbstractArray{<:Quantity}, x::GeoUnit)   = broadcast(+,y, Value(x))
Base.:-(y::AbstractArray{<:Quantity}, x::GeoUnit)   = broadcast(-,y, Value(x))

#for op in (+, *, -, ^,/)

# Broadcasting
Base.broadcasted(::typeof(+), A::GeoUnit,       B::AbstractArray)           = broadcast(+, NumValue(A), B)
Base.broadcasted(::typeof(*), A::GeoUnit,       B::AbstractArray)           = broadcast(*, NumValue(A), B)
Base.broadcasted(::typeof(-), A::GeoUnit,       B::AbstractArray)           = broadcast(-, NumValue(A), B)
Base.broadcasted(::typeof(/), A::GeoUnit,       B::AbstractArray)           = broadcast(/, NumValue(A), B)
Base.broadcasted(::typeof(^), A::GeoUnit,       B::AbstractArray)           = broadcast(^, NumValue(A), B)

Base.broadcasted(::typeof(+), A::AbstractArray, B::GeoUnit)                 = broadcast(+, A, NumValue(B))
Base.broadcasted(::typeof(*), A::AbstractArray, B::GeoUnit)                 = broadcast(*, A, NumValue(B))
Base.broadcasted(::typeof(-), A::AbstractArray, B::GeoUnit)                 = broadcast(-, A, NumValue(B))
Base.broadcasted(::typeof(/), A::AbstractArray, B::GeoUnit)                 = broadcast(/, A, NumValue(B))
Base.broadcasted(::typeof(^), A::AbstractArray, B::GeoUnit)                 = broadcast(^, A, NumValue(B))

Base.broadcasted(::typeof(+), A::GeoUnit,    B::AbstractArray{<:Quantity})  = broadcast(+, Value(A), B)
Base.broadcasted(::typeof(*), A::GeoUnit,    B::AbstractArray{<:Quantity})  = broadcast(*, Value(A), B)
Base.broadcasted(::typeof(-), A::GeoUnit,    B::AbstractArray{<:Quantity})  = broadcast(-, Value(A), B)
Base.broadcasted(::typeof(/), A::GeoUnit,    B::AbstractArray{<:Quantity})  = broadcast(/, Value(A), B)
Base.broadcasted(::typeof(^), A::GeoUnit,    B::AbstractArray{<:Quantity})  = broadcast(^, Value(A), B)

Base.broadcasted(::typeof(+), A::AbstractArray{<:Quantity}, B::GeoUnit)     = broadcast(+, A, Value(B))
Base.broadcasted(::typeof(*), A::AbstractArray{<:Quantity}, B::GeoUnit)     = broadcast(*, A, Value(B))
Base.broadcasted(::typeof(-), A::AbstractArray{<:Quantity}, B::GeoUnit)     = broadcast(-, A, Value(B))
Base.broadcasted(::typeof(/), A::AbstractArray{<:Quantity}, B::GeoUnit)     = broadcast(/, A, Value(B))
Base.broadcasted(::typeof(^), A::AbstractArray{<:Quantity}, B::GeoUnit)     = broadcast(^, A, Value(B))

Base.getindex(x::GeoUnit, i::Int64, j::Int64, k::Int64) = GeoUnit(x.val[i,j,k]*x.unit)
Base.getindex(x::GeoUnit, i::Int64, j::Int64)           = GeoUnit(x.val[i,j]*x.unit)
Base.getindex(x::GeoUnit, i::Int64)                     = GeoUnit(x.val[i]*x.unit)
Base.setindex!(A::GeoUnit{T,U}, val, inds::Vararg{Int,N}) where {T,U,N} = A.val[inds...] = val

"""
    GeoUnits

    Structure that holds parameters used for non-dimensionalization
"""
@with_kw_noshow struct GeoUnits{TYPE}
    # Selectable input parameters
    temperature     =   1               #   Characteristic temperature  [C or K]
    length          =   1               #   Characteristic length unit  [km, m or -]
    stress          =   1               #   Characteristic stress unit  [MPa, Pa or -]
    time            =   1               #   Characteristic time unit    [Myrs, s pr -]
    viscosity       =   1

    # primary characteristic units
    K               =   1;                      # temperature in SI units for material parameter scaling
    s               =   1;                      # time in SI units for material parameter scaling
    m               =   1;                      # length in SI units for material parameter scaling
    Pa              =   1;                      # stress in SI units
    kg              =   upreferred(Pa*m*s^2)    # compute mass from pascal. Note: this may result in very large values

    # Main SI units (used later for nondimensionalization)
    Length          =   m
    Mass            =   kg
    Time            =   s
    Temperature     =   K
    Amount          =   1mol
    Second          =   s

    # Not defined, as they are not common in geodynamics:
    #Current
    #Luminosity

    # Derived units
    N               =   kg*m/s^2        # Newton
    J               =   N*m             # Joule
    W               =   J/s             # Watt
    area            =   m^2             # area   in SI units for material parameter scaling
    volume          =   m^3             # volume
    velocity        =   m/s
    density         =   kg/m^3
    acceleration    =   m/s^2
    force           =   kg*m/s^2
    strainrate      =   1/s
    heatcapacity    =   J/kg/K
    conductivity    =   W/m/K

    # Helpful
    SecYear         =   3600*24*365.25
    Myrs            =   1e6
    cmYear          =   SecYear*100     # to transfer m/s -> cm/yr

end

"""
    GEO_units(;length=1000km, temperature=1000C, stress=10MPa, viscosity=1e20Pas)

Creates a non-dimensionalization object using GEO units.

GEO units implies that upon dimensionalization, `time` will be in `Myrs`, `length` in `km`, stress in `MPa`, etc.
which is more convenient for typical geodynamic simulations than SI units
The characteristic values given as input can be in arbitrary units (`km` or `m`), provided the unit is specified.

# Examples:
```julia-repl
julia> CharUnits = GEO_units()
Employing GEO units
Characteristic values:
         length:      1000 km
         time:        0.3169 Myrs
         stress:      10 MPa
         temperature: 1000.0 °C
julia> CharUnits.velocity
1.0e-7 m s⁻¹
```
If we instead have a crustal-scale simulation, it is likely more appropriate to use a different characteristic `length`:
```julia-repl
julia> CharUnits = GEO_units(length=10km)
Employing GEO units
Characteristic values:
         length:      10 km
         time:        0.3169 Myrs
         stress:      10 MPa
         temperature: 1000.0 °C
```
"""
function GEO_units(;length=1000km, temperature=1000C, stress=10MPa, viscosity=1e20Pas)

    if unit(temperature)==NoUnits;  temperature = temperature*C;        end
    if unit(length)==NoUnits;       length      = length*u"km";         end
    if unit(stress)==NoUnits;       stress      = stress*u"MPa";        end
    if unit(viscosity)==NoUnits;    viscosity   = viscosity*u"Pa*s";    end

    T       =   uconvert(C,      temperature)
    Le      =   uconvert(km,     length);
    Sigma   =   uconvert(MPa,    stress)
    Eta     =   uconvert(Pas,    viscosity)

    T_SI    =   uconvert(K,      T);
    Le_SI   =   uconvert(m,      Le);
    Sigma_SI=   uconvert(Pa,     Sigma)
    Time_SI =   Eta/Sigma_SI;
    t       =   uconvert(Myrs,   Time_SI)

    GeoUnits{GEO}(length=Le, temperature=T,      stress=Sigma,  viscosity=Eta, time=t,
                  m=Le_SI,   K=T_SI,   Pa=Sigma_SI,  s=Time_SI)
end

"""
    CharUnits = SI_units(length=1000m, temperature=1000K, stress=10Pa, viscosity=1e20)

Specify the characteristic values using SI units

# Examples:
```julia-repl
julia> CharUnits = SI_units(length=1000m)
Employing SI units
Characteristic values:
         length:      1000 m
         time:        1.0e19 s
         stress:      10 Pa
         temperature: 1000.0 K
```
Note that the same can be achieved if the input is given in `km`:
```julia-repl
julia> CharUnits = SI_units(length=1km)
```
"""
function SI_units(;length=1000m, temperature=1000K, stress=10Pa, viscosity=1e20Pas)

    if unit(temperature)==NoUnits;  temperature = temperature*K;     end
    if unit(length)==NoUnits;       length      = length*u"m";         end
    if unit(stress)==NoUnits;       stress      = stress*u"Pa";        end
    if unit(viscosity)==NoUnits;    viscosity   = viscosity*u"Pa*s";    end

    T       =   uconvert(K,     temperature)
    Le      =   uconvert(m,     length);
    Sigma   =   uconvert(Pa,    stress)
    Eta     =   uconvert(Pas,   viscosity)

    T_SI    =   uconvert(K,      T);
    Le_SI   =   uconvert(m,      Le);
    Sigma_SI=   uconvert(Pa,     Sigma)
    Time_SI =   Eta/Sigma_SI;
    t       =   uconvert(s,      Time_SI)

    GeoUnits{SI}(length=Le, temperature=T,      stress=Sigma,  viscosity=Eta, time=t,
                 m=Le_SI,   K=T_SI,   Pa=Sigma_SI,  s=Time_SI)
end

"""
    CharUnits = NO_units(length=1, temperature=1, stress=1, viscosity=1)

Specify the characteristic values in non-dimensional units

# Examples:
```julia-repl
julia> using GeoParams;
julia> CharUnits = NO_units()
Employing NONE units
Characteristic values:
         length:      1
         time:        1.0
         stress:      1
         temperature: 1.0
```
"""
function NO_units(;length=1, temperature=1, stress=1, viscosity=1)

    if unit(temperature)!=NoUnits;  error("temperature should not have units")    end
    if unit(length)!=NoUnits;       error("length should not have units")    end
    if unit(stress)!=NoUnits;       error("stress should not have units")    end
    if unit(viscosity)!=NoUnits;    error("viscosity should not have units")    end

    T       =   temperature
    Le      =   length;
    Sigma   =   stress
    Eta     =   viscosity
    Time    =   Eta/Sigma;

    GeoUnits{NONE}(length=Le, temperature=T, stress=Sigma, viscosity=Eta, time=Time,
                 m=Le, K=T, Pa=Sigma, s=Time)
end

"""
    nondimensionalize(param, CharUnits::GeoUnits{TYPE})

Nondimensionalizes `param` using the characteristic values specified in `CharUnits`

# Example 1
```julia-repl
julia> using GeoParams;
julia> CharUnits =   GEO_units();
julia> v         =   3cm/yr
3 cm yr⁻¹
julia> v_ND      =   nondimensionalize(v, CharUnits)
0.009506426344208684
```
# Example 2
In geodynamics one sometimes encounters more funky units
```julia-repl
julia> CharUnits =   GEO_units();
julia> A         =   6.3e-2MPa^-3.05*s^-1
0.063 MPa⁻³·⁰⁵ s⁻¹
julia> A_ND      =   nondimensionalize(A, CharUnits)
7.068716262102384e14
```

In case you are interested to see how the units of `A` look like in different units, use this function from the [Unitful](https://github.com/PainterQubits/Unitful.jl) package:
```julia-repl
julia> uconvert(u"Pa^-3.05*s^-1",A)
3.157479571851836e-20 Pa⁻³·⁰⁵
```
and to see it decomposed in the basic `SI` units of length, mass and time:
```julia-repl
julia> upreferred(A)
3.1574795718518295e-20 m³·⁰⁵ s⁵·¹ kg⁻³·⁰⁵
```
"""
function nondimensionalize(param::GeoUnit{T,U}, g::GeoUnits{TYPE}) where {T,U,TYPE}

    if param.isdimensional
        char_val      = compute_units(param, g);
    #    @show char_val, typeof(char_val)
        val_ND        = upreferred.(param.val*param.unit)/char_val
        val_ND_number = convert(T,val_ND)
        param_ND      = GeoUnit{T,U}(val_ND_number, param.unit, false)          # store new value, but keep original dimensions

    else
        param_ND = param;

    end
    return param_ND
end

# in case the parameter is already non-dimensional:
nondimensionalize(param::String, g::GeoUnits{TYPE}) where {TYPE} = param

# in case it is a unitful quantity
function nondimensionalize(param::Union{Unitful.Quantity{T,K,M}, AbstractArray{<:Quantity{T,K,M}}}, g::GeoUnits{TYPE}) where {TYPE,T,K,M}

    param_Geo   =   GeoUnit(param)
    result      =   nondimensionalize(param_Geo,g)

    return UnitValue(result)

end

# If it is an array, but has no units we cannot know how to nondimensionalize it
nondimensionalize(param::AbstractArray{<:Number}, g::GeoUnits{TYPE}) where {TYPE} = param

# This computes the characteristic value
function compute_units(param::GeoUnit{<:Union{T,AbstractArray{T}},U}, g::GeoUnits{TYPE}) where {T,U,TYPE}
    dim         =   Unitful.dimension(param.unit);              # Basic SI units

    # Note: `char_val`` is type unstable within the routine
    # Yet, the dimensions at the end are known as they must be the same as `U`
    #  (perhaps converted to preferred units such as m instead of km)
    char_val    =   1.0
    foreach((typeof(dim).parameters[1])) do y
        val = upreferred(getproperty(g, Unitful.name(y)))       # Retrieve the characteristic value from structure g
        pow = Float64(y.power)                                  # power by which it should be multiplied
        char_val *= val^pow                                     # multiply characteristic value
    end

    value::T = ustrip(char_val)                                 # numerical value
    char_val_out = upreferred(param.unit)*value;                # this is done for type-stability

    return char_val_out
end

"""
    MatParam_ND = nondimensionalize(MatParam::AbstractMaterialParam, CharUnits::GeoUnits{TYPE})

Non-dimensionalizes a material parameter structure (e.g., Density, CreepLaw)

"""
function nondimensionalize(MatParam::AbstractMaterialParam, g::GeoUnits{TYPE}) where {TYPE}
    for param in fieldnames(typeof(MatParam))
        if isa(getfield(MatParam, param), GeoUnit)
            z=getfield(MatParam, param)

            # non-dimensionalize:
            z = nondimensionalize(z, g)

            # Replace field (using Setfield package):
            MatParam = set(MatParam, Setfield.PropertyLens{param}(), z)

        end
    end
    return MatParam
end

"""
    nondimensionalize(phase_mat::MaterialParams, g::GeoUnits{TYPE})

nondimensionalizes all fields within the Material Parameters structure that contain material parameters
"""
function nondimensionalize(phase_mat::AbstractMaterialParamsStruct, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(phase_mat))
        fld = getfield(phase_mat, param)
        if length(fld)>0
            if typeof(fld[1]) <: AbstractPhaseDiagramsStruct

                # in case we employ a phase diagram
                temp = PerpleX_LaMEM_Diagram(fld[1].Name, CharDim = g)
                fld_new = (temp,)
                #setfield!(phase_mat, param, fld_new)
                phase_mat = set(phase_mat, Setfield.PropertyLens{param}(), fld_new)

            else
                # otherwise non-dimensionalize all fields and create a new tuple
                id = findall(isa.(fld,AbstractMaterialParam))
                if length(id)>0
                    # Create a new tuple with non-dimensionalized fields:
                    fld_new = ntuple(i -> Units.nondimensionalize(fld[id[i]],g), length(id) )
#                    setfield!(phase_mat, param, fld_new)        # to be changed for immutable struct
                    phase_mat = set(phase_mat, Setfield.PropertyLens{param}(), fld_new)

                end
            end
        end
    end
    # phase_mat.Nondimensional = true
    phase_mat = set(phase_mat, Setfield.PropertyLens{:Nondimensional}(), true)
    return phase_mat
end

"""
    dimensionalize(param, param_dim::Unitful.FreeUnits, CharUnits::GeoUnits{TYPE})

Dimensionalizes `param` into the dimensions `param_dim` using the characteristic values specified in `CharUnits`.

# Example
```julia-repl
julia> CharUnits =   GEO_units();
julia> v_ND      =   nondimensionalize(3cm/yr, CharUnits)
0.031688087814028945
julia> v_dim     =   dimensionalize(v_ND, cm/yr, CharUnits)
3.0 cm yr⁻¹
```

"""
function dimensionalize(param_ND, param_dim::Unitful.FreeUnits, g::GeoUnits{TYPE}) where {TYPE}

    char_val    =   compute_units(GeoUnit(1.0*param_dim), g);         # Determine characteristic units
    param       =   uconvert.(param_dim, param_ND*char_val)

    return param
end

function dimensionalize(param_ND::GeoUnit{T,U}, g::GeoUnits{TYPE}) where {T,U,TYPE}

    if isdimensional(param_ND)==false
        char_val    = compute_units(param_ND, g);                       # Determine characteristic units
        val         = uconvert.(param_ND.unit, param_ND.val*char_val)   # dimensionalize
        param       = GeoUnit{T,U}(ustrip.(val), param_ND.unit, true)   # store new value, but keep original dimensions

        return param
    else
        return param_ND
    end
end

"""
    dimensionalize(MatParam::AbstractMaterialParam, CharUnits::GeoUnits{TYPE})

Dimensionalizes a material parameter structure (e.g., Density, CreepLaw)

"""
function dimensionalize(MatParam::AbstractMaterialParam, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(MatParam))
        if isa(getfield(MatParam, param), GeoUnit)
            z       =   getfield(MatParam, param)
            z_dim   =   dimensionalize(z, g)

            # Replace field (using Setfield package):
            MatParam =  set(MatParam, Setfield.PropertyLens{param}(), z_dim)
        end
    end
    return MatParam
end

"""
    Dimensionalize(phase_mat::MaterialParams, g::GeoUnits{TYPE})

Dimensionalizes all fields within the Material Parameters structure that contain material parameters
"""
function dimensionalize(phase_mat::AbstractMaterialParamsStruct, g::GeoUnits{TYPE}) where {TYPE}

    for param in fieldnames(typeof(phase_mat))
        fld = getfield(phase_mat, param)
        if length(fld)>0
            id = findall(isa.(fld,AbstractMaterialParam))
            if length(id)>0
                # Create a new tuple with non-dimensionalized fields:
                fld_new     = ntuple(i -> Units.dimensionalize(fld[id[i]],g), length(id) )

                phase_mat   = set(phase_mat, Setfield.PropertyLens{param}(), fld_new)
            end
        end
    end
    phase_mat = set(phase_mat, Setfield.PropertyLens{:Nondimensional}(), false)

    return phase_mat
end

"""
    isDimensional(MatParam::AbstractMaterialParam)

`true` if MatParam is in dimensional units.
"""
function isDimensional(MatParam::AbstractMaterialParam)
    isDim = false;
    for param in fieldnames(typeof(MatParam))
        if isa(getfield(MatParam, param), GeoUnit)
            z=getfield(MatParam, param)
            if z.isdimensional==true
                isDim=true;
            end
        end
    end
    return isDim
end

# Define a view for the GEO_Units structure
function show(io::IO, g::GeoUnits{TYPE})  where {TYPE}
    print(io, "Employing $TYPE units \n",
              "Characteristic values: \n",
              "         length:      $(g.length)\n",
              "         time:        $(round(ustrip(g.time),digits=4)) $(unit(g.time))\n",
              "         stress:      $(g.stress)\n",
              "         temperature: $(Float64(g.temperature))\n")
end

# This replaces the viewer of the Unitful package, such that the printing of units is done as floats (better)
function Unitful.superscript(i::Rational{Int64}; io=nothing)
    string(superscript(float(i)))
    if io === nothing
        iocontext_value = nothing
    else
        iocontext_value = get(io, :fancy_exponent, nothing)
    end
    if iocontext_value isa Bool
        fancy_exponent = iocontext_value
    else
        v = get(ENV, "UNITFUL_FANCY_EXPONENTS", Sys.isapple() ? "true" : "false")
        t = tryparse(Bool, lowercase(v))
        fancy_exponent = (t === nothing) ? false : t
    end
    if fancy_exponent
        return superscript(float(i))
    else
        return  "^" * string(float(i))
    end

end

Unitful.superscript(i::Float64) = map(repr(i)) do c
    c == '-' ? '\u207b' :
    c == '1' ? '\u00b9' :
    c == '2' ? '\u00b2' :
    c == '3' ? '\u00b3' :
    c == '4' ? '\u2074' :
    c == '5' ? '\u2075' :
    c == '6' ? '\u2076' :
    c == '7' ? '\u2077' :
    c == '8' ? '\u2078' :
    c == '9' ? '\u2079' :
    c == '0' ? '\u2070' :
    c == '.' ? '\u0387' :
    error("unexpected character")
end

end


text:
function __init__()
    register(
        DataDep(
            "BioSemiTestFiles",
            "Manafacturer provided example files",
            ["https://www.biosemi.com/download/BDFtestfiles.zip"];
            post_fetch_method = [file -> run(`unzip $file`)],
        ),
    )

    register(
        DataDep(
            "ExampleSSR",
            "Steady state response data with few channels",
            [
                "https://github.com/rob-luke/Neuroimaging.jl-example-data/archive/refs/heads/master.zip",
            ];
            post_fetch_method = [file -> run(`unzip $file`)],
        ),
    )

    register(
        DataDep(
            "BioSemi64Locations",
            "ELP containing location of 64 channel biosemi sensors",
            [
                "https://raw.githubusercontent.com/rob-luke/NeuroimagingFiles/main/biosemi64.elp",
            ],
        ),
    )
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHopHamiltonianKitaev <: AbstractBondHamiltonian{L,1}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION OF HOP-HAMILTONIAN
#
################################################################################

mutable struct BondHopHamiltonianKitaev{L} <: AbstractBondHamiltonian{L,1}  # L is the label

    # couplings
    Jx :: Float64
    Jy :: Float64
    Jz :: Float64

    # bond labels of these couplings
    x_bonds :: Vector{L}
    y_bonds :: Vector{L}
    z_bonds :: Vector{L}

end

#export the type
export BondHopHamiltonianKitaev

#get the bond term of the hamiltonian

function bondtermhopping(
            h :: BondHopHamiltonianKitaev{L},
            b :: AbstractBond{L,NB}
        ) :: Matrix{Complex} where {L,NB}

    # get the bond label
    l = label(b)

    # define an empty coupling matrix
    coupling_matrix = zeros(1,1)

    # check for all neighbors if the label is in one of the categories
    if l in h.x_bonds
        coupling_matrix[1,1] = h.Jx
        return coupling_matrix
    elseif l in h.y_bonds
        coupling_matrix[1,1] = h.Jy
        return coupling_matrix
    elseif l in h.z_bonds
        coupling_matrix[1,1] = h.Jz
        return coupling_matrix
    end

    # if no bond found, return 0
    return coupling_matrix
end

# FUNCTIONS TO FIND BOND LABELS
function getBondLabelsKitaevX(
        all_couplings :: Vector{L}
    ) :: Vector{L} where {L}

    # all couplings that somehow contain a x (x-bond corresponds to label=1)

    return filter(c->c==1, all_couplings)
end
function getBondLabelsKitaevY(
        all_couplings :: Vector{L}
    ) :: Vector{L} where {L}

    # all couplings that somehow contain a y (y-bond corresponds to label=2)

    return filter(c->c==2, all_couplings)
end
function getBondLabelsKitaevZ(
        all_couplings :: Vector{L}
    ) :: Vector{L} where {L}

    # all couplings that somehow contain a z (z-bond corresponds to label=3)

    return filter(c->c==3, all_couplings)
end

#construction function
function getHopHamiltonianKitaev(
            unitcell :: U
        ) :: BondHopHamiltonianKitaev{L} where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B}}

    # obtain all couplings
    couplings = unique!(label.(bonds(unitcell)))

    # create default couplings
    Jx = 1.0
    Jy = 1.0
    Jz = 1.0

    # create list of labels
    x_bonds = getBondLabelsKitaevX(couplings)
    y_bonds = getBondLabelsKitaevY(couplings)
    z_bonds = getBondLabelsKitaevZ(couplings)

    # create and return a new object
    return BondHopHamiltonianKitaev{L}(Jx,Jy,Jz,x_bonds,y_bonds,z_bonds)
end

# export the construction function
export getHopHamiltonianKitaev


text:
module stubby_corbel_examples
using FinEtools
using FinEtools.AlgoBaseModule: evalconvergencestudy
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule: linearstatics, exportstresselementwise, exportstress
using Statistics: mean
using LinearAlgebra: Symmetric, cholesky

# Isotropic material
E=1000.0;
nu=0.4999; #Taylor data
W=25.0;
H=50.0;
L= 50.0;
htol = minimum([L,H,W])/1000;
uzex =-12.6;
magn = 0.2*uzex/4;
Force =magn*W*H*2;
CTE = 0.0
n = 5 #

function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    copyto!(forceout, [0.0; 0.0; magn])
end

function stubby_corbel_H8_by_hand()
    elementtag = "H8"
    println("""
    Stubby corbel example. Element: $(elementtag)
    """)

    fens,fes = H8block(W, L, H, n, 2*n, 2*n)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sectionL = selectelem(fens, bfes; facing=true, direction = [0.0 +1.0 0.0])
    # 0 cross-section surface  for the reactions
    section0 = selectelem(fens, bfes; facing=true, direction = [0.0 -1.0 0.0])
    # 0 cross-section surface  for the reactions
    sectionlateral = selectelem(fens, bfes; facing=true, direction = [1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    lx0 = connectednodes(subset(bfes, section0))
    setebc!(u,lx0,true,1,0.0)
    setebc!(u,lx0,true,2,0.0)
    setebc!(u,lx0,true,3,0.0)
    lx1 = connectednodes(subset(bfes, sectionlateral))
    setebc!(u,lx1,true,1,0.0)
    applyebc!(u)
    numberdofs!(u)
    println("u.nfreedofs = $(u.nfreedofs)")

    fi = ForceIntensity(FFlt, 3, getfrcL!);
    el2femm = FEMMBase(IntegDomain(subset(bfes, sectionL), GaussRule(2, 2)))
    F2 = distribloads(el2femm, geom, u, fi, 2);
    associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)
    @time K=cholesky(K)
    @time U = K\(F2)
    scattersysvec!(u,U[:])
    @show length(U)
    Tipl = selectnode(fens, box=[0 W L L 0 H], inflate=htol)
    utip = mean(u.values[Tipl, 3], dims=1)
    println("Deflection: $(utip), compared to $(uzex)")

    File =  "stubby_corbel_H8_by_hand.vtk"
    vtkexportmesh(File, fens, fes;  vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
    # modeldata = exportstresselementwise(modeldata)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
    # "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
    # "component"=>collect(1:6))
    # modeldata = exportstresselementwise(modeldata)
    # stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

    true

end # stubby_corbel_H8_by_hand

function allrun()
    println("#####################################################")
    println("# stubby_corbel_H8_by_hand ")
    stubby_corbel_H8_by_hand()
    return true
end # function allrun

end # module stubby_corbel_examples


text:
using LinearAlgebra
using Random
# import Random.rand
# import Random.rand!

function rand_orthomtx!(mtx::AbstractMatrix{T}) where T
    mtx .= qr(rand(T, size(mtx)...)).Q
end

rand(nsolt::AbstractNsolt, args...; kwargs...) = rand!(similar(nsolt), args...; kwargs...)

function rand!(nsolt::CnsoltTypeI{T,D}; isInitMat=true, isPropMat=true, isPropAng=true, isSymmetry=true) where {T,D}
    if isSymmetry
        nsolt.Φ .= Diagonal(cis.(rand(nsolt.nChannels)))
    end

    if isInitMat
        rand_orthomtx!(nsolt.V0)
    end

    for d = 1:D
        if isPropMat
            foreach(nsolt.Wdks[d], nsolt.Udks[d]) do W, U
                rand_orthomtx!(W)
                rand_orthomtx!(U)
            end
        end
        if isPropAng
            foreach(nsolt.θdks[d]) do θ
                θ .= 2pi*(rand(size(θ)...) .- 0.5)
            end
        end
    end
    return nsolt
end

function rand!(nsolt::CnsoltTypeII{T,D}; isInitMat=true, isPropMat=true, isPropAng=true, isSymmetry=true) where {T,D}
    if isSymmetry
        nsolt.Φ .= Diagonal(cis.(rand(nsolt.nChannels)))
    end

    if isInitMat
        rand_orthomtx!(nsolt.V0)
    end

    for d = 1:D
        if isPropMat
            foreach(nsolt.Wdks[d], nsolt.Udks[d], nsolt.Ŵdks[d], nsolt.Ûdks[d]) do W, U, Ŵ, Û
                rand_orthomtx!(W)
                rand_orthomtx!(U)
                rand_orthomtx!(Ŵ)
                rand_orthomtx!(Û)
            end
        end
        if isPropAng
            foreach(nsolt.θ1dks[d], nsolt.θ2dks[d]) do θ1, θ2
                θ1 .= 2pi*(rand(size(θ1)...) .- 0.5)
                θ2 .= 2pi*(rand(size(θ2)...) .- 0.5)
            end
        end
    end
    return nsolt
end

function rand!(nsolt::RnsoltTypeI{T,D}; isInitMat=true, isPropMat=true, isPropAng=true, isSymmetry=true) where {T,D}
    if isInitMat
        rand_orthomtx!(nsolt.W0)
        rand_orthomtx!(nsolt.U0)
    end

    for d = 1:D
        if isPropMat
            foreach(nsolt.Udks[d]) do U
                rand_orthomtx!(U)
            end
        end
    end
    return nsolt
end

function rand!(nsolt::RnsoltTypeII{T,D}; isInitMat=true, isPropMat=true, isPropAng=true, isSymmetry=true) where {T,D}

    if isInitMat
        rand_orthomtx!(nsolt.W0)
        rand_orthomtx!(nsolt.U0)
    end

    for d = 1:D
        if isPropMat
            foreach(nsolt.Wdks[d], nsolt.Udks[d]) do W, U
                rand_orthomtx!(W)
                rand_orthomtx!(U)
            end
        end
    end
    return nsolt
end

function intervals(lns::AbstractVector, offset::Integer=0)
    map(lns, cumsum(lns)) do st, ed
        UnitRange(ed - st + 1 + offset, ed + offset)
    end
end

intervals(lns::Tuple, args...) = (intervals(collect(lns), args...)...,)

function vm1constraint!(ns::RnsoltTypeI)
    θ, μ = mat2rotations(ns.W0)
    θ[1:ns.nChannels[1]-1] .= 0
    ns.W0 .= rotations2mat(θ, μ)

    ns
end

function vm1constraint!(ns::CnsoltTypeI)
    hP = fld(ns.nChannels, 2)
    ws = foldl(ns.Wdks, init = Matrix(I, hP, hP)) do wtmp, mts
        prod([ reverse(mts)..., wtmp ])
    end
    θ, μ = mat2rotations(ns.V0)
    θ[1:ns.nChannels-1] .= 0
    ns.V0 .= cat(ws', Matrix(I, size(ws)...), dims=[1,2]) * rotations2mat(θ, μ)

    ns
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
"""
A graph to test shortest_path functions.

The shortest path of this graph in terms of loss is:
A(1) - M3(5) - B(2)
With costs of loss:10.0, Z:2.0

The shortest path in terms of loss is:
A(1) - M1(3) - B(2)
With costs of loss:2.0, Z:10.0
"""

shortest_path_test = QNetwork()
A = BasicNode("A")
B = BasicNode("B")
M1 = BasicNode("M1")
M2 = BasicNode("M2")
M3 = BasicNode("M3")

AM1 = BasicChannel(A, M1, exp_cost=false)
AM2 = BasicChannel(A, M2, exp_cost=false)
AM3 = BasicChannel(A, M3, exp_cost=false)
M1B = BasicChannel(M1, B, exp_cost=false)
M2B = BasicChannel(M2, B, exp_cost=false)
M3B = BasicChannel(M3, B, exp_cost=false)

AM1.costs = Dict("loss"=>1.0, "Z"=>5.0)
AM2.costs = Dict("loss"=>3.0, "Z"=>3.0)
AM3.costs = Dict("loss"=>5.0, "Z"=>1.0)
M1B.costs = Dict("loss"=>1.0, "Z"=>5.0)
M2B.costs = Dict("loss"=>3.0, "Z"=>3.0)
M3B.costs = Dict("loss"=>5.0, "Z"=>1.0)

for i in [A, B, M1, M2, M3, AM1, AM2, AM3, M1B, M2B, M3B]
    add(shortest_path_test, i)
end

refresh_graph!(shortest_path_test)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
using Optim
import Optim: optimize

# Fitting logistic regression
#
# TODO: add support for dual formulation
struct LogisticOptimizer
    penalty::AbstractPenalty
    dual::Bool
    fit_intercept::Bool
    algo::Optim.AbstractOptimizer
    options::Optim.Options
end

function LogisticOptimizer(; penalty=L2Penalty(0.0),
                             dual=false,
                             fit_intercept=false,
                             algo=Optim.BFGS(),
                             options=Optim.Options())
    return LogisticOptimizer(penalty, dual, fit_intercept, algo, options)
end

function LogisticRegressor(X::AbstractArray{T, 2}, y::AbstractVector{T},
                           logopt::LogisticOptimizer) where T
    Dataset = issparse(X) ? SparseLogitData : LogitData
    data = Dataset(X, y)
    return LogisticRegressor(data, logopt.penalty, NULL_HASH, logopt.fit_intercept, NoLogger())
end

function Optim.optimize(logreg::LogisticRegressor,
                        x0::AbstractVector,
                        algo::Optim.AbstractOptimizer,
                        options::Optim.Options)
    f, g!, _, _ = generate_callbacks(logreg)
    res = Optim.optimize(f, g!, x0, algo, options)
    return res
end

function fit!(logopt::LogisticOptimizer,
              X::AbstractArray{T, 2},
              y::AbstractVector{T},
              x0=zeros(size(X, 2))) where T
    logreg = LogisticRegressor(X, y, logopt)
    res = Optim.optimize(logreg, x0, logopt.algo, logopt.options)
    return res
end

# Fitting optimal regression paramater (ℓ2)
abstract type AbstractL2PenaltyOptimizer end

struct L2PenaltyOptimizer
    kfold::Int
    shuffle::Bool
    fit_intercept::Bool
    outer_algo::Optim.AbstractOptimizer
    inner_algo::Optim.AbstractOptimizer
    options::Optim.Options
end
function L2PenaltyOptimizer( ; kfold=5,
                               shuffle=false,
                               fit_intercept=false,
                               outer_algo=BFGS(),
                               inner_algo=LBFGS(),
                               options=Optim.Options()
                            )
    return L2PenaltyOptimizer(kfold, shuffle, fit_intercept, outer_algo,
                              inner_algo, options)
end

function crossvalid_loss(X, y, γ, penopt::L2PenaltyOptimizer)
    kfold = penopt.kfold
    n, p = size(X)
    # Set-up permutation for cross validation
    if penopt.shuffle
        folds = randperm(n)
    else
        folds = vec(1:n)
    end
    #
    modk = div(n, kfold)
    start = 1
    stop = modk

    # Allocate a Hessian matrix
    ∇²f = zeros(p, p)
    # Get a view on upper-triangular
    h_train = @view ∇²f[triul(p)]
    # Allocate a gradient
    g_test = zeros(p)

    cost = 0.0
    ∇c   = 0.0
    w♯ = zeros(p)

    for i in 1:kfold
        # Split training and testing sets
        ind_train = folds[start:stop]
        ind_test = setdiff(1:n, ind_train)
        start += modk
        # Be careful we stay inside authorized indexes
        stop = min(stop + modk, n)
        Xtrain, ytrain = X[ind_train, :], y[ind_train]
        Xtest, ytest  = X[ind_test, :], y[ind_test]

        # First step: Training logreg
        logtrain = LogisticRegressor(Xtrain, ytrain, fit_intercept=false,
                                     penalty=L2Penalty(γ))
        # Inner optimization
        res = Optim.optimize(logtrain, w♯,
                            penopt.inner_algo,
                            penopt.options)
        p = res.minimizer
        # Compute Hessian
        hessian!(h_train, p, logtrain)

        # Second step: get prediction on testing subset
        # No penalty for testing!
        logtest = LogisticRegressor(Xtest, ytest, fit_intercept=false)
        loss_test = loss(p, logtest)
        # Compute gradient on testing set
        gradient!(g_test, p, logtest)

        # Get gradient w.r.t. current penalty
        ∇θ = - Symmetric(∇²f) \  p
        ∇c += 1 / kfold * dot(g_test, ∇θ)
        # Update cost
        cost += 1 / kfold * loss_test

        # Set hotstart for future computation
        w♯ .= p
    end
    return cost, ∇c
end

function fit!(penopt::L2PenaltyOptimizer,
              X::AbstractArray{T, 2},
              y::AbstractVector{T},
              γ0::Real) where T
    # Compute loss and gradient in a closure for Optim
    function fg!(f, g, x)
        c, ∇c = crossvalid_loss(X, y, x[1], penopt)
        if g != nothing
            g[1] = ∇c
        end
        if f != nothing
            return c
        end
    end
    res = Optim.optimize(Optim.only_fg!(fg!), [γ0], penopt.outer_algo, penopt.options)
    return res.minimum, res.minimizer[1]
end


text:
# Sometimes, it's a good idea to appreciate just how big the ocean is. Using
# the Manhattan distance, how far apart do the scanners get?
#
# In the above example, scanners 2 (1105,-1205,1229) and 3 (-92,-2380,-20) are
# the largest Manhattan distance apart. In total, they are 1197 + 1175 + 1249 =
# 3621 units apart.
#
# What is the largest Manhattan distance between any two scanners?

# read data
iter = eachline("input.txt")
scanners = Vector{Tuple{Int32, Int32, Int32}}[]
while !isnothing(iterate(iter))
  # while condition will consume the "--- scanner X ---" lines
  beacons = map(l -> tuple(parse.(Int32, split(l, ','))...), Iterators.takewhile(l -> l != "", iter))
  push!(scanners, beacons)
end

# Utility iterator to compute all combinations of 2 indices
struct Permutations
  len::Int
end

function Base.iterate(P::Permutations, (i,j)=(1,1))
  j += 1
  if j > P.len
    i += 1
    j = i + 1
    j > P.len && return nothing
  end
  idxs = (i,j)
  return (idxs, idxs)
end

Base.length(P::Permutations) = binomial(P.len, 2)

Base.eltype(P::Permutations) = Tuple{Int, Int}

# Distance is the sqrt of the sum of squares of differences. I'm going to leave
# off the sqrt bit so I'm not beholden to floating point errors.
function distance(scanneridx, i, j)
  (x1,y1,z1) = scanners[scanneridx][i]
  (x2,y2,z2) = scanners[scanneridx][j]
  (widen(x2)-x1)^2 + (widen(y2)-y1)^2 + (widen(z2)-z1)^2
end

distance(scanneridx, (i,j)) = distance(scanneridx, i, j)

distance(scanneridx) = (p) -> distance(scanneridx, p)

# Utility thot finds the distance between all beacons in the a scanner. This is
# used to find points that potentially match between two scanners.
struct RefDistances
  idxs::Vector{Tuple{Int, Int}}
  dists::Vector{Int}
  function RefDistances(baseidx)
    base = scanners[baseidx]
    idxs = collect(Permutations(length(base)))
    dists = map(distance(1), idxs)
    sortedidx = sortperm(dists)
    permute!(idxs, sortedidx)
    permute!(dists, sortedidx)
    new(idxs, dists)
  end
end

Base.getindex(R::RefDistances, i::Int) = R.idxs[i]

function Base.findfirst(R::RefDistances, dist)
  distidx = searchsortedfirst(R.dists, dist)
  distidx <= length(R.dists) && R.dists[distidx] == dist ? distidx : nothing
end

# Functions to re-orient points
const reorientfuncs = [
  identity, ((x,y,z),) -> (-x,-y,z), ((x,y,z),) -> (-x,-z,-y),
  ((x,y,z),) -> (-x,y,-z), ((x,y,z),) -> (-x,z,y), ((x,y,z),) -> (-y,-x,-z),
  ((x,y,z),) -> (-y,-z,x), ((x,y,z),) -> (-y,x,z), ((x,y,z),) -> (-y,z,-x),
  ((x,y,z),) -> (-z,-x,y), ((x,y,z),) -> (-z,-y,-x), ((x,y,z),) -> (-z,x,-y),
  ((x,y,z),) -> (-z,y,x), ((x,y,z),) -> (x,-y,-z), ((x,y,z),) -> (x,-z,y),
  ((x,y,z),) -> (x,z,-y), ((x,y,z),) -> (y,-x,z), ((x,y,z),) -> (y,-z,-x),
  ((x,y,z),) -> (y,x,-z), ((x,y,z),) -> (y,z,x), ((x,y,z),) -> (z,-x,-y),
  ((x,y,z),) -> (z,-y,x), ((x,y,z),) -> (z,x,y), ((x,y,z),) -> (z,y,-x),
]

# Utility that re-orients points and moves them so they are relative to the
# base scanner
struct Positioner
  reorientfunc
  adj::Tuple{Int, Int, Int}
  matched::BitVector
end

(P::Positioner)(t) = P.reorientfunc(t) .+ P.adj

function findpositioner(base, scanner, matches)
  len = length(scanner)
  best = nothing
  bestcnt = 0

  # if there are more than 12 total matches, we can't assume the first match is
  # actually valid
  for tryidx in 1:(length(matches)-11)
    for reorientfunc in reorientfuncs
      # We'll try reorienting the first matched point, then calculate an
      # adjustment. We'll then try the same reorientation and adjustment on the
      # rest of the matches.
      (b2,b1) = matches[tryidx]
      adj = base[b1] .- reorientfunc(scanner[b2])
      matched = falses(len)
      matched[b2] = true
      for k in (tryidx+1):length(matches)
        (b2,b1) = matches[k]
        if isequal(base[b1], reorientfunc(scanner[b2]) .+ adj)
          matched[b2] = true
        end
      end

      cnt = count(matched)
      if cnt > bestcnt && cnt >= 12
        best = Positioner(reorientfunc, adj, matched)
        bestcnt = cnt
      end
    end
  end
  return best
end

# utility to sort a list of numbers by the frequency - also reduces the set to
# unique values
function findmostcommon(v)
  cnts = Dict{eltype(v), Int}()
  common = nothing
  commoncnt = 0
  for e in v
    cnts[e] = get(cnts, e, 0) + 1
    if cnts[e] > commoncnt
      common = e
      commoncnt = cnts[e]
    end
  end
  return common
end

# we need to find the position of all the scanners
scannerpositions = Vector{Tuple{Int,Int,Int}}(undef, length(scanners))
scannerpositions[1] = (0,0,0)

# loop over scanners until we've figured them all out - if the loop ever makes
# no progress, it tries to change the "base" and find other matches that way
finished = falses(length(scanners))
finished[1] = true
baseidx = 1
loopfinishedcnt = 1
refdist = RefDistances(baseidx)
while !all(finished)
  base = scanners[baseidx]
  progress = false
  for i in (baseidx+1):length(scanners)
    finished[i] && continue

    # find matching distances
    scanner = scanners[i]
    matches = Dict{Int, Vector{Int}}()
    for (j,k) in Permutations(length(scanner))
      dist = distance(i, j, k)
      distidx = findfirst(refdist, dist)
      isnothing(distidx) && continue
      (pot1, pot2) = refdist[distidx]
      if haskey(matches, j)
        push!(matches[j], pot1, pot2)
      else
        matches[j] = [pot1, pot2]
      end
      if haskey(matches, k)
        push!(matches[k], pot1, pot2)
      else
        matches[k] = [pot1, pot2]
      end
    end

    # do we have at least 12 matches?
    length(matches) < 12 && continue

    # figure out how this scanner is oriented relative to the first
    positioner = findpositioner(base, scanner, [(k,findmostcommon(v)) for (k,v) in matches])
    if isnothing(positioner)
      println("Had 12+ matches, but couldn't find a positioner for scanner ", i)
      continue
    end

    # add new points
    append!(base, map(positioner, scanner[.!positioner.matched]))
    println("Number of beacons after ", count(finished), " scanners: ", length(base))
    scannerpositions[i] = positioner((0,0,0))

    # rebuild distances and mark this scanner finished
    global refdist = RefDistances(baseidx)
    finished[i] = true
    progress = true
  end

  if !progress
    global baseidx = findnext(i -> !i, finished, baseidx+1)
    if isnothing(baseidx)
      global baseidx = 1
      cnt = count(finished)
      loopfinishedcnt == cnt && error("stuck")
      global loopfinishedcnt = cnt
    end
    println("Switching base to ", baseidx)
    global refdist = RefDistances(baseidx)
  end
end

# find largest Manhattan distance between all scanners
largestdistance = 0
for (i,j) in Permutations(length(scannerpositions))
  dist = sum(abs.(scannerpositions[i] .- scannerpositions[j]))
  if dist > largestdistance
    global largestdistance = dist
  end
end

println("Result: ", largestdistance)


text:
include("../src/SplinePlotting.jl")

function nrbtestcrv()
    pnts = [0.5 1.5 4.5 3.0 7.5 6.0 8.5;
        3.0 5.5 5.5 1.5 1.5 4.0 4.5;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0]
    crv = nrbmak(pnts, [[0, 0, 0, 1/4, 1/2, 3/4, 3/4, 1, 1, 1]])
    return crv
end

function nrbtestcrv3d()
    coefs =[ 6.0  0.0  6.0  1;
        -5.5  0.5  5.5  1;
        -5.0  1.0 -5.0  1;
         4.5  1.5 -4.5  1;
         4.0  2.0  4.0  1;
        -3.5  2.5  3.5  1;
        -3.0  3.0 -3.0  1;
         2.5  3.5 -2.5  1;
         2.0  4.0  2.0  1;
        -1.5  4.5  1.5  1;
        -1.0  5.0 -1.0  1;
         0.5  5.5 -0.5  1;
         0.0  6.0  0.0  1]'
    knots = [0, 0, 0, 0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1, 1, 1, 1]
    crv3d = nrbmak(coefs, [knots])
    return crv3d
end

function nrbtestsrf()
    pnts = zeros(3, 5, 5)

    pnts[:,:,1] = [ 0.0  3.0  5.0  8.0 10.0;     # w*x
                0.0  0.0  0.0  0.0  0.0;         # w*y
                2.0  2.0  7.0  7.0  8.0]         # w*z

    pnts[:,:,2] = [ 0.0  3.0  5.0  8.0 10.0;
                    3.0  3.0  3.0  3.0  3.0;
                    0.0  0.0  5.0  5.0  7.0]

    pnts[:,:,3] = [ 0.0  3.0  5.0  8.0 10.0;
                    5.0  5.0  5.0  5.0  5.0;
                    0.0  0.0  5.0  5.0  7.0]

    pnts[:,:,4] = [ 0.0  3.0  5.0  8.0 10.0;
                    8.0  8.0  8.0  8.0  8.0;
                    5.0  5.0  8.0  8.0 10.0]

    pnts[:,:,5] = [ 0.0  3.0  5.0  8.0 10.0;
                   10.0 10.0 10.0 10.0 10.0;
                    5.0  5.0  8.0  8.0 10.0]

    uKnots = [0, 0, 0, 1/3, 2/3, 1, 1, 1]
    vKnots = [0, 0, 0, 1/3, 2/3, 1, 1, 1]

    srf = nrbmak(pnts, [uKnots, vKnots])
    return srf
end

function nrbtestvol()
    knots = [[0, 0, 0, 1/2, 1, 1, 1], [0, 0, 0, 1, 1, 1],
             [0, 0, 0, 1/6, 2/6, 1/2, 1/2, 4/6, 5/6, 1, 1, 1]];

    coefs = [-1.0000   -0.9734   -0.7071    1.4290    1.0000    3.4172
             0    2.4172         0    0.0148   -2.0000   -1.9734
             0    2.0000    4.9623    9.4508    4.0000    2.0000
        1.0000    1.0000    0.7071    1.0000    1.0000    1.0000
       -0.8536         0   -0.6036    1.9571    1.2071    3.5000
        0.3536    2.5000    0.2500    0.5429   -1.7071   -1.0000
             0    2.0000    4.4900    8.5444    3.4142    2.0000
        0.8536    1.0000    0.6036    1.0000    0.8536    1.0000
       -0.3536   -4.0000   -0.2500   -1.2929    1.7071    1.0000
        0.8536         0    0.6036   -2.7071   -1.2071   -5.0000
             0    2.0000    4.4900   10.0711    3.4142    2.0000
        0.8536    1.0000    0.6036    1.0000    0.8536    1.0000
             0   -4.0000         0    0.7071    2.0000    5.0000
        1.0000    4.0000    0.7071   -0.7071   -1.0000   -5.0000
             0    2.0000    4.9623   14.4142    4.0000    2.0000
        1.0000    1.0000    0.7071    1.0000    1.0000    1.0000
       -2.5000   -4.0000   -1.7678    0.7071    1.0000    5.0000
             0    4.0000         0   -0.7071   -3.5000   -5.0000
             0    2.0000    6.0418   14.4142    4.0000    2.0000
        1.0000    1.0000    0.7071    1.0000    1.0000    1.0000
       -2.4379         0   -1.7238    2.7071    1.9527    5.0000
        0.9527    4.0000    0.6737    1.2929   -3.4379   -1.0000
             0    2.0000    6.6827   10.0711    4.0000    2.0000
        1.0000    1.0000    0.7071    1.0000    1.0000    1.0000
       -0.9734   -1.0000   -0.6883    0.7071    3.4172    1.0000
        2.4172         0    1.7092   -1.4142   -1.9734   -2.0000
             0    4.0000    6.6827    4.9623    4.0000         0
        1.0000    1.0000    0.7071    0.7071    1.0000    1.0000
             0   -0.8536         0    0.8536    3.5000    1.2071
        2.5000    0.3536    1.7678   -1.2071   -1.0000   -1.7071
             0    3.4142    6.0418    4.4900    4.0000         0
        1.0000    0.8536    0.7071    0.6036    1.0000    0.8536
       -4.0000   -0.3536   -2.8284    1.2071    1.0000    1.7071
             0    0.8536         0   -0.8536   -5.0000   -1.2071
             0    3.4142    7.1213    4.4900    4.0000         0
        1.0000    0.8536    0.7071    0.6036    1.0000    0.8536
       -4.0000         0   -2.8284    1.4142    5.0000    2.0000
        4.0000    1.0000    2.8284   -0.7071   -5.0000   -1.0000
             0    4.0000   10.1924    4.9623    4.0000         0
        1.0000    1.0000    0.7071    0.7071    1.0000    1.0000
       -4.0000   -2.5000   -2.8284    0.7071    5.0000    1.0000
        4.0000         0    2.8284   -2.4749   -5.0000   -3.5000
             0    4.0000   10.1924    6.0418    4.0000         0
        1.0000    1.0000    0.7071    0.7071    1.0000    1.0000
             0   -2.4379         0    1.3808    5.0000    1.9527
        4.0000    0.9527    2.8284   -2.4309   -1.0000   -3.4379
             0    4.0000    7.1213    6.6827    4.0000         0
        1.0000    1.0000    0.7071    0.7071    1.0000    1.0000
       -1.0000   -0.9734    0.2071    2.4163    1.0000    3.4172
             0    2.4172   -1.2071   -1.3954   -2.0000   -1.9734
        2.0000    4.0000    7.0178    6.6827    2.0000         0
        1.0000    1.0000    1.0000    0.7071    1.0000    1.0000
       -0.8536         0    0.3536    2.4749    1.2071    3.5000
        0.3536    2.5000   -0.8536   -0.7071   -1.7071   -1.0000
        1.7071    4.0000    6.3498    6.0418    1.7071         0
        0.8536    1.0000    0.8536    0.7071    0.8536    1.0000
       -0.3536   -4.0000    0.8536    0.7071    1.7071    1.0000
        0.8536         0   -0.3536   -3.5355   -1.2071   -5.0000
        1.7071    4.0000    6.3498    7.1213    1.7071         0
        0.8536    1.0000    0.8536    0.7071    0.8536    1.0000
             0   -4.0000    1.2071    3.5355    2.0000    5.0000
        1.0000    4.0000   -0.2071   -3.5355   -1.0000   -5.0000
        2.0000    4.0000    7.0178   10.1924    2.0000         0
        1.0000    1.0000    1.0000    0.7071    1.0000    1.0000
       -2.5000   -4.0000   -0.5429    3.5355    1.0000    5.0000
             0    4.0000   -1.9571   -3.5355   -3.5000   -5.0000
        2.0000    4.0000    8.5444   10.1924    2.0000         0
        1.0000    1.0000    1.0000    0.7071    1.0000    1.0000
       -2.4379         0   -0.0355    3.5355    1.9527    5.0000
        0.9527    4.0000   -1.4497   -0.7071   -3.4379   -1.0000
        2.0000    4.0000    9.4508    7.1213    2.0000         0
        1.0000    1.0000    1.0000    0.7071    1.0000    1.0000]

    coefs = reshape(coefs, 4, 4, 3, 9)
    horseshoe = nrbmak(coefs, knots)

    return horseshoe
end


text:
module LagrangianFESpacesTests

using Test
using Gridap
using GridapGmsh

mshfile = joinpath(@__DIR__,"../demo/","demo.msh")

T = VectorValue{3,Float64}

dirinames = ["boundary1", "boundary2"]
dirimasks = [(true,false,true), (true,true,false)]

ufun(x) = x

a(v,u) = inner(e,e)

# CLagrangianFESpace

fespace = GmshCLagrangianFESpace(T, mshfile, dirinames, dirimasks)

grid = fespace.grid

trian = Triangulation(grid)
quad = CellQuadrature(trian,degree=2)

uh = interpolate(fespace,ufun)

u = CellField(trian,ufun)

e = u - uh

tol = 1.e-8
@test sum( integrate( a(e,e), trian, quad) ) < tol

# DLagrangianFESpace

fespace = GmshDLagrangianFESpace(T, mshfile, dirinames, dirimasks)

grid = fespace.grid

trian = Triangulation(grid)
quad = CellQuadrature(trian,degree=2)

uh = interpolate(fespace,ufun)

u = CellField(trian,ufun)

e = u - uh

tol = 1.e-8
@test sum( integrate( a(e,e), trian, quad) ) < tol

end # module


text:
"""
    Wagner(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)

Two-dimensional aerodynamic model based on Wagner's function with state variables
``\\lambda_1, \\lambda_2``, inputs ``u, v, \\omega``, and parameters ``a, b, a_0,
\\alpha_0, c_{d0}, c_{m0}``
"""
struct Wagner{TF}
    C1::TF
    C2::TF
    eps1::TF
    eps2::TF
end

"""
    Wagner(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)

Initialize a model of type [`Wagner`](@ref)
"""
Wagner(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3) = Wagner(promote(C1, C2, eps1, eps2)...)

# --- Submodel Creation --- #

function Submodel(model::Wagner)

    C1 = model.C1
    C2 = model.C2
    eps1 = model.eps1
    eps2 = model.eps2

    # residual function
    fresid = (dx, x, y, p, t) -> wagner_residual(dx, x, y, p, t; C1, C2, eps1, eps2)

    # number of state, input, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(2)
    ny = Val(3)
    np = Val(6)

    # jacobian definitions
    ratejac = Identity()
    statejac = Linear((dx, x, y, p, t) -> wagner_state_jacobian(dx, x, y, p, t; model.eps1, model.eps2))
    inputjac = Nonlinear((dx, x, y, p, t) -> wagner_input_jacobian(dx, x, y, p, t; model.C1, model.C2, model.eps1, model.eps2))
    paramjac = Nonlinear((dx, x, y, p, t) -> wagner_parameter_jacobian(dx, x, y, p, t; model.C1, model.C2, model.eps1, model.eps2))
    tgrad = Zeros()

    # convenience functions for setting states, inputs, and parameters
    setstate = wagner_set_states!
    setinput = wagner_set_inputs!
    setparam = wagner_set_parameters!

    # convenience functions for separating states, inputs, and parameters
    sepstate = wagner_separate_states
    sepinput = wagner_separate_inputs
    sepparam = wagner_separate_parameters

    # model definition
    return Submodel{false}(fresid, nx, ny, np;
        ratejac = ratejac,
        statejac = statejac,
        inputjac = inputjac,
        paramjac = paramjac,
        tgrad = tgrad,
        setstate = setstate,
        setinput = setinput,
        setparam = setparam,
        sepstate = sepstate,
        sepinput = sepinput,
        sepparam = sepparam,
    )
end

# --- Internal Methods --- #

# residual function
function wagner_residual(dx, x, y, p, t; C1, C2, eps1, eps2)
    # extract state rates
    dλ1, dλ2 = dx
    # extract states
    λ1, λ2 = x
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # return residual
    return SVector(
        dλ1 + eps1*u/b*λ1 - C1*eps1*u/b*(v + (b/2-a*b)*ω - u*α0),
        dλ2 + eps2*u/b*λ2 - C2*eps2*u/b*(v + (b/2-a*b)*ω - u*α0)
    )
end

# state jacobian function
function wagner_state_jacobian(dx, x, y, p, t; eps1, eps2)
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # return jacobian of residual with respect to the states
    return @SMatrix [eps1*u/b 0; 0 eps2*u/b]
end

# input jacobian function
function wagner_input_jacobian(dx, x, y, p, t; C1, C2, eps1, eps2)

    λ1, λ2 = x
    u, v, ω = y
    a, b, a0, α0, cd0, cm0 = p

    tmp1 = v/b + (1/2-a)*ω - 2*α0*u/b
    λ1dot_u = eps1/b*λ1 - C1*eps1*tmp1
    λ2dot_u = eps2/b*λ2 - C2*eps2*tmp1

    tmp2 = u/b
    λ1dot_v = -C1*eps1*tmp2
    λ2dot_v = -C2*eps2*tmp2

    tmp3 = u*(1/2-a)
    λ1dot_ω = -C1*eps1*tmp3
    λ2dot_ω = -C2*eps2*tmp3

    return @SMatrix [λ1dot_u λ1dot_v λ1dot_ω; λ2dot_u λ2dot_v λ2dot_ω]
end

# parameter jacobian function
function wagner_parameter_jacobian(dx, x, y, p, t; C1, C2, eps1, eps2)

    λ1, λ2 = x
    u, v, ω = y
    a, b, a0, α0, cd0, cm0 = p

    λ1dot_a = C1*eps1*u*ω
    λ1dot_b = -eps1*u/b^2*λ1 + C1*eps1/b^2*(u*v - u^2*α0)
    λ1dot_a0 = 0
    λ1dot_α0 = C1*eps1*u^2/b
    λ1dot_cd0 = 0
    λ1dot_cm0 = 0

    λ2dot_a = C2*eps2*u*ω
    λ2dot_b = -eps2*u/b^2*λ2 + C2*eps2/b^2*(u*v - u^2*α0)
    λ2dot_a0 = 0
    λ2dot_α0 = C2*eps2*u^2/b
    λ2dot_cd0 = 0
    λ2dot_cm0 = 0

    return @SMatrix [
        λ1dot_a λ1dot_b λ1dot_a0 λ1dot_α0 λ1dot_cd0 λ1dot_cm0;
        λ2dot_a λ2dot_b λ2dot_a0 λ2dot_α0 λ2dot_cd0 λ2dot_cm0;
        ]
end

# convenience function for defining this model's state vector
wagner_set_states!(x; lambda) = x .= lambda

# convenience function for defining this model's input vector
function wagner_set_inputs!(y; u, v, omega)
    y[1] = u
    y[2] = v
    y[3] = omega
    return y
end

# convenience function for defining this model's parameter vector
function wagner_set_parameters!(p; a, b, a0, alpha0, cd0, cm0)
    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0
    p[5] = cd0
    p[6] = cm0
    return p
end

# convenience function for separating this model's state vector
wagner_separate_states(x) = (lambda = x,)

# convenience function for separating this model's input vector
wagner_separate_inputs(y) = (u=y[1], v=y[2], omega=y[3])

# convenience function for separating this model's parameter vector
wagner_separate_parameters(p) = (a=p[1], b=p[2], a0=p[3], alpha0=p[4], cd0=p[5], cm0=p[6])

# aerodynamic loads per unit span
function wagner_loads(a, b, ρ, c, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # Wagner's function at t = 0.0
    ϕ0 = 1 - C1 - C2
    # Velocity Magnitude (squared)
    V2 = u^2 + v^2
    # Mach Number (squared)
    M2 = V2/c^2
    # Prandtl-Glauert correction factor
    beta = sqrt(1 - min(0.99, M2))
    # normal force at reference point
    N = tmp1*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    # axial force at reference point
    A = -a0*ρ*b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)^2
    # moment at reference point
    M = -tmp2*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N
    # apply compressibility correction
    N = N / beta
    A = A / beta
    M = M / beta
    # add skin friction drag
    A += ρ*b*u^2*cd0
    # return loads
    return SVector(N, A, M)
end


text:
module ImageMagick

using FileIO: DataFormat, @format_str, Stream, File, filename, stream
using InteractiveUtils: subtypes
using ImageCore
using Libdl
using ImageMagick_jll
using Base: convert

Color1{T}           = Color{T,1}
Color2{T,C<:Color1} = TransparentColor{C,T,2}
Color3{T}           = Color{T,3}
Color4{T,C<:Color3} = TransparentColor{C,T,4}

AbstractGray{T}     = Color{T,1}

export readblob, image2wand, magickinfo

include("libmagickwand.jl")

# Image / Video formats

const image_formats = [
    format"BMP",
    format"AVI",
    format"CRW",
    format"CUR",
    format"DCX",
    format"DOT",
    format"EPS",
    format"GIF",
    format"HDR",
    format"ICO",
    format"INFO",
    format"JP2",
    format"JPEG",
    format"PCX",
    format"PDB",
    format"PDF",
    format"PGM",
    format"PNG",
    format"PSD",
    format"RGB",
    format"TIFF",
    format"WMF",
    format"WPG",
    format"TGA"
]

"""
    metadata(file) -> tuple of dimensions

The returned value is (X,Y[,Z]), not (R,C), so the first two dimensions
of the corresponding Array returned by `load` are reversed.
"""
function metadata(filename::AbstractString)
    wand = MagickWand()
    readimage(wand, filename)
    sz, T, _, _ = _metadata(wand)
    return sz, T
end

metadata(imagefile::File{T}) where {T <: DataFormat} = metadata(filename(imagefile))

function _metadata(wand)
    # Determine what we need to know about the image format
    sz = size(wand)
    n = getnumberimages(wand)
    if n > 1
        sz = tuple(sz..., n)
    end

    imtype      = getimagetype(wand)
    havealpha   = getimagealphachannel(wand)
    cs          = getimagecolorspace(wand)
    if imtype == "GrayscaleType" || imtype == "GrayscaleMatteType"
        cs = "Gray"
    end

    depth = getimagedepth(wand)
    # use an even # of fractional bits for depth>8 (see issue 242#issuecomment-68845157)
    evendepth = ((depth+1)>>1)<<1
    if depth <= 8
        if cs == "Gray"
            cdepth = getimagechanneldepth(wand, GrayChannel)
            if n > 1
                for k = 1:n  # while it might seem that this should be 2:n, that doesn't work...
                    nextimage(wand)
                    cdepth = max(cdepth, getimagechanneldepth(wand, GrayChannel))
                end
                resetiterator(wand)
            end
            T = cdepth == 1 ? Bool : N0f8
        else
            T = N0f8     # otherwise use 8 fractional bits
        end
    elseif depth <= 16
        T = Normed{UInt16,evendepth}
    else
        @warn "some versions of ImageMagick give spurious low-order bits for 32-bit TIFFs" maxlog=1
        T = Normed{UInt32,evendepth}
    end

    channelorder = cs
    if havealpha
        if channelorder == "sRGB" || channelorder == "RGB"
            T, channelorder = RGBA{T}, "RGBA"
        elseif channelorder == "Gray"
            T, channelorder = GrayA{T}, "IA"
        else
            error("Cannot parse colorspace $channelorder")
        end
    else
        if channelorder == "sRGB" || channelorder == "RGB"
            T, channelorder = RGB{T}, "RGB"
        elseif channelorder == "Gray"
            T, channelorder = Gray{T}, "I"
        else
            error("Cannot parse colorspace $channelorder")
        end
    end
    sz, T, cs, channelorder
end

load(@nospecialize(imagefile::File{<:DataFormat}), @nospecialize(args...); key_args...) = load_(filename(imagefile), args...; key_args...)
load(filename::AbstractString, @nospecialize(args...); key_args...) = load_(filename, args...; key_args...)
save(@nospecialize(imagefile::File{<:DataFormat}), @nospecialize(args...); key_args...) = save_(filename(imagefile), args...; key_args...)
save(filename::AbstractString, @nospecialize(args...); key_args...) = save_(filename, args...; key_args...)

load(@nospecialize(imgstream::Stream{<:DataFormat}), @nospecialize(args...); key_args...) = load_(stream(imgstream), args...; key_args...)
load(imgstream::IO, args...; key_args...) = load_(imgstream, args...; key_args...)
save(@nospecialize(imgstream::Stream{<:DataFormat}), args...; key_args...) = save_(imgstream, args...; key_args...)

const ufixedtype = Dict(10=>N6f10, 12=>N4f12, 14=>N2f14, 16=>N0f16)

readblob(data::Vector{UInt8}) = load_(data)

function load_(file::Union{AbstractString,IO,Vector{UInt8}}, permute_horizontal=true; ImageType=Array, extraprop="", extrapropertynames=nothing, view=false)
    if ImageType != Array
        error("this function now returns an Array, do not use ImageType keyword.")
    end
    if extraprop != "" || extrapropertynames != nothing
        error("keywords \"extraprop\" and \"extrapropertynames\" no longer work, use magickinfo instead")
    end
    wand = MagickWand()
    readimage(wand, file)
    resetiterator(wand)

    sz, T, cs, channelorder = _metadata(wand)

    # Allocate the buffer and get the pixel data
    buf = Array{T}(undef, sz)
    exportimagepixels!(rawview(channelview(buf)), wand, cs, channelorder)

    orient = getimageproperty(wand, "exif:Orientation", false)
    default = (A,ph) -> ph ? vertical_major(A) : identity(A)
    oriented_buf = get(orientation_dict, orient, default)(buf, permute_horizontal)
    view ? oriented_buf : collect(oriented_buf)
end

function save_(filename::AbstractString, @nospecialize(img), permute_horizontal=true; mapi = identity, quality = nothing, kwargs...)
    wand = image2wand(img, mapi, quality, permute_horizontal; kwargs...)
    writeimage(wand, filename)
end

# This differs from `save_` for files because this is primarily used
# by IJulia, and we want to restrict large images to make display faster.
function save_(s::Stream, img, permute_horizontal=true; mapi = clamp01nan, quality = nothing)
    wand = image2wand(img, mapi, quality, permute_horizontal)
    blob = getblob(wand, formatstring(s))
    write(stream(s), blob)
end

function image2wand(@nospecialize(img), mapi=identity, quality=nothing, permute_horizontal=true; kwargs...)
    local imgw
    try
        imgw = map(x->mapIM(mapi(x)), img)
    catch
        @warn "Mapping to the storage type failed; perhaps your data had out-of-range values?\nTry `map(clamp01nan, img)` to clamp values to a valid range."
        rethrow()
    end
    permute_horizontal && (imgw = collect(vertical_major(imgw)))
    if ndims(imgw) > 3
        error("At most 3 dimensions are supported")
    end
    wand = MagickWand()
    T = eltype(imgw)
    channelorder = T<:Real ? "Gray" : ColorTypes.colorant_string(T)
    if T <: Union{RGB,RGBA,ARGB,BGRA,ABGR}
        # For imagemagick versions <= v6.7.5, this should be "RGB",
        # but since we ship our own version we always use "sRGB"
        cs = "sRGB"
    else
        cs = channelorder
    end
    if channelorder == "Gray"
        channelorder = "I"
    elseif channelorder == "GrayA"
        channelorder = "IA"
    elseif channelorder == "AGray"
        channelorder = "AI"
    end
    tmp = to_explicit(to_contiguous(imgw))
    constituteimage(tmp, wand, cs, channelorder)
    if quality != nothing
        setimagecompressionquality(wand, quality)
    end
    resetiterator(wand)
    setproperties!(wand; kwargs...)
    wand
end

formatstring(s::Stream{DataFormat{S}}) where {S} = string(S)

@doc raw"""
    magickinfo(file::Union{AbstractString,IO})

Reads an image `file` and returns an array of strings describing the property data in the image file.

# Examples
```julia-repl
julia> p = magickinfo("Image.jpeg")
63-element Array{String,1}:
 "date:create"
 "date:modify"
 "exif:ApertureValue"
 "exif:BrightnessValue"
 "exif:ColorSpace"
 "exif:ComponentsConfiguration"
 "exif:DateTime"
 "exif:DateTimeDigitized"
 "exif:DateTimeOriginal"
 "exif:ExifImageLength"
 "exif:ExifImageWidth"
 "exif:ExifOffset"
 "exif:ExifVersion"
 ⋮
 "exif:thumbnail:JPEGInterchangeFormat"
 "exif:thumbnail:JPEGInterchangeFormatLength"
 "exif:thumbnail:ResolutionUnit"
 "exif:thumbnail:XResolution"
 "exif:thumbnail:YResolution"
 "exif:WhiteBalance"
 "exif:XResolution"
 "exif:YCbCrPositioning"
 "exif:YResolution"
 "jpeg:colorspace"
 "jpeg:sampling-factor"
 "unknown"a
```
"""
function magickinfo(file::Union{AbstractString,IO})
    wand = MagickWand()
    readimage(wand, file)
    resetiterator(wand)
    getimageproperties(wand, "*")
end

@doc raw"""
    magickinfo(file::Union{AbstractString,IO}, properties::Union{Tuple,AbstractVector})

Reads an image `file` and returns a Dict containing the values for the requested `properties`.

# Examples
```julia-repl
julia> magickinfo("IMG_6477.jpeg",
           ("exif:DateTime","exif:GPSLatitude","exif:GPSLatitudeRef",
            "exif:GPSLongitude","exif:GPSLongitudeRef","exif:GPSAltitude"))
Dict{String,Any} with 6 entries:
  "exif:GPSLongitudeRef" => "W"
  "exif:GPSLatitude"     => "44/1, 22/1, 2326/100"
  "exif:GPSLatitudeRef"  => "N"
  "exif:GPSLongitude"    => "71/1, 13/1, 5301/100"
  "exif:DateTime"        => "2020:01:22 13:17:41"
  "exif:GPSAltitude"     => "261189/757"
```
"""
function magickinfo(file::Union{AbstractString,IO}, properties::Union{Tuple,AbstractVector})
    wand = MagickWand()
    readimage(wand, file)
    resetiterator(wand)

    props = Dict{String,Any}()
    for p in properties
        props[p] = getimageproperty(wand, p)
    end
    props
end
magickinfo(file::Union{AbstractString,IO}, properties...) = magickinfo(file, properties)

function setproperties!(wand; fps=nothing)
    if fps != nothing
        setimagedelay(wand, fps)
    end
    wand
end

# ImageMagick element-mapping function. Converts to RGB/RGBA and uses
# N0f8 "inner" element type.
mapIM(c::Color1) = mapIM(convert(Gray, c))
mapIM(c::Gray{T}) where {T} = convert(Gray{N0f8}, c)
mapIM(c::Gray{T}) where {T<:Normed} = c

mapIM(c::Color2) = mapIM(convert(GrayA, c))
mapIM(c::GrayA{T}) where {T} = convert(GrayA{N0f8}, c)
mapIM(c::GrayA{T}) where {T<:Normed} = c

mapIM(c::Color3) = mapIM(convert(RGB, c))
mapIM(c::RGB{T}) where {T} = convert(RGB{N0f8}, c)
mapIM(c::RGB{T}) where {T<:Normed} = c

mapIM(c::Color4) = mapIM(convert(RGBA, c))
mapIM(c::RGBA{T}) where {T} = convert(RGBA{N0f8}, c)
mapIM(c::RGBA{T}) where {T<:Normed} = c

mapIM(x::UInt8) = reinterpret(N0f8, x)
mapIM(x::UInt16) = reinterpret(N0f16, x)
mapIM(x::UInt32) = reinterpret(N0f32, x)
mapIM(x::Bool) = x
mapIM(x::AbstractFloat) = convert(N0f8, x)
mapIM(x::Normed) = x

# Make the data contiguous in memory, this is necessary for
# imagemagick since it doesn't handle stride.
to_contiguous(A::Array) = A
to_contiguous(A::AbstractArray) = collect(A)
to_contiguous(A::BitArray) = convert(Array{Bool}, A)
if isdefined(ImageCore, :ColorView)
    to_contiguous(A::ColorView) = to_contiguous(channelview(A))
end

to_explicit(A::Array{C}) where {C<:Colorant} = to_explicit(channelview(A))
function to_explicit(A::AbstractArray)
    As = similar(A)
    As .= A
    to_explicit(As)
end
to_explicit(A::Array{Bool}) = A
to_explicit(A::Array{T}) where {T<:Normed} = rawview(A)
to_explicit(A::Array{T}) where {T<:AbstractFloat} = to_explicit(convert(Array{N0f8}, A))

include("precompile.jl")
_precompile_()

end # module


text:
function append_choice!(r::Real, x1::Vector, x1_idx::Int, x2::Vector, gp_state::Dict)
    # This part is pretty ugly, but can't help it
    data_x = gp_state[:data_x]
    data_c = gp_state[:data_c]

    choices  = nothing
    xc = r*x1 + (1-r)*x2
    xm = (x1 + x2) / 2

    start_idx = size(data_x, 2) + 1
    data_new, choices = begin
        if(abs(r - 1.0) < 0.02)
            # Corner solution: x1 is favored
            # [x1 x_mid x2]
            x_new   = hcat(xm, x2)
            choices = [x1_idx, start_idx, start_idx+1]
            x_new, choices
        elseif(abs(r - 0.0) < 0.02)
            # Corner solution: x2 is favored
            # [x2 x_mid x1]
            x_new   = hcat(xm, x2)
            choices = [start_idx+1, start_idx, x1_idx]
            x_new, choices
        else
            # [x_c x1 x2]
            x_new   = hcat(x2, xc)
            choices = [start_idx+1, x1_idx, start_idx]
            x_new, choices
        end
    end
    data_x = hcat(data_x, data_new)
    data_c = vcat(data_c, reshape(choices, (1, 3)))
    gp_state[:data_x] = data_x
    gp_state[:data_c] = data_c
end

function append_choice!(r::Real, x1::Vector, x2::Vector, gp_state::Dict)
    # This part is pretty ugly, but can't help it
    data_x = gp_state[:data_x]
    data_c = gp_state[:data_c]

    choices  = nothing
    xc = r*x1 + (1-r)*x2
    xm = (x1 + x2) / 2

    start_idx = size(data_x, 2) + 1
    data_new, choices = begin
        if(abs(r - 1.0) < 0.02)
            # Corner solution: x1 is favored
            x_new   = hcat(xm, x2, x1)
            choices = [start_idx+2, start_idx+1, start_idx]
            x_new, choices
        elseif(abs(r - 0.0) < 0.02)
            # Corner solution: x2 is favored
            x_new   = hcat(xm, x1, x2)
            choices = [start_idx+2, start_idx+1, start_idx]
            x_new, choices
        else
            # [x_c x1 x2]
            x_new   = hcat(x1, x2, xc)
            choices = [start_idx+2, start_idx+1, start_idx]
            x_new, choices
        end
    end
    data_x = hcat(data_x, data_new)
    data_c = vcat(data_c, reshape(choices, (1, 3)))
    gp_state[:data_x] = data_x
    gp_state[:data_c] = data_c
end

function train_gp!(prng,
                   gp_state::Dict,
                   mcmc_samples::Int,
                   mcmc_burnin::Int,
                   mcmc_thin::Int,
                   scale::Real,
                   marginalize::Bool)

    dims   = size(gp_state[:data_x], 1)
    priors = Product(
        vcat([Normal(0, 2), Normal(-2, 2)],
             [Normal(log(0.5), log(2)) for i = 1:dims]))

    if(marginalize)
        θ_init = begin
            if(haskey(gp_state, :θ))
                mean(gp_state[:θ], dims=2)[:,1]
            else
                rand(prng, priors)
            end
        end

        latent_init = zeros(size(gp_state[:data_x], 2))
        θ, f, a, K = pm_ess(prng,
                            mcmc_samples,
                            mcmc_burnin,
                            mcmc_thin,
                            θ_init,
                            latent_init,
                            priors,
                            scale,
                            gp_state[:data_x],
                            gp_state[:data_c])
        k = precompute_lgp(θ)
        gp_state[:θ] = θ
        gp_state[:f] = f
        gp_state[:a] = a
        gp_state[:K] = K
        gp_state[:k] = k
    else
        θ_init = begin
            if(haskey(gp_state, :θ))
                θ_init
            else
                rand(prng, priors)
            end
        end

        θ, f, _, a, _, K = map_laplace(gp_state[:data_x],
                                       gp_state[:data_c],
                                       θ_init,
                                       scale,
                                       priors;
                                       verbose=true)
        k = construct_kernel(θ...)
        gp_state[:θ] = θ
        gp_state[:f] = f
        gp_state[:a] = a
        gp_state[:K] = K
        gp_state[:k] = k
    end
end

function prefbo_next_query(prng,
                           gp_state::Dict,
                           search_budget::Real;
                           verbose::Bool=false)
    dims = size(gp_state[:data_x],1)
    K    = gp_state[:K]
    a    = gp_state[:a]
    k    = gp_state[:k]
    X    = gp_state[:data_x]

    x_opt, y_opt, opt_idx = optimize_mean(dims, search_budget, X, K, a, k)
    x_sol, y_sol          = optimize_acquisition(
        dims, search_budget, x_opt, y_opt, X, K, a, k;
        verbose=verbose, prng=prng)

    x1 = x_opt
    x2 = x_sol
    x1, opt_idx, x2
end