text:
using Manifolds
using ApproxManifoldProducts
using TensorCast
using LinearAlgebra
using Test

##

@testset "test calcProductGaussians" begin
##

M = TranslationGroup(2)
u = [[1;1.0], [0.0;0]]
c = [([1.0;1]), ([1.0;1])]

u_, C_ = calcProductGaussians(M, u, c)
@test isapprox( u_, [0.5, 0.5])
@test isapprox( C_, [0.5 0.0; 0.0 0.5])

##
end

@testset "comparison test with basic product" begin
## simply multiply two beliefs, sim2

N = 10
d = 2
M = TranslationGroup(d)

#densities to multiply
pts1 = [randn(d) for _ in 1:N]
P1 = manikde!(M, pts1, bw=[1;1.0])

pts2 = [randn(d)  for _ in 1:N]
P2 = manikde!(M, pts2, bw=[1;1.0])

##

sl = Vector{Vector{Int}}()

P12 = manifoldProduct([P1;P2], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P12)
@test P12._partial === nothing
# @test isapprox( mean(P12)[1], 0, atol=1 )
# @test isapprox( mean(P12)[2], 0, atol=1 )

(x->println()).(1:5)
@show sl;

P12

## validate selected labels are working properly, with addEntropy=false

for sidx in 1:N

  bw1 = getBW(P1)[:,1] .^ 2
  bw2 = getBW(P2)[:,1] .^ 2

  u1 = pts1[sl[sidx][1]]
  u2 = pts2[sl[sidx][2]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);

  @test isapprox( u12, getPoints(P12)[sidx] )

end

##
end

@testset "test dim=2 product with one partial/marginal" begin
## basic test one full with one partial

N = 50
d = 2
M = TranslationGroup(d)

#densities to multiply
pts1 = [randn(d) for _ in 1:N]
P1 = manikde!(M, pts1)

pts2 = [randn(d) for _ in 1:N]
(x->(x[2]+=100)).(pts2)
P2_ = manikde!(M, pts2, partial=[1;])

##

sl = Vector{Vector{Int}}()

P12_ = manifoldProduct([P1;P2_], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P12_)

@test isapprox( mean(P12_)[1], 0, atol=1 )
@test isapprox( mean(P12_)[2], 0, atol=1 )

(x->println()).(1:5)
# @show sl

P12_

##

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw2 = getBW(P2_, false)[:,1] .^2

  u1 = pts1[sl[sidx][1]]
  u2 = pts2[sl[sidx][2]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);

  @test isapprox( u12[1], getPoints(P12_)[sidx][1] , atol=0.1)
  @test isapprox( pts1[sl[sidx][1]][2], getPoints(P12_)[sidx][2] )

end

##
end

@testset "test dim=2 product with one full and two similar marginals" begin
## basic test one full with one partial

N = 50
d = 2
M = TranslationGroup(d)

#densities to multiply
pts1 = [randn(d) for _ in 1:N]
P1 = manikde!(M, pts1)

pts2 = [randn(d) for _ in 1:N]
(x->(x[2]+=100)).(pts2)
P2_ = manikde!(M, pts2, partial=[1;])

pts3 = [randn(d) for _ in 1:N]
(x->(x[2]+=100)).(pts3)
P3_ = manikde!(M, pts3, partial=[1;])

##

sl = Vector{Vector{Int}}()

P123_ = manifoldProduct([P1;P2_;P3_], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P123_)

@test isapprox( mean(P123_)[1], 0, atol=1 )
@test isapprox( mean(P123_)[2], 0, atol=1 )

(x->println()).(1:5)
# @show sl

P123_

##

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw2 = getBW(P2_, false)[:,1] .^2
  bw3 = getBW(P3_, false)[:,1] .^2

  u1 = pts1[sl[sidx][1]]
  u2 = pts2[sl[sidx][2]]
  u3 = pts3[sl[sidx][3]]

  u123, = calcProductGaussians(M, [u1,u2,u3], [bw1,bw2,bw3]);

  @test isapprox( u123[1], getPoints(P123_)[sidx][1] , atol=0.1)
  @test isapprox( pts1[sl[sidx][1]][2], getPoints(P123_)[sidx][2] )

end

##
end

@testset "test dim=2 product with one full and two different marginals" begin
## random data

N = 50
M = TranslationGroup(2)

pts1 = [randn(2) .- 10.0 for _ in 1:N]
pts2 = [randn(2) for _ in 1:N]
pts3 = [randn(2) .+ 10.0 for _ in 1:N]

# get different marginals

P1 = marginal(manikde!(M, pts1), [1;])
P2 =          manikde!(M, pts2)
P3 = marginal(manikde!(M, pts3), [2;])

##

sl = Vector{Vector{Int}}()
P = manifoldProduct([P2;P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P)

(x->println()).(1:5)
# @show sl;
P

## check the results

pts = getPoints(P)
@cast pGM[i,j] := pts[j][i]

@test 0.7*N < sum(-10 .< pGM[1,:] .< 0)
@test 0.7*N < sum(0 .< pGM[2,:] .< 10 )

## check the selection of labels and resulting Gaussian products are correct

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw2 = getBW(P2, false)[:,1] .^2
  bw3 = getBW(P3, false)[:,1] .^2

  # full density first
  u2 = pts2[sl[sidx][1]]
  u1 = pts1[sl[sidx][2]]
  u3 = pts3[sl[sidx][3]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);
  u23, = calcProductGaussians(M, [u2,u3], [bw2,bw3]);

  @test isapprox( u12[1], getPoints(P)[sidx][1])
  @test isapprox( u23[2], getPoints(P)[sidx][2])

end

##
end

@testset "product of only one marginal per each of two dimensions" begin
## random data

N = 50
M = TranslationGroup(2)

pts1 = [randn(2) .- 10.0 for _ in 1:N]
pts3 = [randn(2) .+ 10.0 for _ in 1:N]

# different marginals

P1 = marginal(manikde!(M, pts1), [1;])
P3 = marginal(manikde!(M, pts3), [2;])

##

sl = Vector{Vector{Int}}()

P_ = manifoldProduct([P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P_)

(x->println()).(1:5)
# @show sl

##

pts = getPoints(P_)
@cast pGM[i,j] := pts[j][i]

@test 0.7*N < sum(-13 .< pGM[1,:] .< -7)
@test 0.7*N < sum(  7 .< pGM[2,:] .< 13 )

## check the selection of labels and resulting Gaussian products are correct

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw3 = getBW(P3, false)[:,1] .^2

  u1 = pts1[sl[sidx][1]]
  u3 = pts3[sl[sidx][2]]

  @test isapprox( u1[1], pts[sidx][1])
  @test isapprox( u3[2], pts[sidx][2])

end

##
end

@testset "test dim=2 product of only marginals, two per dimension, 4 factors total" begin
##

N = 50
d = 2
M = TranslationGroup(d)

pts4 = [randn(d) .- 10.0 for _ in 1:N]
(x->x[2]-=90.0).(pts4)
pts5 = [randn(d) .+ 10.0 for _ in 1:N]
(x->x[1]+=90.0).(pts5)

P4 = marginal(manikde!(M, pts4), [1;])
P5 = marginal(manikde!(M, pts5), [d;])

# test duplication
pts4_ = [randn(d) .- 10.0 for _ in 1:N]
(x->x[2]-=90.0).(pts4_)
pts5_ = [randn(d) .+ 10.0 for _ in 1:N]
(x->x[1]+=90.0).(pts5_)

P4_ = marginal(manikde!(M, pts4_), [1;])
P5_ = marginal(manikde!(M, pts5_), [d;])

##

sl = Vector{Vector{Int}}()

P45__ = manifoldProduct([P4;P4_;P5;P5_], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P45__)

(x->println()).(1:5)
# @show sl;

P45__

## check the selection of labels and resulting Gaussian products are correct

# println("getPoints(P45__) = ")
# getPoints(P45__) .|> println
# println()

# sidx = 1
for sidx in 1:N

  bw1 = getBW(P4, false)[:,1] .^2
  bw2 = getBW(P4_, false)[:,1] .^2
  bw3 = getBW(P5, false)[:,1] .^2
  bw4 = getBW(P5_, false)[:,1] .^2

  u1 = pts4[ sl[sidx][1]]
  u2 = pts4_[sl[sidx][2]]
  u3 = pts5[ sl[sidx][3]]
  u4 = pts5_[sl[sidx][4]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);
  u34, = calcProductGaussians(M, [u3,u4], [bw3,bw4]);

  # @info "DEBU" string(u1) string(u2) string(u3) string(u4) u12[1], u34[2]
  @test isapprox( u12[1], getPoints(P45__)[sidx][1])
  @test isapprox( u34[2], getPoints(P45__)[sidx][2])

end

##
end

@testset "test dim=3 product with one full and two different marginals" begin
## random data

d = 3
N = 50
M = TranslationGroup(3)

pts1 = [randn(d) .- 10.0 for _ in 1:N]
pts2 = [randn(d) for _ in 1:N]
pts3 = [randn(d) .+ 10.0 for _ in 1:N]

## get different marginals

P1 = marginal(manikde!(M, pts1), [1;])
P2 =          manikde!(M, pts2)
P3 = marginal(manikde!(M, pts3), [d;])

##

# @test_broken false
# # P = manifoldProduct([P1;P2;P3])
# @error "manifoldProduct needs complete point type for keyword `oldPoints`, current tests assume no marginal beliefs at front of product array"
sl = Vector{Vector{Int}}()
P = manifoldProduct([P2;P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P)

(x->println()).(1:5)
# @show sl;
P

## check the results

pts = getPoints(P)
@cast pGM[i,j] := pts[j][i]

@test 0.6*N < sum(-10 .< pGM[1,:] .< 0)
@test 0.6*N < sum(  0 .< pGM[3,:] .< 10 )

## check the selection of labels and resulting Gaussian products are correct

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw2 = getBW(P2, false)[:,1] .^2
  bw3 = getBW(P3, false)[:,1] .^2

  # full density first
  u2 = pts2[sl[sidx][1]]
  u1 = pts1[sl[sidx][2]]
  u3 = pts3[sl[sidx][3]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);
  u23, = calcProductGaussians(M, [u2,u3], [bw2,bw3]);

  @test isapprox( u12[1], getPoints(P)[sidx][1])
  @test isapprox(  u2[2], getPoints(P)[sidx][2] )
  @test isapprox( u23[3], getPoints(P)[sidx][3])

end

##
end

@testset "test dim=3 product with one full and two different marginals (marginal first in product)" begin
## random data

d = 3
N = 50
M = TranslationGroup(3)

pts1 = [randn(d) .- 10.0 for _ in 1:N]
pts2 = [randn(d) for _ in 1:N]
pts3 = [randn(d) .+ 10.0 for _ in 1:N]

## get different marginals

P1 = marginal(manikde!(M, pts1), [1;])
P2 =          manikde!(M, pts2)
P3 = marginal(manikde!(M, pts3), [d;])

##

sl = Vector{Vector{Int}}()
P = manifoldProduct([P1;P2;P3], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test !isPartial(P)

(x->println()).(1:5)
# @show sl;
P

## check the results

pts = getPoints(P)
@cast pGM[i,j] := pts[j][i]

@test 0.8*N < sum(-10 .< pGM[1,:] .<  0 )
@test 0.8*N < sum(  0 .< pGM[3,:] .< 10 )

## check the selection of labels and resulting Gaussian products are correct

for sidx in 1:N

  bw1 = getBW(P1, false)[:,1] .^2
  bw2 = getBW(P2, false)[:,1] .^2
  bw3 = getBW(P3, false)[:,1] .^2

  # full density first
  u1 = pts1[sl[sidx][1]]
  u2 = pts2[sl[sidx][2]]
  u3 = pts3[sl[sidx][3]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);
  u23, = calcProductGaussians(M, [u2,u3], [bw2,bw3]);

  @test isapprox( u12[1], getPoints(P)[sidx][1])
  @test isapprox(  u2[2], getPoints(P)[sidx][2])
  @test isapprox( u23[3], getPoints(P)[sidx][3])

end

##
end

@testset "test dim=3 product with two different marginals and one open dimension" begin
## random data

d = 3
N = 50
M = TranslationGroup(3)

pts1 = [randn(d) .- 10.0 for _ in 1:N]
pts3 = [randn(d) .+ 10.0 for _ in 1:N]

## get different marginals

P1 = marginal(manikde!(M, pts1), [1;])
P3 = marginal(manikde!(M, pts3), [d;])

##

sl = Vector{Vector{Int}}()
P = manifoldProduct([P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=false)

@test isPartial(P)
@test P._partial == [1;3]

(x->println()).(1:5)
# @show sl;
P

## check the results

pts = getPoints(P, false)
@cast pGM[i,j] := pts[j][i]

@test 0.7*N < sum(-13 .< pGM[1,:] .< -7 )
@test 0.7*N < sum(  7 .< pGM[3,:] .< 13 )

## check the selection of labels and resulting Gaussian products are correct

for sidx in 1:N

  bw1 = getBW(P1)[:,1] .^2
  bw3 = getBW(P3)[:,1] .^2

  # full density first
  u1 = pts1[sl[sidx][1]]
  u3 = pts3[sl[sidx][2]]

  @test isapprox( u1[1], getPoints(P, false)[sidx][1])
  @test isapprox( u3[3], getPoints(P, false)[sidx][3])

end

##
end

#


text:
include("../PlotUtils/PlotUtils.jl")
include("../PlotUtils/imports.jl")

# Reproduce sims metrics

# Name of output file
OUTPUT_FILE = "output.bson"

println("Producing metrics ...")

results_dirs = ["results/sim-runs-$(dsize)/mm0/"
                for dsize in ("small", "large")]

for d in results_dirs
  println(d)
  metrics = PlotUtils.make_metrics(d, OUTPUT_FILE, thresh=.01)
end

# Missmech vs LPML, DIC
function mm_compare(datsize)
  if datsize == "small"
    results_dir = "results/sim-runs-small"
    output_paths = ["$(results_dir)/mm$(mm)/Kmcmc_5/output.bson" for mm in 0:2]
  else
    results_dir = "results/sim-runs-large"
    output_paths = ["$(results_dir)/mm$(mm)/Kmcmc_10/output.bson" for mm in 0:2]
  end

  LPML = Float64[]
  DIC = Float64[]
  q = [(0., 0.25, 0.5),  # mm0
       (0, 0.2, 0.4),  # mm1
       (0, 0.15, 0.3)]  # mm2
  rho = (.05, .8, .05)

  for p in output_paths
    metrics = BSON.load(p)[:metrics]
    append!(LPML, metrics[:LPML])
    append!(DIC, metrics[:DIC])
  end

  # Data frame comparing missing mechanisms and LPML / DIC
  mm_compare_df = DataFrames.DataFrame(missing_mech=0:2, q=q,
                                       rho=rho, LPML=LPML, DIC=DIC)

  # Write data frame. (separated by semi-colon)
  CSV.write("$(results_dir)/mm_compare.csv", mm_compare_df, delim=';')
end

foreach(mm_compare, ["small", "large"])

println("DONE!")

ext: jl

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
@testset "world" begin
    straight_agent = Agent(0.2, 0.0)
    circling_agent = Agent(0.2, 10.0 / 180.0 * pi)
    landmarks =
        [Landmark([2.0, -2.0], 0), Landmark([-1.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    m = Map()
    camera1 = IdealCamera(landmarks)
    camera2 = IdealCamera(landmarks)
    push!(m, landmarks[1])
    push!(m, landmarks[2])
    push!(m, landmarks[3])

    robot1 = IdealRobot([2.0, 3.0, pi / 6], straight_agent, camera1, 0.05, "blue")
    robot2 = IdealRobot([-2.0, -1.0, pi / 5 * 6], circling_agent, camera2, 0.05, "orange")
    xlim = [-5.5, 10]
    ylim = [-5.5, 10]
    world = World(xlim, ylim)
    push!(world, robot1)
    push!(world, robot2)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        obsv1 = observations(camera1, robot1.pose_)
        obsv2 = observations(camera2, robot2.pose_)
        p = draw(world, nothing)
        plot(p)
        v, ω = decision(straight_agent, obsv1)
        state_transition(robot1, v, ω, dt)
        v, ω = decision(circling_agent, obsv2)
        state_transition(robot2, v, ω, dt)
    end
    if GUI
        gif(anim, "ch3_robot11.gif", fps = 20)
    end
end

@testset "ch03_robot11" begin
    straight_agent = Agent(0.2, 0.0)
    circling_agent = Agent(0.2, 10.0 / 180.0 * pi)
    landmarks =
        [Landmark([2.0, -2.0], 0), Landmark([-1.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    m = Map()
    camera1 = IdealCamera(landmarks)
    camera2 = IdealCamera(landmarks)
    push!(m, landmarks)
    robot1 = IdealRobot([2.0, 3.0, pi / 6], straight_agent, camera1, 0.05, "blue")
    robot2 = IdealRobot([-2.0, -1.0, pi / 5 * 6], circling_agent, camera2, 0.05, "orange")
    xlim = [-5.5, 10]
    ylim = [-5.5, 10]
    world = World(xlim, ylim)
    push!(world, robot1)
    push!(world, robot2)
    push!(world, m)
    dt = 0.1
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        obsv1 = observations(camera1, robot1.pose_)
        obsv2 = observations(camera2, robot2.pose_)
        p = draw(world, annota)
        plot(p)
        v, ω = decision(straight_agent, obsv1)
        state_transition(robot1, v, ω, dt)
        v, ω = decision(circling_agent, obsv2)
        state_transition(robot2, v, ω, dt)
    end
    if GUI
        gif(anim, "ch3_robot11.gif", fps = 20)
    end
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text: Dict{String, Any}("a" => Any[Dict{String, Any}("b" => Any[Dict{String, Any}("c" => Dict{String, Any}("d" => Dict{String, Any}("value" => "val0", "type" => "string"))), Dict{String, Any}("c" => Dict{String, Any}("d" => Dict{String, Any}("value" => "val1", "type" => "string")))])])


text:
using StrainRecon
using Base.Test

if StrainRecon.hasGurobi
	@testset "getMAPConvex" begin

	M  = [1 1 1; 1 1 0; 1 0 1; 1 0 0; 0 1 1; 0 1 0; 0 0 1;]
	(m,n) = size(M)
	R     = StrainRecon.sdiag(ones(m)*100);
	atrue = [.6; 0.3; 0.1]
	Data  = M*atrue

	param = getStrainReconParam(Data,n,R=R)

	Mmap,amap = solveConvexMAP(param)
	@test norm(amap - atrue)/norm(atrue) < 1e-8
	@test norm(Mmap*amap-Data) < 1e-7
	end
end


text:
module SLPA

using StatsBase
using Graphs

export
    slpa,
    post_processing!,
    getgrp,
    nmi,
    msp2grp,
    grp2msp,
    readgrp,
    readmsp,
    writegrp,
    writemsp

include("slpafunc.jl")
include("utils.jl")
include("nmi.jl")

end # module


text:
using MathOptInterface, GalacticOptim, Optim, Test, Random, Flux, ForwardDiff, Zygote
using Nonconvex
Nonconvex.@load MMA
Nonconvex.@load Ipopt
Nonconvex.@load NLopt
Nonconvex.@load BayesOpt
Nonconvex.@load Juniper
Nonconvex.@load Pavito
Nonconvex.@load Hyperopt
Nonconvex.@load MTS

rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p  = [1.0, 100.0]

f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
l1 = rosenbrock(x0, _p)
prob = OptimizationProblem(f, x0, _p)
Random.seed!(1234)
sol = solve(prob, SimulatedAnnealing())
@test 10*sol.minimum < l1

Random.seed!(1234)
prob = OptimizationProblem(f, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, SAMIN())
@test 10*sol.minimum < l1

using CMAEvolutionStrategy
sol = solve(prob, CMAEvolutionStrategyOpt())
@test 10*sol.minimum < l1

rosenbrock(x, p) =  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
prob = OptimizationProblem(rosenbrock, x0)
sol = solve(prob, Optim.NelderMead())
@test 10*sol.minimum < l1

cons= (x,p) -> [x[1]^2 + x[2]^2]
optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons)

prob = OptimizationProblem(optprob, x0)

sol = solve(prob, Flux.ADAM(0.1), maxiters = 1000)
@test 10*sol.minimum < l1

sol = solve(prob, Optim.BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [Inf])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-5.0], ucons = [10.0])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [Inf], lb = [-500.0,-500.0], ub=[50.0,50.0])
sol = solve(prob, IPNewton())
@test sol.minimum < l1

function con2_c(x,p)
    [x[1]^2 + x[2]^2, x[2]*sin(x[1])-x[1]]
end

optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= con2_c)
prob = OptimizationProblem(optprob, x0, lcons = [-Inf,-Inf], ucons = [Inf,Inf])
sol = solve(prob, IPNewton())
@test 10*sol.minimum < l1

cons_circ = (x,p) -> [x[1]^2 + x[2]^2]
optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons_circ)
prob = OptimizationProblem(optprob, x0, lcons = [-Inf], ucons = [0.25^2])
sol = solve(prob, IPNewton())
@test sqrt(cons(sol.u,nothing)[1]) ≈ 0.25 rtol = 1e-6

optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, Flux.ADAM(), maxiters = 1000, progress = false)
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Optim.Fminbox())
@test 10*sol.minimum < l1

Random.seed!(1234)
prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Optim.SAMIN())
@test 10*sol.minimum < l1

optprob = OptimizationFunction((x,p) -> -rosenbrock(x,p), GalacticOptim.AutoZygote())
prob = OptimizationProblem(optprob, x0; sense = GalacticOptim.MaxSense)

import Ipopt
sol = solve(prob, Ipopt.Optimizer())
@test 10*sol.minimum < l1

optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
prob = OptimizationProblem(optprob, x0, _p; sense = GalacticOptim.MinSense)

sol = solve(prob, Ipopt.Optimizer())
@test 10*sol.minimum < l1

sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(Ipopt.Optimizer, "max_cpu_time" => 60.0))
@test 10*sol.minimum < l1

import NLopt
prob = OptimizationProblem(optprob, x0)
sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
@test 10*sol.minimum < l1

sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LN_BOBYQA))
@test 10*sol.minimum < l1

sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LD_LBFGS))
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LD_LBFGS))
@test 10*sol.minimum < l1

sol = solve(prob, NLopt.Opt(:G_MLSL_LDS, 2), local_method = NLopt.Opt(:LD_LBFGS, 2), maxiters=10000)
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, NLopt.LN_BOBYQA())
@test 10*sol.minimum < l1

sol = solve(prob, NLopt.LD_LBFGS())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, NLopt.LD_LBFGS())
@test 10*sol.minimum < l1

sol = solve(prob, NLopt.G_MLSL_LDS(), local_method = NLopt.LD_LBFGS(), local_maxiters=10000, maxiters=10000, population=10)
@test 10*sol.minimum < l1

# using MultistartOptimization
# sol = solve(prob, MultistartOptimization.TikTak(100), local_method = NLopt.LD_LBFGS())
# @test 10*sol.minimum < l1

# using QuadDIRECT
# sol = solve(prob, QuadDirect(); splits = ([-0.5, 0.0, 0.5],[-0.5, 0.0, 0.5]))
# @test 10*sol.minimum < l1

using Evolutionary
prob = GalacticOptim.OptimizationProblem(optprob, x0)
sol = solve(prob, CMAES(μ =40 , λ = 100),abstol=1e-15)
@test 10*sol.minimum < l1

using BlackBoxOptim
prob = GalacticOptim.OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited())
@test 10*sol.minimum < l1

using Metaheuristics
prob = GalacticOptim.OptimizationProblem(f, x0, lb=[-1.0, -1.0], ub=[1.5, 1.5])
sol = solve(prob, ECA())
@test 10*sol.minimum < l1

sol = solve(prob, Metaheuristics.DE())
@test 10*sol.minimum < l1

sol = solve(prob, PSO())
@test 10*sol.minimum < l1

sol = solve(prob, ABC())
@test 10*sol.minimum < l1

sol = solve(prob, CGSA())
@test 10*sol.minimum < l1

sol = solve(prob, SA())
@test 10*sol.minimum < l1

sol = solve(prob, WOA())
@test 10*sol.minimum < l1

sol = solve(prob, ECA())
@test 10*sol.minimum < l1

sol = solve(prob, Metaheuristics.DE(), use_initial=true)
@test 10*sol.minimum < l1

sol = solve(prob, PSO(), use_initial=true)
@test 10*sol.minimum < l1

sol = solve(prob, ABC(), use_initial=true)
@test 10*sol.minimum < l1

sol = solve(prob, CGSA(), use_initial=true)
@test 10*sol.minimum < l1

sol = solve(prob, SA(), use_initial=true)
@test 10*sol.minimum < l1

sol = solve(prob, WOA(), use_initial=true)
@test 10*sol.minimum < l1

using ModelingToolkit
rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p = [1.0, 100.0]

f = OptimizationFunction(rosenbrock,ModelingToolkit.AutoModelingToolkit(),x0,_p,grad=true,hess=true)
prob = OptimizationProblem(f,x0,_p)
sol = solve(prob,Optim.Newton())

### Nonconvex test
prob = GalacticOptim.OptimizationProblem(f, x0, _p, lb = [-1.0,-1.0], ub = [1.5,1.5])

sol = solve(prob, MMA02())
@test 10*sol.minimum < l1

sol = solve(prob, GCMMA())
@test 10*sol.minimum < l1

sol = solve(prob, IpoptAlg())
@test 10*sol.minimum < l1

sol = solve(prob, NLoptAlg(:LN_NELDERMEAD))
@test 10*sol.minimum < l1

sol = solve(prob, NLoptAlg(:LD_LBFGS))
@test 10*sol.minimum < l1

sol = solve(prob, MTSAlg())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler())
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler())
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

sol = solve(prob, MMA02(), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, GCMMA(), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, IpoptAlg(), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, NLoptAlg(:LN_NELDERMEAD), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, NLoptAlg(:LD_LBFGS), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, MTSAlg(), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
@test 10*sol.minimum < l1

### suboptions
sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1

sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
@test 10*sol.minimum < l1


text:
using Bridge, StaticArrays, Distributions, PyPlot
using Test
import Base.Math.gamma
#import Bridge: b, σ, a, transitionprob
const percentile = 3.0
const SV = SVector{2,Float64}
const SM = SMatrix{2,2,Float64,4}
kernel(x, a=0.001) = exp(Bridge.logpdfnormal(x, a*I))

TEST = false
CLASSIC = false

@inline _traceB(t, K, P) = trace(Bridge.B(t, P))

traceB(tt, u::T, P) where {T} = solve(Bridge.R3(), _traceB, tt, u, P)

runmean(x, cx = cumsum(x)) = [cx[n]/n for n in 1:length(x)]

using Bridge: outer

# Define a diffusion process
if ! @isdefined(Target)
struct Target  <: ContinuousTimeProcess{SV}
    c::Float64
    κ::Float64
end
end

if ! @isdefined(Linear)
struct Linear  <: ContinuousTimeProcess{SV}
    T::Float64
    v::SV
    b11::Float64
    b21::Float64
    b12::Float64
    b22::Float64
end
end

g(t, x) = sin(x)
gamma(t, x) = 1.2 - sech(x)/2

# define drift and sigma of Target

Bridge.b(t, x, P::Target) = SV(P.κ*x[2] - P.c*x[1],  -P.c*x[2] + g(t, x[2]))::SV
Bridge.σ(t, x, P::Target) =  SM(0.5, 0.0, 0.0, gamma(t, x[2]))
Bridge.a(t, x, P::Target) = SM(0.25, 0, 0, outer(gamma(t, x[2])))
Bridge.constdiff(::Target) = false

# define drift and sigma of Linear approximation

Bridge.b(t, x, P::Linear) = SV(P.b11*x[1] + P.b12*x[2], P.b21*x[1] + P.b22*x[2] + g(P.T, P.v[2]))
Bridge.B(t, P::Linear) = SM(P.b11, P.b21, P.b12, P.b22)

Bridge.β(t, P::Linear) = SV(0, g(P.T, P.v[2]))

Bridge.σ(t, x, P::Linear) = SM(0.5, 0, 0, gamma(P.T, P.v[2]))
Bridge.a(t, x, P::Linear) = SM(0.25, 0, 0, outer(gamma(P.T, P.v[2])))
Bridge.a(t, P::Linear) = SM(0.25, 0, 0, outer(gamma(P.T, P.v[2])))
Bridge.constdiff(::Linear) = false

c = 0.0
κ = 3.0

t = 1.0
T = 1.5
n = 401
dt = (T-t)/(n-1)
tt = t:dt:T
m = 200_000

u = @SVector [0.1, 0.1]
v = @SVector [0.3, -0.6]

P = Target(c, κ)
Pt = Linear(T, v, -c-0.1, -0.1, κ-0.1, -c/2)

B = Bridge.B(0, Pt)
β = Bridge.β(0, Pt)
a = Bridge.a(0, Pt)
σ = sqrt(Bridge.a(0, Pt))

Phi = exp(B*(T-t))
Λ = lyap(B, a)

if CLASSIC
    Pt2 = LinPro(B, -B\β, σ)
    lpt2 = lp(t, u, T, v, Pt2)
end

GP = Bridge.GuidedBridge(tt, (u,v), P, Pt)
lpt = Bridge.logpdfnormal(v - Bridge.gpmu(tt, u, Pt), Hermitian(Bridge.gpK(tt, zero(SM), Pt)))

if TEST
    @test norm(Bridge.a(0,0, Pt) - Bridge.a(0,0, Pt2)) < sqrt(eps())

    @test norm(Bridge.K(t, T, Pt2) - Bridge.gpK(tt, zero(SM), Pt)) < 1e-6
    # norm(Bridge.gpmu(tt, u, Pt) - Bridge.mu(t, u, T,  Pt2))

    @test norm(Phi*u + sum(exp(B*(T-t))*β*dt for t in tt) - Bridge.gpmu(tt, u, Pt)) < 5e-2

    @test norm(Bridge.gpmu(tt, u, Pt) - ( Phi*u + (Phi-I)*(B\β))) < 1e-6

    @test norm(Bridge.mu(t, u, T,  Pt2) - ( Phi*u + (Phi-I)*(B\β))) < 1e-6
end

W = sample(tt, Wiener{SV}())

# Xtilde

Yt = SV[]
Xt = SamplePath(tt, zeros(SV, length(tt)))
Xts = SamplePath(tt, zeros(SV, length(tt)))
best = Inf
for i in 1:m
    W = sample!(W, Wiener{SV}())
    Bridge.solve!(Euler(), Xt, u, W, Pt)
    push!(Yt, Xt.yy[end])
    nrm = norm(v-Xt.yy[end])

    if nrm < best
        best = nrm
        Xts.yy .= Xt.yy
    end
end

# Target

Y = SV[]
X = SamplePath(tt, zeros(SV, length(tt)))
Xs = SamplePath(tt, zeros(SV, length(tt)))
best = Inf
for i in 1:m
    W = sample!(W, Wiener{SV}())
    Bridge.solve!(Euler(), X, u, W, P)
    push!(Y, X.yy[end])
    nrm = norm(v-X.yy[end])

    if nrm < best
        best = nrm
        Xs.yy .= X.yy
    end
end

# Proposal
lpthat = log(mean(kernel.(collect(y - v for y in Yt))))

Z = Float64[]
Xo = SamplePath(tt, zeros(SV, length(tt)))
@time for i in 1:m
    W = sample!(W, Wiener{SV}())
    Bridge.bridge!(Bridge.Euler(), Xo, W, GP)
    z = llikelihood(LeftRule(), Xo, GP) + lpt
    push!(Z, z)
end

# compare
if CLASSIC

    GP2 = GuidedProp(P, tt[1], u, tt[end], v, Pt2)

    Bridge.H(t, x, P::GuidedProp) = Bridge.H(t, P.t1, Pt)

    Z2 = Float64[]
    Xo2 = SamplePath(tt, zeros(SV, length(tt)))
    @time for i in 1:m
        W = sample!(W, Wiener{SV}())
        Bridge.bridge!(Xo2, W, GP2)
        z = llikelihood(Xo2, GP2) + lpt
        push!(Z2, z)
    end

end

if TEST
    @test norm(Bridge.bridge!(Xo2, W, GP2).yy  - Bridge.bridge!(Bridge.Euler(), Xo, W, GP).yy) < 1e-3
    @test norm(llikelihood(Xo2, GP2) - llikelihood(LeftRule(), Xo, GP)) < 1e-3

    # Some tests

    @test norm(lpt2 - lpt) < 10*sqrt(eps())/dt

    @test norm(mean(Yt)[2] - Bridge.mu(t, u, T, Pt2)[2]) < 1.9*std(last.(Yt))/sqrt(m)
end

Ytv = collect(y - v for y in Yt)
Yv = collect(y - v for y in Y)

@show mean(kernel.(Ytv))
@show exp(lpt)
@show std(kernel.(Ytv))

@show mean(kernel.(Yv))
@show std(kernel.(Yv))

@show mean(exp.(Z))
@show std(exp.(Z))

# Target

figure()
subplot(211)
X = SamplePath(tt, zeros(SV, length(tt)))
for i in 1:10
    W = sample!(W, Wiener{SV}())
    Bridge.solve!(Euler(), X, u, W, P)
    display(plot(X.tt, X.yy))
end
subplot(212)
plot(first.(Y[1:1000]), last.(Y[1:1000]), ".")
plot(v[1], v[2], "o")

figure()
subplot(411)
plot(Xs.tt, Xs.yy, label="X*")
legend()

subplot(412)
plot(Xts.tt, Xts.yy, label="Xt*")
legend()

subplot(413)
plot(Xo.tt, Xo.yy, label="Xo")
legend()

subplot(414)
step = 10
plot(runmean(exp.(Z))[1:step:end], label="Xo")
plot(runmean(kernel.(Yv))[1:step:end], label="X")
plot(runmean(kernel.(Ytv))[1:step:end], label="Xt")
legend()
axis([1, div(m,step), 0, 2*exp(lpt)])

r = Bridge.ri(n-1,Xo.yy[end-1], GP)
println( (Bridge.b(Xo[end-1]..., P)-Bridge.b(Xo[end-1]..., Pt))'*r)
println(trace((Bridge.a(Xo[end-1]..., P)-a)*inv(GP.K[end-1])))
println(r'*(Bridge.a(Xo[end-1]..., P)-a)*r)

ex = Dict(
"u" => u,
"v" => v,
"Xo" => Xo,
"Xs" => Xs,
"Xts" => Xts,
"Yt" => Yt,
"Y" => Y,
"Z" => Z,
"lpt" => lpt
)


text:
module JungleHelperDragonfly

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Dragonfly" Dragonfly(x::Integer, y::Integer, wingsColor::String="FFFFFF", sprite::String="")

const placements = Ahorn.PlacementDict(
    "Dragonfly (Jungle Helper)" => Ahorn.EntityPlacement(
        Dragonfly
    )
)

function Ahorn.selection(entity::Dragonfly)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x, y, 32, 16)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Dragonfly, room::Maple.Room)
    color = Ahorn.argb32ToRGBATuple(parse(Int, entity.wingsColor, base=16))[1:3] ./ 255
	color = (color..., 1.0)
    Ahorn.drawSprite(ctx, "JungleHelper/Dragonfly/dragonfly_body00", 0, 0, jx=0, jy=0)
    Ahorn.drawSprite(ctx, "JungleHelper/Dragonfly/dragonfly_wings00", 0, 0, jx=0, jy=0, tint=color)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module MultiFieldFEFunctionsTests

using BlockArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using Gridap.CellData
using Gridap.MultiField
using Gridap.ReferenceFEs
using Test

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)
degree = order
quad = CellQuadrature(trian,degree)

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:H1)
Q = TestFESpace(model,ReferenceFE(lagrangian,Float64,order-1),conformity=:L2)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

free_values = rand(num_free_dofs(X))
xh = FEFunction(X,free_values)
test_fe_function(xh)
uh, ph = xh

cell_values = get_cell_dof_values(xh)
@test isa(cell_values[1],BlockArrayCoo)

end # module


text:
# -------------- #
# B-Spline Basis #
# -------------- #

# constructor to accept spline params arguments, do some pre-processing
function Basis(::Spline, breaks::Vector, evennum::Int, k::Int)
    # error handling
    k < 0 && error("spline order must be positive")
    length(breaks) < 2 && error("Must have at least two breakpoints")
    any(diff(breaks) .< 0) && error("Breakpoints must be non-decreasing")

    if evennum == 0  # 43
        if length(breaks) == 2  # 44
            evennum = 2
        end
    else
        if length(breaks) == 2
            breaks = linspace(breaks[1], breaks[2], evennum)
        else
            error("Breakpoint squence must contain 2 values when evennum > 0")
        end
    end

    n = length(breaks) + k - 1
    a = breaks[1]
    b = breaks[end]
    Basis(Spline(), n, a, b, SplineParams(breaks, evennum, k))
end

# define methods for SplineParams type
Basis(p::SplineParams) = Basis(Spline(), p.breaks, p.evennum, p.k)

"""
Construct interpolation nodes, given SplineParams

Note that `p.k - 1` additional nodes will be inserted

##### Arguments

- `p::SplineParams`: `SplineParams` instance

##### Returns

- `x::Vector`: The Vector of 1d interpolation nodes.

"""
function nodes(p::SplineParams)
    breaks, evennum, k = p.breaks, p.evennum, p.k
    a = breaks[1]  # 20
    b = breaks[end]  # 21
    n = length(breaks) + k - 1  # 22
    x = cumsum(vcat(fill(a, k), breaks, fill(b, k)))  # 23
    x = (x[1+k:n+k] - x[1:n]) / k  # 24
    x[1] = a  # 25
    x[end] = b  # 26
    x
end

function derivative_op(p::SplineParams, order=1)
    breaks, evennum, k = p.breaks, p.evennum, p.k

    any(order .> k) && error("Order of differentiation can't be greater than k")

    # 38-40
    n = length(breaks) + k - 1
    kk = max(k - 1, k - order - 1)
    augbreaks = vcat(fill(breaks[1], kk), breaks, fill(breaks[end], kk))

    D = Array(SparseMatrixCSC{Float64,Int64}, abs(order), 1)

    if order > 0  # derivative
        temp = k ./ (augbreaks[k+1:n+k-1] - augbreaks[1:n-1])
        D[1] = spdiagm((-temp, temp), 0:1, n-1, n)

        for i=2:order
            temp = (k+1-i) ./ (augbreaks[k+1:n+k-i] - augbreaks[i:n-1])
            D[i] = spdiagm((-temp, temp), 0:1, n-i, n+1-i)*D[i-1]
        end
    else
        error("not implemented")
    end

    D, SplineParams(breaks, evennum, k-order)
end

function _chk_evalbase(p::SplineParams, x, order)
    breaks, evennum, k = p.breaks, p.evennum, p.k

    # error handling
    k < 0 && error("spline order must be positive")
    !(issorted(breaks)) && error("Breakpoints must be non-decreasing")
    any(order .>= k) && error("Order of differentiation must be less than k")
    size(x, 2) > 1 && error("x must be a column vector")

    m = size(x, 1)  # 54
    minorder = minimum(order)  # 55

    # Augment the breakpoint sequence 57-59
    n = length(breaks)+k-1
    a = breaks[1]
    b = breaks[end]
    augbreaks = vcat(fill(a, k-minorder), breaks, fill(b, k-minorder))

    ind = lookup(augbreaks, x, 3)  # 69

    n, m, minorder, augbreaks, ind
end

evalbase(p::SplineParams, x=nodes(p), order::Int=0) = evalbase(p, x, [order])[1]

"""
Evaluate spline basis matrices for a certain order derivative at x

##### Arguments

- `p::SplineParams`: A `SplineParams` summarizing spline properties
- `x(nodes(p))` : the nodes at which to evaluate the basis matrices
- `order(0)` : The order(s) of derivative for which to evaluate the basis
matrices. `order=0` corresponds to the function itself, negative numbers
correspond to integrals.

##### Returns

- `B::SparseMatrixCSC` : Matrix containing the evaluation of basis functions
at each point in `x`. Each column represents a basis function.
- `x`: Points at which the functions were evaluated

"""
function evalbase(p::SplineParams, x, order::AbstractVector{Int})
    n, m, minorder, augbreaks, ind = _chk_evalbase(p, x, order)

    bas = zeros(m, p.k-minorder+1)  # 73
    bas[:, 1] = 1.0  # 74
    B = Array(SparseMatrixCSC{Float64,Int}, length(order))  # 75

    # 76
    if maximum(order) > 0
        D = derivative_op(p, maximum(order))[1]
    end

    # 77
    if minorder < 0
        I = derivative_op(p, minorder)[1]
    end

    for j=1:p.k-minorder  # 78
        for jj=j:-1:1  # 79
            for ix in eachindex(ind)
                b0 = augbreaks[ind[ix]+jj-j]  # 80
                b1 = augbreaks[ind[ix]+jj]  # 81
                temp = bas[ix, jj] / (b1 - b0)  # 82
                bas[ix, jj+1] = (x[ix] - b0) * temp + bas[ix, jj+1]  # 83
                bas[ix, jj] = (b1-x[ix]) * temp  # 84
            end
        end

        # bas now contains the order `j` spline basis
        ii = findfirst(order, p.k-j)  # 87
        if ii > 0
            # Put values in appropriate columns of a sparse matrix
            r = collect(1:m)''  # 91
            r = r[:, fill(1, p.k-order[ii]+1)]  # 91
            c = collect((order[ii] - p.k:0) - (order[ii] - minorder))'  # 92
            c = c[fill(1, m), :] + ind[:, fill(1, p.k-order[ii] + 1)]  # 93
            B[ii] = sparse(vec(r), vec(c), vec(bas[:, 1:p.k-order[ii]+1]),
                           m, n-order[ii])

            # 96-100
            if order[ii] > 0
                B[ii] = B[ii] * D[order[ii]]
            elseif order[ii] < 0
                B[ii] = B[ii] * I[-order[ii]]
            end
        end
    end

    B
end


text:
module Data

using PyCall: pyimport
using Dates: DateTime, DateFormat
using DataFrames: DataFrame, outerjoin
using Serialization: Serialization
using Scratch: @get_scratch!
using CSV: CSV

# Raw data is expected to be in the rawdata folder with the structure

# rawdata
# └───ESGI_Observations
#     ├───Precipitation_Observations
#     │   ├───Spain
#     │   └───United Kingdom
#     └───Temperature_Observations
#         ├───Spain
#         └───United Kingdom

const OBS_RAWDATA = joinpath(@__DIR__, "..", "rawdata", "ESGI_Observations")
const OBS_PRECIP = joinpath(OBS_RAWDATA, "Precipitation_Observations")
const OBS_TEMP = joinpath(OBS_RAWDATA, "Temperature_Observations")

"""
    convert_obs_dir(dir)

Load observation data in NumPy format and return as a named tuple with fields `(country,
area, datetime, value)`.
"""
function convert_obs_dir(dir)
    NP = pyimport("numpy")
    pattern = r"\\(\w*).hourly_(.......)\.npy$"
    df = DateFormat("uuuyyyy")
    countries = String[]
    areas = String[]
    dates = DateTime[]
    values = Union{Missing,Float64}[]
    for country in readdir(dir)
        countrypath = joinpath(dir, country)
        if !isdir(countrypath)
            continue
        end
        for filename in readdir(countrypath; join=true)
            res = match(pattern, filename)
            if res === nothing
                @warn "File name did not match expected pattern" filename
                continue
            else
                (area, month) = res.captures
            end
            data = replace(NP.load(filename), NaN => missing)'  # NaNs represent missing values
            basedate = DateTime(month, df)
            @assert size(data, 1) == 24  # 24 hours in a day
            interval = range(basedate; step=Hour(1), stop=basedate + Month(1) - Hour(1))
            n = length(interval)
            # Append
            append!(countries, fill(country, n))
            append!(areas, fill(area, n))
            append!(dates, interval)
            append!(values, vec(data)[1:n])  # the data is always stored as 31 days but in months with 30 days the last day is a duplicate of the 31st day of the previous month
        end
    end
    return (country=countries, area=areas, date=dates, value=values)
end

"""
    convert_obs()

Convert temperature and precipitation data from NumPy format into a `DataFrame`.
"""
function convert_obs()
    p = convert_obs_dir(OBS_PRECIP)
    precip = DataFrame((country=p.country, area=p.area, date=p.date, precipitation=p.value))
    t = convert_obs_dir(OBS_TEMP)
    temp = DataFrame((country=t.country, area=t.area, date=t.date, temperature=t.value))
    return outerjoin(precip, temp; on=[:country, :area, :date])
end

"""
    get_obs()

Return temperature and precipitation data as a `DataFrame`. Data is cached in a scratch
space for speed.
"""
function get_obs()
    path = @get_scratch!("data")
    obsfile = joinpath(path, "observations.jls")
    if isfile(obsfile)
        return Serialization.deserialize(obsfile)
    else
        obs = convert_obs()
        Serialization.serialize(obsfile, obs)
        return obs
    end
end

function convert_prices()
    spinach = CSV.File("rawdata/spinach-uk.csv"; dateformat="dd/mm/yyyy")
    return DataFrame((crop="spinach", date=c.date, price=c.price) for c in spinach)
end

function get_prices()
    path = @get_scratch!("data")
    pricefile = joinpath(path, "prices.jls")
    if isfile(pricefile)
        return Serialization.deserialize(pricefile)
    else
        prices = convert_prices()
        Serialization.serialize(pricefile, prices)
        return prices
    end
end

end


text:
module MCMC

export runGibbs

using DataFrames
using CategoricalArrays
using StatsModels
using MixedModels
using Distributions,LinearAlgebra #samplers
using StatsBase
using Printf

include("equations.jl")
include("runTime.jl")
include("samplers.jl")

runGibbs = function(formula,userHints,userData,nChain,nBurn,nThin;blockThese=[],VCV=[],userPedData=[],map=[],genotypes...)
	idY,A,yVec,FE,RE,ME,regionSizes = equations.mme(formula,userHints,userData,blocks=blockThese,path2ped=userPedData,paths2geno=genotypes)
        samplers.runSampler(idY,A,yVec,FE,RE,nChain,nBurn,nThin,VCV,ME,map,regionSizes)
        return(yVec,FE,RE,ME)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 12bae146-12e7-4a78-afac-f5c2d6b86b66
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 7be5652e-3e5a-4199-a099-40f2da28053c
using LinearAlgebra, Arpack, LinearMaps, SparseArrays

# ╔═╡ 4b63c1d9-f043-448f-9e05-911f52d4227d
begin
	using Random
	Random.seed!(421)
	n=100
	A=Matrix(Symmetric(rand(n,n)))
	# Or: A = rand(5,5) |> t -> t + t'
	x=rand(n)
	k=10
end

# ╔═╡ f647d0dd-1fe4-42cf-b55c-38baa12f2db8
md"""
# Symmetric Eigenvalue Decomposition - Lanczos Method

If the matrix $A$ is large and sparse and/or if only some eigenvalues and their eigenvectors are desired, iterative methods are the methods of choice. For example, the power method can be useful to compute the eigenvalue with the largest modulus. The basic operation in the power method is matrix-vector multiplication, and this can be performed very fast if $A$ is sparse. Moreover, $A$ need not be stored in the computer -- the input for the algorithm can be just a function which, given some vector $x$, returns the product $Ax$.

An _improved_ version of the power method, which efficiently computes some eigenvalues (either largest in modulus or near some target value $\mu$) and the corresponding eigenvectors, is the Lanczos method.

For more details, see [I. Slapničar, Symmetric Matrix Eigenvalue Techniques, pp. 55.1-55.25](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

__Prerequisites__

The reader should be familiar with concepts of eigenvalues and eigenvectors, related perturbation theory, and algorithms.

__Competences__

The reader should be able to recognise matrices which warrant use uf Lanczos method, to apply the method and to assess the accuracy of the solution.
"""

# ╔═╡ d5f270bd-94c1-4da6-a8e6-a53337488020
md"""
# Lanczos method

 $A$ is a real symmetric matrix of order $n$.

## Definitions

Given a nonzero vector $x$ and an index $k<n$, the __Krylov matrix__ is defined as

$$
K_k=\begin{bmatrix} x & Ax & A^2 x &\cdots & A^{k-1}x \end{bmatrix}.$$

__Krilov subspace__ is the subspace spanned by the columns of $K_k$.

## Facts

1. The Lanczos method is based on the following observation. If $K_k=XR$ is the $QR$ factorization of the matrix $K_k$, then the $k\times k$ matrix $T=X^T A X$ is tridiagonal. The matrices $X$ and $T$ can be computed by using only matrix-vector products in $O(kn)$ operations.

2. Let $T=Q\Lambda Q^T$ be the EVD of $T$. Then $\lambda_i$ approximate well some of the largest and smallest eigenvalues of $A$, and the columns of the matrix $U=XQ$ approximate the corresponding eigenvectors.

3. As $k$ increases, the largest (smallest) eigenvalues of the matrix $T_{1:k,1:k}$ converge towards some of the largest (smallest) eigenvalues of $A$ (due to the Cauchy interlace property). The algorithm can be redesigned to compute only largest or smallest eigenvalues. Also, by using shift and invert strategy, the method can be used to compute eigenvalues near some specified value. In order to obtain better approximations, $k$ should be greater than the number of required eigenvalues. On the other side, in order to obtain better accuracy and efficacy, $k$ should be as small as possible.

4. The last computed element, $\mu=T_{k+1,k}$, provides information about accuracy:

$$
\begin{aligned}
\|AU-U\Lambda\|_2&=\mu, \\
\|AU_{:,i}-\lambda_i U_{:,i}\|_2&=\mu |Q_{ki}|, \quad  i=1,\ldots,k.
\end{aligned}$$

Further, there are $k$ eigenvalues $\tilde\lambda_1,\ldots,\tilde\lambda_k$ of $A$ such that $|\lambda_i-\tilde\lambda_i|\leq \mu$, and for the corresponding eigenvectors, we have

$$\sin2\Theta(U_{:,i},\tilde U_{:,i}) \leq \frac{2\mu}{\min_{j\neq i} |\lambda_i-\tilde \lambda_j|}.$$

5. In practical implementations, $\mu$ is usually used to determine the index $k$.

6. The Lanczos method has inherent numerical instability in the floating-point arithmetic: since the Krylov vectors are, in fact, generated by the power method, they converge towards an eigenvector of $A$. Thus, as $k$ increases, the Krylov vectors become more and more parallel, and the recursion in the function `Lanczos()` becomes numerically unstable and the computed columns of $X$ cease to be sufficiently orthogonal. This affects both the convergence and the accuracy of the algorithm. For example, several eigenvalues of $T$ may converge towards a simple eigenvalue of $A$ (the, so called, __ghost eigenvalues__).

7. The loss of orthogonality is dealt with by using the __full reorthogonalization__ procedure: in each step, the new ${\bf z}$ is orthogonalized against all previous columns of $X$, that is, in function `Lanczos()`, the formula
```
z=z-Tr.dv[i]*X[:,i]-Tr.ev[i-1]*X[:,i-1]
```
is replaced by
```
z=z-sum(dot(z,Tr.dv[i])*X[:,i]-Tr.ev[i-1]*X[:,i-1]
```

To obtain better orthogonality, the latter formula is usually executed twice. The full reorthogonalization raises the operation count to $O(k^2n)$.

8. The __selective reorthogonalization__ is the procedure in which the current $z$ is orthogonalized against some selected columns of $X$, in order to attain sufficient numerical stability and not increase the operation count too much. The details are very subtle and can be found in the references.

9. The Lanczos method is usually used for sparse matrices. Sparse matrix $A$ is stored in the sparse format in which only values and indices of nonzero elements are stored. The number of operations required to multiply some vector by $A$ is also proportional to the number of nonzero elements.

10. The function `eigs()` implements Lanczos method real for symmetric matrices and more general Arnoldi method for general matrices.
"""

# ╔═╡ 97e0dbf8-b9be-4503-af5d-4cf6e86311eb
function Lanczos(A::Array{T}, x::Vector{T}, k::Int) where T
    n=size(A,1)
    X=Array{T}(undef,n,k)
    dv=Array{T}(undef,k)
    ev=Array{T}(undef,k-1)
    X[:,1]=x/norm(x)
    for i=1:k-1
        z=A*X[:,i]
        dv[i]=X[:,i]⋅z
        # Three-term recursion
        if i==1
            z=z-dv[i]*X[:,i]
        else
            # z=z-dv[i]*X[:,i]-ev[i-1]*X[:,i-1]
            # Full reorthogonalization - once or even twice
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
        end
        μ=norm(z)
        if μ==0
            Tr=SymTridiagonal(dv[1:i-1],ev[1:i-2])
            return eigvals(Tr), X[:,1:i-1]*eigvecs(Tr), X[:,1:i-1], μ
        else
            ev[i]=μ
            X[:,i+1]=z/μ
        end
    end
    # Last step
    z=A*X[:,end]
    dv[end]=X[:,end]⋅z
    z=z-dv[end]*X[:,end]-ev[end]*X[:,end-1]
    μ=norm(z)
    Tr=SymTridiagonal(dv,ev)
    eigvals(Tr), X*eigvecs(Tr), X, μ
end

# ╔═╡ 219ce78b-8bd3-4df3-93df-c08fad30e33f
λ,U,X,μ=Lanczos(A,x,10)

# ╔═╡ 1ef82905-e422-4644-ad00-26c448cb0e3a
# Orthogonality of X
norm(X'*X-I)

# ╔═╡ 595b22a5-4456-41fb-b4d6-861833bc6d47
# Tridiagonalization
X'*A*X

# ╔═╡ b2e762d4-650c-4846-ae30-32614b516fe3
# Residual
norm(A*U-U*Diagonal(λ)), μ

# ╔═╡ bd2b8d59-0525-4db5-99cd-e2f98f735eda
U'*A*U

# ╔═╡ 6791903f-0cf7-4bbc-adfb-ff9e80b373dc
# Orthogonality of U
norm(U'*U-I)

# ╔═╡ 218a2a21-e391-4a90-a96a-71bbb0d2f895
# Full eigenvalue decomposition
λeigen,Ueigen=eigen(A);

# ╔═╡ e23c1bc6-c4d8-4c5d-9a83-2b00c5d93b25
# ?eigs

# ╔═╡ fc00c272-d66f-42be-9964-7a934b32015c
# Lanczos method from Arpack.jl
λeigs,Ueigs=eigs(A; nev=k, which=:LM, ritzvec=true, v0=x)

# ╔═╡ 91eec1a8-413c-4960-bb52-8da2435e9e4a
[λ λeigs λeigen[sortperm(abs.(λeigen),rev=true)[1:k]] ]

# ╔═╡ aa20cf19-2bc8-4831-ae26-fb883614c63f
md"""
We see that `eigs()` computes `k` eigenvalues with largest modulus. What eigenvalues did `Lanczos()` compute?
"""

# ╔═╡ 9ce7a012-723d-4184-8ef3-00f468e61281
for i=1:k
    println(minimum(abs,λeigen.-λ[i]))
end

# ╔═╡ ed7d1400-e6d9-44e2-a633-7f6b3df74272
md"""
Conslusion is that the naive implementation of Lanczos is not enough. However, it is fine, when all eigenvalues are computed. Why?
"""

# ╔═╡ d2f31110-5587-4dd3-a5cf-cc3e046e1ee0
λall,Uall,Xall,μall=Lanczos(A,x,100)

# ╔═╡ 04763c96-3a2b-4092-8969-99c18cc8fd54
# Residual and relative errors
norm(A*Uall-Uall*Diagonal(λall)), norm((λeigen-λall)./λeigen)

# ╔═╡ 95ef94d2-129f-4dbb-b705-9a2c8660e22e
md"""
# Operator version

We can use Lanczos method with operator which, given vector `x`, returns the product `A*x`. We use the function `LinearMap()` from the package [LinearMaps.jl](https://github.com/Jutho/LinearMaps.jl)
"""

# ╔═╡ b8716ffe-8405-4018-bdb9-29c8cd2243dc
# ?LinearMap

# ╔═╡ 3bd1283d-91fa-4fdd-bcf8-946ac83b848c
# Operator from the matrix
C=LinearMap(A)

# ╔═╡ 0f8596c1-ee32-4b0d-b13b-d36fb611c99e
begin
	λC,UC=eigs(C; nev=k, which=:LM, ritzvec=true, v0=x)
	λeigs-λC
end

# ╔═╡ 05cd7b4e-18cf-485c-8113-4855f22ac8a0
md"""
Here is an example of `LinearMap()` with the function.
"""

# ╔═╡ a8aeb3ce-5a67-48c7-9c15-8cc11e7ec39b
f(x)=A*x

# ╔═╡ 86bcfb51-6508-4644-ab3d-ee5d9675e1d6
D=LinearMap(f,n,issymmetric=true)

# ╔═╡ 5c08501f-010a-4d17-96b7-e339b0299262
begin
	λD,UD=eigs(D, nev=k, which=:LM, ritzvec=true, v0=x)
	λeigs-λD
end

# ╔═╡ 8a8d59cb-0a69-412e-948e-8110f04419fe
md"""
# Sparse matrices
"""

# ╔═╡ f63b3be6-f5d4-49c7-b8c7-a03826f21ad4
# ?sprand

# ╔═╡ 7beaf848-ad66-47ff-9411-b4ea94476f38
# Generate a sparse symmetric matrix
C₁=sprand(n,n,0.05) |> t -> t+t'

# ╔═╡ d468bf54-54ad-4f07-87d7-ce3d666471aa
issymmetric(C₁)

# ╔═╡ 0f9a4d39-c00c-44b2-8831-eeca1b34e79c
eigs(C₁; nev=k, which=:LM, ritzvec=true, v0=x)

# ╔═╡ a2384a29-511e-43b5-831d-b5b9c2e85ef0

# ╔═╡ Cell order:
# ╟─12bae146-12e7-4a78-afac-f5c2d6b86b66
# ╟─f647d0dd-1fe4-42cf-b55c-38baa12f2db8
# ╟─d5f270bd-94c1-4da6-a8e6-a53337488020
# ╠═7be5652e-3e5a-4199-a099-40f2da28053c
# ╠═97e0dbf8-b9be-4503-af5d-4cf6e86311eb
# ╠═4b63c1d9-f043-448f-9e05-911f52d4227d
# ╠═219ce78b-8bd3-4df3-93df-c08fad30e33f
# ╠═1ef82905-e422-4644-ad00-26c448cb0e3a
# ╠═595b22a5-4456-41fb-b4d6-861833bc6d47
# ╠═b2e762d4-650c-4846-ae30-32614b516fe3
# ╠═bd2b8d59-0525-4db5-99cd-e2f98f735eda
# ╠═6791903f-0cf7-4bbc-adfb-ff9e80b373dc
# ╠═218a2a21-e391-4a90-a96a-71bbb0d2f895
# ╠═e23c1bc6-c4d8-4c5d-9a83-2b00c5d93b25
# ╠═fc00c272-d66f-42be-9964-7a934b32015c
# ╠═91eec1a8-413c-4960-bb52-8da2435e9e4a
# ╟─aa20cf19-2bc8-4831-ae26-fb883614c63f
# ╠═9ce7a012-723d-4184-8ef3-00f468e61281
# ╟─ed7d1400-e6d9-44e2-a633-7f6b3df74272
# ╠═d2f31110-5587-4dd3-a5cf-cc3e046e1ee0
# ╠═04763c96-3a2b-4092-8969-99c18cc8fd54
# ╟─95ef94d2-129f-4dbb-b705-9a2c8660e22e
# ╠═b8716ffe-8405-4018-bdb9-29c8cd2243dc
# ╠═3bd1283d-91fa-4fdd-bcf8-946ac83b848c
# ╠═0f8596c1-ee32-4b0d-b13b-d36fb611c99e
# ╟─05cd7b4e-18cf-485c-8113-4855f22ac8a0
# ╠═a8aeb3ce-5a67-48c7-9c15-8cc11e7ec39b
# ╠═86bcfb51-6508-4644-ab3d-ee5d9675e1d6
# ╠═5c08501f-010a-4d17-96b7-e339b0299262
# ╟─8a8d59cb-0a69-412e-948e-8110f04419fe
# ╠═f63b3be6-f5d4-49c7-b8c7-a03826f21ad4
# ╠═7beaf848-ad66-47ff-9411-b4ea94476f38
# ╠═d468bf54-54ad-4f07-87d7-ce3d666471aa
# ╠═0f9a4d39-c00c-44b2-8831-eeca1b34e79c
# ╠═a2384a29-511e-43b5-831d-b5b9c2e85ef0

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L4d Symmetric Eigenvalue Decomposition - Lanczos Method.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using Iconv
using Test

@test iconv("笨熊", "UTF-8", "GBK") == togbk("笨熊") == g"笨熊"
@test iconv(g"笨熊", "GBK", "UTF-8") == toutf8(g"笨熊") == b"笨熊"


text:
"""
    ustrip(x::MetadataArray)

Strips unit from `x.parent` but stores it in `x.metadata` for safekeeping.
"""
Unitful.ustrip(x::MetadataArray) = MetadataArray(ustrip.(x.parent), (x.metadata..., unit=unit(eltype(x))))
"""
    upreferred(x::MetadataArray)

Converts `x.parent` to SI unit but keeps `x.metadata` for safekeeping.
"""
function Unitful.upreferred(x::MetadataArray)
    parent = upreferred.(x.parent)
    processing_str = "converted from $(unit(eltype(x.parent))) to preferred $(unit(eltype(parent)))"
    if haskey(x.metadata, :processing)
        metadata = x.metadata
        push!(metadata.processing, processing_str)
    else
        metadata = (x.metadata..., processing=[processing_str])
    end
    return MetadataArray(parent, metadata)
end
"""
    *(x::MetadataArray, q::Quantity)

    Preserve metadata (and history of unit conversions) when multiplying `x` by a quantity.
"""
function Base.:*(x::MetadataArray, q::Quantity)
    parent = x.parent * q
    processing_str = "× $q"
    if haskey(x.metadata, :processing)
        metadata = x.metadata
        push!(metadata.processing, processing_str)
    else
        metadata = (x.metadata..., processing=[processing_str])
    end
    return MetadataArray(parent, metadata)
end

"""
    onlykeep(x::MetadataArray, idx)

Returns `x[idx]` but also applies the index to the metadata that is originally of the same length as `x`.
"""
function onlykeep(x::MetadataArray, idx)
    xout = x.parent[idx]
    mout = (;((v isa Vector && length(v) == length(x)) ? (k,v[idx]) : (k,v) for (k,v) in pairs(x.metadata))...)
    return MetadataArray(xout, mout)
end
export onlykeep

"""
    AgedJacobianFactors

Type containing the Jacobian Factors and the age of the Jacobian.
This allows for the Shamanskii method to not update the Jacobian at each iterate.
"""
mutable struct AgedJacobianFactors{T}
    fac::T     # Jacobian factors — can be Real, Complex, Dual, or HyperDual
    age::Int   # age of the Jacobian
end

"""
    \\(Jf::JacobianFactors, y)

Dispatches backslash to work with all `JacobianFactors` subtypes.
"""
function Base.:\(Jf::AgedJacobianFactors, y)
    return Jf.fac \ y
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

import .Serializer: known_object_data, object_number, serialize_cycle, deserialize_cycle, writetag,
                      __deserialized_types__, serialize_typename, deserialize_typename,
                      TYPENAME_TAG, object_numbers

type ClusterSerializer{I<:IO} <: AbstractSerializer
    io::I
    counter::Int
    table::ObjectIdDict

    sent_objects::Set{UInt64} # used by serialize (track objects sent)

    ClusterSerializer(io::I) = new(io, 0, ObjectIdDict(), Set{UInt64}())
end
ClusterSerializer(io::IO) = ClusterSerializer{typeof(io)}(io)

function deserialize(s::ClusterSerializer, ::Type{TypeName})
    full_body_sent = deserialize(s)
    number = read(s.io, UInt64)
    if !full_body_sent
        tn = get(known_object_data, number, nothing)::TypeName
        if !haskey(object_numbers, tn)
            # set up reverse mapping for serialize
            object_numbers[tn] = number
        end
        deserialize_cycle(s, tn)
    else
        tn = deserialize_typename(s, number)
    end
    return tn
end

function serialize(s::ClusterSerializer, t::TypeName)
    serialize_cycle(s, t) && return
    writetag(s.io, TYPENAME_TAG)

    identifier = object_number(t)
    send_whole = !(identifier in s.sent_objects)
    serialize(s, send_whole)
    write(s.io, identifier)
    if send_whole
        serialize_typename(s, t)
        push!(s.sent_objects, identifier)
    end
#   println(t.module, ":", t.name, ", id:", identifier, send_whole ? " sent" : " NOT sent")
    nothing
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module ParameterLearning

using DataFrames
using LPSjl.BayesianNetworks
include("BaseParameterLearning.jl")

"""
    calculateEmpiricalDistribution(::AbstractDataFrame)

Returns the Empirical Distribution of a dataFrame
"""
function calculateEmpiricalDistribution(df::AbstractDataFrame)
    combinedDf = combine(groupby(df, names(df)), nrow => :p)
    combinedDf.p = combinedDf.p / size(df)[1]

    return sort(combinedDf, names(df))
end

export
    calculateEmpiricalDistribution,
    MaximumLikelihoodEstimator,
    estimate

end  # module ParameterLearning


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Resolve

using ..Versions

using Printf
using Random
using UUIDs

export resolve, sanity_check, Graph, pkgID

####################
# Requires / Fixed #
####################
const Requires = Dict{UUID,VersionSpec}

struct Fixed
    version::VersionNumber
    requires::Requires
end
Fixed(v::VersionNumber) = Fixed(v, Requires())

Base.:(==)(a::Fixed, b::Fixed) = a.version == b.version && a.requires == b.requires
Base.hash(f::Fixed, h::UInt) = hash((f.version, f.requires), h + (0x68628b809fd417ca % UInt))

Base.show(io::IO, f::Fixed) = isempty(f.requires) ?
    print(io, "Fixed(", repr(f.version), ")") :
    print(io, "Fixed(", repr(f.version), ",", f.requires, ")")

struct ResolverError <: Exception
    msg::AbstractString
    ex::Union{Exception,Nothing}
end
ResolverError(msg::AbstractString) = ResolverError(msg, nothing)

function Base.showerror(io::IO, pkgerr::ResolverError)
    print(io, pkgerr.msg)
    if pkgerr.ex !== nothing
        pkgex = pkgerr.ex
        if isa(pkgex, CompositeException)
            for cex in pkgex
                print(io, "\n=> ")
                showerror(io, cex)
            end
        else
            print(io, "\n")
            showerror(io, pkgex)
        end
    end
end

const uuid_julia = UUID("1222c4b2-2114-5bfd-aeef-88e4692bbb3e")

include("versionweights.jl")
include("fieldvalues.jl")
include("graphtype.jl")
include("maxsum.jl")

"Resolve package dependencies."
function resolve(graph::Graph)
    # attempt trivial solution first
    greedy_ok, sol = greedysolver(graph)

    greedy_ok && @goto solved

    log_event_global!(graph, "greedy solver failed")

    # trivial solution failed, use maxsum solver
    maxsum_ok, sol, staged = maxsum(graph)

    maxsum_ok && @goto solved

    log_event_global!(graph, "maxsum solver failed")

    # the problem is unsat, force-trigger a failure
    # in order to produce a log - this will contain
    # information about the best that the solver could
    # achieve
    trigger_failure!(graph, sol, staged)

    @label solved

    # verify solution (debug code) and enforce its optimality
    @assert verify_solution(sol, graph)
    greedy_ok || enforce_optimality!(sol, graph)

    log_event_global!(graph, "the solver found $(greedy_ok ? "an optimal" : "a feasible") configuration")

    # return the solution as a Dict mapping UUID => VersionNumber
    return compute_output_dict(sol, graph)
end

"""
Scan the graph for (explicit or implicit) contradictions. Returns a list of problematic
(package,version) combinations.
"""
function sanity_check(graph::Graph, sources::Set{UUID} = Set{UUID}(), verbose::Bool = true)
    req_inds = graph.req_inds
    fix_inds = graph.fix_inds

    id(p) = pkgID(p, graph)

    isempty(req_inds) || @warn("sanity check called on a graph with non-empty requirements")
    if !any(is_julia(graph, fp0) for fp0 in fix_inds)
        @warn("sanity check called on a graph without julia requirement, adding it")
        add_fixed!(graph, Dict(uuid_julia=>Fixed(VERSION)))
    end
    if length(fix_inds) ≠ 1
        @warn("sanity check called on a graph with extra fixed requirements (besides julia)")
    end

    isources = isempty(sources) ?
        Set{Int}(1:graph.np) :
        Set{Int}(graph.data.pdict[p] for p in sources)

    simplify_graph!(graph, isources)

    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    data = graph.data
    pkgs = data.pkgs
    pdict = data.pdict
    pvers = data.pvers
    eq_classes = data.eq_classes

    problematic = Tuple{String,VersionNumber}[]

    np == 0 && return problematic

    vers = [(pkgs[p0],pvers[p0][v0]) for p0 = 1:np for v0 = 1:(spp[p0]-1)]
    sort!(vers, by=pv->(-length(gadj[pdict[pv[1]]])))

    nv = length(vers)

    svdict = Dict{Tuple{UUID,VersionNumber},Int}(vers[i] => i for i = 1:nv)

    checked = falses(nv)

    last_str_len = 0

    for (i,(p,vn)) in enumerate(vers)
        if verbose
            frac_compl = i / nv
            print("\r", " "^last_str_len, "\r")
            progr_msg = @sprintf("%.3i/%.3i (%i%%) — problematic so far: %i", i, nv, round(Int, 100 * frac_compl), length(problematic))
            print(progr_msg)
            last_str_len = length(progr_msg)
        end

        length(gadj[pdict[p]]) == 0 && break
        checked[i] && continue

        push_snapshot!(graph)

        # enforce package version
        # TODO: use add_reqs! instead...
        p0 = graph.data.pdict[p]
        v0 = graph.data.vdict[p0][vn]
        fill!(graph.gconstr[p0], false)
        graph.gconstr[p0][v0] = true
        push!(graph.req_inds, p0)

        ok = false
        try
            simplify_graph_soft!(graph, Set{Int}([p0]), log_events = false)
        catch err
            isa(err, ResolverError) || rethrow()
            @goto done
        end

        ok, sol = greedysolver(graph)
        ok && @goto done
        ok, sol = maxsum(graph)

        @label done

        if !ok
            for vneq in eq_classes[p][vn]
                push!(problematic, (id(p), vneq))
            end
        else
            @assert verify_solution(sol, graph)
            sol_dict = compute_output_dict(sol, graph)
            for (sp,svn) in sol_dict
                j = svdict[sp,svn]
                checked[j] = true
            end
        end

        # state reset
        empty!(graph.req_inds)
        pop_snapshot!(graph)
    end
    if verbose
        print("\r", " "^last_str_len, "\r")
        println("found $(length(problematic)) problematic versions")
    end
    return sort!(problematic)
end

"""
Translate the solver output (a Vector{Int} of package states) into a Dict which
associates a VersionNumber to each installed package UUID.
"""
function compute_output_dict(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    fix_inds = graph.fix_inds
    pkgs = graph.data.pkgs
    pvers = graph.data.pvers
    pruned = graph.data.pruned

    want = Dict{UUID,VersionNumber}()
    for p0 = 1:np
        p0 ∈ fix_inds && continue
        p = pkgs[p0]
        s0 = sol[p0]
        s0 == spp[p0] && continue
        vn = pvers[p0][s0]
        want[p] = vn
    end
    for (p,vn) in pruned
        @assert !haskey(want, p)
        want[p] = vn
    end

    return want
end

"""
Preliminary solver attempt: tries to maximize each version; bails out as soon as
some non-trivial requirement is detected.
"""
function greedysolver(graph::Graph)
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    np = graph.np

    push_snapshot!(graph)
    gconstr = graph.gconstr

    # initialize solution: all uninstalled
    sol = [spp[p0] for p0 = 1:np]

    # packages which are not allowed to be uninstalled
    # (NOTE: this is potentially a superset of graph.req_inds,
    #        since it may include implicit requirements)
    req_inds = Set{Int}(p0 for p0 = 1:np if !gconstr[p0][end])

    # set up required packages to their highest allowed versions
    for rp0 in req_inds
        # look for the highest version which satisfies the requirements
        rv0 = findlast(gconstr[rp0])
        @assert rv0 ≢ nothing && rv0 ≠ spp[rp0]
        sol[rp0] = rv0
        fill!(gconstr[rp0], false)
        gconstr[rp0][rv0] = true
    end

    # propagate the requirements
    try
        simplify_graph_soft!(graph, req_inds, log_events = false)
    catch err
        err isa ResolverError || rethrow()
        pop_snapshot!(graph)
        return (false, Int[])
    end

    # we start from required packages and explore the graph
    # following dependencies
    staged = req_inds
    seen = copy(staged)

    while !isempty(staged)
        staged_next = Set{Int}()
        for p0 in staged
            s0 = sol[p0]
            @assert s0 < spp[p0]

            # scan dependencies
            for (j1,p1) in enumerate(gadj[p0])
                msk = gmsk[p0][j1]
                # look for the highest version which satisfies the requirements
                v1 = findlast(msk[:,s0] .& gconstr[p1])
                v1 == spp[p1] && continue # p1 is not required by p0's current version
                # if we found a version, and the package was uninstalled
                # or the same version was already selected, we're ok;
                # otherwise we can't be sure what the optimal configuration is
                # and we bail out
                old_v1 = sol[p1]
                if v1 ≡ nothing || (old_v1 ≠ v1 && old_v1 ≠ spp[p1])
                    pop_snapshot!(graph)
                    return (false, Int[])
                elseif old_v1 == spp[p1]
                    sol[p1] = v1
                    push!(staged_next, p1)
                end
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    pop_snapshot!(graph)

    for p0 = 1:np
        log_event_greedysolved!(graph, p0, sol[p0])
    end

    return true, sol
end

"""
Verifies that the solver solution fulfills all hard constraints
(requirements and dependencies). This is intended as debug code.
"""
function verify_solution(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    gconstr = graph.gconstr

    @assert length(sol) == np
    @assert all(sol .> 0)

    # verify constraints and dependencies
    for p0 = 1:np
        s0 = sol[p0]
        gconstr[p0][s0] || (@warn("gconstr[$p0][$s0] fail"); return false)
        for (j1,p1) in enumerate(gadj[p0])
            msk = gmsk[p0][j1]
            s1 = sol[p1]
            msk[s1,s0] || (@warn("gmsk[$p0][$p1][$s1,$s0] fail"); return false)
        end
    end
    return true
end

"""
Uninstall unreachable packages:
start from the required ones and keep only the packages reachable from them along the graph.
"""
function _uninstall_unreachable!(sol::Vector{Int}, why::Vector{Union{Symbol,Int}}, graph::Graph)
    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    gconstr = graph.gconstr

    uninst = trues(np)
    staged = Set{Int}(p0 for p0 = 1:np if !gconstr[p0][end])
    seen = copy(staged) ∪ Set{Int}(p0 for p0 = 1:np if sol[p0] == spp[p0]) # we'll skip uninstalled packages

    while !isempty(staged)
        staged_next = Set{Int}()
        for p0 in staged
            s0 = sol[p0]
            @assert s0 ≠ spp[p0]
            uninst[p0] = false
            for (j1,p1) in enumerate(gadj[p0])
                p1 ∈ seen && continue            # we've already seen the package, or it is uninstalled
                gmsk[p0][j1][end,s0] && continue # the package is not required by p0 at version s0
                push!(staged_next, p1)
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    for p0 in findall(uninst)
        sol[p0] = spp[p0]
        why[p0] = :uninst
    end
end

"""
Push the given solution to a local optimum if needed: keeps increasing
the states of the given solution as long as no constraints are violated.
It might also install additional packages, if needed to bump the ones already
installed.
It also removes unnecessary parts of the solution which are unconnected
to the required packages.
"""
function enforce_optimality!(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    gconstr = graph.gconstr
    pkgs = graph.data.pkgs

    # keep a track for the log
    why = Union{Symbol,Int}[0 for p0 = 1:np]

    # Strategy:
    # There's a cycle in which first the unnecessary (unconnected) packages are removed,
    # then we make a pass over the whole packages trying to bump each of them.
    # We repeat the above two steps until no further action is allowed.
    # When attempting to bump a package, we may attempt to bump or install other packages
    # if needed. Except if the bump would uninstall a package, in which cases we don't
    # touch anything else: we do it only if it has no consequence at all. This strategy
    # favors installing packages as needed.
    # During the bumping pass, we keep an upper and lower bound for each package, which
    # progressively shrink. These are used when adjusting for the effect of an attempted bump.
    # The way it's written should ensure that no package is ever downgraded (unless it was
    # originally unneeded, and then got removed, and later reinstalled to a lower version as
    # a consequence of a bump of some other package).

    # move_up is used to keep track of which packages can move up
    # (they start installed and can be bumped) and which down (they start uninstalled and
    # can be installed)
    move_up = BitVector(undef, length(sol))
    # lower and upper bounds on the valid range of each package
    upperbound = similar(spp)
    lowerbound = similar(spp)
    # backup space for restoring the state if an attempted bump fails
    bk_sol = similar(sol)
    bk_lowerbound = similar(lowerbound)
    bk_upperbound = similar(upperbound)

    # auxiliary sets to perform breadth-first search on the graph
    staged = Set{Int}()
    staged_next = Set{Int}()

    old_sol = similar(sol)       # to detect if we made any changes
    allsols = Set{Vector{Int}}() # used to make 100% sure we avoid infinite loops

    while true
        copy!(old_sol, sol)
        push!(allsols, copy(sol))

        # step 1: uninstall unneded packages
        _uninstall_unreachable!(sol, why, graph)

        # setp 2: try to bump each installed package in turn

        move_up .= sol .≠ spp
        copy!(upperbound, spp)
        let move_up = move_up
            lowerbound .= [move_up[p0] ? sol[p0] : 1 for p0 = 1:np]
        end

        for p0 = 1:np
            s0 = sol[p0]
            s0 == spp[p0] && (why[p0] = :uninst; continue) # the package is not installed
            move_up[p0] || continue # the package is only installed as a result of a previous bump, skip it

            @assert upperbound[p0] == spp[p0]

            # pick the next version that doesn't violate a constraint (if any)
            bump_range = collect(s0+1:spp[p0])
            bump = let gconstr = gconstr
                findfirst(v0->gconstr[p0][v0], bump_range)
            end

            # no such version was found, skip this package
            bump ≡ nothing && (why[p0] = :constr; continue)

            # assume that we will succeed in bumping the version (otherwise we
            # roll-back at the end)

            new_s0 = bump_range[bump]
            try_uninstall = new_s0 == spp[p0] # are we trying to uninstall a package?

            copy!(bk_sol, sol)
            copy!(bk_lowerbound, lowerbound)
            copy!(bk_upperbound, upperbound)
            sol[p0] = new_s0

            # if we're trying to uninstall, the bump is "soft": we don't update the
            # lower bound so that the package can be reinstalled later in the pass
            # if needed by another package
            try_uninstall || (lowerbound[p0] = new_s0) # note that we're in the move_up case

            empty!(staged)
            empty!(staged_next)
            push!(staged, p0)

            while !isempty(staged)
                for f0 in staged
                    for (j1,f1) in enumerate(gadj[f0])
                        s1 = sol[f1]
                        msk = gmsk[f0][j1]
                        if f1 == p0 || try_uninstall
                            # when uninstalling or looking at p0, no further changes are allowed
                            bump_range = [s1]
                        else
                            lb1 = lowerbound[f1]
                            ub1 = upperbound[f1]
                            @assert lb1 ≤ s1 ≤ ub1
                            if move_up[f1]
                                s1 > lb1 && @assert s1 == spp[f1]
                                # the arrangement of the range gives precedence to improving the
                                # current situation, but allows reinstalling a package if needed
                                bump_range = vcat(s1:ub1, s1-1:-1:lb1)
                            else
                                bump_range = collect(ub1:-1:lb1)
                            end
                        end
                        bump = let gconstr = gconstr
                            findfirst(v1->(gconstr[f1][v1] && msk[v1, sol[f0]]), bump_range)
                        end
                        if bump ≡ nothing
                            why[p0] = f1 # TODO: improve this? (ideally we might want the path from p0 to f1)
                            @goto abort
                        end
                        new_s1 = bump_range[bump]
                        sol[f1] = new_s1
                        new_s1 == s1 && continue
                        push!(staged_next, f1)
                        if move_up[f1]
                            lowerbound[f1] = new_s1
                        else
                            upperbound[f1] = new_s1
                        end
                    end
                end
                staged, staged_next = staged_next, staged
                empty!(staged_next)
            end

            # if we're here the bump was successful, there's nothing more to do
            continue

            ## abort the bumping: restore the solution
            @label abort

            copy!(sol, bk_sol)
            copy!(lowerbound, bk_lowerbound)
            copy!(upperbound, bk_upperbound)
        end
        sol ≠ old_sol || break
        # It might be possible in principle to contrive a situation in which
        # the solutions oscillate
        sol ∈ allsols && break
    end

    @assert verify_solution(sol, graph)

    for p0 = 1:np
        log_event_maxsumsolved!(graph, p0, sol[p0], why[p0])
    end
end

function apply_maxsum_trace!(graph::Graph, sol::Vector{Int})
    gconstr = graph.gconstr

    for (p0,s0) in enumerate(sol)
        s0 == 0 && continue
        gconstr0 = gconstr[p0]
        old_constr = copy(gconstr0)
        @assert old_constr[s0]
        fill!(gconstr0, false)
        gconstr0[s0] = true
        gconstr0 ≠ old_constr && log_event_maxsumtrace!(graph, p0, s0)
    end
end

function trigger_failure!(graph::Graph, sol::Vector{Int}, staged::Tuple{Int,Int})
    apply_maxsum_trace!(graph, sol)
    simplify_graph_soft!(graph, Set(findall(sol .> 0)), log_events = true) # this may throw an error...

    np = graph.np
    gconstr = graph.gconstr
    p0, v0 = staged

    @assert gconstr[p0][v0]
    fill!(gconstr[p0], false)
    gconstr[p0][v0] = true
    log_event_maxsumtrace!(graph, p0, v0)
    simplify_graph!(graph) # this may throw an error...
    outdict = resolve(graph) # ...otherwise, this MUST throw an error
    open(io->showlog(io, graph, view=:chronological), "logchrono.errresolve.txt", "w")
    error("this is not supposed to happen... $(Dict(pkgID(p, graph) => vn for (p,vn) in outdict))")
end

end # module


text:
# TODO: find a way to make this a parameter.
_STRICTLY_FLOAT_EPSILON = 1.0e-9

"""
Bridges `CP.Strictly` to linear constraints.
"""
struct Strictly2LPBridge{T} <: MOIBC.AbstractBridge
    con::Union{
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}},
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}},
    }
end

function MOIBC.bridge_constraint(
    ::Type{Strictly2LPBridge{T}},
    model,
    f::MOI.VariableIndex,
    s::Union{CP.Strictly{MOI.LessThan{T}, T}, CP.Strictly{MOI.GreaterThan{T}, T}},
) where {T}
    return MOIBC.bridge_constraint(
        Strictly2LPBridge{T},
        model,
        MOI.ScalarAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Strictly2LPBridge{T}},
    model,
    f::MOI.ScalarAffineFunction{T},
    s::CP.Strictly{MOI.LessThan{T}, T},
) where {T <: Real}
    con = MOI.add_constraint(
        model,
        f,
        MOI.LessThan(s.set.upper - _STRICTLY_FLOAT_EPSILON)
    )

    return Strictly2LPBridge(con)
end

function MOIBC.bridge_constraint(
    ::Type{Strictly2LPBridge{T}},
    model,
    f::MOI.ScalarAffineFunction{T},
    s::CP.Strictly{MOI.LessThan{T}, T},
) where {T <: Integer}
    con = MOI.add_constraint(
        model,
        f,
        MOI.LessThan(s.set.upper - one(T))
    )

    return Strictly2LPBridge(con)
end

function MOIBC.bridge_constraint(
    ::Type{Strictly2LPBridge{T}},
    model,
    f::MOI.ScalarAffineFunction{T},
    s::CP.Strictly{MOI.GreaterThan{T}, T},
) where {T <: Real}
    con = MOI.add_constraint(
        model,
        f,
        MOI.GreaterThan(s.set.lower + _STRICTLY_FLOAT_EPSILON)
    )

    return Strictly2LPBridge(con)
end

function MOIBC.bridge_constraint(
    ::Type{Strictly2LPBridge{T}},
    model,
    f::MOI.ScalarAffineFunction{T},
    s::CP.Strictly{MOI.GreaterThan{T}, T},
) where {T <: Integer}
    con = MOI.add_constraint(
        model,
        f,
        MOI.GreaterThan(s.set.lower + one(T))
    )

    return Strictly2LPBridge(con)
end

function MOI.supports_constraint(
    ::Type{Strictly2LPBridge{T}},
    ::Union{Type{MOI.VariableIndex}, Type{MOI.ScalarAffineFunction{T}}},
    ::Union{Type{CP.Strictly{MOI.LessThan{T}, T}}, Type{CP.Strictly{MOI.GreaterThan{T}, T}}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{Strictly2LPBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{Strictly2LPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
    ]
end

function MOI.get(
    b::Strictly2LPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return (b.con isa MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}) ? 1 : 0
end

function MOI.get(
    b::Strictly2LPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return (b.con isa MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}) ? 1 : 0
end

function MOI.get(
    b::Strictly2LPBridge{T},
    ::Union{
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
        },
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
        },
    }
) where {T}
    return [b.con]
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()
using PackageCompiler
create_sysimage(
  :MocosSimLauncher,
  sysimage_path=joinpath(@__DIR__, "..", "sysimage.img"),
  precompile_execution_file=joinpath(@__DIR__, "run_example.jl"),
  #precompile_statements_file=joinpath(@__DIR__, "precompile_statements.jl"),
  project=joinpath(@__DIR__, ".."),
  filter_stdlibs=true,
  incremental=false,
  )


text:
using NNPACK, Test, Random

Random.seed!(144)

@testset "activation" begin
    xs = rand(5,5)

    @test all(sum(softmax(xs), dims = 1) .≈ Float32(1))

    @test sum(softmax(vec(xs))) ≈ Float32(1)
end

@testset "conv2d" begin

    x = reshape(Float32[1:16;], 4, 4, 1, 1)
    w = reshape(Float32[1:9;], 3, 3, 1, 1)

    @test dropdims(conv(x, w, pad=1), dims=(3,4)) ≈ Float32.([
        29   99  207  263
        62  192  372  446
        83  237  417  485
        75  198  330  365])

    x = reshape(Float64[1:20;], 5, 4, 1, 1)
    w = reshape(Float64[1:4;], 2, 2, 1, 1)

    @test dropdims(conv(x, w), dims = (3,4)) ≈ Float32.([
        29  79 129;
        39  89 139;
        49  99 149;
        59 109 159
    ])

    @test_throws Exception dropdims(conv(x, w; stride=2), dims = (3,4)) ≈ Float32.([
        29 49;
        129 149
    ])

    @test dropdims(conv(x, w; pad=1), dims = (3,4)) ≈ Float32.([
        1    9   29   49   48;
        4   29   79  129  115;
        7   39   89  139  122;
        10  49   99  149  129;
        13  59  109  159  136;
        10  40   70  100   80
    ])

    @test size(∇conv_filter(reshape(rand(4,3), 4, 3, 1, 1), x, w)) == size(w)
    @test size(∇conv_data(reshape(rand(4,3), 4, 3, 1, 1), x, w)) == size(x)

end

@testset "maxpool2d" begin

	x = reshape(Float32[1:16;], 4, 4, 1, 1)

	@test dropdims(maxpool(x, (2,2)), dims = (3,4)) ≈ Float32.([
        6.0 14.0;
        8.0 16.0
    ])

    x = reshape(Float64[1:20;], 5, 4, 1, 1)

    @test_throws Exception dropdims(maxpool(x, (2,2)), dims = (3,4)) ≈ Float32.([7 17; 9 19])
    @test_throws Exception dropdims(maxpool(x, (2,2); stride=(2,2)), dims = (3,4)) ≈ Float32.([7 17; 9 19])
    @test_throws Exception dropdims(maxpool(x, (2,2); pad=(1,1)), dims = (3,4)) ≈ Float32.([
        1  11  16;
        3  13  18;
        5  15  20;
    ])

    x = reshape(Float64[1:16;], 4, 4, 1, 1)

    @test dropdims(maxpool(x, (2,2)), dims = (3,4)) == [6 14; 8 16]
    @test dropdims(maxpool(x, (2,2); stride=(2,2)), dims = (3,4)) == [6 14; 8 16]
    @test dropdims(maxpool(x, (3,3); stride=(1,1), pad=(1,1)), dims = (3,4)) == [
         6.0  10.0  14.0  14.0;
         7.0  11.0  15.0  15.0;
         8.0  12.0  16.0  16.0;
         8.0  12.0  16.0  16.0;
    ]

end


text:
# The 5-digit number, 16807=7^5, is also a fifth power. Similarly, the 9-digit
# number, 134217728=8^9, is a ninth power.
#
# How many n-digit positive integers exist which are also an nth power?

using ProjectEulerSolutions

# Brute for approach, cycling through nth powers on single digits until we see
# that no more signal digits have the right number of digits (using the built-in
# ndigits function).
function p063solution_bruteforce(max_n::Integer=10)::Integer

    digs = BigInt.(1:9)
    count = 0
    for n = 1:max_n
        c = sum(ndigits.(digs .^ n) .== n)
        if c == 0
           break
        end
        count += c
    end
    return count
end

# Recognize that a single digit number k will have n digits as long as
# 10^(n-1) <= k^n < 10^n.  We know k^n < 10^n when k is 1:9.  Solve for k
# 10^((n-1) / n) <= k and cycle through n values to see the number of values for
# n where the relation holds.  Vectorize and try with k = 1:9 to solve the
# problem in a single line.
function p063solution_explicit_powers(max_n::Integer=10)::Integer
    return Integer(sum(9 .- floor.(10 .^((0:(max_n-1))./(1:max_n)))) + 1)
end

# Alternatively, solve directly for n: 10^(n-1) = k^n, (n-1) log(10) = n log(k),
# n = log(10) / (log(10) - log(k)), where we pass values for k = 1:9.
function p063solution_explicit_logs(max_n::Integer=10)::Integer
    return Integer(sum(floor.(log.(10)./(log.(10).-log.(1:9)))))
end

p063 = Problems.Problem(Dict("Brute-force" => p063solution_bruteforce,
                             "Explicit powers" => p063solution_explicit_powers,
                             "Explicit logs" => p063solution_explicit_logs))

Problems.benchmark(p063, 30)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using KitBase
using KitBase.Plots
using KitBase.SpecialFunctions
using FluxRC

cd(@__DIR__)
ps = UnstructFRPSpace("naca0012.msh", 3)

scatter(ps.rs[:, 1], ps.rs[:, 2], ratio = 1)
i = 6;
scatter(ps.xp[i, :, 1], ps.xp[i, :, 2], ratio = 1);

using PyCall
const itp = pyimport("scipy.interpolate")
const np = pyimport("numpy")

# Lagrange interpolation: hand-written vs. scipy
## value
xGauss = legendre_point(2)
yGauss = exp.(xGauss)
ll = lagrange_point(xGauss, -1.0)
yGauss .* ll |> sum

poly = itp.lagrange(xGauss, yGauss)
np.polyval(poly, -1.0)

## derivative
lpdm = ∂lagrange(xGauss)
yGauss[1] * lpdm[1, 1] + yGauss[2] * lpdm[1, 2] + yGauss[3] * lpdm[1, 3]

p1 = np.polyder(poly, 1)
np.polyval(p1, xGauss[1])

# 2D Lagrange interpolation
x = legendre_point(2)
y = [-0.9, -0.45, 0.0, 0.45, 0.9]

lx = lagrange_point(x, 0.77)
ly = lagrange_point(y, 0.9)

coords = cat(meshgrid(x, y)[1] |> permutedims, meshgrid(x, y)[2] |> permutedims, dims = 3)
val = exp.(coords[:, :, 1] + coords[:, :, 2] .^ 2)

v = 0.0
for i in axes(weights, 1), j in axes(weights, 2)
    v += val[i, j] * lx[i] * ly[j]
end

lagrange_point([0, 0.5, 0.5, 1], 0.77)

legendre_point(2) .- -1


text:
function check_branches!(branches::Array{<:Branch,1})
    checkanglelimits!(branches)
    check_ascending_order([b.connectionpoints.from.number for b in branches], "Branch")
end

# function check_angle_limits(anglelimits::(max::Float64, min::Float64))
function checkanglelimits!(branches::Array{<:Branch,1})
    for (ix,l) in enumerate(branches)
        if isa(l,Union{MonitoredLine,Line})
            orderedlimits(l.anglelimits, "Angles")

            hist = (false,l.anglelimits)
            new_max = 1.57
            new_min = -1.57

            if (l.anglelimits.max/1.57 > 3) | (-1*l.anglelimits.min/1.57 > 3)

                @warn "The angle limits provided is larger than 3π/2 radians.\n PowerSystems inferred the data provided in degrees and will transform it to radians" maxlog=PS_MAX_LOG

                (l.anglelimits.max/1.57 >= 0.99) ? new_max = min(l.anglelimits.max*(π/180),new_max) : new_max = min(l.anglelimits.max,new_max)

                (-1*l.anglelimits.min/1.57 > 0.99) ? new_min = max(l.anglelimits.min*(π/180),new_min) : new_min = max(l.anglelimits.min,new_min)

                hist = (true,(min = new_min, max = new_max))

            else

                (l.anglelimits.max >= 1.57 && l.anglelimits.min <= -1.57) ? hist = (true,(min = -1.57,max = 1.57)) : true
                (l.anglelimits.max >= 1.57 && l.anglelimits.min >= -1.57) ? hist =(true, (min = l.anglelimits.min,max = 1.57)) : true
                (l.anglelimits.max <= 1.57 && l.anglelimits.min <= -1.57) ? hist = (true,(min = -1.57,max = l.anglelimits.max)) : true
                (l.anglelimits.max == 0.0 && l.anglelimits.min == 0.0) ? hist = (true,(min = -1.57,max = 1.57)) : true

            end

            if hist[1]

                branches[ix] = Line(deepcopy(l.name),deepcopy(l.available),
                                    deepcopy(l.connectionpoints),deepcopy(l.r),
                                    deepcopy(l.x),deepcopy(l.b),deepcopy(l.rate),
                                    hist[2])
            end
        end
    end
end

function linerate_calculation(l::Line)
    theta_max = max(abs(l.anglelimits.min), abs(l.anglelimits.max))
    g =  l.r / (l.r^2 + l.x^2)
    b = -l.x / (l.r^2 + l.x^2)
    y_mag = sqrt(g^2 + b^2)
    fr_vmax = l.connectionpoints.from.voltagelimits.max
    to_vmax =  l.connectionpoints.to.voltagelimits.max

    if isa(fr_vmax,Nothing) || isa(to_vmax,Nothing)
        fr_vmax = 1.0
        to_vmax = 0.9
        diff_angle = abs(l.connectionpoints.from.angle -l.connectionpoints.to.angle)
        new_rate = y_mag*fr_vmax*to_vmax*cos(theta_max)
    else

        m_vmax = max(fr_vmax, to_vmax)
        c_max = sqrt(fr_vmax^2 + to_vmax^2 - 2*fr_vmax*to_vmax*cos(theta_max))
        new_rate = y_mag*m_vmax*c_max

    end

    return new_rate

end

function calculatethermallimits!(branches::Array{<:Branch,1},basemva::Float64)
    for (ix,l) in enumerate(branches)

        if isa(l,Line)

            flag = false

            #This is the same check as implemented in PowerModels
            if l.rate <= 0.0
                (flag, rate) = (true,linerate_calculation(l))
            elseif l.rate > linerate_calculation(l)
                    (flag, rate) = (true,linerate_calculation(l))
            else
                rating= l.rate
            end

            if (l.rate/basemva) > 20
                @warn "Data for line rating is 20 times larger than the base MVA for the system\n. Power Systems inferred the Data Provided is in MVA and will transform it using a base of $("basemva")" maxlog=PS_MAX_LOG
                (flag, rate) = (true,l.rate/basemva)
            end

            if flag
                branches[ix] = Line(deepcopy(l.name),deepcopy(l.available),
                                    deepcopy(l.connectionpoints),deepcopy(l.r),
                                    deepcopy(l.x),deepcopy(l.b),
                                    rate,deepcopy(l.anglelimits))
            end
        end
    end
end


text:
"""
    State(scans::Vector{PageScan})

State of the scans.
"""
struct State
    scans::Dict{String,PageScan}

    function State(scans::Vector{PageScan})
        urls = [scan.page.url for scan in scans]
        if !(length(scans) == length(unique(urls)))
            msg = """
                Not all URLs were unique.
                Cannot store multiple scans for one page.
                """
            error(msg)
        end
        dic = Dict(zip(urls, scans))
        return new(dic)
    end
end

"""
    Repo

Abstract type for repositories.
These Git repositories are used to store the state.
Of course, a database would be better but repositories are available when running in a CI job and databases are not.
"""
abstract type Repo end

"""
    MockRepo(state=State(PageScan[])) <: Repo

Mock repository used for testing.
"""
struct MockRepo <: Repo
    state::State

    function MockRepo(state=State(PageScan[]))
        return new(state)
    end
end

struct MockFileRepo <: Repo
    dir::String

    MockFileRepo(dir=mktempdir()) = new(dir)
end

"""
    GitHubRepo(;
        dir=mktempdir(),
        user=ENV["GITHUB_ACTOR"],
        token=ENV["GITHUB_TOKEN"],
        repo=ENV["GITHUB_REPOSITORY"],
        branch="skan"
    )

Assuming that we run on a GitHub Runner, then no extra information is needed.
The token and repo can be obtained via `GITHUB_TOKEN` and `GITHUB_REPOSITORY`.
For public repositories, set `token` to the empty string `""`.
`repo` is of the form `octocat/Hello-World`.
"""
struct GitHubRepo <: Repo
    dir::String
    user::String # octocat
    token::String
    repo::String # octocat/Hello-World
    branch::String

    function GitHubRepo(;
            dir=mktempdir(),
            user=ENV["GITHUB_ACTOR"],
            token=ENV["GITHUB_TOKEN"],
            repo=ENV["GITHUB_REPOSITORY"],
            branch="skan"
        )
        return new(
            dir,
            string(user)::String,
            string(token)::String,
            string(repo)::String,
            string(branch)::String
        )
    end
end

store!(repo::MockRepo, state::State) = MockRepo(state; branch=repo.branch)

state_path(repo::Repo) = joinpath(repo.dir, "state.toml")

"""
    toml(scans::Dict{String,PageScan})::String

Return `scans` as a valid TOML string.
Unlike Julia's default TOML printer, this method uses TOML's multiline strings and the ordering of the output is always the same (sorted alphabetically).
"""
function toml(scans::Dict{String,PageScan})::String
    sorted_keys = sort(collect(keys(scans)))
    items = map(sorted_keys) do key
        content = scans[key].content
        content = replace(content, TRIPLEQUOTE => "”””")

        # Avoid invalid escape characters.
        content = replace(content, '\\' => "\\\\")

        "\"$key\" = $TRIPLEQUOTE\n$content$TRIPLEQUOTE"
    end
    return join(items, '\n')
end

function store!(repo::Repo, state::State)
    path = state_path(repo)
    text = toml(state.scans)
    write(path, text)
    return repo.dir
end

retrieve(repo::MockRepo)::State = repo.state

"""
    retrieve(repo::Repo) -> State

Return state from `repo`.
"""
function retrieve(repo::Repo)
    path = state_path(repo)
    if isfile(path)
        data = read(path, String)
        dic = tomlparse(data)
        scans = [PageScan(WebPage(k), dic[k]) for k in keys(dic)]
        return State(scans)
    else
        return State(PageScan[])
    end
end

function retrieve(state::State, page::Page)::Union{PageScan,Nothing}
    # State could be a dict to speed things up but this is good enough for now.
    dic = state.scans
    key = page.url
    if key in keys(dic)
        return dic[key]
    else
        return nothing
    end
end

function update!(repo::MockRepo, state::State, scan::PageScan)
    key = scan.page.url
    repo.state.scans[key] = scan
    return repo
end

function update!(repo::MockRepo, state::State, scans::Vector)
    filter!(!isnothing, scans)
    foreach(scan -> update!(repo, state, scan), scans)
    return repo
end

function update!(repo::Repo, state::State, scans::Vector)
    if isempty(scans)
        return repo
    end
    for scan in scans
        key = scan.page.url
        state.scans[key] = scan
    end
    store!(repo, state)
    commit!(repo)
    return repo
end