text:
using JuMP,CPLEX,Statistics,MAT,Redis,Printf

#connect to redi
conn = RedisConnection()

function min_(X,model)
    if(length(X)!=2)
        error("supperted only two sided minimum")
    end
    d=@variable(model,binary = true)
    dmin=@variable(model)
    @constraints(model, begin
           dmin<=X[1]
           dmin<=X[2]
           dmin>=X[1]-100*(d)
           dmin>=X[2]-100*(1-d)
    end)
    return dmin
end

function rate(F,T,h)
    return [(1-F)*T[1,h],
        (1-F)*T[3,h],
        (F)*T[2,h],
        ];
end

function buildCtrl(jump,delta,dt,dt2,iH,oH,H,maxS,MU,X0,tgt,N,Nr)

    model = Model(CPLEX.Optimizer)
    set_optimizer_attribute(model, "CPX_PARAM_SCRIND", 0)
    set_optimizer_attribute(model, "CPX_PARAM_THREADS", 1)

    @variable(model,X[i=1:N,j=0:H]>=0)
    @variable(model,T[i=1:Nr,j=0:H]>=0)
    @variable(model,E_abs[i=1:H]>=0)
    #@variable(model,E_abs>=0)
    @variable(model,avgX>=0)
    @variable(model,S[i=1:2]>=0)

    @constraint(model,S[1]<=maxS[1])
    @constraint(model,S[1]>=1)
    @constraint(model,S[2]<=maxS[2])
    @constraint(model,X.<=sum(X0))

    @variable(model,D[i=1:2,j=0:H]>=0)
    @constraint(model,[h=0:H],D[1,h]<=X[2,h])
    @constraint(model,[h=0:H],D[1,h]<=S[1]-X[3,h])
    @constraint(model,[h=0:H],D[2,h]<=X[3,h])
    @constraint(model,[h=0:H],D[2,h]<=S[2])

    #INIT STATE constraint
    @constraint(model,X0c[i=1:N],X[i,0]==X0[i])
    #Constraints per lo sforzo di controllo
    @constraint(model,SC_u[i=1:2],S[i]<=maxS[i]);
    @constraint(model,SC_l[i=1:2],S[i]>=1);

    for h=0:H
        @constraint(model,T[1,h]==MU[1]*X[1,h])
        if(h<=H)
            @constraint(model,T[2,h]==delta*min_([X[2,h],S[1]-X[3,h]],model))
            @constraint(model,T[3,h]==MU[2]*min_([X[3,h],S[2]],model))
        else
            @constraint(model,T[2,h]==delta*D[1,h])
            @constraint(model,T[3,h]==MU[2]*D[2,h])
        end
    end

    local h=0
    for oIdx=1:(oH)
        for iIdx=1:iH
            dy=(jump'*rate(1,T,h))'
            @constraint(model,[i=1:N],X[i,h+1]==dy[i]*dt2+X[i,h])
            h=h+1;
        end
        #integrate the slow system for one step
        dy2=(jump'*rate(0,T,h))'
        @constraint(model,[i=1:N],X[i,h+1]==dy2[i]*dt+X[i,h])
        h=h+1;
    end

    #avg constraint
    @constraint(model,avgX==(1.0/H)*sum(X[1,:]))

    #add absolute value constraint
    for i=1:H
        @constraints(model, begin
               E_abs[i]>=(X[1,i]-tgt)
               E_abs[i]>=-(X[1,i]-tgt)
        end)
    end

    #@constraints(model, begin
    #    E_abs>=(avgX-tgt)
    #    E_abs>=-(avgX-tgt)
    #end)

    for i=0:H
        @constraints(model, begin
               X[3,i]<=S[1]
        end)
    end

    @objective(model,Min,sum(E_abs))
    return [model,SC_l,SC_u,X0c,S]
end

jump=[[-1,1,0],
      [1,0,-1],
      [0,-1,1]];
delta=10^2
dt=0.3*10^-2  #time step lento
dt2=0.5*dt #time step veloce
iH=1
oH=parse(Int64, ARGS[3])
maxS=[40,40];
H=oH*(1+iH)
MU=[parse(Float64, ARGS[1]),parse(Float64, ARGS[2])]
X0=[20,0,0]
alfa=parse(Float64, ARGS[4])
tgt=round((sum(X0)*((1/MU[1])/(1/MU[1]+1/MU[2]))),digits=4)
tgt=tgt*alfa
N=3
Nr=3

isStarted=false
objOpt=buildCtrl(jump,delta,dt,dt2,iH,oH,H,maxS,MU,X0,tgt,N,Nr)

model=objOpt[1]
SC_l=objOpt[2]
SC_u=objOpt[3]
X0c=objOpt[4]
S=objOpt[5]

X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
prevS=[0,0]
while(true)
    isAppRunning=get(conn, "app_running")
    if(isAppRunning=="1")
        if(isStarted==true)
            set_normalized_rhs(SC_l[1],max(1,prevS[1]-2))
            set_normalized_rhs(SC_l[2],max(1,prevS[2]-10))
            #set_normalized_rhs(SC_u[1],min(maxS[1],prevS[1]+10)
            #set_normalized_rhs(SC_u[2],min(maxS[2],prevS[2]+1))
        end

        if(isStarted==false)
            global isStarted=true
            set(conn, "ctrl","1")
        end
        for j=1:N
            set_normalized_rhs(X0c[j],X0[j])
        end

        JuMP.optimize!(model)

        if(termination_status(model)==MOI.OPTIMAL)
            local soltime=MOI.get(model, MOI.SolveTime())

            Sopt=Int.(round.(value.(S)))
            #Sopt=[Int(round(value(S[1]))),value(S[2])]

            println(Sopt,X0,soltime,tgt)

            global prevS=Sopt

            set(conn, "swthread",@sprintf("%d",Sopt[1]))
            set(conn, "hwthread",@sprintf("%.5f",Sopt[2]))
        else
            println(termination_status(model), X0)
        end

        global X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
    else
        println("App stopped")
        global X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
        global isStarted=false
        global prevS=[0,0]
        global objOpt=buildCtrl(jump,delta,dt,dt2,iH,oH,H,maxS,MU,X0,tgt,N,Nr)
        global model=objOpt[1]
        global SC_l=objOpt[2]
        global SC_u=objOpt[3]
        global X0c=objOpt[4]
        global S=objOpt[5]
    end
end

#Controlled layered three-tier system by varying #threads HW/SW at each tier. Simulation (i.e., 100 runs) Vs Target reguirement.
#println(mean(mean(XS[:,:,:],dims=3)[end-100:end,1])," ",abs(mean(mean(XS[:,:,:],dims=3)[end-100:end,1])-tgt)*100/tgt," ",mean(stimes)," ",maximum(stimes))
#println(mean(mean(Shw,dims=3),dims=1))
# plot(t,mean(XS[:,:,:],dims=3)[:,:],label = "Controlled-simulation",legend=false)
# hline!([tgt],label = "Target",reuse=true,legend=false)
# xlabel!("Time(s)")
# ylabel!("QueueLength")

#p=plot(t,mean(Shw,dims=3)[:,:],reuse = false,title="SW/HW Core")

# matwrite("LP_LQNMPC.mat", Dict(
# 	"XS" => XS,
# 	"NT" => Tsw,
#     "NC" => Shw,
#     "Time"=>stimes
# ); compress = false)

max_stars_repo_name: giulio-garbi/fitMS

max_issues_repo_name: giulio-garbi/fitMS

max_forks_repo_name: giulio-garbi/fitMS


text:
using Destruct
using BenchmarkTools

function unpack_bc(w::Array{<:Tuple})
    Tuple((v->v[i]).(w) for i=1:length(w[1]))
end

# some transform
f(a, b) = a+b, a*b, a-b
# some transform with different return types
g(a, b) = a+b+1im*(a-b), a*b, convert(Int, round(a-b))

shape = (10,10)
a_1 = rand(shape)
a_2 = rand(shape)
a_3 = rand(shape)
unpacked = (a_1, a_2, a_3)
packed = collect(zip(a_1, a_2, a_3))

@assert unpack_bc(packed) == unpacked
@assert destruct(packed) == unpacked

println("Array of NTuple")
println("---------------")
for sz=[10 100 200]
    packed = f.(rand(sz,sz,1), rand(1,1,sz))
    shape = size(packed)
    println("shape: $shape")
    println("broadcast")
    @btime unpack_bc($packed)
    println("cartesian")
    @btime destruct($packed)
    println("")
end
println("Array of Tuple")
println("--------------")
for sz=[10 100 200]
    packed = g.(rand(sz,sz,1), rand(1,1,sz))
    shape = size(packed)
    println("shape: $shape")
    println("broadcast")
    @btime unpack_bc($packed)
    println("cartesian")
    @btime destruct($packed)
    println("")
end

ext: jl


text:
"""
    block(i, j)

Return the linear index of the `[i,j]` position ("block") in the row-major packed lower triangle.

Use the row-major ordering in this case because the result depends only on `i`
and `j`, not on the overall size of the array.

When `i == j` the value is the same as `kp1choose2(i)`.
"""
function block(i::Integer, j::Integer)
    0 < j ≤ i || throw(ArgumentError("[i,j] = [$i,$j] must be in the lower triangle"))
    return kchoose2(i) + j
end

"""
    kchoose2(k)

The binomial coefficient `k` choose `2` which is the number of elements
in the packed form of the strict lower triangle of a matrix.
"""
function kchoose2(k)      # will be inlined
    return (k * (k - 1)) >> 1
end

"""
    kp1choose2(k)

The binomial coefficient `k+1` choose `2` which is the number of elements
in the packed form of the lower triangle of a matrix.
"""
function kp1choose2(k)
    return (k * (k + 1)) >> 1
end

"""
    ltriindprs

A row-major order `Vector{NTuple{2,Int}}` of indices in the strict lower triangle.
"""
const ltriindprs = NTuple{2,Int}[]

function checkindprsk(k::Integer)
    kc2 = kchoose2(k)
    if length(ltriindprs) < kc2
        sizehint!(empty!(ltriindprs), kc2)
        for i in 1:k, j in 1:(i - 1)
            push!(ltriindprs, (i, j))
        end
    end
    return ltriindprs
end


text:
using Cascadia, Gumbo, HTTP, AbstractTrees

pages = Dict(
  :en => ["Wikipedia", "Osama_bin_Laden_(elephant)", "List_of_lists_of_lists", "Josephine_Butler", "Canadian_football", "Judaism"],
  :it => ["Wikipedia", "Ludovico_Einaudi", "Filosofia_della_scienza", "Pizza", "Effie_Gray", "Galeazzo_Maria_Sforza", "Ebraismo"],
  :fr => ["Wikipedia", "Philosophie_des_sciences", "Seconde_Guerre_mondiale", "Eric_Hakonsson"],
  :es => ["Wikipedia", "Chorizo", "Historia_de_Barcelona", "Espana", "Las_Vegas_Strip", "Judaismo"],
  :da => ["Wikipedia", "H.C._Andersen", "L.A._Ring", "Jiangxi", "NATO", "Thomas_Edison", "Bangladesh"])

rawpage(url) = parsehtml(String(HTTP.get(url).body)).root

function innerText(dom)
  text = IOBuffer()
  for elem in PreOrderDFS(dom)
    elem isa HTMLText && print(text, elem.text)
  end
  return String(text)
end

content(url) = join(innerText.(matchall(sel".mw-parser-output > p", rawpage(url))), "\n")

cd(@__DIR__)
mkpath("corpus")

for (lang, ps) in pages
  open("corpus/$lang.txt", "w") do io
    for p in ps
      write(io, content("https://$lang.wikipedia.org/wiki/$p"))
    end
  end
end

max_stars_repo_name: shashikdm/model-zoo

max_issues_repo_name: shashikdm/model-zoo

max_forks_repo_name: shashikdm/model-zoo


text:
#=
    Used from https://github.com/lanl-ansi/Juniper.jl
=#

function get_default_options()
    nl_solver                           = nothing
    log_levels                          = [:Options,:Table,:Info]
    silent                              = false
    atol                                = 1e-6
    # Obj cuts
    incumbent_constr                    = false
    obj_epsilon                         = 0
    # :UserLimit
    time_limit                          = Inf
    mip_gap                             = 1e-4
    best_obj_stop                       = NaN
    solution_limit                      = 0
    all_solutions                       = false
    list_of_solutions                   = false

    mip_solver                          = nothing
    allow_almost_solved                 = true
    allow_almost_solved_integral        = true
    registered_functions                = nothing
    # Parallel
    processors                          = 1

    return SolverOptions(nl_solver,log_levels,silent,atol,incumbent_constr,obj_epsilon,time_limit,mip_gap,best_obj_stop,solution_limit,all_solutions,
        list_of_solutions, mip_solver, allow_almost_solved, allow_almost_solved_integral, registered_functions, processors)
end

function combine_options(options)
    options_dict = Dict{Symbol,Any}()
    for kv in options
        if !in(kv[1], fieldnames(SolverOptions))

        else
            options_dict[kv[1]] = kv[2]
        end
    end
    if haskey(options_dict, :log_levels)
        if length(options_dict[:log_levels]) == 0
            options_dict[:log_levels] = Symbol[]
        end
    end

    defaults = get_default_options()

    return defaults
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
using FactorioUtils
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
test_blueprint_string = "0eNqV0V0KgzAMAOC75LmCfzjpVcYYVcMIaJQ2jon07msVxkAG7qkkTb6UdIWmn3GyxAJ6BWpHdqCvKzh6sOljTpYJQQMJDqCAzRAjsYbdNFpJGuwFvALiDl+gM39TgCwkhLu0Bcud56FBGwo+BrFDKyGnYBpdaBg5zgtIEoqWeAS3I4vtfpd7deDyX086oOfN4rSZnjbLP83Ux0VuS9dff6SgN6EtMuhi+ETrtt6qvuRZXVRpUXr/BsQ6nMQ="
test_json_string = """{"blueprint":{"icons":[{"signal":{"type":"item","name":"transport-belt"},"index":1}],"entities":[{"entity_number":1,"name":"inserter","position":{"x":-1,"y":-1},"direction":2},{"entity_number":2,"name":"transport-belt","position":{"x":1,"y":-1},"direction":2},{"entity_number":3,"name":"transport-belt","position":{"x":0,"y":-1},"direction":2},{"entity_number":4,"name":"transport-belt","position":{"x":0,"y":0}}],"item":"blueprint","label":"test","version":68721836034}}}"""

@test decode_blueprint(test_blueprint_string)== test_json_string
@test decode_blueprint(IOBuffer(test_blueprint_string)) == test_json_string
@test decode_blueprint(convert(Array{UInt8, 1}, test_blueprint_string)) == test_json_string


text:
using ParamLevelSet
using jInv.Mesh
using LinearAlgebra

println("Testing the rotation itself (extended RBFs):")

n = [64,64,64];
Mesh = getRegularMesh([0.0;5.0;0.0;5.0;0.0;5.0],n);
midmesh = [(Mesh.domain[2]+Mesh.domain[1])/2 ; (Mesh.domain[4]+Mesh.domain[3])/2; (Mesh.domain[5]+Mesh.domain[6])/2];
lenmesh = [(Mesh.domain[2]-Mesh.domain[1]) ; (Mesh.domain[4]-Mesh.domain[3]); (Mesh.domain[6]-Mesh.domain[5])];
numParamOfRBF = 10;
nRBF = 5
m = zeros(nRBF*numParamOfRBF);
for k=1:nRBF
	offset = (k-1)*numParamOfRBF + 1;
	m[offset] = randn();
	#A = randn(3,3); A = 0.3*A'*A + 0.5*eye(3);
	A = diagm(rand(3).+1.0);
	# L = chol(A)';
	# m[(offset+1):(offset+6)] = L[L.!=0.0];
	m[(offset+1):(offset+6)] = A[idx_trilA()];
	m[(offset+7):(offset+9)] = midmesh .+ 0.05*lenmesh.*randn(size(lenmesh));
end
u, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0,numParamOfRBF=numParamOfRBF);

theta_phi = deg2rad.([37.0 24.0 ; 27.0 14.0]);
b = [1.0 5.0 7.0 ; -2.3 5.2 3.5]*sum((Mesh.h))/3.0;

Xc = convert(Array{Float32,2},getCellCenteredGrid(Mesh));

N = prod(size(u));
ur1 = zeros(N ,2);
ur2 = copy(ur1);

u = reshape(u,tuple(n...));

@elapsed ur1[:,1] = reshape(rotateAndMove3D(u,theta_phi[1,:],b[1,:]./Mesh.h)[1],N,1);
@elapsed ur1[:,2] = reshape(rotateAndMove3D(u,theta_phi[2,:],b[2,:]./Mesh.h)[1],N,1);;

mr2,Jrot = rotateAndMoveRBFsimple(m,Mesh,theta_phi,b;computeJacobian = 0,numParamOfRBF=numParamOfRBF);
@elapsed ur2[:,1] = ParamLevelSetModelFunc(Mesh,mr2[:,1];computeJacobian=0,Xc = Xc,numParamOfRBF=numParamOfRBF)[1];
@elapsed ur2[:,2] = ParamLevelSetModelFunc(Mesh,mr2[:,2];computeJacobian=0,Xc = Xc,numParamOfRBF=numParamOfRBF)[1];

err = norm(ur1-ur2)./(2*norm(u));
if err < 2e-1
	println("RBF Rotation test succeeded: ",err);
else
	error("RBF Rotation test not succeeded: ",err);
end

println("Testing the sensitivity of the rotation (simple version - fixed theta/phi and b):");
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
mr2,Jrot = rotateAndMoveRBFsimple(m,Mesh,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
m0 = m;
dm = 0.01*randn(length(m));

# for k=1:nRBF
	# offset = (k-1)*numParamOfRBF + 1;
	# dm[offset] = 0.01*randn();
	# dm[(offset+1):(offset+6)] = 0.01*randn(6);
	# dm[(offset+7):(offset+9)] = 0.01*randn(3);
# end

hh = 1.0;
for k=0:5
	hhh = (0.5^k)*hh;
	mt = rotateAndMoveRBFsimple(m + hhh*dm,Mesh,theta_phi,b;computeJacobian = 0,numParamOfRBF=numParamOfRBF)[1];
	mtt= mr2 + reshape(Jrot*hhh*dm,size(mr2));
	println("norm(mt-mr2): ",norm(mt[:]-mr2[:]),", norm(mt - mr2 - Jrot*dm): ",norm(mt[:] - mr2[:] - Jrot*hhh*dm));
end

println("Testing the sensitivity of the rotation (with varying theta/phi and b):");
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
nRBF = div(length(m),numParamOfRBF);
nMoves = size(theta_phi,1);
m_wrapped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b);
m_t,theta_phi_t,b_t = splitRBFparamAndRotationsTranslations(m_wrapped,nRBF,nMoves,numParamOfRBF);

if norm(m_t - m) > 1e-5 || norm(theta_phi_t - theta_phi)>1e-5 || norm(b_t - b)>1e-5
	error("wrap/split RBFparamAndRotationsTranslations failed");
end

m0 = wrapRBFparamAndRotationsTranslations(m,theta_phi,b);
dm = 0.01*randn(length(m0));

dm,dtheta_phi,db = splitRBFparamAndRotationsTranslations(m_wrapped,nRBF,nMoves,numParamOfRBF);
dm = wrapRBFparamAndRotationsTranslations(dm,dtheta_phi,db);

mr2,Jrot = rotateAndMoveRBF(m,Mesh,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);

hh = 1.0;
for k=0:10
	hhh = (0.5^k)*hh;
	(mt,theta_phi_t,b_t) = splitRBFparamAndRotationsTranslations(m0+hhh*dm,nRBF,nMoves,numParamOfRBF);
	mt2 = rotateAndMoveRBF(mt,Mesh,theta_phi_t,b_t;computeJacobian = 0,numParamOfRBF=numParamOfRBF)[1];
	println("norm(mt2-mr2): ",norm(mt2[:]-mr2[:]),", norm(mt2 - mr2 - Jrot*dm): ",norm(mt2[:] - mr2[:] - Jrot*hhh*dm));
end


text:
module EventTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Event" => Ahorn.EntityPlacement(
        Maple.EventTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.EventTrigger)
    return Dict{String, Any}(
        "event" => Maple.event_trigger_events
    )
end

end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the Scalar type.
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test
import Random.rand

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Constructor tests

@testset "Scalar: inner constructor" begin
    #=
      Notes
      -----
      * Test value of constructed instance
    =#

    # --- Scalar{T}(value::AbstractFloat)

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    for precision_type in subtypes(AbstractFloat)
        for value_type in subtypes(AbstractFloat)
            converted_test_value = value_type(test_value)

            # value != 0, value != 1
            S = Scalar{precision_type}(converted_test_value)
            @test S isa Scalar{precision_type}
            @test S.value isa precision_type
            @test S.value == precision_type(converted_test_value)

            # value = 0
            S = Scalar{precision_type}(value_type(0))
            @test S isa Zero{precision_type}

            # value = 1
            S = Scalar{precision_type}(value_type(1))
            @test S isa One{precision_type}
        end
    end

    # --- Scalar{T}(value::Integer)

    test_value = (rand() > 0.5) ? 10 : -10

    for precision_type in subtypes(AbstractFloat)
        # subtypes(Signed)
        for value_type in subtypes(Signed)
            converted_test_value = value_type(test_value)

            # value != 0, value != 1
            S = Scalar{precision_type}(converted_test_value)
            S isa Scalar{precision_type}
            @test S.value isa precision_type
            @test S.value == precision_type(converted_test_value)

            # value = 0
            S = Scalar{precision_type}(value_type(0))
            S isa Zero{precision_type}

            # value = 1
            S = Scalar{precision_type}(value_type(1))
            S isa One{precision_type}
        end

        # subtypes(Unsigned)
        for value_type in subtypes(Unsigned)
            converted_test_value = value_type(abs(test_value))

            # value != 0, value != 1
            S = Scalar{precision_type}(converted_test_value)
            S isa Scalar{precision_type}
            @test S.value isa precision_type
            @test S.value == precision_type(converted_test_value)

            # value = 0
            S = Scalar{precision_type}(value_type(0))
            @test S isa Zero{precision_type}

            # value = 1
            S = Scalar{precision_type}(value_type(1))
            S isa One{precision_type}
        end

        # Bool
        S = Scalar{precision_type}(true)
        @test S isa One{precision_type}

        S = Scalar{precision_type}(false)
        @test S isa Zero{precision_type}
    end
end

@testset "Scalar: outer constructor - basic constructors" begin
    #=
      Notes
      -----
      * Test type of constructed instances. Correct construction of instances
        is tested by the inner constructor tests.
    =#

    # --- Scalar(value::AbstractFloat)

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    for precision_type in subtypes(AbstractFloat)
        # value != 0, value != 1
        converted_test_value = precision_type(test_value)
        S = Scalar(converted_test_value)
        @test S isa Scalar{precision_type}
        @test S.value == precision_type(converted_test_value)

        # value = 0
        converted_test_value = precision_type(0)
        S = Scalar(converted_test_value)
        @test S isa Zero{precision_type}

        # value = 1
        converted_test_value = precision_type(1)
        S = Scalar(converted_test_value)
        @test S isa One{precision_type}
    end

    # --- Scalar(value::AbstractIrrational)

    # value = π
    test_value = π
    S = Scalar(test_value)
    @test S isa Scalar{Float64}
    @test S.value == Float64(test_value)

    # --- Scalar(value::Integer)

    test_value = (rand() > 0.5) ? 10 : -10

    # subtypes(Signed)
    for value_type in subtypes(Signed)
        # value != 0, value != 1
        S = Scalar(value_type(test_value))
        @test S isa Scalar{Float64}
        @test S.value == test_value

        # value = 0
        converted_test_value = value_type(0)
        S = Scalar(converted_test_value)
        @test S isa Zero{Float64}

        # value = 1
        converted_test_value = value_type(1)
        S = Scalar(converted_test_value)
        @test S isa One{Float64}
    end

    # subtypes(Unsigned)
    for value_type in subtypes(Unsigned)
        # value != 0, value != 1
        S = Scalar(value_type(abs(test_value)))
        @test S isa Scalar{Float64}
        @test S.value == abs(test_value)

        # value = 0
        converted_test_value = value_type(0)
        S = Scalar(converted_test_value)
        @test S isa Zero{Float64}

        # value = 1
        converted_test_value = value_type(1)
        S = Scalar(converted_test_value)
        @test S isa One{Float64}
    end

    # Bool
    S = Scalar(true)
    @test S isa One{Float64}

    S = Scalar(false)
    @test S isa Zero{Float64}

    # --- Scalar(value::Rational)

    # value != 0, value != 1
    test_value = 3//5
    S = Scalar(test_value)
    @test S isa Scalar{Float64}
    @test S.value == Float64(test_value)

    # value = 0
    S = Scalar(Rational(0))
    @test S isa Zero{Float64}

    # value = 1
    S = Scalar(Rational(1))
    @test S isa One{Float64}
end

# --- Test attribute methods

@testset "Scalar: AbstractMultivector attribute functions" begin
    # --- Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # --- Tests

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        converted_test_value = precision_type(test_value)

        # value > 0
        positive_test_value =
            converted_test_value > 0 ? converted_test_value : -converted_test_value
        S = Scalar{precision_type}(positive_test_value)
        @test dim(S) == 0
        @test grades(S) == [0]
        @test blades(S) == [S]
        @test norm(S) isa precision_type
        @test norm(S) == abs(converted_test_value)

        @test S[0] == [S]
        @test S[1] == []

        # value < 0
        negative_test_value =
            converted_test_value < 0 ? converted_test_value : -converted_test_value
        S = Scalar{precision_type}(negative_test_value)
        @test norm(S) isa precision_type
        @test norm(S) == abs(converted_test_value)

        # value = Inf
        S = Scalar(precision_type(Inf))
        @test norm(S) isa precision_type
        @test norm(S) == Inf

        # value = -Inf
        S = Scalar(precision_type(-Inf))
        @test norm(S) isa precision_type
        @test norm(S) == Inf
    end
end

@testset "Scalar: AbstractBlade attribute functions" begin
    # --- Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # --- Tests

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        converted_test_value = precision_type(test_value)

        # value > 0
        positive_test_value =
            converted_test_value > 0 ? converted_test_value : -converted_test_value
        S = Scalar(positive_test_value)
        @test grade(S) == 0
        @test basis(S) == 1
        @test volume(S) isa precision_type
        @test volume(S) == positive_test_value
        @test sign(S) == 1

        # value < 0
        negative_test_value =
            converted_test_value > 0 ? -converted_test_value : converted_test_value
        S = Scalar(negative_test_value)
        @test volume(S) isa precision_type
        @test volume(S) == negative_test_value
        @test sign(S) == -1

        # value = Inf
        S = Scalar(precision_type(Inf))
        @test volume(S) isa precision_type
        @test volume(S) == precision_type(Inf)
        @test sign(S) == 1

        # value = -Inf
        S = Scalar(precision_type(-Inf))
        @test volume(S) isa precision_type
        @test volume(S) == precision_type(-Inf)
        @test sign(S) == -1
    end
end

@testset "Scalar: AbstractScalar attribute functions" begin
    # Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # Exercise functionality and check results
    for precision_type in subtypes(AbstractFloat)
        # Preparations
        converted_test_value = precision_type(test_value)

        # value > 0
        positive_test_value =
            converted_test_value > 0 ? converted_test_value : -converted_test_value
        S = Scalar(positive_test_value)
        @test value(S) isa precision_type
        @test value(S) == positive_test_value

        # value < 0
        negative_test_value =
            converted_test_value > 0 ? -converted_test_value : converted_test_value
        S = Scalar(negative_test_value)
        @test value(S) isa precision_type
        @test value(S) == negative_test_value

        # value = Inf
        S = Scalar(precision_type(Inf))
        @test value(S) isa precision_type
        @test value(S) == precision_type(Inf)

        # value = -Inf
        S = Scalar(precision_type(-Inf))
        @test value(S) isa precision_type
        @test value(S) == precision_type(-Inf)
    end
end

# --- Tests for AbstractMultivector interface functions

@testset "Scalar: -(B)" begin
    # Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # Tests
    for precision_type in subtypes(AbstractFloat)
        B = Scalar{precision_type}(test_value)

        negative_B = -B
        @test negative_B isa Scalar{precision_type}
        @test negative_B == Scalar{precision_type}(-test_value)
    end
end

@testset "Scalar: reverse(B)" begin
    # --- Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # --- Exercise functionality and check results

    # value != 0
    B = Scalar(test_value)
    reverse_B = reverse(B)
    @test reverse_B === B
    @test B * reverse_B ≈ norm(B)^2

    # value = Inf
    B = Scalar(Inf)
    @test reverse(B) === B

    # value = -Inf
    B = Scalar(-Inf)
    @test reverse(B) === B
end

@testset "Scalar: dual(B)" begin
    # Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # Tests
    for precision_type in subtypes(AbstractFloat)
        B = Scalar{precision_type}(test_value)

        for test_dim in 5:8
            dual_B = dual(B, test_dim)

            expected_result = if mod(test_dim, 4) < 2
                Pseudoscalar{precision_type}(test_dim, precision_type(test_value))
            else
                Pseudoscalar{precision_type}(test_dim, precision_type(-test_value))
            end
            @test dual_B isa Pseudoscalar{precision_type}
            @test dual_B == expected_result
        end

        expected_message = "The dual of a scalar is not well-defined if `dim` is not specified"
        @test_throws ErrorException(expected_message) dual(B)
    end
end

@testset "Scalar: convert(S)" begin
    # Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # Exercise functionality and check results
    for precision_type_converted in subtypes(AbstractFloat)
        for precision_type_src in subtypes(AbstractFloat)
            converted_test_value = precision_type_src(test_value)
            S = Scalar{precision_type_src}(converted_test_value)
            S_converted = convert(precision_type_converted, S)
            @test S_converted isa Scalar{precision_type_converted}
            if precision_type_src == precision_type_converted
                @test S_converted === S
            else
                @test S_converted !== S
                @test S_converted ≈ S
            end
        end
    end
end

# --- Tests for AbstractBlade interface functions

@testset "Scalar: inv(B)" begin
    # Preparations

    test_value = get_random_value(2) # add 2 to avoid 0 and 1

    # Tests
    for precision_type in subtypes(AbstractFloat)
        converted_value = precision_type(test_value)

        # value > 0
        B = Scalar(abs(converted_value))
        inverse_B = inv(B)
        @test inverse_B isa Scalar{precision_type}
        @test inverse_B == Scalar(1 / abs(converted_value))
        @test B * inverse_B ≈ 1

        # value < 0
        negative_value = -(abs(converted_value))
        B = Scalar(negative_value)
        inverse_B = inv(B)
        @test inverse_B isa Scalar{precision_type}
        @test inverse_B == Scalar(1 / negative_value)
        @test B * inverse_B ≈ 1

        # value = Inf
        B = Scalar(precision_type(Inf))
        @test inv(B) isa Zero{precision_type}

        # value = -Inf
        B = Scalar(precision_type(-Inf))
        @test inv(B) isa Zero{precision_type}
    end
end


text:
"""
# `IPFitting.jl`

Package for defining and fitting interatomic potentials based on the
N-Body expansion (ANOVA, HDMR, ...). The main exported type is
`NBodyIP` which is a `JuLIP` calculator.

See `?...` on how to
* `?IPFitting.Fitting` : fit an `NBodyIP`
* `?IPFitting.Data` : load data sets
* `?IPFitting.IO` : write and read an `NBodyIP`
"""
__precompile__()
module IPFitting

using Reexport

include("tools.jl")

include("prototypes.jl")

include("obsiter.jl")

include("datatypes.jl")

# loading data
include("data.jl")
@reexport using IPFitting.Data

include("lsq_db.jl")
@reexport using IPFitting.DB

include("lsqerrors.jl")
@reexport using IPFitting.Errors

include("lsq.jl")
@reexport using IPFitting.Lsq

include("preconditioners.jl")
@reexport using IPFitting.Preconditioners

include("auxiliary.jl")

end # module


text:
@everywhere module TestDecisionObjective

using StochasticPrograms
using Test

function test_SingleDecision_objective(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Min, x)
    end
    @second_stage sp = begin
        @recourse(model, y)
        @objective(model, Min, y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionRef
    @test JuMP.objective_function(sp, 1) == x
    @test x == @inferred JuMP.objective_function(sp, 1, DecisionRef)
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2, JuMP.objective_function(sp))
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.objective_function(sp) == x
        @test x == @inferred JuMP.objective_function(sp, DecisionRef)
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        @test JuMP.isequal_canonical(x1 + y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(x2 + y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Max, x)
    end
    @second_stage sp = begin
        @recourse(model, y)
        @objective(model, Max, y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionRef
    @test JuMP.objective_function(sp, 1) == x
    @test x == @inferred JuMP.objective_function(sp, 1, DecisionRef)
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x + 0.5*y1 + 0.5*y2)
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 1) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test JuMP.objective_function(sp) == x
        @test x == @inferred JuMP.objective_function(sp, DecisionRef)
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        @test JuMP.isequal_canonical(x1 + y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(x2 + y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_DecisionAffExpr_objective(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Min, 2x)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, y)
        @objective(model, Min, a*y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionAffExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 2x)
    @test JuMP.isequal_canonical(
        2x, @inferred JuMP.objective_function(sp, 1, DecisionAffExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
    @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        @test JuMP.isequal_canonical(
            2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x + y1 + 2y2)
        @test JuMP.isequal_canonical(2x + y1 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        @test JuMP.isequal_canonical(
            2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x)
        @test JuMP.isequal_canonical(2x,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x1 + 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x2 + 4y2)
        @test JuMP.isequal_canonical(
            2x1 + 2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            2x2 + 4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Max, x + 3x + 1)
    end
    @second_stage sp = begin
        @known x
        @uncertain a
        @recourse(model, y)
        @objective(model, Max, x + a*x + 1 + y + 3y + 1)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionAffExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 4x + 1)
    @test JuMP.isequal_canonical(
        4x + 1, @inferred JuMP.objective_function(sp, 1, DecisionAffExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
    @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 3x + 4y1 + 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 5x + 4y2 + 2)
        @test JuMP.isequal_canonical(
            3x + 4y1 + 2, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            5x + 4y2 + 2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 8x + 2y1 + 2y2 + 3)
        @test JuMP.isequal_canonical(8x + 2y1 + 2y2 + 3,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 3x + 4y1 + 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 5x + 4y2 + 2)
        @test JuMP.isequal_canonical(
            3x + 4y1 + 2, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            5x + 4y2 + 2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 4x + 1)
        @test JuMP.isequal_canonical(4x + 1,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 7x1 + 4y1 + 3)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 9x2 + 4y2 + 3)
        @test JuMP.isequal_canonical(
            7x1 + 4y1 + 3, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            9x2 + 4y2 + 3, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_DecisionQuadExpr_objective(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Min, x^2 + 2x)
    end
    @second_stage sp = begin
        @known x
        @uncertain a
        @recourse(model, y)
        @objective(model, Min, y^2 + a*y + x^2 + 2x)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionQuadExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), x^2 + 2x)
    @test JuMP.isequal_canonical(
        x^2 + 2x, @inferred JuMP.objective_function(sp, 1, DecisionQuadExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionQuadExpr{Float64}
    # Structure specific
    x = DecisionRef(sp[1,:x])
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), y1^2 + 2y1 + x^2 + 2x)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), y2^2 + 4y2 + x^2 + 2x)
        @test JuMP.isequal_canonical(
            y1^2 + 2y1 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            y2^2 + 4y2 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x^2 + 4x + 0.5*y1^2 + y1 + 0.5*y2^2 + 2y2)
        @test JuMP.isequal_canonical(2x^2 + 4x + 0.5*y1^2 + y1 + 0.5*y2^2 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), y1^2 + 2y1 + x^2 + 2x)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), y2^2 + 4y2 + x^2 + 2x)
        @test JuMP.isequal_canonical(
            y1^2 + 2y1 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            y2^2 + 4y2 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x^2 + 2x)
        @test JuMP.isequal_canonical(x^2 + 2x,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x1^2 + 4x1 + y1^2 + 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x2^2 + 4x2 + y2^2 + 4y2)
        @test JuMP.isequal_canonical(
            2x1^2 + 4x1 + y1^2 + 2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            2x2^2 + 4x2 + y2^2 + 4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_decision_objective_coefficient_modifiction(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Min, x)
    end
    @second_stage sp = begin
        @known x
        @recourse(model, y)
        @objective(model, Min, y)
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_objective_coefficient(sp, x, 1, 4.0)
    x = DecisionRef(x)
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 4x)
    # Second-stage
    x = sp[1,:x]
    y = sp[2,:y]
    JuMP.set_objective_coefficient(sp, y, 2, 1, 4.0)
    JuMP.set_objective_coefficient(sp, y, 2, 2, 4.0)
    y1 = DecisionRef(y, 1)
    y2 = DecisionRef(y, 2)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 7x + 2y1 + 2y2)
        @test JuMP.isequal_canonical(7x + 2y1 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 4x)
        @test JuMP.isequal_canonical(4x,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        @test JuMP.isequal_canonical(4x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(model, x)
        @objective(model, Min, x^2 + x)
    end
    @second_stage sp = begin
        @known x
        @recourse(model, y)
        @objective(model, Min, x^2 + y^2 + y)
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_objective_coefficient(sp, x, 1, 4.0)
    x = DecisionRef(x)
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), x^2 + 4x)
    # Second-stage
    x = sp[1,:x]
    y = sp[2,:y]
    JuMP.set_objective_coefficient(sp, y, 2, 1, 4.0)
    JuMP.set_objective_coefficient(sp, y, 2, 2, 4.0)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4y2)
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x^2 + 7x + 0.5*y1^2 + 2y1 + 0.5*y2^2 + 2y2)
        @test JuMP.isequal_canonical(2x^2 + 7x + 0.5*y1^2 + 2y1 + 0.5*y2^2 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.VerticalStructure
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4y2)
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x^2 + 4x)
        @test JuMP.isequal_canonical(x^2 + 4x,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.HorizontalStructure
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x^2 + y1^2 + 4x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x^2 + y2^2 + 4x + 4y2)
        @test JuMP.isequal_canonical(2x^2 + y1^2 + 4x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(2x^2 + y2^2 + 4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x^2 + y2^2 + 4x + 4y2)
        @test JuMP.isequal_canonical(2x^2 + y1^2 + 2x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(2x^2 + y2^2 + 4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function runtests()
    @testset "DecisionObjective" begin
        for structure in [(Deterministic(),),
                          (Vertical(),),
                          (Horizontal(),),
                          (Deterministic(), () -> MOIU.MockOptimizer(MOIU.Model{Float64}()))]
            name = length(structure) == 1 ? "$(structure[1])" : "$(structure[1]) with decision bridges"
            @testset "$name" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

function run_dtests()
    @testset "DecisionObjective" begin
        for structure in [(DistributedVertical(),), (DistributedHorizontal(),)]
            @testset "$(structure[1])" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

end


text:
using DataDeps

register(DataDep(
    "WikiText-103",
    """
    WikiText Long Term Dependency Language Modeling Dataset
    This is a language modelling dataset under Creative Commons Attribution-ShareAlike License. I contains over 100
    million tokens.
    """,
    "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-103-v1.zip",
    post_fetch_method = function (file)
        unpack(file)
        dir = "wikitext-103"
        files = readdir(dir)
        mv.(joinpath.(dir, files), files)
        rm(dir)

        corpus = read(open("wiki.train.tokens", "r"), String)
        corpus = split(lowercase(corpus), '\n')
        deleteat!(corpus, findall(x -> isequal(x, "")||isequal(x, " ")||(isequal(x[1:2], " =")&&isequal(x[prevind(x, lastindex(x), 1):end], "= ")), corpus))
        corpus = strip.(corpus)
        corpus .*= " <eos>"
        corpus = join(corpus, ' ')
        ftrain = open("wiki.train.tokens", "w")
        write(ftrain, corpus)
        close(ftrain)
    end
))


text:
using Base.Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Moving Averages" begin

@testset "Moving Averages on TimeArrays" begin
    @testset "sma" begin
        ta = sma(cl, 10)
        @test isapprox(ta.values[1]  , 98.782, atol=.01)  # TTR value 98.782
        @test isapprox(ta.values[2]  , 97.982, atol=.01)  # TTR value 97.982
        @test isapprox(ta.values[3]  , 98.388, atol=.01)  # TTR value 98.388
        @test isapprox(ta.values[490], 21.266, atol=.01)  # TTR value 21.266
        @test isapprox(ta.values[491], 21.417, atol=.01)  # TTR value 21.417
        @test ta.timestamp[491] == Date(2001,12,31)

        ta = sma(ohlc, 10)
        @test isapprox(ta.values[1,:]  , [100.692, 103.981, 96.001, 98.782] , atol=.01)
        @test isapprox(ta.values[2,:]  , [100.304, 103.331, 95.876, 97.982], atol=.01)
        @test isapprox(ta.values[3,:]  , [100.041, 103.144, 96.095,98.388], atol=.01)
        @test isapprox(ta.values[490,:], [21.081, 21.685, 20.797, 21.266], atol=.01)
        @test isapprox(ta.values[491,:], [21.259, 21.868, 20.971, 21.417], atol=.01)
        @test ta.timestamp[491] == Date(2001,12,31)
    end

    @testset "ema" begin
        ta = ema(cl, 10)
        @test isapprox(ta.values[1]  , 98.782, atol=.01)   # TTR value 98.78200
        @test isapprox(ta.values[2]  , 99.719, atol=.01)   # TTR value 99.71982
        @test isapprox(ta.values[3]  , 100.963, atol=.01)  # TTR value 100.96349
        @test isapprox(ta.values[490], 21.585, atol=.01)   # TTR value 21.58580
        @test isapprox(ta.values[491], 21.642, atol=.01)   # TTR value 21.64293
        @test ta.timestamp[491] == Date(2001,12,31)

        ta = ema(ohlc, 10)
        @test isapprox(ta.values[1,:]  , [100.692, 103.981, 96.001, 98.782], atol=.01)
        @test isapprox(ta.values[2,:]  , [100.748, 104.348, 96.808, 99.719], atol=.01)
        @test isapprox(ta.values[3,:]  , [101.634, 105.148, 98.003, 100.963], atol=.01)
        @test isapprox(ta.values[490,:], [21.369, 22.030, 21.125, 21.585], atol=.01)
        @test isapprox(ta.values[491,:], [21.576, 22.145, 21.253, 21.642], atol=.01)
        @test ta.timestamp[491] == Date(2001,12,31)
    end

    @testset "ema wilder true" begin
        ta = ema(cl, 10, wilder=true)
        @test isapprox(ta.values[1]  , 98.782, atol=.01)  # TTR value 98.7820
        @test isapprox(ta.values[2]  , 99.297, atol=.01)  # TTR value 99.2978
        @test isapprox(ta.values[3]  , 100.024, atol=.01)  # TTR value 100.0240
        @test isapprox(ta.values[490], 21.345, atol=.01)  # TTR value 21.34556
        @test isapprox(ta.values[491], 21.401, atol=.01)  # TTR value 21.40100
        @test ta.timestamp[491] == Date(2001,12,31)

        ta = ema(ohlc, 10, wilder=true)
        @test isapprox(ta.values[1,:]  , [100.692, 103.981, 96.001, 98.782], atol=.01)
        @test isapprox(ta.values[2,:]  , [100.723, 104.183, 96.444, 99.297], atol=.01)
        @test isapprox(ta.values[3,:]  , [101.213, 104.64, 97.138, 100.024], atol=.01)
        @test isapprox(ta.values[490,:], [21.184, 21.776, 20.847, 21.345], atol=.01)
        @test isapprox(ta.values[491,:], [21.317, 21.865, 20.945, 21.401], atol=.01)
        @test ta.timestamp[491] == Date(2001,12,31)
    end

    @testset "kama" begin
        ta = kama(cl)
        @test isapprox(ta.values[1], 98.9052, atol=.01)
        @test isapprox(ta.values[2], 99.0098, atol=.01)
        @test isapprox(ta.values[3], 99.4499, atol=.01)
        @test ta.timestamp[1] == Date(2000, 1, 18)
        @test ta.colnames     == ["kama"]

        ta = kama(ohlc)
        @test length(ta.colnames) == 4
        @test ta.timestamp[1]     == Date(2000, 1, 18)

        ta = kama(TimeArray(collect(Date(2011, 1, 1):Date(2011, 1, 20)), 1:20))
        @test ta.timestamp[end] == Date(2011, 1, 20)

        ta = kama(TimeArray(collect(Date(2011, 1, 1):Date(2011, 1, 20)), fill(42, 20)))
        @test ta.values == fill(42, length(ta))
    end

    @testset "env" begin
        ta = env(cl, 10)
        @test isapprox(ta.values[1,2]  , 108.66, atol=.01)
        @test isapprox(ta.values[2,2]  , 107.78, atol=.01)
        @test isapprox(ta.values[490,2]  , 23.3926, atol=.01)
        @test isapprox(ta.values[491,2]  , 23.5587, atol=.01)
        @test isapprox(ta.values[1,1]  , 88.9038, atol=.01)
        @test isapprox(ta.values[2,1]  , 88.1838, atol=.01)
        @test isapprox(ta.values[490,1]  , 19.1394, atol=.01)
        @test isapprox(ta.values[491,1]  , 19.2753, atol=.01)

        ta = env(ohlc, 10)
        @test isapprox(ta.values[1,5:end]  , [110.761, 114.379, 105.601, 108.66] , atol=.01)
        @test isapprox(ta.values[2,5:end]  , [110.334, 113.664, 105.464, 107.78], atol=.01)
        @test isapprox(ta.values[490,5:end]  , [23.1891, 23.8535, 22.8767, 23.3926], atol=.01)
        @test isapprox(ta.values[491,5:end]  , [23.3849, 24.0548, 23.0681, 23.5587], atol=.01)
        @test isapprox(ta.values[1,1:4]  , [90.6228, 93.5829, 86.4009, 88.9038] , atol=.01)
        @test isapprox(ta.values[2,1:4]  , [90.2736, 92.9979, 86.2884, 88.1838], atol=.01)
        @test isapprox(ta.values[490,1:4]  , [18.9729, 19.5165, 18.7173, 19.1394], atol=.01)
        @test isapprox(ta.values[491,1:4]  , [19.1331, 19.6812, 18.8739, 19.2753], atol=.01)
    end
end

@testset "Moving Averages on arrays" begin
    @testset "Array dispatch on sma" begin
        arr = sma(cl.values, 10)
        # same values as above TimeArray dispatch
        @test isapprox(arr[1]  , 98.782, atol=.01)
        @test isapprox(arr[2]  , 97.982, atol=.01)
        @test isapprox(arr[3]  , 98.388, atol=.01)
        @test isapprox(arr[490], 21.266, atol=.01)
        @test isapprox(arr[491], 21.417, atol=.01)

        arr = sma(ohlc.values, 10)
        @test isapprox(arr[1,:]  , [100.692, 103.981, 96.001, 98.782] , atol=.01)
        @test isapprox(arr[2,:]  , [100.304, 103.331, 95.876, 97.982], atol=.01)
        @test isapprox(arr[3,:]  , [100.041, 103.144, 96.095,98.388], atol=.01)
        @test isapprox(arr[490,:], [21.081, 21.685, 20.797, 21.266], atol=.01)
        @test isapprox(arr[491,:], [21.259, 21.868, 20.971, 21.417], atol=.01)
    end

    @testset "Array dispatch on ema" begin
        arr = ema(cl.values, 10)
        @test isapprox(arr[1]  , 98.782, atol=.01)
        @test isapprox(arr[2]  , 99.719, atol=.01)
        @test isapprox(arr[3]  , 100.963, atol=.01)
        @test isapprox(arr[490], 21.585, atol=.01)
        @test isapprox(arr[491], 21.642, atol=.01)

        arr = ema(ohlc.values, 10)
        @test isapprox(arr[1,:]  , [100.692, 103.981, 96.001, 98.782], atol=.01)
        @test isapprox(arr[2,:]  , [100.748, 104.348, 96.808, 99.719], atol=.01)
        @test isapprox(arr[3,:]  , [101.634, 105.148, 98.003, 100.963], atol=.01)
        @test isapprox(arr[490,:], [21.369, 22.030, 21.125, 21.585], atol=.01)
        @test isapprox(arr[491,:], [21.576, 22.145, 21.253, 21.642], atol=.01)
        @test ema(ohlc, 10).timestamp[491] == Date(2001,12,31)
    end

    @testset "Array dispatch on ema wilder true" begin
        arr = ema(cl.values, 10, wilder=true)
        @test isapprox(arr[1]  , 98.782, atol=.01)
        @test isapprox(arr[2]  , 99.297, atol=.01)
        @test isapprox(arr[3]  , 100.024, atol=.01)
        @test isapprox(arr[490], 21.345, atol=.01)
        @test isapprox(arr[491], 21.401, atol=.01)

        arr = ema(ohlc.values, 10, wilder=true)
        @test isapprox(arr[1,:]  , [100.692, 103.981, 96.001, 98.782], atol=.01)
        @test isapprox(arr[2,:]  , [100.723, 104.183, 96.444, 99.297], atol=.01)
        @test isapprox(arr[3,:]  , [101.213, 104.64, 97.138, 100.024], atol=.01)
        @test isapprox(arr[490,:], [21.184, 21.776, 20.847, 21.345], atol=.01)
        @test isapprox(arr[491,:], [21.317, 21.865, 20.945, 21.401], atol=.01)
    end

    @testset "Array dispatch on env" begin
        ta = env(cl.values, 10)
        @test isapprox(ta[1,:]  , [88.9038, 108.66], atol=.01)
        @test isapprox(ta[2,:]  , [88.1838, 107.78], atol=.01)
        @test isapprox(ta[490,:]  , [19.1394, 23.3926], atol=.01)
        @test isapprox(ta[491,:]  , [19.2753, 23.5587], atol=.01)

        ta = env(ohlc.values, 10)
        @test isapprox(ta[1,:]  , [90.6228, 93.5829, 86.4009, 88.9038, 110.761, 114.379, 105.601, 108.66] , atol=.01)
        @test isapprox(ta[2,:]  , [90.2736, 92.9979, 86.2884, 88.1838, 110.334, 113.664, 105.464, 107.78], atol=.01)
        @test isapprox(ta[490,:]  , [18.9729, 19.5165, 18.7173, 19.1394, 23.1891, 23.8535, 22.8767, 23.3926], atol=.01)
        @test isapprox(ta[491,:]  , [19.1331, 19.6812, 18.8739, 19.2753, 23.3849, 24.0548, 23.0681, 23.5587], atol=.01)
    end
end

end  # @testset "Moving Averages"


text:
"""
*TAX*

`TAX` constructs a model object for transfer of attention exchange.

- `δ`: transfer of attention parameter
- `γ`: probability weighting parameter
- `β`: utility curvature

Constructor
````julia
TAX(;δ=.80, β=.3, γ=.70)
````
*References*

Birnbaum, M. H., & Chavez, A. (1997). Tests of theories of decision making: Violations of branch independence and distribution independence. Organizational Behavior and human decision Processes, 71(2), 161-194.
Birnbaum, M. H. (2008). New paradoxes of risky decision making. Psychological review, 115(2), 463.

"""
mutable struct TAX{T1,T2,T3} <:UtilityModel
    δ::T1
    γ::T2
    β::T3
end

function TAX(;δ=-1.0, β=1.0, γ=.70)
    return TAX(δ, γ, β)
end

"""
*compute_utility*

`compute_utility` computes utility of gamble outcomes according to TAX

- `model`: a model object for TAX
- `gamble`: a gamble object

Function Signature
````julia
compute_utility(model::TAX, gamble::Gamble)
````
"""
function compute_utility(model::TAX, gamble)
    @unpack β = model
    @unpack v = gamble
    utility = @. sign(v)*abs(v)^β
    return utility
end

tax_weight(p, γ) = (p^γ)

function ω(p, pk, n, δ, γ)
    if δ > 0
        return δ*tax_weight(pk, γ)/(n+1)
    end
    return δ*tax_weight(p, γ)/(n+1)
end

function sort!(model::TAX, gamble)
    @unpack p,v = gamble
    i = sortperm(v)
    p .= p[i]; v .= v[i]
    return nothing
end

"""
*mean*

`mean` a method for computing the mean for the TAX model

- `model`: a model M <: UtilityModel
- `gamble`: a gamble object

Function Signature
````julia
mean(model::TAX, gamble::Gamble)
````
"""
function mean(model::TAX, gamble::Gamble)
    @unpack p,v = gamble
    @unpack γ,δ = model
    n = length(p)
    sort!(model, gamble)
    utility = compute_utility(model, gamble)
    eu = 0.0
    sum_weight = 0.0
    for i in 1:n
        eu += tax_weight(p[i], γ)*utility[i]
        sum_weight += tax_weight(p[i], γ)
        for k in 1:(i-1)
            eu += (utility[i] - utility[k])*ω(p[i], p[k], n, δ, γ)
        end
    end
    return eu/sum_weight
end

function var(model::TAX, gamble::Gamble)
    error("var not implimented for TAX")
    return -100.0
end


text:
using Bio.Seq
using BioBenchmarks

@benchmark collect(open(datafile("chr1.fa"), FASTA))


text:
"""
    findmax(C::ClimGrid; skipnan::Bool=false)

Return the maximum element of ClimGrid C and its index. If there are multiple maximal elements, then the first one will be returned. If any data element is NaN, this element is returned. The result is in line with max. As climate data can often have NaN values (irregular polygons, missing weather data), the option to skip NaNs is provided as a keyword argument.

"""
function Base.findmax(C::ClimGrid; skipnan::Bool=false)
    # Get data
    data = C[1].data

    if skipnan
        idx = findall(data .== NaNMath.maximum(data))[1]
        val = data[idx]

    else
        idx = argmax(data)
        val = data[idx]
    end

    return val, idx

end

"""
    findmin(C::ClimGrid; skipnan::Bool=false)

Return the minimum element of ClimGrid C and its index. If there are multiple minimal elements, then the first one will be returned. If any data element is NaN, this element is returned. The result is in line with min. As climate data can often have NaN values (irregular polygons, missing weather data), the option to skip NaNs is provided as a keyword argument.

"""
function Base.findmin(C::ClimGrid; skipnan::Bool=false)
    # Get data
    data = C[1].data

    if skipnan
        idx = findall(data .== NaNMath.minimum(data))[1]
        val = data[idx]

    else
        idx = argmin(data)
        val = data[idx]
    end

    return val, idx

end

# """
#     dashboard(C::ClimGrid)
#
# This function returns the summary of ClimGrid C, such as maps of mean, maximum and minimum values in ClimGrid C. The annual cycle, the histogram and estimated probability density function and time series.
# """
# function dashboard(C::ClimGrid; poly = ([]))
#
#     # MAPS, WHILE NICE TO LOOK AT, ARE ALREADY COVERED BY THE FUNCTION "mapclimgrid". Should focus here on annual cycles, and other diagnostics such as annualmax and annualmin
#     # # Mean, maximum and minimum values for each grid points
#     # mean_map = mean(annualmean(C)[1].data, 1)[1, :, :]
#     # FD = AxisArray(mean_map, Axis{:lon}(C[1][Axis{:lon}][:]), Axis{:lat}(C[1][Axis{:lat}][:]))
#     # meanmap = ClimGrid(FD, model = C.model, experiment = C.experiment, run = C.run, filename = C.filename, dataunits = C.dataunits, latunits = C.latunits, lonunits = C.lonunits, variable = string(C.variable, " - temporal mean"), typeofvar = C.typeofvar, typeofcal = C.typeofcal)
#     #
#     # max_map = maximum(annualmax(C)[1].data, 1)[1, :, :]
#     # FD = AxisArray(max_map, Axis{:lon}(C[1][Axis{:lon}][:]), Axis{:lat}(C[1][Axis{:lat}][:]))
#     # maxmap = ClimGrid(FD, model = C.model, experiment = C.experiment, run = C.run, filename = C.filename, dataunits = C.dataunits, latunits = C.latunits, lonunits = C.lonunits, variable = string(C.variable, " - temporal maximum"), typeofvar = C.typeofvar, typeofcal = C.typeofcal)
#     #
#     # min_map = minimum(annualmin(C)[1].data, 1)[1, :, :]
#     # FD = AxisArray(min_map, Axis{:lon}(C[1][Axis{:lon}][:]), Axis{:lat}(C[1][Axis{:lat}][:]))
#     # minmap = ClimGrid(FD, model = C.model, experiment = C.experiment, run = C.run, filename = C.filename, dataunits = C.dataunits, latunits = C.latunits, lonunits = C.lonunits, variable = string(C.variable, " - temporal minimum"), typeofvar = C.typeofvar, typeofcal = C.typeofcal)
#
#     # Spatial-temporal histogram and PDF
#
#     # Annual cycles
#
#
#     # Timeseries (annual mean, max and min)
#     mean_ts = mean(annualmean(C)[1].data,2:3)[:, 1, 1]
#     max_ts = mean(annualmax(C)[1].data,2:3)[:, 1, 1]
#     min_ts = mean(annualmin(C)[1].data,2:3)[:, 1, 1]
#
#
#     # Annual cycle
#
#
#     # Plot everything
#     fig = figure("pyplot_subplot_column",figsize=(10,10))
#     subplot(231) # Create the 1st axis of a 3x1 array of axes
#     title("Mean")
#     subplot(232) # Create the 2nd axis of a 3x1 arrax of axes
#     ax = gca() # Get the handle of the current axis
#     #ax[:set_yscale]("log") # Set the y axis to a logarithmic scale
#     grid("on")
#     ylabel("Log Scale")
#     title("Minimum")
#     subplot(233) # Create the 3rd axis of a 3x1 array of axes
#     ax = gca()
#     ax[:set_xscale]("log") # Set the x axis to a logarithmic scale
#     #xlabel("Log Scale")
#     title("Maximum")
#     fig[:canvas][:draw]() # Update the figure
#     suptitle("3x1 Subplot")
#
#     return true
#
#
# end


text:
function local_updates(mc::AbstractDQMC)
  p = mc.p
  s = mc.s
  l = mc.l
  a = mc.a

  acc_rat = 0.0
  @inbounds for i in 1:l.sites
    @views new_op = p.hsfield[:,i,s.current_slice] + randuniform(p.box, p.opdim)
    exp_delta_S_boson = exp(- calc_boson_action_diff(mc, i, s.current_slice, new_op) )
    detratio = calc_detratio(mc,i,new_op)

    if p.opdim == 3
      p_acc_fermion = real(detratio)
    elseif p.opdim < 3
      p_acc_fermion = real(detratio * conj(detratio))
    end

    if p.opdim == 3 && abs(imag(detratio)/real(detratio)) > 1e-4
      @printf("%d, %d \t Determinant ratio isn't real. \t abs imag: %.3e \t relative: %.1f%%\n", s.current_slice, i, abs.(imag(detratio)), abs.(imag(detratio))/abs.(real(detratio))*100)
    # elseif real(detratio) < 0
    #   println("Negative fermion weight.")
    # elseif detratio == 0
    #   println("Encountered non-invertible M with det = 0.")
    elseif p.opdim < 3 && p_acc_fermion < 0
      println("Negative fermion weight!")
    end

    p_acc = exp_delta_S_boson * p_acc_fermion

    if p_acc > 1.0 || rand() < p_acc
      acc_rat += 1
      @views p.hsfield[:,i,s.current_slice] = new_op
      p.boson_action += -log(exp_delta_S_boson)
      update_greens!(mc,i)
    end
  end
  return acc_rat / l.sites
end

@inline function calc_detratio(mc::AbstractDQMC, i::Int, new_op::Vector{Float64})
  p = mc.p
  s = mc.s
  l = mc.l

  @mytimeit mc.a.to "calc_detratio" begin

  interaction_matrix_exp_op!(mc,p.hsfield[:,i,s.current_slice],-1.,s.eVop1) #V1i
  interaction_matrix_exp_op!(mc,new_op,1.,s.eVop2) #V2i
  mul!(s.eVop1eVop2, s.eVop1, s.eVop2)
  s.delta_i .= s.eVop1eVop2 .- s.eye_flv
  s.Mtmp .= s.eye_flv .- s.greens[i:l.sites:end,i:l.sites:end]
  mul!(s.Mtmp2, s.delta_i, s.Mtmp)
  s.M .= s.eye_flv .+ s.Mtmp2

  end #timeit
  return det(s.M)
end

@inline function update_greens!(mc::AbstractDQMC, i::Int)
  p = mc.p
  s = mc.s
  l = mc.l
  g = mc.s.greens
  ab = mc.s.AB

  @mytimeit mc.a.to "update_greens!" begin

  s.A = s.greens[:,i:l.sites:end]

  if p.opdim == 3
    @simd for k in 0:3
        s.A[i+k*l.sites,k+1] -= 1.
    end
  else
    @simd for k in 0:1
        s.A[i+k*l.sites,k+1] -= 1.
    end
  end

  s.A *= inv(s.M)
  mul!(s.B, s.delta_i, s.greens[i:l.sites:end,:])

  # benchmark: most of time (>90% is spend below)
  mul!(s.AB, s.A, s.B)

  # more explicit way of doing s.greens .+= s.AB
  @inbounds @simd for i in eachindex(g)
    g[i] += ab[i]
  end

  end #timeit
  nothing
end

# function get_full_delta(mc,i)
#   delta_i = mc.s.delta_i
#   N = mc.l.sites
#   Delta = zeros(mc.s.greens)
#   Delta[i:N:end, i:N:end] = delta_i
#   return Delta
# end

# function cols2mat(mc,cols,i)
#   N = mc.l.sites
#   mat = zeros(mc.s.greens)
#   mat[:,i:N:end] .= cols
#   return mat
# end

# function shouldbe(mc,i)
#   Delta = get_full_delta(mc,i)
#   D = sparse(Delta)
#   A = (mc.s.eye_full - mc.s.greens)
#   return A*D #[:,i:mc.l.sites:end]
# end

# function shouldbe2(mc,i)
#   S = shouldbe(mc,i)
#   return inv(one(S) + S) #[:,i:mc.l.sites:end]
# end

# function update_greens_max!(mc::AbstractDQMC, i::Int, cols, colsTimesG)
#   G = geltype(mc)
#   N = mc.l.sites
#   flv = mc.p.flv
#   g = mc.s.greens
#   delta_i = mc.s.delta_i

#   site = i

#   # Calculate (1-G)*Delta
#   # The result is a (sparse) matrix with just a bunch of non-zero columns.
#   # Hence, we only calculate those columns.
#   # [O(N)]
#   # cols = zeros(G, flv*N, flv)
#   @inbounds for c in 1:flv
#     for r in 1:flv*N
#       cols[r,c] = - g[r,site] * delta_i[1, c]
#     end
#     cols[site, c] += delta_i[1,c]
#     for dr in 2:flv
#       for r in 1:flv*N # max had 4 instead of flv. mistake!?
#         cols[r,c] += - g[r, site + (dr-1)*N] * delta_i[dr,c]
#       end
#       cols[site + (dr-1)*N, c] += delta_i[dr,c]
#     end
#   end

#   # return cols2mat(mc,cols,site)
#   # return cols

#   # Now, calculate [1 + (1-G)*Delta]^(-1) based on (1-G)*Delta result above.
#   # The result is, again, a sparse matrix with only few (p.flv) non-zero columns.
#   # However, there are ones on the diagonal as well.
#   determinant = 1.0 + 0.0im
#   @inbounds for j in 1:flv
#     col = cols[:,j] #view?
#     for k in (j-1):-1:1
#       col[site + (k-1)*N] .= 0
#     end
#     for k in (j-1):-1:1
#       col .+= cols[site + (k-1)*N, j] * cols[:,k]
#     end

#     # @views begin
#       divisor = 1 + col[site + (j-1)*N]
#       cols[:,j] .= -1/divisor * col
#       cols[site + (j-1)*N, j] += 1
#       for k in (j-1):-1:1
#         cols[:,k] .-= (cols[site + (j-1)*N, k] / divisor) * col;
#       end
#       determinant *= divisor
#     # end
#   end
#   # @show determinant
#   # return cols

#   # Compensate for already included diagonal entries of I in cols
#   cols[site,1] -= 1
#   cols[site+N,2] -= 1
#   if p.opdim == 3
#     cols[site+2*N,3] -= 1
#     cols[site+3*N,4] -= 1
#   end

#   perform_update_greens!(mc, i, cols, colsTimesG)
#   nothing

#   # Based on the following by Max:
#   # //****
#   # //Compute the determinant and inverse of I + Delta*(I - G)
#   # //based on Sherman-Morrison formula / Matrix-Determinant lemma
#   # //****

#   # //Delta*(I - G) is a sparse matrix containing just 4 rows:
#   # //site, site+N, site+2N, site+3N
#   # //Compute the values of these rows [O(N)]:
#   # for (uint32_t r = 0; r < MatrixSizeFactor; ++r) {
#   #     //TODO: Here are some unnecessary operations:
#   #     //delta_forsite contains many repeated elements, and even
#   #     //some zeros
#   #     rows[r] = VecData(MatrixSizeFactor*pars.N);
#   #     for (uint32_t col = 0; col < MatrixSizeFactor*pars.N; ++col) {
#   #         rows[r][col] = -delta_forsite(r,0) * g.col(col)[site];
#   #     }
#   #     rows[r][site] += delta_forsite(r,0);
#   #     for (uint32_t dc = 1; dc < MatrixSizeFactor; ++dc) {
#   #         for (uint32_t col = 0; col < 4*pars.N; ++col) {
#   #             rows[r][col] += -delta_forsite(r,dc) * g.col(col)[site + dc*pars.N];
#   #         }
#   #         rows[r][site + dc*pars.N] += delta_forsite(r,dc);
#   #     }
#   # }

#   # // [I + Delta*(I - G)]^(-1) again is a sparse matrix
#   # // with two (four) rows site, site+N(, site+2N, site+3N).
#   # // compute them iteratively, together with the determinant of
#   # // I + Delta*(I - G)
#   # // Apart from these rows, the remaining diagonal entries of
#   # // [I + Delta*(I - G)]^(-1) are 1
#   # //
#   # // before this loop rows[] holds the entries of Delta*(I - G),
#   # // after the loop rows[] holds the corresponding rows of [I + Delta*(I - G)]^(-1)
#   # DataType det = 1;
#   # for (uint32_t l = 0; l < MatrixSizeFactor; ++l) {
#   #     VecData row = rows[l];
#   #     for (int k = int(l)-1; k >= 0; --k) {
#   #         row[site + unsigned(k)*pars.N] = 0;
#   #     }
#   #     for (int k = int(l)-1; k >= 0; --k) {
#   #         row += rows[l][site + unsigned(k)*pars.N] * rows[(unsigned)k];
#   #     }
#   #     DataType divisor = DataType(1) + row[site + l*pars.N];
#   #     rows[l] = (-1.0/divisor) * row;
#   #     rows[l][site + l*pars.N] += 1;
#   #     for (int k = int(l) - 1; k >= 0; --k) {
#   #         rows[(unsigned)k] -= (rows[unsigned(k)][site + l*pars.N] / divisor) * row;
#   #     }
#   #     det *= divisor;
#   # }

#   # //compensate for already included diagonal entries of I in invRows
#   # rows[0][site] -= 1;
#   # rows[1][site + pars.N] -= 1;
#   # if (OPDIM == 3) {
#   #     rows[2][site + 2*pars.N] -= 1;
#   #     rows[3][site + 3*pars.N] -= 1;
#   # }
#   # //compute G' = G * [I + Delta*(I - G)]^(-1) = G * [I + invRows]
#   # // [O(N^2)]
#   # MatData gTimesInvRows(MatrixSizeFactor*pars.N,
#   #                       MatrixSizeFactor*pars.N);
#   # auto& G = g;
#   # for (uint32_t col = 0; col < MatrixSizeFactor*pars.N; ++col) {
#   #     for (uint32_t row = 0; row < MatrixSizeFactor*pars.N; ++row) {
#   #         gTimesInvRows(row, col) =
#   #             G(row, site)            * rows[0][col]
#   #             + G(row, site + pars.N)   * rows[1][col];
#   #         if (OPDIM == 3) {
#   #             gTimesInvRows(row, col) +=
#   #                 G(row, site + 2*pars.N) * rows[2][col]
#   #                 + G(row, site + 3*pars.N) * rows[3][col];
#   #         }
#   #     }
#   # }
#   # g += gTimesInvRows

# end

# function perform_update_greens!(mc,i,cols,colsTimesG)
#   N = mc.l.sites
#   flv = mc.p.flv
#   g = mc.s.greens
#   site = i

#   # Compute G' = [I + (I - G)*Delta]^(-1) * G = [I + cols] * G
#   # [O(N^2)]
#   # colsTimesG = zeros(g)
#   @inbounds for c in 1:flv*N
#     @simd for r in 1:flv*N
#       colsTimesG[r,c] = cols[r,1] * g[site, c] + cols[r,2] * g[site + N, c]
#       if mc.p.opdim == 3
#         colsTimesG[r,c] += cols[r,3] * g[site + 2*N, c] + cols[r,4] * g[site + 3*N, c]
#       end
#     end
#   end
#   # TODO how to make this fast???

#   # g.+=colsTimesG;
#   @inbounds @simd for i in eachindex(g)
#     g[i] += colsTimesG[i]
#   end

#   # @show maximum(abs.(colsTimesG))
#   # return colsTimesG
#   nothing
# end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
# # Holomorphic integration on the unit disk
# In this example, we explore integration of a harmonic function:
# ```math
#   f(x,y) = \frac{x^2-y^2+1}{(x^2-y^2+1)^2+(2xy+1)^2},
# ```
# over the unit disk. In this case, we know from complex analysis that the
# integral of a holomorphic function is equal to $\pi \times f(0,0)$.
# We analyze the function on an $N\times M$ tensor product grid defined by:
# ```math
# \begin{aligned}
# r_n & = \cos\left[(n+\tfrac{1}{2})\pi/2N\right],\quad{\rm for}\quad 0\le n < N,\quad{\rm and}\\
# \theta_m & = 2\pi m/M,\quad{\rm for}\quad 0\le m < M;
# \end{aligned}
# ```
# we convert the function samples to Chebyshev×Fourier coefficients using
# `plan_disk_analysis`; and finally, we transform the Chebyshev×Fourier
# coefficients to Zernike polynomial coefficients using `plan_disk2cxf`.
#
# For the storage pattern of the arrays, please consult the
# [documentation](https://MikaelSlevinsky.github.io/FastTransforms).

using FastTransforms, LinearAlgebra, Plots
const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated")
!isdir(GENFIGS) && mkdir(GENFIGS)
plotlyjs()

# Our function $f$ on the disk:
f = (x,y) -> (x^2-y^2+1)/((x^2-y^2+1)^2+(2x*y+1)^2)

# The Zernike polynomial degree:
N = 15
M = 4N-3

# The radial grid:
r = [sinpi((N-n-0.5)/(2N)) for n in 0:N-1]

# The angular grid (mod $\pi$):
θ = (0:M-1)*2/M

# On the mapped tensor product grid, our function samples are:
F = [f(r*cospi(θ), r*sinpi(θ)) for r in r, θ in θ]

# We superpose a surface plot of $f$ on top of the grid:
X = [r*cospi(θ) for r in r, θ in θ]
Y = [r*sinpi(θ) for r in r, θ in θ]
scatter3d(vec(X), vec(Y), vec(0F); markersize=0.75, markercolor=:red)
surface!(X, Y, F; legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "zernike.html"))
###```@raw html
###<object type="text/html" data="../zernike.html" style="width:100%;height:400px;"></object>
###```

# We precompute a (generalized) Zernike--Chebyshev×Fourier plan:
α, β = 0, 0
P = plan_disk2cxf(F, α, β)

# And an FFTW Chebyshev×Fourier analysis plan on the disk:
PA = plan_disk_analysis(F)

# Its Zernike coefficients are:
U = P\(PA*F)

# The Zernike coefficients are useful for integration. The integral of $f(x,y)$
# over the disk should be $\pi/2$ by harmonicity. The coefficient of $Z_{0,0}$
# multiplied by `√π` is:
U[1, 1]*sqrt(π)

# Using an orthonormal basis, the integral of $[f(x,y)]^2$ over the disk is
# approximately the square of the 2-norm of the coefficients:
norm(U)^2, π/(2*sqrt(2))*log1p(sqrt(2))

# But there's more! Next, we repeat the experiment using the Dunkl-Xu
# orthonormal polynomials supported on the rectangularized disk.
N = 2N
M = N

# We analyze the function on an $N\times M$ mapped tensor product $xy$-grid defined by:
# ```math
# \begin{aligned}
# x_n & = \cos\left(\frac{2n+1}{2N}\pi\right) = \sin\left(\frac{N-2n-1}{2N}\pi\right),\quad {\rm for} \quad 0 \le n < N,\quad{\rm and}\\
# z_m & = \cos\left(\frac{2m+1}{2M}\pi\right) = \sin\left(\frac{M-2m-1}{2M}\pi\right),\quad {\rm for} \quad 0 \le m < M,\\
# y_{n,m} & = \sqrt{1-x_n^2}z_m.
# \end{aligned}
# ```
# Slightly more accuracy can be expected by using an auxiliary array:
# ```math
#   w_n = \sin\left(\frac{2n+1}{2N}\pi\right),\quad {\rm for} \quad 0 \le n < N,
# ```
# so that $y_{n,m} = w_nz_m$.
#
# The x grid
w = [sinpi((n+0.5)/N) for n in 0:N-1]
x = [sinpi((N-2n-1)/(2N)) for n in 0:N-1]

# The z grid
z = [sinpi((M-2m-1)/(2M)) for m in 0:M-1]

# On the mapped tensor product grid, our function samples are:
F = [f(x[n], w[n]*z) for n in 1:N, z in z]

# We superpose a surface plot of $f$ on top of the grid:
X = [x for x in x, z in z]
Y = [w*z for w in w, z in z]
scatter3d(vec(X), vec(Y), vec(0F); markersize=0.75, markercolor=:green)
surface!(X, Y, F; legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "dunklxu.html"))
###```@raw html
###<object type="text/html" data="../dunklxu.html" style="width:100%;height:400px;"></object>
###```

# We precompute a Dunkl-Xu--Chebyshev plan:
P = plan_rectdisk2cheb(F, β)

# And an FFTW Chebyshev² analysis plan on the rectangularized disk:
PA = plan_rectdisk_analysis(F)

# Its Dunkl-Xu coefficients are:
U = P\(PA*F)

# The Dunkl-Xu coefficients are useful for integration. The integral of $f(x,y)$
# over the disk should be $\pi/2$ by harmonicity. The coefficient of $P_{0,0}$
# multiplied by `√π` is:
U[1, 1]*sqrt(π)

# Using an orthonormal basis, the integral of $[f(x,y)]^2$ over the disk is
# approximately the square of the 2-norm of the coefficients:
norm(U)^2, π/(2*sqrt(2))*log1p(sqrt(2))


text:
using OffsetArrays: OffsetArray
using Oceananigans.Architectures: device_event

#####
##### General halo filling functions
#####

fill_halo_regions!(::Nothing, args...) = nothing

"""
    fill_halo_regions!(fields::Union{Tuple, NamedTuple}, arch, args...)

Fill halo regions for each field in the tuple `fields` according to their boundary
conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples.
"""
function fill_halo_regions!(fields::Union{Tuple, NamedTuple}, arch, args...)

    for field in fields
        fill_halo_regions!(field, arch, args...)
    end

    return nothing
end

# Some fields have `nothing` boundary conditions, such as `FunctionField` and `ZeroField`.
fill_halo_regions!(c::OffsetArray, ::Nothing, args...; kwargs...) = nothing

"Fill halo regions in x, y, and z for a given field's data."
function fill_halo_regions!(c::OffsetArray, field_bcs, arch, grid, args...; kwargs...)

    fill_halos! = [
        fill_west_and_east_halo!,
        fill_south_and_north_halo!,
        fill_bottom_and_top_halo!,
    ]

    field_bcs_array_left = [
        field_bcs.west,
        field_bcs.south,
        field_bcs.bottom,
    ]

    field_bcs_array_right = [
        field_bcs.east,
        field_bcs.north,
        field_bcs.top,
    ]

    perm = sortperm(field_bcs_array_left, lt=fill_first)
    fill_halos! = fill_halos![perm]
    field_bcs_array_left  = field_bcs_array_left[perm]
    field_bcs_array_right = field_bcs_array_right[perm]

    for task = 1:3

        barrier = device_event(arch)

        fill_halo!  = fill_halos![task]
        bc_left     = field_bcs_array_left[task]
        bc_right    = field_bcs_array_right[task]

        events      = fill_halo!(c, bc_left, bc_right, arch, barrier, grid, args...; kwargs...)

        wait(device(arch), events)

    end

    return nothing
end

@inline validate_event(::Nothing) = NoneEvent()
@inline validate_event(event) = event

# Fallbacks split into two calls
function fill_west_and_east_halo!(c, west_bc, east_bc, args...; kwargs...)
     west_event = validate_event(fill_west_halo!(c, west_bc, args...; kwargs...))
     east_event = validate_event(fill_east_halo!(c, east_bc, args...; kwargs...))
    multi_event = MultiEvent((west_event, east_event))
    return multi_event
end

function fill_south_and_north_halo!(c, south_bc, north_bc, args...; kwargs...)
    south_event = validate_event(fill_south_halo!(c, south_bc, args...; kwargs...))
    north_event = validate_event(fill_north_halo!(c, north_bc, args...; kwargs...))
    multi_event = MultiEvent((south_event, north_event))
    return multi_event
end

function fill_bottom_and_top_halo!(c, bottom_bc, top_bc, args...; kwargs...)
    bottom_event = validate_event(fill_bottom_halo!(c, bottom_bc, args...; kwargs...))
       top_event = validate_event(fill_top_halo!(c, top_bc, args...; kwargs...))
     multi_event = MultiEvent((bottom_event, top_event))
     return multi_event
end

#####
##### Halo filling order
#####

fill_first(bc1::PBC, bc2)      = false
fill_first(bc1, bc2::PBC)      = true
fill_first(bc1::PBC, bc2::PBC) = true
fill_first(bc1, bc2)           = true


text:
using AbstractPlotting:
    NoConversion,
    convert_arguments,
    conversion_trait

using IntervalSets

@testset "Conversions" begin

    # NoConversion
    struct NoConversionTestType end
    conversion_trait(::NoConversionTestType) = NoConversion()

    let nctt = NoConversionTestType(),
        ncttt = conversion_trait(nctt)
        @test convert_arguments(ncttt, 1, 2, 3) == (1, 2, 3)
    end

end

@testset "functions" begin
    x = -pi..pi
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end

    x = range(-pi, stop=pi, length=100)
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end
end