text:
"""
A struct to contain all elements required to simulate a fishery.
"""
struct SpatQSimSetup{P,F,M,R,D,S}
    init_pop::P
    fleet::F
    movement::M
    pop_dynamics::R
    domain::D
    n_years::Int
    spec::S

    function SpatQSimSetup(init_pop::P,
                           fleet::F,
                           movement::M,
                           pop_dynamics::R,
                           domain::D,
                           n_years::Int,
                           spec::S) where {
                               P<:PopState,
                               F<:Fleet,
                               M<:MovementModel,
                               R<:PopulationDynamicsModel,
                               D<:AbstractFisheryDomain,
                               S<:SpatQSimSpec}
        new{P,F,M,R,D,S}(init_pop, fleet, movement, pop_dynamics, domain, n_years, spec)
    end
end

function SpatQSimSetup(spec::SpatQSimSpec,
                       prep::SpatQSimPrep)
    SpatQSimSetup(init_pop(prep),
                  fleet(spec, prep),
                  movement(prep),
                  pop_dynamics(spec),
                  domain(spec),
                  n_years(spec),
                  spec)
end

function SpatQSimSetup(spec::SpatQSimSpec)
    prep = SpatQSimPrep(spec)
    SpatQSimSetup(spec, prep)
end

function simulate(setup::SpatQSimSetup)
    P, C = simulate(setup.init_pop,
                    setup.fleet,
                    setup.movement,
                    setup.pop_dynamics,
                    setup.domain,
                    setup.n_years)
    SpatQSimResult(P, C, setup.spec)
end

function SpatQSimSetup(spec::HabQSpec, prep::HabQPrep)
    fleet = Fleet([survey_vessel(spec, prep),
                   comm_vessel(spec, prep)],
                  [length(survey_targeting(spec)), 40_000],
                  [1, 2])

    SpatQSimSetup(init_pop(prep),
                  fleet,
                  movement(prep),
                  pop_dynamics(spec),
                  domain(prep),
                  15,
                  spec)
end


text:
using Documenter, MyPackage

makedocs(;
    modules=[MyPackage],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.umn.edu/wulab/MyPackage.jl/blob/{commit}{path}#L{line}",
    sitename="MyPackage.jl",
    authors="Xiaoping Wu",
    assets=String[],
)

deploydocs(;
    repo="github.umn.edu/wulab/MyPackage.jl",
)

ext: jl


text:
"""
    parameters2hdf5(p::Params, filename)

Save `p` to HDF5 file (e.g. `.out.h5`).
"""
function parameters2hdf5(p::Params, filename::String)
  isfile(filename) ? f = h5open(filename, "r+") : f = h5open(filename, "w")
  for i in 1:length(fieldnames(Params))
    field = fieldname(Params, i)
    field_name = string(field)
    field_value = getfield(p, field)
    field_type = typeof(field_value)

    if field_type == Bool
      # HDF5 1.8.x doesn't support Bool fields, use Int instead.
      field_value = Int(field_value)
    elseif field_type == Dates.DateTime
      field_value = string(field_value)
    elseif field_type == Set{Symbol}
      field_value = string.(field_value)
    end

    try
      if HDF5.exists(f, "params/" * field_name)
        if read(f["params/"*field_name]) != field_value
          close(f)
          error(field_name, " exists but differs from current ")
        end
      else
        f["params/" * field_name] = field_value
      end
    catch e
      close(f)
      @warn "Error in dumping parameters to hdf5: $e"
    end
  end
  close(f)
end

"""
    hdf52parameters!(p::Params, input_h5)

Load `p` from HDF5 file (e.g. `.out.h5`).
"""
function hdf52parameters!(p::Params, input_h5::String)
  fields = fieldnames(Params)
  HDF5.h5open(input_h5, "r") do f
    try
      for field_name in names(f["params"])
        field = Symbol(field_name)
        if field in fields

          value = read(f["params/$(field_name)"])
          if field_name in ["global_updates", "chkr", "Bfield", "all_checks", "resume", "edrun"] # handle Bools
            value = Bool(value)
          elseif field_name in ["walltimelimit"] # handle DateTimes
            value = Dates.DateTime(value)
          elseif field_name == "obs"
            value = Set(Symbol.(value))
          end

          setfield!(p, field, value)

        else
          @warn "HDF5 contains \"params/$(field_name)\" which is not a field of type Params! Maybe old file? Try `hdf52parameters!_old`."
        end
      end
    catch e
      error("Error in loading Params object from HDF5: ", e)
    end
  end

  p.lattice_file = "" # let deduce_remaining_parameters handle this

  deduce_remaining_parameters(p)

  nothing
end

"""
    hdf52parameters!_old(p::Params, filename)

Deprecated version of `hdf52parameters!`. Use it only for old data (created before 28.11.2017).
"""
function hdf52parameters!_old(p::Params, input_h5::String)
  @warn "DEPRECATED: Only use for old data (where we dumped `params::Dict`). Should now use `hdf52parameters!` instead."
  # READ Params from h5 file
  if input_h5[end-2:end] == "jld"
    f = jldopen(input_h5, "r")
  else
    f = HDF5.h5open(input_h5, "r")
  end

  params = Dict{Any, Any}()

  try
    for e in names(f["params"])
           params[e] = read(f["params/$e"])
    end
  catch e
    println(e)
  end

  set_parameters(p, params)
  close(f)

  return params
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
@testset "BondGraph Properties" begin
    bg = BondGraph("newBG")
    @test type(bg) == :BG
    @test name(bg) == :newBG
    @test isempty(vertices(bg))

    @test eltype(BondGraph) == AbstractNode
    @test eltype(bg) == AbstractNode

    @test isempty(edges(bg))
    @test edgetype(BondGraph) == Graphs.AbstractSimpleEdge{Integer}
    @test edgetype(bg) == Graphs.AbstractSimpleEdge{Integer}
    @test is_directed(bg)

    @test size(zero(BondGraph)) == size(BondGraph())
    @test size(zero(bg)) == size(BondGraph())
end

@testset "Adding and removing elements" begin
    c = Component(:C, :C1)
    r = Component(:R, :C1, numports = 1)
    j = EqualEffort()

    bg = BondGraph()

    @test add_vertex!(bg, c)
    b = add_edge!(bg, c, j)

    @test src(b) == vertex(c)
    @test dst(b) == vertex(j)

    add_vertex!(bg, r)
    add_vertex!(bg, j)

    @test ne(bg) == 1
    @test has_edge(bg, vertex(c), vertex(j))
    @test !has_edge(bg, vertex(c), vertex(r))

    @test nv(bg) == 3
    @test has_vertex(bg, j)
    @test has_vertex(bg, 3)
    @test !has_vertex(bg, 0)

    @test components(bg) == [c, r]
    @test junctions(bg) == [j]

    @test inneighbors(bg, vertex(c)) == []
    @test outneighbors(bg, vertex(c)) == [3]

    @test rem_edge!(bg, c, j) == b
    @test ne(bg) == 0
    @test rem_vertex!(bg, r)
    @test nv(bg) == 2
end

@testset "BondGraphNode" begin
    bg = BondGraph(:RCI)
    bgn = BondGraphNode(bg)

    @test bgn.type == :BG
    @test bgn.name == :RCI
    @test bgn.freeports == Bool[]
end

@testset "Printing" begin
    C = Component(:C)
    SS = Component(:SS, :Source)
    J0 = EqualEffort(name = :J)
    b1 = Bond(C, J0)
    b2 = Bond(J0, SS)
    port = b1.srcport
    bg = BondGraph(:newbg)
    bgn = BondGraphNode(bg)

    # repr returns the output of the 'show' function
    @test repr(C) == "C:C"
    @test repr(SS) == "SS:Source"
    @test repr(b1) == "Bond C:C ⇀ J"
    @test repr(b2) == "Bond J ⇀ SS:Source"
    @test repr(port) == "Port C:C (1)"
    @test repr(bg) == "BondGraph BG:newbg (0 Nodes, 0 Bonds)"
    @test repr(bgn) == "BG:newbg"

    add_vertex!(bg, C)
    add_vertex!(bg, SS)
    add_vertex!(bg, J0)
    add_edge!(bg, C, J0)
    add_edge!(bg, J0, SS)
    @test repr(bg) == "BondGraph BG:newbg (3 Nodes, 2 Bonds)"
end

@testset "Graphs.jl Extra Functions" begin
    c1 = Component(:C)
    c2 = Component(:R)
    c3 = Component(:I)
    j = EqualFlow()

    bg = BondGraph()

    add_vertex!(bg, c1)
    add_vertex!(bg, c2)
    add_vertex!(bg, c3)
    add_vertex!(bg, j)

    add_edge!(bg, c1, j)
    add_edge!(bg, j, c2)
    add_edge!(bg, j, c3)

    # Testing on a selection of common graph functions
    @test Δ(bg) == 3
    @test Graphs.density(bg) == 0.25
    @test Array(adjacency_matrix(bg)) == [0 0 0 1; 0 0 0 0; 0 0 0 0; 0 1 1 0]
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
#!/usr/bin/env julia
#
a=[]

try
    push!(a,1)
    println(a)
catch err
    println("Error caught.")
    showerror(STDOUT, err, backtrace());println()
end

println("Continuing after error")

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
include("../src/util.jl")
include("../src/LUMND.jl")

using .LUMND
using Test
using Random
using SHA
using MatrixDepot

Random.seed!(0)

@testset "Partition" begin
    for i = 1:50
        d = rand(2:3)
        n = d == 2 ? rand(5:20) : rand(5:10)
        A = Ad(n,d)
        N = size(A,1)
        maxLvl = rand(1:10)
        (subseps, hrch, dofs) = LUMND.mnd(A, maxLvl)
        # Check sizes
        @test length(dofs) == maxLvl
        @test length(hrch) == maxLvl
        @test all([length(dofs[i]) == 2^(maxLvl-i) for i = 1:maxLvl])
        @test all([length(hrch[i]) == 2^(maxLvl-i) for i = 1:maxLvl])
        @test length(subseps) == N
        # Check ND consistency
		for (l1, s1) in LUMND.BinTreeIt(maxLvl)
			for (l2, s2) in LUMND.BinTreeIt(maxLvl)
                d1, d2 = dofs[l1][s1], dofs[l2][s2]
                if ! are_connected(l1, s1, l2, s2, maxLvl)
                    @test norm(A[d1,d2]) == 0.0
                    @test norm(A[d2,d1]) == 0.0
                end
            end
		end
        # Check hierarchy consistency
		for (l, s) in LUMND.BinTreeIt(maxLvl)
            h = hrch[l][s]
            # One root
            @test length(h.nit[1]) == 1
            # Right depth
            @test length(h.nit) == l
            ss = subseps[dofs[l][s]]
            @test all(sort(unique(ss)) .== 1:length(h.nit[end])) # Will pass if ss = []
            for l_ = 1:length(h.nit)
            	# At least one node per level, all numbered 1 ... n
                l__ = length(h.nit[l_])
                @assert l__ > 0
                @test all(1:l__ .== [n.i for n in h.nit[l_]])
            	# At least one children, except at leaves
                for n in h.nit[l_]
                    if l_ < length(h.nit)
                        @assert length(n.c) > 0
                    else
                        @test length(n.c) == 0
                    end
                end
            end
        end
    end
end

function checkAccuracy(A, lvl)
    b = rand(size(A, 1))
    tree = LUMND.factorize(A, lvl)
    x = LUMND.solve(tree, b)
    xref = A\b
    if issymmetric(A)
        k = cond(A, 1)
    else # Crazy slow otherwise
        k = 100
    end
    resref = norm(A*xref-b)/norm(b)
    @test(norm(A*x-b)/norm(b)     < 10 * resref)
    @test(norm(x-xref)/norm(xref) < k  * 1e-12 )
    return b, k, x, xref
end

@testset "Solve-Smallworld-Gen" begin
    for n in 2:20:500
        A = matrixdepot("smallworld", n)
        leaf = rand([1,2,5])
        N = size(A,1)
        A = A + sparse(UniformScaling(10.0), N, N) # Make it non-singular
        (I,J,V) = findnz(A)
        sig = rand([1e-1, 1e-2, 1e-3, 1e-4]) # Make it non-symmetric
        V += sig .* randn(size(V)) .* V
        A = sparse(I,J,V,N,N)
        lvl = Int64(round(log(N / (leaf)) / log(2) + 1))
        (b, k, x, xref) = checkAccuracy(A, lvl)
        @printf("Smallworld | %3d -> %4.2e %4.2e\n", N, norm(A*x-b)/norm(b), norm(x-xref)/norm(xref))
    end
end

@testset "Solve-Wathen-SPD" begin
    for n in 1:5:30
        A = matrixdepot("wathen", n)
        leaf = rand([1,2,5])
        N = size(A,1)
        lvl = Int64(round(log(N / (leaf)) / log(2) + 1))
        (b, k, x, xref) = checkAccuracy(A, lvl)
        @printf("Wathen | %5d -> %4.2e %4.2e\n", N, norm(A*x-b)/norm(b), norm(x-xref)/norm(xref))
    end
end

@testset "Solve-Lapl-SPD" begin
    @printf("  # |     N lvl         k ->   |Ax-b|   |x-x*| hash\n")
    for i = 1:100
        d = rand(2:3)
        n = d == 2 ? rand(5:30) : rand(5:10)
        leaf = rand([1,2,5,10,15,20])
        lvl = Int64(round(log(n^d / (leaf)) / log(2) + 1))
        A = Ad(n,d)
        (b, k, x, xref) = checkAccuracy(A, lvl)
        @printf("%3d | %5d  %2d  %2.2e -> %4.2e %4.2e %d\n", i, size(A,1), lvl, k, norm(A*x-b)/norm(b), norm(xref - x)/norm(xref), hash(x))
    end
end

@testset "Solve-Lapl-Gen" begin
    @printf("  # |     N lvl         k ->   |Ax-b|   |x-x*| hash\n")
    for i = 1:100
        d = rand(2:3)
        n = d == 2 ? rand(5:30) : rand(5:10)
        leaf = rand([1,2,5,10,15,20])
        lvl = Int64(ceil(log(n^d / (leaf)) / log(2)))
        A = Ad(n,d)
        N = size(A,1)
        (I,J,V) = findnz(A)
        sig = rand([1e-1, 1e-2, 1e-3, 1e-4]) # Make it non-symmetric
        V += sig .* randn(size(V)) .* V
        A = sparse(I,J,V,N,N)
        (b, k, x, xref) = checkAccuracy(A, lvl)
        @printf("%3d | %5d  %2d  %2.2e -> %4.2e %4.2e %d\n", i, size(A,1), lvl, k, norm(A*x-b)/norm(b), norm(xref - x)/norm(xref), hash(x))
    end
end

@testset "Solve-Blur-SPD" begin
    for n in 3:5:53
        A = matrixdepot("blur", n)
        leaf = rand([1,2,5])
        N = size(A,1)
        lvl = Int64(round(log(N / (leaf)) / log(2)))
        (b, k, x, xref) = checkAccuracy(A, lvl)
        @printf("Blur | %2d -> %4.2e %4.2e\n", n, norm(A*x-b)/norm(b), norm(x-xref)/norm(xref))
    end
end


text:
@testset "Move" begin
    g = Grid("maps/empty.yaml")
    r = Robot()
    r.x = 3
    r.y = 3
    move!(r, g)
    @test r.x == 3
    @test r.y == 2
    move!(r, g; forward=false)
    @test r.x == 3
    @test r.y == 3
end

@testset "Rotate right" begin
    g = Grid("maps/empty.yaml")
    r = Robot()
    r.x = 2
    r.y = 2
    turn_right!(r)
    move!(r, g)
    @test r.x == 3
    @test r.y == 2
    turn_right!(r)
    move!(r, g)
    @test r.x == 3
    @test r.y == 3
    turn_right!(r)
    move!(r, g)
    @test r.x == 2
    @test r.y == 3
    turn_right!(r)
    move!(r, g)
    @test r.x == 2
    @test r.y == 2
end

@testset "Rotate left" begin
    g = Grid("maps/empty.yaml")
    r = Robot()
    r.x = 2
    r.y = 2
    turn_left!(r)
    move!(r, g)
    @test r.x == 1
    @test r.y == 2
    turn_left!(r)
    move!(r, g)
    @test r.x == 1
    @test r.y == 3
    turn_left!(r)
    move!(r, g)
    @test r.x == 2
    @test r.y == 3
    turn_left!(r)
    move!(r, g)
    @test r.x == 2
    @test r.y == 2
end

@testset "Boundaries" begin
    g = Grid("maps/empty.yaml")
    r = Robot()
    move!(r, g)
    @test r.x == 1
    @test r.y == 1
    turn_right!(r)
    r.x = size(g.cells, 2)
    move!(r, g)
    @test r.x == size(g.cells, 2)
    turn_right!(r)
    r.y = size(g.cells, 1)
    move!(r, g)
    @test r.y == size(g.cells, 1)
    turn_right!(r)
    r.x = 1
    move!(r, g)
    @test r.x == 1
end

@testset "Walls" begin
    g = Grid("maps/empty.yaml")
    r = Robot()
    r.x = 1
    r.y = 2
    g.cells[1, 1].obj = 1.0
    move!(r, g)
    @test r.x == 1
    @test r.y == 2
end


text:
include(joinpath(dirname(@__FILE__),"..", "src/WebBLAST", "WebBLAST.jl"))
using WebBLAST
using BioSeq

seq = "QIKDLLVSSSTDLDTTLVLVNAIYFKGMWKTAFNAEDTREMPFHVTKQESKPVQMMCMNNSFNVATLPAE
KMKILELPFASGDLSMLVLLPDEVSDLERIEKTINFEKLTEWTNPNTMEKRRVKVYLPQMKIEEKYNLTS
VLMALGMTDLFIPSANLTGISSAESLKISQAVHGAFMELSEDGIEMAGSTGVIEDIKHSPESEQFRADHP
FLFLIKHNPTNTIVYFGRYWSP"

expected_first_result = "gi|400977284|pdb|3VVJ|A"

results = webblast(seq)

@test results[1].id == expected_first_result

results = webblast(aminoacid(seq))

@test results[1].id == expected_first_result


text:
using SimulatedAnnealing
using Test

@testset "SimulatedAnnealing.jl" begin
    # Write your tests here.
end


text:
"""
isLessOrEqual(b::AbstractBoolVar, x::AbstractIntVar, y::AbstractIntVar, trailer::Trailer)

Equivalence between a boolean variable and the inequality between variables, states that `b ⟺ x ≤ y`
"""
struct isLessOrEqual <: Constraint
    b       ::AbstractBoolVar
    x       ::AbstractIntVar
    y       ::AbstractIntVar
    active  ::StateObject{Bool}

    function isLessOrEqual(b, x, y, trailer)
        constraint = new(b, x, y, StateObject(true, trailer))
        addOnDomainChange!(b, constraint)
        addOnDomainChange!(x, constraint)
        addOnDomainChange!(y, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::isLessOrEqual, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`isLessOrEqual` propagation function.
"""
function propagate!(constraint::isLessOrEqual, toPropagate::Set{Constraint}, prunedDomains::CPModification)

    if !isbound(constraint.b)
        if maximum(constraint.x.domain) <= minimum(constraint.y.domain)
            prunedB = remove!(constraint.b.domain, false)
            if !isempty(prunedB)
                addToPrunedDomains!(prunedDomains, constraint.b, prunedB)
                triggerDomainChange!(toPropagate, constraint.b)
            end
            setValue!(constraint.active, false)
        elseif minimum(constraint.x.domain) > maximum(constraint.y.domain)
            prunedB = remove!(constraint.b.domain, true)
            if !isempty(prunedB)
                addToPrunedDomains!(prunedDomains, constraint.b, prunedB)
                triggerDomainChange!(toPropagate, constraint.b)
            end
            setValue!(constraint.active, false)
        end
        return true
    else
        if assignedValue(constraint.b)
            # propagate constraint x <= y
            if maximum(constraint.x.domain) <= minimum(constraint.y.domain)
                setValue!(constraint.active, false)
            end

            prunedX = removeAbove!(constraint.x.domain, maximum(constraint.y.domain))
            if !isempty(prunedX)
                addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
                triggerDomainChange!(toPropagate, constraint.x)
            end

            prunedY = removeBelow!(constraint.y.domain, minimum(constraint.x.domain))
            if !isempty(prunedY)
                addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
                triggerDomainChange!(toPropagate, constraint.y)
            end
            return !isempty(constraint.x.domain) && !isempty(constraint.y.domain)
        else
            # propagate constraint x > y which is y <= x -1
            if maximum(constraint.y.domain) <= minimum(constraint.x.domain) - 1
                setValue!(constraint.active, false)
            end

            prunedY = removeAbove!(constraint.y.domain, maximum(constraint.x.domain) - 1)
            if !isempty(prunedY)
                addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
                triggerDomainChange!(toPropagate, constraint.y)
            end

            prunedX = removeBelow!(constraint.x.domain, minimum(constraint.y.domain) + 1)
            if !isempty(prunedX)
                addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
                triggerDomainChange!(toPropagate, constraint.x)
            end
            return !isempty(constraint.x.domain) && !isempty(constraint.y.domain)
        end
    end
    # should not get there
    return false
end

variablesArray(constraint::isLessOrEqual) = [constraint.b, constraint.x, constraint.y]

function Base.show(io::IO, ::MIME"text/plain", con::isLessOrEqual)
    println(io, typeof(con), ": ", con.b.id, " ≡ ", con.x.id, " ≤ ", con.y.id, ", active = ", con.active)
    println(io, "   ", con.b)
    println(io, "   ", con.x)
    println(io, "   ", con.y)
end

function Base.show(io::IO, con::isLessOrEqual)
    print(io, typeof(con), ": ", con.b.id, " ≡ ", con.x.id, " ≤ ", con.y.id)
end


text:
# Operate fields at local level

"""
"""
function operate_fields(op::Function,args...)
  k = FieldOpKernel(op)
  apply_kernel_to_field(k,args...)
end

"""
    field_operation_axes(axs::Tuple...)

Returns the axes after applying a FieldOpKernel to some arrays with axes `axs...`
"""
function field_operation_axes(axs::Tuple...)
  n = maximum( map(length,axs) )
  rs = []
  for i in 1:n
    r = 1:0
    for ax in axs
      if length(ax) >= i
        if length(r) > 1
          @assert length(ax[i]) == 1 || length(ax[i]) == length(r)
        end
        if length(ax[i]) > length(r)
          r = ax[i]
        end
      end
    end
    push!(rs,r)
  end
  Tuple(rs)
end

function field_operation_metasize(axs::Tuple...)
  n = maximum( map(length,axs) )
  rs = []
  for i in 1:n
    r = 1
    for ax in axs
      if length(ax) >= i
        if ax[i] == (:)
          r = ax[i]
        end
      end
    end
    push!(rs,r)
  end
  Tuple(rs)
end

# Operate fields at global level

"""
"""
function operate_arrays_of_fields(op::Function,args...)
  k = FieldOpKernel(op)
  apply_to_field_array(k,args...)
end

function operate_arrays_of_fields(::Type{T},op::Function,args...) where T
  k = FieldOpKernel(op)
  apply_to_field_array(T,k,args...)
end

# Pre-define some operations
# some of them only make sense for fields, not for arrays of fields.

for op in (:+,:-,:tr, :transpose, :adjoint, :symmetric_part)
  @eval begin

    function ($op)(f::Field)
      operate_fields($op,f)
    end

    #function ($op)(f::AbstractArray{<:Field})
    #  operate_arrays_of_fields($op,f)
    #end

  end
end

for op in (:+,:-,:*,:inner,:outer,:dot)
  @eval begin

    function ($op)(f::Field,g::Field)
      operate_fields($op,f,g)
    end

    #function ($op)(f::AbstractArray{<:Field},g::AbstractArray{<:Field})
    #  operate_arrays_of_fields($op,f,g)
    #end

  end
end

# Operations between field values

struct FieldOpKernel{T} <: Kernel
  op::T
end

# FieldOpKernel does broadcast by default. It will work always assuming that the trial bases have shape (np,1,ndof)
# but perhaps inefficient for blocked matrices.
# In any case, optimizations for block matrices will be done at the global level (for all cells)
# instead of at the cell level in this kernel.
# In other words, we can assume that this kernel receives standard non-blocked arrays in practice.

function kernel_cache(k::FieldOpKernel,args...)
  bk = bcast(k.op)
  kernel_cache(bk,args...)
end

@inline function apply_kernel!(cache,k::FieldOpKernel,args...)
  bk = bcast(k.op)
  apply_kernel!(cache,bk,args...)
end

# Define gradients at local and global level

function apply_kernel_gradient(k::FieldOpKernel,args...)
  @notimplemented "The gradient of the result of operation $(k.op) is not yet implemented."
end

function apply_gradient(k::Valued{<:FieldOpKernel},args...)
  @notimplemented "The gradient of the result of operation $(k.op) is not yet implemented."
end

for op in (:+,:-)
  @eval begin

    # Local level

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},a,b)
      ga = field_gradient(a)
      gb = field_gradient(b)
      apply_kernel_to_field(k,ga,gb)
    end

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},a)
      ga = field_gradient(a)
      apply_kernel_to_field(k,ga)
    end

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},a::Number,b)
      gb = field_gradient(b)
      apply_kernel_to_field(k,gb)
    end

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},b,a::Number)
      gb = field_gradient(b)
      gb
    end

    # Global level

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},a,b)
      ga = field_array_gradient(a)
      gb = field_array_gradient(b)
      apply(k,ga,gb)
    end

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},a)
      ga = field_array_gradient(a)
      apply(k,ga)
    end

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},a::Number,b)
      gb = field_array_gradient(b)
      apply(k,gb)
    end

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},b,a::Number)
      gb = field_array_gradient(b)
      gb
    end

  end
end

for op in (:*,⋅)
  @eval begin

    # Local level

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},a::Number,b)
      gb = field_gradient(b)
      apply_kernel_to_field(k,a,gb)
    end

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},b,a::Number)
      gb = field_gradient(b)
      apply_kernel_to_field(k,gb,a)
    end

    function apply_kernel_gradient(k::FieldOpKernel{typeof($op)},a,b)
      ga = field_gradient(a)
      gb = field_gradient(b)
      f1 = apply_kernel_to_field(k,ga,b)
      f2 = apply_kernel_to_field(k,a,gb)
      operate_fields(+,f1,f2)
    end

    # Global level

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},a::AbstractArray{<:Number},b)
      gb = field_array_gradient(b)
      apply(k,a,gb)
    end

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},b,a::AbstractArray{<:Number})
      gb = field_array_gradient(b)
      apply(k,gb,a)
    end

    function apply_gradient(k::Valued{FieldOpKernel{typeof($op)}},a,b)
      ga = field_array_gradient(a)
      gb = field_array_gradient(b)
      f1 = apply(k,ga,b)
      f2 = apply(k,a,gb)
      operate_arrays_of_fields(+,f1,f2)
    end

  end
end

# Local level

function apply_kernel_gradient(k::FieldOpKernel{typeof(/)},b,a::Number)
  gb = field_gradient(b)
  apply_kernel_to_field(k,gb,a)
end

# Global level

function apply_gradient(k::Valued{FieldOpKernel{typeof(/)}},b,a::Number)
  gb = field_array_gradient(b)
  apply(k,gb,a)
end

# Move the value of a test basis into "trial" state

"""
"""
function trialize_basis(f)
  apply_kernel_to_field(trialize_basis_value,f)
end

"""
"""
function trialize_array_of_bases(af)
  apply_to_field_array(trialize_basis_value,af)
end

@inline function trialize_basis_value(a::AbstractMatrix)
  TrializedMatrix(a)
end

function apply_kernel_gradient(::typeof(trialize_basis_value),a)
  g = field_gradient(a)
  trialize_basis(g)
end

struct TrializedMatrix{T,A} <: AbstractArray{T,3}
  matrix::A
  @inline function TrializedMatrix(matrix::AbstractMatrix{T}) where T
    A = typeof(matrix)
    new{T,A}(matrix)
  end
end

TrializedMatrix{T,A}(u::UndefInitializer,s::Tuple) where {T,A} = TrializedMatrix(A(u,(s[1],s[3])))

Base.size(a::TrializedMatrix) = (size(a.matrix,1),1,size(a.matrix,2))

Base.IndexStyle(::Type{<:TrializedMatrix{T,A}}) where {T,A} = IndexStyle(A)

@inline Base.getindex(a::TrializedMatrix,i::Integer,j::Integer,k::Integer) = a.matrix[i,k]

@inline Base.getindex(a::TrializedMatrix,i::Integer) = a.matrix[i]

@inline Base.setindex!(a::TrializedMatrix,v,i::Integer,j::Integer,k::Integer) = (a.matrix[i,k] = v)

@inline Base.setindex!(a::TrializedMatrix,v,i::Integer) = (a.matrix[i] = v)

# Optimizations for block matrices at global level (for all cells)

function apply(f::typeof(trialize_basis_value),a::VectorOfBlockMatrixCoo)
  blocks = map(b->apply(f,b),a.blocks)
  blockids = broadcast(ij->(ij[1],1,ij[2]),a.blockids)
  axs = apply( ax -> (ax[1],blockedrange([1]),ax[2]), a.axes)
  VectorOfBlockArrayCoo(blocks,blockids,axs)
end

# Unary operations
# Assumption: op is linear wrt a
function apply(k::FieldOpKernel,a::VectorOfBlockArrayCoo)
  blocks = map(b->apply(k,b), a.blocks)
  VectorOfBlockArrayCoo(blocks,a.blockids,a.axes,a.ptrs)
end

# Binary test/field or trial/field
# Assumption: op is linear wrt a
function apply(k::FieldOpKernel,a::VectorOfBlockArrayCoo,f::AbstractArray{<:AbstractVector})
  blocks = map(b->apply(k,b,f), a.blocks)
  VectorOfBlockArrayCoo(blocks,a.blockids,a.axes,a.ptrs)
end

function apply(k::FieldOpKernel,a::VectorOfBlockArrayCoo,f::AbstractArray{<:Number})
  blocks = map(b->apply(k,b,f), a.blocks)
  VectorOfBlockArrayCoo(blocks,a.blockids,a.axes,a.ptrs)
end

# Binary field/test or field/trial
# Assumption: op is linear wrt a
function apply(k::FieldOpKernel,f::AbstractArray{<:AbstractVector},a::VectorOfBlockArrayCoo)
  blocks = map(b->apply(k,f,b), a.blocks)
  VectorOfBlockArrayCoo(blocks,a.blockids,a.axes,a.ptrs)
end

function apply(k::FieldOpKernel,f::AbstractArray{<:Number},a::VectorOfBlockArrayCoo)
  blocks = map(b->apply(k,f,b), a.blocks)
  VectorOfBlockArrayCoo(blocks,a.blockids,a.axes,a.ptrs)
end

# Binary test/test or trial/trial
# Assumption: op is a linear combination of a and b
function apply(
  k::FieldOpKernel,a::VectorOfBlockArrayCoo{Ta,N} where Ta,b::VectorOfBlockArrayCoo{Tb,N} where Tb) where N
  @assert size(a.ptrs) == size(b.ptrs)
  blocks = []
  blockids = NTuple{N,Int}[]
  for (I,aI) in enumerateblocks(a)
    bI = b[I]
    if is_nonzero_block(a,I) || is_nonzero_block(b,I)
      block = apply(k,aI,bI)
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,a.axes)
end

# Binary + test/test or trial/trial
function apply(
  k::FieldOpKernel{typeof(+)},
  a::VectorOfBlockArrayCoo{Ta,N} where Ta,
  b::VectorOfBlockArrayCoo{Tb,N} where Tb) where N

  @assert size(a.ptrs) == size(b.ptrs)
  blocks = []
  blockids = NTuple{N,Int}[]
  for (I,aI) in enumerateblocks(a)
    bI = b[I]
    if is_nonzero_block(a,I) && is_nonzero_block(b,I)
      block = apply(k,aI,bI)
      push!(blocks,block)
      push!(blockids,I.n)
    elseif is_nonzero_block(a,I)
      block = aI
      push!(blocks,block)
      push!(blockids,I.n)
    elseif is_nonzero_block(b,I)
      block = bI
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,a.axes)
end

# Binary - test/test or trial/trial
function apply(
  k::FieldOpKernel{typeof(-)},
  a::VectorOfBlockArrayCoo{Ta,N} where Ta,
  b::VectorOfBlockArrayCoo{Tb,N} where Tb) where N

  @assert size(a.ptrs) == size(b.ptrs)
  blocks = []
  blockids = NTuple{N,Int}[]
  for (I,aI) in enumerateblocks(a)
    bI = b[I]
    if is_nonzero_block(a,I) && is_nonzero_block(b,I)
      block = apply(k,aI,bI)
      push!(blocks,block)
      push!(blockids,I.n)
    elseif is_nonzero_block(a,I)
      block = aI
      push!(blocks,block)
      push!(blockids,I.n)
    elseif is_nonzero_block(b,I)
      block = apply(k,bI)
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,a.axes)
end

# Binary test/trial
# Assumption: op is a product of a and b
function apply(
  k::FieldOpKernel,a::VectorOfBlockMatrixCoo,b::VectorOfBlockArrayCoo{Tb,3} where Tb)
  axs = apply( (a1,a2) -> (a1[1],a1[2],a2[3]) ,a.axes,b.axes)
  blocks = []
  blockids = NTuple{3,Int}[]
  nfield1 = size(a.ptrs,2)
  nfield2 = size(b.ptrs,3)
  for f1 in 1:nfield1
    I1 = Block(1,f1)
    for f2 in 1:nfield2
      I2 = Block(1,1,f2)
      if is_nonzero_block(a,I1) && is_nonzero_block(b,I2)
        block = apply(k,a[I1],b[I2])
        push!(blocks,block)
        push!(blockids,(1,f1,f2))
      end
    end
  end
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,axs)
end

# Binary trial/test
# Assumption: op is a product of a and b
function apply(
  k::FieldOpKernel,a::VectorOfBlockArrayCoo{Tb,3} where Tb,b::VectorOfBlockMatrixCoo)
  axs = apply( (a1,a2) -> (a1[1],a2[2],a1[3]) ,a.axes,b.axes)
  blocks = []
  blockids = NTuple{3,Int}[]
  nfield1 = size(b.ptrs,2)
  nfield2 = size(a.ptrs,3)
  for f1 in 1:nfield1
    I1 = Block(1,f1)
    for f2 in 1:nfield2
      I2 = Block(1,1,f2)
      if is_nonzero_block(b,I1) && is_nonzero_block(a,I2)
        block = apply(k,b[I1],a[I2])
        push!(blocks,block)
        push!(blockids,(1,f1,f2))
      end
    end
  end
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,axs)
end

# General operation
# Assumption only one block per direction. This assumption is OK since it is the only
# efficient situation.
function apply(
  k::FieldOpKernel,
  b::VectorOfBlockArrayCoo,
  c::VectorOfBlockArrayCoo,
  d::Union{VectorOfBlockArrayCoo,AbstractArray{<:AbstractVector}}...)
  a = (b,c,d...)
  I = _general_op_block_id(a...)
  aI = map(_get_first_block,a)
  blocks = (apply(k,aI...),)
  blockids = [I,]
  axs = _general_op_axes(a...)
  VectorOfBlockArrayCoo(blocks,blockids,axs)
end

function _general_op_block_id(a...)
  I = zeros(Int,3)
  I[1] = 1
  for ai in a
    i = _get_first_block_id(ai)
    if I[length(i)] == 0
        I[length(i)] = i[end]
      else
      @notimplementedif I[length(i)] != i[end]
    end
  end
  if I[3] == 0
    (I[1],I[2])
  else
    if I[2] == 0
      I[2] = 1
    end
    Tuple(I)
  end
end

_get_first_block_id(a) = (1,)

function _get_first_block_id(a::VectorOfBlockArrayCoo)
  @notimplementedif length(a.blocks) != 1
  first(a.blockids)
end

_get_first_block(a) = a
_get_first_block(a::VectorOfBlockArrayCoo) = first(a.blocks)

function  _general_op_axes(a...)
  test = nothing
  for ai in a
    if isa(ai,VectorOfBlockArrayCoo{T,2} where T)
      test = ai
      break
    end
  end
  trial = nothing
  for ai in a
    if isa(ai,VectorOfBlockArrayCoo{T,3} where T)
      trial = ai
      break
    end
  end
  if test != nothing && trial != nothing
    axs = apply( (a1,a2) -> (a1[1],a1[2],a2[3]),test.axes,trial.axes)
  elseif test == nothing && trial != nothing
    axs = trial.axes
  elseif test != nothing && trial == nothing
    axs = test.axes
  else
    @unreachable
  end
  axs
end

# Integration of elem vectors
function apply(k::IntKernel,f::VectorOfBlockArrayCoo{T,2} where T,w::AbstractArray,j::AbstractArray)
  ax = apply(a->(a[2],),f.axes)
  blocks = map(block->apply(k,block,w,j),f.blocks)
  blockids = [ (ids[2],) for ids in f.blockids ]
  VectorOfBlockArrayCoo(blocks,blockids,ax)
end

# Integration of elem matrices
function apply(k::IntKernel,f::VectorOfBlockArrayCoo{T,3} where T,w::AbstractArray,j::AbstractArray)
  ax = apply(a->(a[2],a[3]),f.axes)
  blocks = map(block->apply(k,block,w,j),f.blocks)
  blockids = [ (ids[2], ids[3]) for ids in f.blockids ]
  VectorOfBlockArrayCoo(blocks,blockids,ax)
end


text:
module CUTENSOR

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cudaDataType
using ..CUDA: libcutensor,  @retry_reclaim

using CEnum

using Memoize

using DataStructures

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

# cache for created, but unused handles
const handle_cache_lock = ReentrantLock()
const idle_handles = DefaultDict{CuContext,Vector{Base.RefValue{cutensorHandle_t}}}(()->Base.RefValue{cutensorHandle_t}[])

function handle()
    ctx = context()
    get!(task_local_storage(), (:CUTENSOR, ctx)) do
        handle = @lock handle_cache_lock begin
            if isempty(idle_handles[ctx])
                handle = Ref{cutensorHandle_t}()
                cutensorInit(handle)
                handle
            else
                pop!(idle_handles[ctx])
            end
        end

        finalizer(current_task()) do task
            @spinlock handle_cache_lock begin
                push!(idle_handles[ctx], handle)
            end
        end
        # TODO: destroy to preserve memory, or at exit?

        handle
    end::Base.RefValue{cutensorHandle_t}
end

end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

const AGF_DIR = joinpath(@__DIR__, "downloads", "glasscat") # contains SCHOTT.agf, Sumita.agf, etc.
const GLASSCAT_DIR = joinpath(@__DIR__, "..", "src", "GlassCat") # contains GlassCat.jl (pre-existing)
const JL_DIR = joinpath(GLASSCAT_DIR, "data") # contains AGFGlasscat.jl, SCHOTT.jl, etc.

const SOURCES_PATH = joinpath(@__DIR__, "sources.txt")
const AGFGLASSCAT_PATH = joinpath(JL_DIR, "AGFGlassCat.jl")

include(joinpath(GLASSCAT_DIR, "GlassTypes.jl"))
include("sources.jl")
include("generate.jl")

mkpath(AGF_DIR)
mkpath(JL_DIR)

# Build/verify a source directory using information from sources.txt
sources = [split(line, " ") for line in readlines(SOURCES_PATH)]
verify_sources!(sources, AGF_DIR)
verified_source_names = [source[1] for source in sources]

# Use verified sources to generate required .jl files
@info "$(isfile(AGFGLASSCAT_PATH) ? "Re-g" : "G")enerating $AGFGLASSCAT_PATH"
@info "Using sources: $(join(verified_source_names, ", ", " and "))"
generate_jls(verified_source_names, AGFGLASSCAT_PATH, JL_DIR, AGF_DIR)


text:
using RoME
using Test
using Statistics

## MWE Pose2Point2 from #388
@testset "basic Pose2Point2 test" begin

    fg = initfg()

    addVariable!(fg, :x1, Pose2)
    addVariable!(fg, :l1, Point2)

    addFactor!(fg, [:x1], PriorPose2(MvNormal([0.,0, 0], [0.01, 0.01, 0.01])))

    addFactor!(fg, [:x1; :l1], Pose2Point2(MvNormal([0.0,-1], [0.1,0.1])))

    ensureAllInitialized!(fg)

    tree, smt, hist = solveTree!(fg)

    @test isapprox(mean(getVal(fg, :x1),dims=2), [0,0,0], atol = 0.05)
    @test isapprox(mean(getVal(fg, :l1),dims=2), [0,-1], atol = 0.05)
end


text:
struct ProblemParSetter{NS,NP,NO,IT}
    optinfo::NTuple{NO, Tuple{Symbol,Symbol,IT}}
    statemap::NTuple{NS, IT} # position_in_u0 -> posopt, 0: not optimized
    parmap::NTuple{NP, IT} # position_in_p -> posopt, 0: not optimized
    statesyms::NTuple{NS, Symbol}
    parsyms::NTuple{NP, Symbol}
    # isuopt::NTuple{NP, Bool}
    # ispopt::NTuple{NP, Bool}
end

"""
    ProblemParSetter(state_names,par_names,popt_names)
    ProblemParSetter(sys::ODESystem, popt_names; strip=false)

Helps keeping track of a subset of initial states and paraemters to be optimized.

# Arguments
- `state_names`: Tuple or AbstractVector of all the initial states of the problem
- `par_names`: all the parameters of the problem
- `popt_names`: the parameter/initial states to be optimized.

The states and parameters can be extracted from an `ModelingToolkit.ODESystem`.
If `strip=true`, then namespaces of parameteres of a composed system are removed,
e.g. `subcomp₊p` becomes `p`.
"""
function ProblemParSetter(
    state_syms::NTuple{NS,Symbol},par_syms::NTuple{NP,Symbol},popt_syms::NTuple{NO,Symbol},
    ::Val{allow_missing_popt}=Val(true); it::Type{IT}=Int64) where {NO, NS,NP,IT,allow_missing_popt}
    posu_orig = Dict(parname => convert(IT,pos) for (pos, parname) in enumerate(state_syms))
    posp_orig = Dict(parname => convert(IT,pos) for (pos, parname) in enumerate(par_syms))
    #
    optinfo_u = Tuple((:state, statesym, posu_orig[statesym]) for
        statesym in popt_syms if haskey(posu_orig, statesym))
    optinfo_p = Tuple((:par, parsym, posp_orig[parsym]) for
        parsym in popt_syms if haskey(posp_orig, parsym))
    optinfo = (optinfo_u..., optinfo_p...)
    #
    map_posu_to_posopt = Dict(pos => posopt for
        (posopt, (source, parname, pos)) in enumerate(optinfo) if source == :state)
    map_posp_to_posopt = Dict(pos => posopt for
        (posopt, (source, parname, pos)) in enumerate(optinfo) if source == :par)
        #
    statemap = NTuple{NS}((get(map_posu_to_posopt,posu,0) for posu in one(IT):NS))
    parmap = NTuple{NP}((get(map_posp_to_posopt,posp,0) for posp in one(IT):NP))
    # isuopt = NTuple{NS}((i != 0 for i in statemap))
    # ispopt = NTuple{NP}((i != 0 for i in parmap))
    # @show NS, NP, NO
    # @show optinfo
    # @show statemap
    # @show parmap
    popt_syms_in = popt_syms
    if allow_missing_popt && (NO != length(optinfo))
        optinfo_psyms = getindex.(optinfo,2)
        miss = setdiff(popt_syms, optinfo_psyms)
        @warn(
            "missing optimization parameters in system: " * join(miss,", "))
        popt_syms_in = filter((x -> !isnothing(findfirst(==(x), optinfo_psyms))), popt_syms)
        ps = ProblemParSetter{NS,NP,length(optinfo),IT}(
            optinfo, statemap, parmap, state_syms, par_syms,)
    else
        # note using NO here for type stability
        ps = ProblemParSetter{NS,NP,NO,IT}(
            optinfo, statemap, parmap, state_syms, par_syms,)
    end
    # make sure that states go first: smalles p position > highest u position
    paroptsyms(ps) == popt_syms_in || error(
        "Expected states to optimize before parameters to optimize. But got $popt_syms_in.")
    ps
end

function ProblemParSetter(state_names,par_names,popt_names)
    # not type-stable
    ProblemParSetter(
        Tuple(i for i in symbol.(state_names)),
        Tuple(i for i in symbol.(par_names)),
        Tuple(i for i in symbol.(popt_names))
    )
end

function ProblemParSetter(sys::ODESystem,popt_names; strip=false)
    state_names = strip ? strip_namespace.(statesyms(sys)) : statesyms(sys)
    par_names = strip ? strip_namespace.(parsyms(sys)) : parsyms(sys)
    ProblemParSetter(state_names, par_names, popt_names)
end

"""
    count_states(::ProblemParSetter)
    count_par(::ProblemParSetter)
    count_paropt(::ProblemParSetter)

Report the number of problem states, problem parameters and optimized parameters
respectively.
"""
function count_states(::ProblemParSetter{NS}) where {NS}; NS; end,
function count_par(::ProblemParSetter{NS,NP}) where {NS,NP}; NP; end,
function count_paropt(::ProblemParSetter{NS,NP,NO}) where {NS,NP,NO}; NO; end

"""
    statesyms(ps::ProblemParSetter)
    parsyms(ps::ProblemParSetter)
    paroptsyms(ps::ProblemParSetter)

Report the names, i.e. symbols of problem states, problem parameters and
optimized parameters respectively.
"""
function statesyms(ps::ProblemParSetter); ps.statesyms; end,
function parsyms(ps::ProblemParSetter); ps.parsyms; end,
function paroptsyms(ps::ProblemParSetter)
    first.(Base.Iterators.drop.(ps.optinfo,1))
end

"""
    prob_new = update_statepar(ps::ProblemParSetter, popt, prob::ODEProblem)
    u0new, pnew = update_statepar(ps::ProblemParSetter, popt, u0, p)

Return an updated problem or updates states and parameters where
values corresponding to positions in `popt` are set.
"""
function update_statepar(ps::ProblemParSetter, popt, prob::ODEProblem)
    u0,p = update_statepar(ps, popt, prob.u0, prob.p)
    remake(prob; u0, p)
end

function update_statepar(ps::ProblemParSetter, popt, u0::TU, p::TP) where {TU,TP}
    @assert length(popt) == count_paropt(ps) "expected $(count_paropt(ps)) parameters for "
    # @show TU
    # u2 = TU((ps.statemap[i] == 0 ? u0[i] : popt[ps.statemap[i]] for i in 1:length(u0)))
    # p2 = TP((ps.parmap[i] == 0 ? p[i] : popt[ps.parmap[i]] for i in 1:length(p)))
    # (u2,p2)
    # take care that popt might be of different type, like FowwardDiff
    # need to convert to new type
    u00 = map(x -> x * zero(eltype(popt)), u0)
    u0g = (ps.statemap[i] == 0 ? u0[i] : popt[ps.statemap[i]] for i in 1:length(u0))
    #u0new = typed_from_generator(TU, u00, u0g)
    u0new = typed_from_generator(u00, u0g)
    #u0new = typeof(u00)(u0g)
    #
    p0 = map(x -> x * zero(eltype(popt)), p)
    # pg = (ps.parmap[i] == 0 ? p[i] : popt[ps.parmap[i]] for i in 1:length(p))
    # pnew = convert(typeof(p0), p0 .+ pg)::typeof(p0)
    pg = (ps.parmap[i] == 0 ? p[i] : popt[ps.parmap[i]] for i in 1:length(p))
    #pnew = typed_from_generator(TP, p0, pg)
    pnew = typed_from_generator(p0, pg)
    #pnew = typeof(p0)(pg)
    #
    (u0new, pnew)
end

# function typed_from_generator(type::Type, v0, vgen) where T
#     if type <: AbstractArray
#     end
#     typeof(v0)(vgen)
# end
# # AbstractVector{T} is not more specific than ::Type, need to support all used concrete
# #typed_from_generator(::Type{AbstractVector{T}}, v0, vgen) where T = convert(typeof(v0), v0 .+ vgen)::typeof(v0)
# typed_from_generator(::Type{Vector{T}}, v0, vgen) where T = convert(typeof(v0), v0 .+ vgen)::typeof(v0)

typed_from_generator(v0, vgen) = typeof(v0)(vgen)
function typed_from_generator(v0::AbstractVector, vgen)
    # convert to std vector, because typeof(v0) does not contain all entries
    T = Vector{eltype(v0)}
    convert(T, v0 .+ vgen)::T
end

"""
    get_paropt(ps::ProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt(ps::ProblemParSetter, u0, p)

    get_paropt_labeled(ps::ProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt_labeled(ps::ProblemParSetter, u0, p)

Extract the initial states and parameters corresponding to the positions
that are optimized.
If both u0 and p are AbstractVectors, the result is a Vector, otherwise the result is a Tuple.

The _lebeled versions additionally call `label_paropt` (see [`label_state`](@ref))
on the return value.
"""
function get_paropt(ps::ProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt(ps, prob.u0, prob.p; kwargs...)
end
function get_paropt_labeled(ps::ProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt_labeled(ps, prob.u0, prob.p; kwargs...)
end

function get_paropt(ps::ProblemParSetter, u0::AbstractVector, p::AbstractVector)
    v = [(first(t) == :par) ? p[last(t)] : u0[last(t)] for t in ps.optinfo]
end

function get_paropt_labeled(ps::ProblemParSetter, u0::AbstractVector, p::AbstractVector)
    v = get_paropt(ps, u0, p)
    label_paropt(ps,v)
end

function get_paropt(ps::ProblemParSetter{NS,NP,NO}, u0, p) where {NS,NP,NO}
    t0 = NTuple{NO}(((first(t) == :par) ? p[last(t)] : u0[last(t)]
        for t in ps.optinfo))
    # need to explicitly assure full type for julia 1.6 for type stability
    t1 = t0::NTuple{NO,eltype(t0)}
end

function get_paropt_labeled(ps::ProblemParSetter{NS,NP,NO}, u0, p) where {NS,NP,NO}
    t1 = get_paropt(ps,u0,p)
    label_paropt(ps,t1)
end

"""
    label_state(ps, u::AbstractVector) = LArray{statesyms(ps)}(u)
    label_par(ps, par::AbstractVector) = LArray{parsyms(ps)}(par)
    label_paropt(ps, popt::AbstractVector) = LArray{paroptsyms(ps)}(popt)

Produce a labeled version of a sequence of initial states, parameters, or
optimized parameters respectively.
The return type differs given the input
- SVector -> SLVector
- NTuple -> NamedTuple
- AbstractVector -> LArray
"""
function label_state(ps, u::AbstractVector); LArray{statesyms(ps)}(u); end,
function label_par(ps, par::AbstractVector); LArray{parsyms(ps)}(par); end,
function label_paropt(ps, popt::AbstractVector); LArray{paroptsyms(ps)}(popt); end

label_state(ps, u::SVector) = SLVector(label_state(ps, Tuple(u)))
label_state(ps, u::NTuple) = NamedTuple{statesyms(ps)}(u)

label_par(ps, par::SVector) = SLVector(label_par(ps, Tuple(par)))
label_par(ps, par::NTuple) = NamedTuple{parsyms(ps)}(par)

label_paropt(ps, popt::SVector) = SLVector(label_paropt(ps, Tuple(popt)))
label_paropt(ps, popt::NTuple) = NamedTuple{paroptsyms(ps)}(popt)

# extends Base.merge to work on SVector
merge(x::T, y::NamedTuple) where T<:SLArray = T(merge(NamedTuple(x),y)...)

# and on Labelled Arrays
function merge(x::T, y::NamedTuple) where T<:LArray
    xnew = deepcopy(x)
    for (key,val) in pairs(y)
        xnew[key] = val
    end
    xnew
end

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 34beb7d6-fba9-11eb-15b8-e34afccc9f88
begin
	using Pkg
    Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
	using CSV, DataFrames, Distributions, Plots, MLJ, MLJLinearModels, Random, OpenML
end

# ╔═╡ 4e30f307-7b07-4cb7-a34e-ada7e2f6946d
using PlutoUI; PlutoUI.TableOfContents()

# ╔═╡ e688a9de-2dba-4fab-b4e6-9803c5361a62
begin
	using MLCourse
	MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ 4c536f7c-dc4b-4cf7-aefa-694b830a8e6a
md"# Transformations of the Input

## Splines"

# ╔═╡ 026de033-b9b2-4165-80e7-8603a386381c
function spline_features(data, degree, knots)
	q = length(knots) + degree
	names = [Symbol("H", i) for i in 1:q] # create the feature names
	features = [data.x .^ d for d in 1:degree] # compute the first features
	append!(features, [max.(0, (data.x .- k).^degree) for k in knots]) # remaining features
	DataFrame(features, names)
end

# ╔═╡ 69b9c8ed-acf1-43bf-8478-3a9bc57bc36a
md"We will use an OpenML dataset that was collected in 1985 in the US to determine the factors influencing the wage of different individuals. Here we will just fit the wage as a function of the age, i.e. we will transform the age into a spline representation and perform linear regression on that. Wages are given in USD/hour (see `OpenML.describe_dataset(534)` for more details)."

# ╔═╡ 6532d7c4-2385-401a-9e99-8dc86655687f
wage = OpenML.load(534) |> DataFrame

# ╔═╡ 001c88b3-c0db-4ad4-9c0f-04a2ea7b090d
md"spline degree = $(@bind spline_degree Slider(1:5, show_value = true))

knot 1 = $(@bind knot1 Slider(10:50, default = 30, show_value = true))

knot 2 = $(@bind knot2 Slider(20:60, default = 40, show_value = true))

knot 3 = $(@bind knot3 Slider(30:70, default = 50, show_value = true))

knot 4 = $(@bind knot4 Slider(40:80, default = 60, show_value = true))"

# ╔═╡ 2e1c3fa4-e424-4844-82c1-53d02271e7a6
begin
	# Here we compute the spline representation of the AGE.
	# the values knot1, knot2, knot3, knot4 and spline_degree
	# are defined by the sliders below.
	knots = [knot1, knot2, knot3, knot4]
    spline_input = spline_features(DataFrame(x = wage.AGE), spline_degree, knots)
end

# ╔═╡ d35acc45-7404-4b47-98fe-6e13db8b657e
spline_fit = fit!(machine(LinearRegressor(), spline_input, wage.WAGE));

# ╔═╡ 5774b82f-a09e-46f8-a458-0f85fc4c53d8
fitted_params(spline_fit)

# ╔═╡ 0815ba19-5c7d-40d6-b948-8edccfb5c386
begin
	scatter(wage.AGE, wage.WAGE)
	test_data = spline_features(DataFrame(x = minimum(wage.AGE):maximum(wage.AGE)),
		                        spline_degree, knots)
	predicted_avg_wage = predict(spline_fit, test_data)
	plot!(minimum(wage.AGE):maximum(wage.AGE), predicted_avg_wage, w = 3)
	vline!(knots, linestyle = :dash, legend = false)
end

# ╔═╡ 004958cf-0abb-4d8e-90b1-5966412cba91
md"## Categorical Predictors"

# ╔═╡ b43365c6-384e-4fde-95f3-b2ad1ebc421a
cdata = DataFrame(gender = ["male", "male", "female", "female", "female"],
                  treatment = [1, 2, 2, 1, 3])

# ╔═╡ a91c6298-98c4-4dc6-9450-279da1c34146
coerce!(cdata, :gender => Multiclass, :treatment => Multiclass)

# ╔═╡ 2755577e-8b4c-4cdc-902e-fe661aa91731
MLJ.transform(fit!(machine(OneHotEncoder(), cdata)), cdata)

# ╔═╡ 81fc52cd-e67b-4188-a028-79abcc77cb92
MLJ.transform(fit!(machine(OneHotEncoder(drop_last = true), cdata)), cdata)

# ╔═╡ e928e4f7-0ebc-4db5-b83e-0c5d22c0ff3c
md"We will now apply this transformation to fit the wage data."

# ╔═╡ 1bf42b12-5fd1-4b4c-94c6-14a52317b15f
begin
	preprocessor = machine(OneHotEncoder(drop_last = true), select(wage, Not(:WAGE))) |> fit!
	wage_input = MLJ.transform(preprocessor, select(wage, Not(:WAGE)))
end

# ╔═╡ 0a22b68c-1811-4de1-b3f8-735ee50299d2
wage_mach = machine(LinearRegressor(), wage_input, wage.WAGE) |> fit!;

# ╔═╡ 6323a352-a5f3-4ba7-93bd-903fb1bf22a0
fitted_params(wage_mach)

# ╔═╡ 85c9650d-037a-4e9f-a15c-464ad479ff44
md"From the fitted parameters we can see, for example, that years of education correlate positively with wage, women earn on average ≈2 USD per hour less than men (keeping all other factors fixed) or persons in a management position earn ≈4 USD per hour more than those in a sales position (keeping all other factors fixed)."

# ╔═╡ aa002263-4dc8-4238-ab27-a9094947600c
md"## Dealing with Missing Data

In the following artififial dataset the age and the gender is missing for some of the subjects.
"

# ╔═╡ c7dab115-b990-4760-85d0-4214d33ada5d
datam = DataFrame(age = [12, missing, 41, missing, 38, 35], gender = categorical([missing, "female", "male", "male", "male", "female"]))

# ╔═╡ e88f9e54-f86c-4d3a-80d7-1a9c5006476a
md"Removing rows with missing values is a safe way to deal with such datasets."

# ╔═╡ 89cbf454-8d1b-4c82-ae39-0672a225e7dd
dropmissing(datam)

# ╔═╡ 315f536b-4f10-4471-8e22-18c385ca80f2
md"Alternatively one can fill in the missing values with some standard values. For our artificial dataset the default `FillImputer` fills in the median of all age values rounded to the next integer and gender \"male\" because this dataset contains more \"male\"s than \"female\"s."

# ╔═╡ 69bc41cb-862c-4334-8b71-7f9974fedfe2
MLJ.transform(fit!(machine(FillImputer(), datam)), datam)

# ╔═╡ 76155286-978b-4a05-b428-4be4089d8b9a
md"## Standardization

Some methods (like k-nearest neighbors) are sensitive to rescaling of data. For example, in the following artificial dataset the height of people is measured in cm and their weights in kg, but one could have just as well used meters and grams or inches and pounds to report the height and the weight. By standardizing the data these arbitrary choices of the units do not matter anymore."

# ╔═╡ b2f03b4a-b200-4b30-8793-6edc02f8136d
height_weight = DataFrame(height = [165., 175, 183, 152, 171],
	                      weight = [60., 71, 89, 47, 70])

# ╔═╡ 794de3b3-df22-45c1-950d-a515c4f41409
height_weight_mach = fit!(machine(Standardizer(), height_weight));

# ╔═╡ adca81d5-2453-4954-91ba-7cc24d45e8ef
md"In the fitted parameters of the standardizing machine we see the mean and the standard deviation of the different columns of our dataset."

# ╔═╡ bf879ed5-4241-4abe-a42c-a562d49b21d2
fitted_params(height_weight_mach)

# ╔═╡ 2492bd04-cf24-4252-a37f-05fce8c1b87c
data_height_weight_s = MLJ.transform(height_weight_mach, height_weight)

# ╔═╡ 32020e6d-5ab2-4832-b565-4264ebcc82eb
md"The transformed data has now mean 0 and standard deviation 1 in every column."

# ╔═╡ fd19d841-f99e-46f4-a213-46be7a3d598d
combine(data_height_weight_s,
	    [:height, :weight] .=> mean,
	    [:height, :weight] .=> std)

# ╔═╡ 8a1a61a0-d657-4858-b8a7-d4e7d9bbf306
md"Note that standardization does not change the shape of the distribution of the data; it just shifts and scales the data, as we can see in the example below."

# ╔═╡ ee2f2a32-737e-49b5-885b-e7467a7b93f6
data = DataFrame(x = 10*rand(Beta(2, 8), 10^4) .+ 16);

# ╔═╡ cc550977-1f37-4ad1-8170-d6150cc01a9c
plot(histogram(data.x, nbins = 100))

# ╔═╡ f099be37-82d2-4031-955c-6a797c8400f2
begin
	st_mach = fit!(machine(Standardizer(), data))
	st_data = MLJ.transform(st_mach, data)
	plot(histogram(st_data.x, nbins = 100))
end

# ╔═╡ 75dda04c-c0af-4a05-b840-91f5c1aea53c
(mean(st_data.x), std(st_data.x))

# ╔═╡ a10b3fa1-358b-4a94-b927-5e0f71edd1f3
md"## Vector Features

In this section we look at a version of the famous xor problem. It is a binary classification problem with 2-dimensional input. It is called xor problem, because the class label is true if and only if the ``X_1`` and ``X_2`` coordinate of a point have different sign, which is reminiscent of the [logical exclusive or operation](https://en.wikipedia.org/wiki/Exclusive_or). The decision boundary is non-linear. By taking the scalar product with 4 different 2-dimensional vectors and setting to 0 all scalar product that would be negative, we find a 4-dimensional feature representation for which linear logistic regression can classify each point correctly."

# ╔═╡ b66c7efc-8fd8-4152-9b0f-9332f760b51b
function xor_generator(; n = 200)
	x = 2 * rand(n, 2) .- 1
	DataFrame(X1 = x[:, 1], X2 = x[:, 2],
		      y = coerce((x[:, 1] .> 0) .⊻ (x[:, 2] .> 0), Binary))
end

# ╔═╡ ebe38042-5160-469b-8c07-484b7c019063
xor_data = xor_generator()

# ╔═╡ fa7cf2ae-9777-4ba9-9fe0-df545b94c5d8
begin
	scatter(xor_data.X1, xor_data.X2, c = int.(xor_data.y) .+ 1,
            label = nothing, xlabel = "X₁", ylabel = "X₂")
	hline!([0], label = nothing, c = :black)
	vline!([0], label = nothing, c = :black, size = (400, 300))
end

# ╔═╡ 99d3ae30-6bc2-47bf-adc4-1cc1d3ff178d
begin
	lin_mach = fit!(machine(LogisticClassifier(penalty = :none),
			                select(xor_data, Not(:y)),
			                xor_data.y))
	lin_pred = predict_mode(lin_mach, select(xor_data, Not(:y)))
	(training_misclassification_rate = mean(lin_pred .!= xor_data.y),)
end

# ╔═╡ a7f12901-e155-4f61-b252-c5af6669bfee
function vector_features(data, vectors...)
	q = length(vectors) # number of vector features
	names = [Symbol("H", i) for i in 1:q]
	feature_representation = [max.(0, data.X1 * v[1] .+ data.X2 * v[2])
		                      for v in vectors]
	DataFrame(feature_representation, names)
end

# ╔═╡ ad2cd0c5-e308-4911-aba9-3b56eda8a2e2
transformed_xor_input = vector_features(xor_data, [1, 1], [1, -1], [-1, 1], [-1, -1])

# ╔═╡ 1b8239fb-3894-47d6-87bf-4592b2fd078d
begin
	lin_mach_f = fit!(machine(LogisticClassifier(penalty = :none),
			                  transformed_xor_input,
			                  xor_data.y))
	lin_pred_f = predict_mode(lin_mach_f, transformed_xor_input)
	(training_misclassification_rate = mean(lin_pred_f .!= xor_data.y),)
end

# ╔═╡ cfcad2d7-7dd9-43fe-9633-7ce9e51e1e27
md"# Transformations of the Output

Sometimes it is reasonable to transform the output data. For example, if a method works best for normally distributed output data, it may be a good idea to transform the output data such that its distribution is closer to normal.
"

# ╔═╡ 5e6d6155-c254-4ae9-a0e1-366fc6ce6403
weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"), DataFrame);

# ╔═╡ ed239411-185a-4d9f-b0ec-360400f24dc7
normal_fit = Distributions.fit(Normal, weather.LUZ_wind_peak)

# ╔═╡ 7c04a091-93af-476e-8769-b14c08d3ae10
begin
	histogram(weather.LUZ_wind_peak, normalize = :pdf,
		      xlabel = "LUZ_wind_peak", label = nothing)
	plot!(x -> pdf(normal_fit, x), w = 2, label = "fitted Normal distribution")
end

# ╔═╡ 3c3aa56f-cb19-4115-9268-b73337fb6a49
lognormal_fit = Distributions.fit(Normal, log.(weather.LUZ_wind_peak))

# ╔═╡ 15faf5e7-fa72-4519-a83b-3007486450dd
begin
	histogram(log.(weather.LUZ_wind_peak), normalize = :pdf,
		      xlabel = "log(LUZ_wind_peak)", label = nothing)
	plot!(x -> pdf(lognormal_fit, x), w = 2, label = "fitted Normal distribution")
end

# ╔═╡ 1e4ba3ab-1bf1-495a-b9c4-b9ea657fc91c
md"Instead of transforming the output variable one can also assume that it has a specific distribution, e.g. a Gamma distribution."

# ╔═╡ 74f5178c-a308-4184-bc94-4a04f6f9ffdc
gamma_fit = Distributions.fit_mle(Gamma, weather.LUZ_wind_peak)

# ╔═╡ 62491145-48b9-4ac5-8e7b-0676e9616fe9
begin
	histogram(weather.LUZ_wind_peak, normalize = :pdf,
		      xlabel = "LUZ_wind_peak", label = nothing)
	plot!(x -> pdf(gamma_fit, x), w = 2, label = "fitted Gamma distribution")
end

# ╔═╡ 04397866-4e00-46d2-bda4-0a61cf34e780
md"# Exercises

## Conceptual
1. Suppose, for a regression problem with one-dimensional input ``X`` and one-dimensional output ``Y`` we use the feature functions ``h_1(X) = X``, ``h_2(X) = (X - 1)^2I(X\ge 1)`` (where ``I`` is the indicator function, being 1 when the condition inside the brackets is true and 0 otherwise), fit the linear regression model ``Y = \beta_0 + \beta_1h_1(X) + \beta_2h_2(X) + \varepsilon,`` and obtain coefficient estimates ``\hat{\beta}_0 = 1``, ``\hat{\beta}_1 = 1``, ``\hat{\beta}_2 = -2``. Draw the estimated curve between ``X = -2`` and ``X = 2``.
2. The following table shows one-hot encoded input for columns `mean_of_transport` (car, train, airplane, ship) and `energy_source` (petrol, electric). Indicate for each row the mean of transport and the energy source.
| car | train | airplane | petrol|
|-----|-------|----------|-------|
| 1   | 0     | 0        |  0    |
| 0   | 0     | 0        |  1    |
| 0   | 1     | 0        |  1    |
| 0   | 0     | 0        |  0    |

3. In the \"Vector Features\" section we said that the xor problem has a non-linear decision boundary in the original ``X_1, X_2`` coordinates, but a linear decision boundary in the ``H_1, H_2, H_3, H_4`` coordinates. Here we prove that ``-H_1 + H_2 + H_3 - H_4 = 0`` defines indeed a linear decision boundary that solves the xor-problem.
    * Show that ``-H_1 + H_2 + H_3 - H_4  < 0`` for all points with ``X_1 > 0`` and ``X_2 > 0``.
    * Show that ``-H_1 + H_2 + H_3 - H_4  < 0`` for all points with ``X_1 < 0`` and ``X_2 < 0``.
    * Show that ``-H_1 + H_2 + H_3 - H_4  > 0`` for all other points.

## Applied
1. Load the mushroom data set `OpenML.load(24)` and determine if a linear logistic regression can correctly classify all datapoints if
   * rows with missing values are dropped.
   * rows with missing values are imputed with the most common class.
2. You are given the following artificial dataset. The data is not linearly separable, meaning that there is no linear decision boundary that would perfectly seperate the blue from the red data points.
    * Find a 2-dimensional feature representation ``H_1 = f_1(X_1, X_2), H_2 = f_2(X_1, X_2)``, such that the transformed data is linearly separable.
    * Generate a similar plot as the one below for the transformed data to visualize that the decision boundary has become linear for the transformed data.
"

# ╔═╡ 9a0e0bf7-44e6-4385-ac46-9a6a8e4245bb
begin
	cldata = DataFrame(2*rand(400, 2).-1, :auto)
	cldata.y = cldata.x1.^2 .> cldata.x2
	cldata
end

# ╔═╡ 6f4427a0-18eb-4637-a19a-ec9aa7b6fda8
scatter(cldata.x1, cldata.x2,
        c = Int.(cldata.y) .+ 1, legend = false,
        xlabel = "X₁", ylabel = "X₂")

# ╔═╡ 951f6957-0ff7-4b2e-ba91-1b69122bbe47
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─4e30f307-7b07-4cb7-a34e-ada7e2f6946d
# ╠═34beb7d6-fba9-11eb-15b8-e34afccc9f88
# ╟─4c536f7c-dc4b-4cf7-aefa-694b830a8e6a
# ╠═026de033-b9b2-4165-80e7-8603a386381c
# ╟─69b9c8ed-acf1-43bf-8478-3a9bc57bc36a
# ╠═6532d7c4-2385-401a-9e99-8dc86655687f
# ╠═2e1c3fa4-e424-4844-82c1-53d02271e7a6
# ╠═d35acc45-7404-4b47-98fe-6e13db8b657e
# ╠═5774b82f-a09e-46f8-a458-0f85fc4c53d8
# ╟─001c88b3-c0db-4ad4-9c0f-04a2ea7b090d
# ╠═0815ba19-5c7d-40d6-b948-8edccfb5c386
# ╟─004958cf-0abb-4d8e-90b1-5966412cba91
# ╠═b43365c6-384e-4fde-95f3-b2ad1ebc421a
# ╠═a91c6298-98c4-4dc6-9450-279da1c34146
# ╠═2755577e-8b4c-4cdc-902e-fe661aa91731
# ╠═81fc52cd-e67b-4188-a028-79abcc77cb92
# ╟─e928e4f7-0ebc-4db5-b83e-0c5d22c0ff3c
# ╠═1bf42b12-5fd1-4b4c-94c6-14a52317b15f
# ╠═0a22b68c-1811-4de1-b3f8-735ee50299d2
# ╠═6323a352-a5f3-4ba7-93bd-903fb1bf22a0
# ╟─85c9650d-037a-4e9f-a15c-464ad479ff44
# ╟─aa002263-4dc8-4238-ab27-a9094947600c
# ╠═c7dab115-b990-4760-85d0-4214d33ada5d
# ╟─e88f9e54-f86c-4d3a-80d7-1a9c5006476a
# ╠═89cbf454-8d1b-4c82-ae39-0672a225e7dd
# ╟─315f536b-4f10-4471-8e22-18c385ca80f2
# ╠═69bc41cb-862c-4334-8b71-7f9974fedfe2
# ╟─76155286-978b-4a05-b428-4be4089d8b9a
# ╠═b2f03b4a-b200-4b30-8793-6edc02f8136d
# ╠═794de3b3-df22-45c1-950d-a515c4f41409
# ╟─adca81d5-2453-4954-91ba-7cc24d45e8ef
# ╠═bf879ed5-4241-4abe-a42c-a562d49b21d2
# ╠═2492bd04-cf24-4252-a37f-05fce8c1b87c
# ╟─32020e6d-5ab2-4832-b565-4264ebcc82eb
# ╠═fd19d841-f99e-46f4-a213-46be7a3d598d
# ╟─8a1a61a0-d657-4858-b8a7-d4e7d9bbf306
# ╠═ee2f2a32-737e-49b5-885b-e7467a7b93f6
# ╠═cc550977-1f37-4ad1-8170-d6150cc01a9c
# ╠═f099be37-82d2-4031-955c-6a797c8400f2
# ╠═75dda04c-c0af-4a05-b840-91f5c1aea53c
# ╟─a10b3fa1-358b-4a94-b927-5e0f71edd1f3
# ╠═b66c7efc-8fd8-4152-9b0f-9332f760b51b
# ╠═ebe38042-5160-469b-8c07-484b7c019063
# ╟─fa7cf2ae-9777-4ba9-9fe0-df545b94c5d8
# ╠═99d3ae30-6bc2-47bf-adc4-1cc1d3ff178d
# ╠═a7f12901-e155-4f61-b252-c5af6669bfee
# ╠═ad2cd0c5-e308-4911-aba9-3b56eda8a2e2
# ╠═1b8239fb-3894-47d6-87bf-4592b2fd078d
# ╟─cfcad2d7-7dd9-43fe-9633-7ce9e51e1e27
# ╠═5e6d6155-c254-4ae9-a0e1-366fc6ce6403
# ╠═ed239411-185a-4d9f-b0ec-360400f24dc7
# ╟─7c04a091-93af-476e-8769-b14c08d3ae10
# ╠═3c3aa56f-cb19-4115-9268-b73337fb6a49
# ╟─15faf5e7-fa72-4519-a83b-3007486450dd
# ╟─1e4ba3ab-1bf1-495a-b9c4-b9ea657fc91c
# ╠═74f5178c-a308-4184-bc94-4a04f6f9ffdc
# ╟─62491145-48b9-4ac5-8e7b-0676e9616fe9
# ╟─04397866-4e00-46d2-bda4-0a61cf34e780
# ╠═9a0e0bf7-44e6-4385-ac46-9a6a8e4245bb
# ╠═6f4427a0-18eb-4637-a19a-ec9aa7b6fda8
# ╟─e688a9de-2dba-4fab-b4e6-9803c5361a62
# ╟─951f6957-0ff7-4b2e-ba91-1b69122bbe47

max_stars_repo_name: jbrea/MLCourse

max_issues_repo_name: davibarreira/MLCourse

max_forks_repo_name: davibarreira/MLCourse


text:
nodes2edges(walk::Array{Int64}) = [Edge(o,d) for (o,d) in zip(walk[1:end-1], walk[2:end])]

walk2prodres(walk::Array{<:Edge}, E::Array{<:Edge}) = sum([map(x-> x==e, E) for e in walk])

function get_prodres(g::SimpleDiGraph, NP::Int64, rnd::AbstractRNG, seed::Int64; max_attempts=10_000)
    walks = Set{Array}()
    i = 1
    while length(walks)<NP && i < max_attempts
        w = self_avoiding_walk(g, rand(rnd, 1:nv(g)),rand(rnd, 1:ne(g)); seed = rand(rnd, 1:(seed+i)))
        if length(w)>1 push!(walks, w) end # Only walks that have some edges
        i+=1
    end
    collect(walks)
    walks_e = [nodes2edges(w) for w in collect(walks)]

    E = collect(edges(g))
    prodress = Array{Bool,1}[walk2prodres(w,E) for w in walks_e]
end

"""
Pricing problem with directed graph product structure, edges are resources
 - seed = seed for generating graph
 - c = initial resource capacity
 - T = same selling period end for all resurces and products
 - expected_res = expected number of consumed resources if all requests were satisfied, used to generate random demand
 - res_budget_μ = mean budget per resource, used to set the customer budgets.
"""

function graph_pp(;NV::Int64=5, NE::Int64=8, NP::Int64=20,
                seed::Int64=2,
                c::Int64 = 5,
                T::Int64 = 20,
                expected_res::Float64 = 20.,
                res_budget_μ::Float64 = 5.,
                objective::Symbol = :revenue)

    """
    Graph
    """
    g = SimpleDiGraph(NV, NE, seed=seed)

    # generate connected graph
    attempts = 0
    N_GRAPH_ATTEMPTS = 1000
    while length(connected_components(g))>1 && attempts <= N_GRAPH_ATTEMPTS
        seed = seed*7879+1
        g = SimpleDiGraph(NV, NE, seed=seed)
        attempts+=1
    end
    if attempts >= N_GRAPH_ATTEMPTS
        throw(ErrorException("Connected graph couldnt be found."))
    end

    graph_pp(g, NP; seed = seed, c=c, T=T,
             expected_res=expected_res,
             res_budget_μ=res_budget_μ)
end

function graph_pp(g::SimpleDiGraph, NP::Int64;
                seed::Int64 = 2,
                c::Int64 = 5,
                T::Int64 = 20,
                expected_res::Float64 = 20.,
                res_budget_μ::Float64 = 5.,
                objective::Symbol = :revenue)
    """
    Capacity
    """
    nᵣ = ne(g)
    c₀ = [c for e in 1:nᵣ]

    """
    Products
    """
    rnd = Xorshift1024Plus(seed)
    P = [PMDPs.Product(prodres, T) for prodres in  get_prodres(g, NP, rnd, seed)]

    D = PMDPs.rand_demand(P, expected_res, rnd)

    B = PMDPs.normal_budgets_per_resource(P, res_budget_μ, res_budget_μ/2)

    A = PMDPs.action_space(P, res_budget_μ)

    PMDPs.PMDPProblem(P, c₀, D, B, A, objective; info=(;g=g))
end


text:
import SparseGrids.grow!

"""
    Grow state space grid from node M.G[i].
"""
function grow!(M::Model, id, bounds = ones(Int, size(M.SP, 2)) * 10)

    bounds[map(v -> isa(v, Stochastic), State(M.variables))] = M.G.L[map(v -> isa(v, Stochastic), State(M.variables))]

    G = deepcopy(M.G)
    grow!(G, id, bounds)
    X = M.G(M.X, values(G))
    G, S, SP, ProbWeights = initS(State(M.variables), typeof(M.G).parameters[2], G)

    M.G     = G
    M.SP    = SP
    M.ProbWeights = ProbWeights
    M.X     = X
    M.XP    = zeros(size(SP, 1), size(M.XP, 2))
    M.Fval  = zeros(length(G), size(M.Fval, 2))
    initD(M)
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
using Revise
using LPDM
using Distributions
using POMDPs
import POMDPs.pdf

const LD1Obs    = Float64

d = Distributions.Normal(10.0,2.0)
rng = RNGVector(1,UInt32(42))
a = rand(rng,d)

POMDPs.pdf(d::Distributions.Normal, o::LD1Obs) = Distributions.pdf(d,o)
POMDPs.pdf(d,11.0)

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
using JuLIP, NBodyIPs, LinearAlgebra

path = homedir() * "/Dropbox/PIBmat/Ti_DFTB_Data/Ti_N54_vartemp.xyz"

# Loading the training data
data = NBodyIPs.Data.read(path)
data_sm = data[1:10]

println("generating basis functions...")
r0 = rnn(:W)
rcuts = [5.2, 5.2, 4.5]  # 2-body, 3-body, 4-body
degrees = [10, 8, 6]
TRANSFORM = "@analytic r -> ($r0/r)^3"

# generate everything else
CUTOFFS = ["(:cos, $(0.66*rcut), $rcut)" for rcut in rcuts]
dicts = Dictionary.(TRANSFORM, CUTOFFS)
B = vcat(NBody(1.0), gen_basis.([2,3,4], dicts, degrees, purify=false)...)
@show length(B)

B3 = B[ findall(bodyorder.(B) .< 4) ]
# normalize_basis!(B3, data_sm)

Φ, _, _ = NBodyIPs.assemble_lsq(B, data_sm; verbose=true, nforces=10)

Q, R = qr(Φ)
A = R' * R

@show norm(A - Φ'*Φ, Inf)

display(A[1:13,1:13])

idx = 12  # this should be a 2-body
idx = 20  # this should be a pure 3-body
interesting_b = B3[idx]
# TODO how do we figure out that this is almost in the 2-body set
#  Φ[:,idx] almost in the span of Φ[:, 1:11]
Φ2 = Φ[:,1:11]
Q2, R2 = qr(Φ2)
for idx = 12:length(B3)
   c = Q2' * Φ[:,idx]
   err = Q2 * c - Φ[:,idx]
   @show B3[idx].t
   @show norm(err, Inf) / norm(Φ[:,idx], Inf)
end

# # [3] do the regression and construct and IP
# c = regression(B, train_data, nforces = Inf, stab = 0.0)
# @show norm(c, Inf)
# IP = NBodyIP(B, c)
# # check error => the normalisation is w.r.t. natoms, not a genuine
# # relative error; we can discuss
# rE, rF, mE, mF = fiterrors(IP, test_data)
# println("   E-rms, E-mae on testset = ", rE, ", ", mE)
# println("   F-rms, F-mae on testset = ", rF, ", ", mF)
#
# # [4] save the IP
# NBodyIPs.IO.write(@__DIR__() * "/ip.jld", IP)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

include("../model.jl")
include("../model_for_mock.jl")

MOIU.@model(ModelForCachingOptimizer,
            (MOI.ZeroOne, MOI.Integer),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.ExponentialCone, MOI.DualExponentialCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.RootDetConeTriangle,
             MOI.LogDetConeTriangle),
            (),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

@testset "Test default attributes" begin
    # Without an optimizer attached (i.e., `MOI.state(model) == NO_OPTIMIZER`) we
    # need throw nice errors for attributes that are based on the optimizer. For
    # `AbstractModelAttribute`s that `is_set_by_optimize` returns `true` for, we
    # overload `TerminationStatus`, `PrimalStatus`, or `DualStatus` to return
    # sane default values. Otherwise we throw a nice error.
    model = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), MOIU.MANUAL)
    @test MOIU.state(model) == MOIU.NO_OPTIMIZER

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
    x = MOI.add_variables(model, 2)
    attr = MOI.VariablePrimal()
    exception = ErrorException(
        "Cannot query $(attr) from caching optimizer because no optimizer" *
        " is attached.")
    @test_throws exception MOI.get(model, MOI.VariablePrimal(), x[1])
    @test_throws exception MOI.get(model, MOI.VariablePrimal(), x)

    attr = MOI.SolverName()
    exception = ErrorException(
        "Cannot query $(attr) from caching optimizer because no optimizer" *
        " is attached.")
    @test_throws exception MOI.get(model, attr)
    attr = MOI.Silent()
    exception = ErrorException(
        "Cannot query $(attr) from caching optimizer because no optimizer" *
        " is attached.")
    @test_throws exception MOI.get(model, attr)

    attr = MOI.ResultCount()
    exception = ErrorException(
        "Cannot query $(attr) from caching optimizer because no optimizer" *
        " is attached.")
    @test_throws exception MOI.get(model, attr)
end

@testset "Copyable solver attributes" begin
    cache = MOIU.UniversalFallback(ModelForCachingOptimizer{Float64}())
    cached = MOIU.CachingOptimizer(cache, MOIU.MANUAL)
    MOI.set(cached, MOI.Silent(), true)
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(ModelForMock{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
    MOI.set(cached, MOI.Silent(), false)
    @test !MOI.get(mock, MOI.Silent())
    @test !MOI.get(cached, MOI.Silent())
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(ModelForMock{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test !MOI.get(mock, MOI.Silent())
    @test !MOI.get(cached, MOI.Silent())
    MOI.set(cached, MOI.Silent(), true)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(ModelForMock{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
end

@testset "CachingOptimizer MANUAL mode" begin
    m = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), MOIU.MANUAL)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
    @test MOI.is_empty(s)
    MOIU.reset_optimizer(m, s)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

    v = MOI.add_variable(m)
    x = MOI.add_variables(m, 2)
    saf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], [v; x]), 0.0)
    @test MOI.supports(m, MOI.ObjectiveFunction{typeof(saf)}())
    MOI.set(m, MOI.ObjectiveFunction{typeof(saf)}(), saf)
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveFunction{typeof(saf)}())) ≈ saf
    @test MOI.get(m, MOI.ObjectiveFunction{typeof(saf)}()) ≈ saf

    @test_throws AssertionError MOI.optimize!(m)

    MOIU.attach_optimizer(m)
    @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.ObjectiveFunction{typeof(saf)}())) ≈ saf

    @test MOI.supports(m, MOI.ObjectiveSense())
    MOI.set(m, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(m, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveSense())) == MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.ObjectiveSense())) == MOI.MAX_SENSE

    @test MOI.supports(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()))
    MOI.set(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()), 10)
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute())) == 10

    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.ResultCount()), 1)
    @test MOI.supports(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), typeof(v))
    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v, 3.0)

    MOI.optimize!(m)

    @test MOI.get(m, MOI.VariablePrimal(), v) == 3.0
    @test MOI.get(m, MOI.VariablePrimal(), [v]) == [3.0]
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v) == 3.0

    # ModelForMock doesn't support RotatedSecondOrderCone
    @test !MOI.supports_constraint(m, MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)

    @test MOI.supports_constraint(m.model_cache, MOI.SingleVariable, MOI.LessThan{Float64})
    @test MOI.supports_constraint(m.optimizer.inner_model, MOI.SingleVariable, MOI.LessThan{Float64})
    @test MOI.supports_constraint(m.optimizer, MOI.SingleVariable, MOI.LessThan{Float64})
    @test MOI.supports_constraint(m, MOI.SingleVariable, MOI.LessThan{Float64})
    lb = MOI.add_constraint(m, MOI.SingleVariable(v), MOI.LessThan(10.0))
    MOI.set(m, MOI.ConstraintSet(), lb, MOI.LessThan(11.0))
    @test MOI.get(m, MOI.ConstraintSet(), lb) == MOI.LessThan(11.0)
    @test MOI.get(m, MOI.ConstraintFunction(), lb) == MOI.SingleVariable(v)

    MOIU.drop_optimizer(m)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    MOI.set(m, MOI.ConstraintSet(), lb, MOI.LessThan(12.0))
    @test MOI.get(m, MOI.ConstraintSet(), lb) == MOI.LessThan(12.0)

    MOI.delete(m, x[2])
    @test_throws MOI.InvalidIndex{typeof(x[2])} MOI.delete(m, x[2])
    @test !MOI.is_valid(m, x[2])

    # TODO: test more constraint modifications

end

@testset "CachingOptimizer AUTOMATIC mode" begin
    m = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), MOIU.AUTOMATIC)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    v = MOI.add_variable(m)
    @test MOI.supports(m, MOI.VariableName(), typeof(v))
    MOI.set(m, MOI.VariableName(), v, "v")
    @test MOI.get(m, MOI.VariableName(), v) == "v"

    s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
    @test MOI.is_empty(s)
    MOIU.reset_optimizer(m, s)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

    saf = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, v)], 0.0)
    @test MOI.supports(m, MOI.ObjectiveFunction{typeof(saf)}())
    MOI.set(m, MOI.ObjectiveFunction{typeof(saf)}(), saf)
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveFunction{typeof(saf)}())) ≈ saf

    MOI.optimize!(m)
    @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.TerminationStatus())) == MOI.OPTIMIZE_NOT_CALLED

    @test MOI.get(m, MOI.VariableName(), v) == "v"
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.VariableName()), v) == "v"
    # The caching optimizer does not copy names
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariableName()), v) == ""

    @test MOI.supports(m, MOI.ObjectiveSense())
    MOI.set(m, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveSense())) == MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.ObjectiveSense())) == MOI.MAX_SENSE

    @test MOI.supports(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()))
    MOI.set(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()), 10)
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute())) == 10

    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.ResultCount()), 1)
    @test MOI.supports(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), typeof(v))
    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v, 3.0)

    MOI.optimize!(m)

    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v) == 3.0

    @testset "Modify not allowed" begin
        s.modify_allowed = false
        MOI.modify(m, MOI.ObjectiveFunction{typeof(saf)}(),
                    MOI.ScalarConstantChange(1.0))
        s.modify_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
        ci = MOI.add_constraint(m, saf, MOI.EqualTo(0.0))
        s.modify_allowed = false
        MOI.modify(m, ci, MOI.ScalarCoefficientChange(v, 1.0))
        s.modify_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    end

    # Simulate that constraints cannot be added
    @testset "Add constraint not allowed" begin
        s.add_con_allowed = false
        MOI.add_constraint(m, MOI.VectorOfVariables([v]), MOI.SecondOrderCone(1))
        s.add_con_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOI.empty!(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    end

    @testset "Add variable not allowed" begin
        s.add_var_allowed = false # Simulate optimizer that cannot add variables incrementally
        MOI.add_variable(m)
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        @test_throws MOI.AddVariableNotAllowed MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

        s.add_var_allowed = true
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
        s.add_var_allowed = false
        MOI.add_variables(m, 2)
        s.add_var_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    end

    @testset "Delete not allowed" begin
        vi = MOI.add_variable(m)
        s.delete_allowed = false # Simulate optimizer that cannot delete variable
        MOI.delete(m, vi)
        s.delete_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER

        vi = MOI.add_variable(m)
        ci = MOI.add_constraint(m, MOI.SingleVariable(vi), MOI.EqualTo(0.0))
        s.delete_allowed = false # Simulate optimizer that cannot delete constraint
        MOI.delete(m, ci)
        s.delete_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    end

end

@testset "Constructor with optimizer" begin
    @testset "Empty model and optimizer" begin
        s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
        model = ModelForCachingOptimizer{Float64}()
        m = MOIU.CachingOptimizer(model, s)
        @test m isa MOIU.CachingOptimizer{typeof(s), typeof(model)}
        @test MOI.is_empty(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
        @test MOIU.mode(m) == MOIU.AUTOMATIC
        @test MOI.get(m, MOI.SolverName()) == "Mock"
    end
    @testset "Non-empty optimizer" begin
        s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
        MOI.add_variable(s)
        model = ModelForCachingOptimizer{Float64}()
        @test MOI.is_empty(model)
        @test !MOI.is_empty(s)
        @test_throws AssertionError MOIU.CachingOptimizer(model, s)
    end
    @testset "Non-empty model" begin
        s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
        model = ModelForCachingOptimizer{Float64}()
        MOI.add_variable(model)
        @test !MOI.is_empty(model)
        @test MOI.is_empty(s)
        @test_throws AssertionError MOIU.CachingOptimizer(model, s)
    end
end

for state in (MOIU.NO_OPTIMIZER, MOIU.EMPTY_OPTIMIZER, MOIU.ATTACHED_OPTIMIZER)
    @testset "Optimization tests in state $state and mode $mode" for mode in (MOIU.MANUAL, MOIU.AUTOMATIC)
        m = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), mode)
        if state != MOIU.NO_OPTIMIZER
            s = MOIU.MockOptimizer(ModelForMock{Float64}(), supports_names=false)
            MOIU.reset_optimizer(m, s)
            if state == MOIU.ATTACHED_OPTIMIZER
                MOIU.attach_optimizer(m)
            end
        end
        @test MOIU.state(m) == state
        @test MOIU.mode(m) == mode

        @testset "Name test" begin
            MOIT.nametest(m)
        end

        @testset "Copy test" begin
            MOIT.failcopytestc(m)
            MOIT.failcopytestia(m)
            MOIT.failcopytestva(m)
            MOIT.failcopytestca(m)
            MOIT.copytest(m, Model{Float64}())
        end

        config = MOIT.TestConfig(solve=false)
        @testset "Unit" begin
            MOIT.unittest(m, config)
        end
        @testset "Continuous Linear" begin
            exclude = ["partial_start"] # VariablePrimalStart not supported.
            MOIT.contlineartest(m, config, exclude)
        end
    end
end


text:
"""
    AngularCentralGaussian(M; params...)

The Angular Central Gaussian (ACG) distribution on the manifold ``M``.

Accepted manifolds are `Sphere`, `ProjectiveSpace`, `Stiefel`, `Grassmann`, `Rotations`
and `SpecialOrthogonal`.

For manifolds with matrix points, this is also called the Matrix Angular Central Gaussian
distribution.

# Constructors

    AngularCentralGaussian(M; P)
    AngularCentralGaussian(M; L)

For a manifold ``M ⊂ 𝔽^{n × k}``, construct the ACG distribution parameterized either by
the inverse ``P = Σ^{-1}`` of an ``n × n`` positive definite matrix ``Σ`` or by the lower
Cholesky factor ``L`` of ``Σ``, such that ``Σ = L L^\\mathrm{H}``.
"""
struct AngularCentralGaussian{M,N,T} <: ParameterizedMeasure{N}
    manifold::M
    par::NamedTuple{N,T}
end
AngularCentralGaussian(M; params...) = AngularCentralGaussian(M, NamedTuple(params))

function Base.show(io::IO, mime::MIME"text/plain", μ::AngularCentralGaussian)
    return show_manifold_measure(io, mime, μ)
end

Manifolds.base_manifold(d::AngularCentralGaussian) = getfield(d, :manifold)

function MeasureTheory.basemeasure(μ::AngularCentralGaussian)
    return normalize(Hausdorff(base_manifold(μ)))
end

# Chikuse, 2003 Eq. 2.4.3
# extended to complex Stiefel using Eq. 6 of
# Wróblewska J. A note on some extensions of the matrix angular central Gaussian distribution. 2020.
# https://arxiv.org/abs/2010.03243
function MeasureTheory.logdensity(μ::AngularCentralGaussian{<:Any,(:P,)}, x)
    M = base_manifold(μ)
    d = real_dimension(number_system(M))
    s = representation_size(M)
    n, k = length(s) > 1 ? s : (first(s), 1)
    P = μ.P
    logdetx′Px = isone(k) ? real(dot(x, P, x)) : real(logdet(x' * P * x))
    return d * (k * real(logdet(P)) - n * logdetx′Px) / 2
end
function MeasureTheory.logdensity(μ::AngularCentralGaussian{<:Any,(:L,)}, x)
    M = base_manifold(μ)
    d = real_dimension(number_system(M))
    s = representation_size(M)
    n, k = length(s) > 1 ? s : (first(s), 1)
    L = LowerTriangular(μ.L)
    z = L \ x
    logdetx′Px = isone(k) ? log(real(norm(z))) : real(logdet(z'z)) / 2
    return -d * (k * real(logdet(L)) + n * logdetx′Px)
end

function Base.rand(rng::AbstractRNG, T::Type, d::AngularCentralGaussian)
    p = default_point(d, T)
    return Random.rand!(rng, p, d)
end

function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::AngularCentralGaussian{M,(:L,)}
) where {M}
    z = randn!(rng, p)
    return project!(base_manifold(d), p, d.L * z)
end


text:
"""
    QDLDL inplace functionality
"""
mutable struct LDLSolver{Tf<:AbstractFloat,Ti<:Integer} <: LinearSolver
    # QDLDL Factorization
    F::QDLDL.QDLDLFactorisation{Tf,Ti}
    # Allocate memory
    Pr::Vector{Ti}
    Pc::Vector{Ti}
    Pv::Vector{Tf}
    num_entries::Vector{Ti}
end

function LDLSolver(A::SparseMatrixCSC{Tv,Ti}, F::QDLDL.QDLDLFactorisation{Tv,Ti}) where {Tv<:AbstractFloat,Ti<:Integer}
    Pr = zeros(Ti, nnz(A))
    Pc = zeros(Ti, size(A, 1) + 1)
    Pv = zeros(Tv, nnz(A))
    num_entries = zeros(Ti, size(A, 2))
    return LDLSolver{Tv,Ti}(F, Pr, Pc, Pv, num_entries)
end

function factorize!(s::LDLSolver{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}) where {Tv<:AbstractFloat, Ti<:Integer}
    # Reset the pre-allocated fields
    s.Pr .= 0
    s.Pc .= 0
    s.Pv .= 0.0
    s.num_entries .= 0

    # Triangularize the matrix with the allocation-free method.
    A = permute_symmetricAF(A, s.F.iperm, s.Pr, s.Pc, s.Pv, s.num_entries)  #returns an upper triangular matrix

    # Update the workspace, triuA is the only field we need to update
    s.F.workspace.triuA.nzval .= A.nzval

    #factor the matrix
    QDLDL.factor!(s.F.workspace, s.F.logical)

    return nothing
end

function permute_symmetricAF(A::SparseMatrixCSC{Tv, Ti}, iperm::AbstractVector{Ti},
    Pr::AbstractVector{Ti}, Pc::AbstractVector{Ti}, Pv::AbstractVector{Tv},
    num_entries::AbstractVector{Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    # 1. count number of entries that each column of P will have
    n = size(A, 2)
    # num_entries = zeros(Ti, n)
    Ar = A.rowval
    Ac = A.colptr
    Av = A.nzval
    # count the number of upper-triangle entries in columns of P, keeping in mind the row permutation
    for colA = 1:n
        colP = iperm[colA]
        # loop over entries of A in column A...
        for row_idx = Ac[colA]:Ac[colA+1]-1
            rowA = Ar[row_idx]
            rowP = iperm[rowA]
            # ...and check if entry is upper triangular
            if rowA <= colA
                # determine to which column the entry belongs after permutation
                col_idx = max(rowP, colP)
                num_entries[col_idx] += one(Ti)
            end
        end
    end
    # 2. calculate permuted Pc = P.colptr from number of entries
    Pc[1] = one(Ti)
    @inbounds for k = 1:n
        Pc[k + 1] = Pc[k] + num_entries[k]

        # reuse this vector memory to keep track of free entries in rowval
        num_entries[k] = Pc[k]
    end
    # use alias
    row_starts = num_entries

    # 3. permute the row entries and position of corresponding nzval
    for colA = 1:n
        colP = iperm[colA]
        # loop over rows of A and determine where each row entry of A should be stored
        for rowA_idx = Ac[colA]:Ac[colA+1]-1
            rowA = Ar[rowA_idx]
            # check if upper triangular
            if rowA <= colA
                rowP = iperm[rowA]
                # determine column to store the entry
                col_idx = max(colP, rowP)

                # find next free location in rowval (this results in unordered columns in the rowval)
                rowP_idx = row_starts[col_idx]

                # store rowval and nzval
                Pr[rowP_idx] = min(colP, rowP)
                Pv[rowP_idx] = Av[rowA_idx]

                # increment next free location
                row_starts[col_idx] += 1
            end
        end
    end
    P = SparseMatrixCSC{Tv, Ti}(n, n, Pc, Pr, Pv)
    # order row indices within P.rowcal[P.colptr[k]:P.colptr[k+1]-1]
    return (P')'
end

"""
    LDL solver
"""
function ldl_solver(A::SparseMatrixCSC{T,Int}) where T
    LDLSolver(A, qdldl(A))
end

ldl_solver(A::Array{T, 2}) where T = ldl_solver(sparse(A))

function linear_solve!(solver::LDLSolver{Tv,Ti}, x::Vector{Tv}, A::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv}) where {Tv<:AbstractFloat,Ti<:Integer}
    factorize!(solver, A) # factorize
    x .= b
    QDLDL.solve!(solver.F, x) # solve
end

function linear_solve!(solver::LDLSolver{Tv,Ti}, X::Matrix{Tv}, A::AbstractMatrix{Tv}, B::AbstractMatrix{Tv}) where {Tv<:AbstractFloat,Ti<:Integer}
    X .= A \ B # TODO: fix
end