text:
using DAI
using Base.Test

x = Var(0,2)
y = Var(1,2)

println("label: ", label(x))
println("states: ", states(x))

vs0 = VarSet(x) - x
@test labels(vs0) == Int[]

vs1 = VarSet(x,y)
vs2 = VarSet(Var(3,3))

@test x in vs1
@test y in vs1
@test !(x in vs2)
@test !(y in vs2)

vs3 = vs1 + vs2

@test vs1 < vs3
@test !(vs1 > vs3)
@test !(vs1 >= vs3)
@test vs1 <= vs3
@test vs2 < vs3
@test vs3 >= vs2

println("length: ", length(vs3))
println("nrstates: ", nrStates(vs3))

for i=1:nrStates(vs3)
  println("#################")
  states = calcState(vs3, i)
  print("state $i -> states: $(states')")
  state = calcLinearState(vs3, states)
  println("-> state: $state")
  @test i==state
end

vs = VarSet(Var(0,2), Var(2,2))
v = Var(1,2)
for i=1:nrStates(vs)
  i1,i2 = conditionalStateBoth(v,vs,i)
  for j=1:2
    state = conditionalState(v,vs,j,i)
    state_pris = DAI.pris_conditionalState(v,vs,j,i)
    @test state==state_pris
    if j==1
      @test i1 == state
    else
      @test i2 == state
    end
  end
end

vs = VarSet(Var(0,2), Var(2,2))
v = Var(3,2)
for i=1:nrStates(vs), j=1:2
  state = conditionalState(v,vs,j,i)
  state_pris = DAI.pris_conditionalState(v,vs,j,i)
  @test state==state_pris
end

vs = VarSet(Var(1,2), Var(2,2))
v = Var(0,2)
for i=1:nrStates(vs), j=1:2
  state = conditionalState(v,vs,j,i)
  state_pris = DAI.pris_conditionalState(v,vs,j,i)
  @test state==state_pris
end

# Memory allocation tests
vs = VarSet(Var(1,2), Var(2,2))
for i=1:100000
  vs = deepcopy(vs)
end
gc()
allvars = vars(vs)
@test allvars[1].label == 1
@test allvars[2].label == 2
@test allvars[1].states == 2
@test allvars[2].states == 2

#testvals = zeros(int(nrStates(vs3)))
#testvals[1] = 3.0
#testvals[2] = 3.0
#testvals[12] = 3.0

testvals = rand(int(nrStates(vs3)))

# Memory allocation tests
vs = VarSet(Var(1,2), Var(2,2))
fac = Factor(vs3, testvals)
for i=1:10000
  fac = deepcopy(fac)
end
gc()
@test vars(fac) == vs3
@test p(fac) == testvals

fac = Factor(vs3, testvals)

println("normalize factor: ", normalize!(fac))
tot = 0.0
for i=1:nrStates(fac)
  val = fac[i]
  println("val: $val")
  tot += val
end
@test_approx_eq abs(tot) 1.0

println("\nNow testing a marginal factor:")
fac2 = marginal(fac, vs2)
tot = 0.0
for i=1:nrStates(fac2)
  val = fac2[i]
  println("val: $val")
  tot += val
end
@test_approx_eq abs(tot) 1.0
fac2[1]=22.3

fac3 = Factor(Var(5,2) + x)
fac3[1] = 10.0
fac3[2] = 5.0
fac3[3] = 3.0

println("")
println("Testing FactorGraph")
fg = FactorGraph(fac, fac2, fac3)
println("Number of edges: $(numEdges(fg))")
println("Number of factors: $(numFactors(fg))")
println("Number of vars: $(numVars(fg))")

# Memory allocation tests
fg = FactorGraph(fac, fac2, fac3)
for i=1:10000
  fg = deepcopy(fg)
end
gc()
@test fg[1] == fac
@test fg[2] == fac2
@test fg[3] == fac3

ia = InfAlg(fg)
ia2 = deepcopy(ia)
@test ia2.hdl != ia.hdl

init!(ia)
run!(ia)
println("iterations: $(iterations(ia))")
println(properties(ia))
marg = marginal(ia, vs1)
bel = belief(ia, vs1)

@test_approx_eq sum(p(marg)) 1
@test_approx_eq sum(p(bel)) 1

@test fg[1] == fac
@test fg[2] == fac2
@test fg[3] == fac3

## Test factor copying
facc = copy(fac)
@test facc == fac
x = p(facc)
x[1] += 2.0
@test facc != fac
println("Factor copying checks out")

#Test fg and infalg copying
fgc = copy(fg)
@test fgc == fg
iac = copy(ia)
#@test iac == ia # doesn't exist yet

## Test fg equality
setBackedFactor!(fgc,1,facc)
@test fgc != fg
restoreFactors!(fgc)
@test fgc == fg

## Test fg backups
fgc[1] = facc
@test fgc != fg
restoreFactors!(fgc)
@test fgc != fg

## Test unsafe vars
facc = copy(fac)
@test facc == fac
vs = vars(fac)
vsc1 = vars(fac)
vsc2 = copy(vars(fac))
temp = vars(copy(fac))
@test vsc2 == temp

@test vs == vsc1
@test vs == vsc2
@test vsc1 == vsc2
vs2 = vs - y
@test vs == vsc1
@test vs == vsc2
@test vsc1 == vsc2

@test vs2 != vsc1
@test vs2 != vsc2
@test vs2 != vs

@test facc == fac

erase!(vs, y)

@test facc != fac #uhh.. fac is no longer a proper factor here
@show fac # because vs = vars(fac) is tampering directly with internals

@test vs == vsc1
@test vs != vsc2
@test vsc1 != vsc2

@test vs2 == vsc1
@test vs2 != vsc2
@test vs2 == vs

insert!(vs, y)

@test vs == vsc1
@test vs == vsc2
@test vsc1 == vsc2

@test vs2 != vsc1
@test vs2 != vsc2
@test vs2 != vs

@test facc == fac

@test labels(fac) == [0,1,3]
@test labels(vs1) == [0,1]

############
vs = VarSet(Var(1,2), Var(2,2))

fac = Factor(vs)
fac2 = embed(fac, Var(3,2))
fac3 = embed(fac, vs+Var(3,2))
@test fac2 == fac3
#

#y = vars(fg)
#x = vars(vs1)

#println(x)
#println(vs1)
#println(fac)
#println(fg)

v1 = Var(1,2)
v2 = Var(2,2)
v3 = Var(3,2)

f1 = Factor(v1)
f2 = Factor(VarSet(v1,v2))
f3 = Factor(VarSet(v1,v3))

f1[1] = 0.0
f1[2] = 1.0

f2[1] = 0.6
f2[2] = 0.2
f2[3] = 0.8
f2[4] = 0.2

f3[1] = 0.6
f3[2] = 0.2
f3[3] = 0.8
f3[4] = 0.2

vs = VarSet(v1,v2)
fg = FactorGraph(f1, f2, f3)

println("")
println("#########################################")
println("Testing Marginalization")
println("#########################################")
println("")
ia = InfAlg(fg)
ia2 = deepcopy(ia)
@test ia2.hdl != ia.hdl

init!(ia)
run!(ia)
println("iterations: $(iterations(ia))")
println(properties(ia))
marginal(ia, vs)
@time marg = marginal(ia, vs)
belief(ia, vs)
@time bel = belief(ia, vs)

@show p(marg)
@show p(bel)

@test_approx_eq sum(p(marg)) 1
@test_approx_eq sum(p(bel)) 1
println("")
println("#########################################")
println("")


text:
using Distributed
using DataFrames
using CSV

nprocs() == 1 && addprocs(; exeflags = "--project")
@everywhere include("twoblocks.jl")
@everywhere using Main.TwoBlocksPipeline

dataset = getiris()
X = dataset[:,1:4]
Y = dataset[:,5] |> collect
X = dataset[:,2:end-1]

results=TwoBlocksPipeline.twoblockspipelinesearch(X,Y)
println(results)
println("best pipeline is: ", results[1,1],", with mean-sd ",results[1,2], " ± ",results[1,3])

ext: jl


text:
import DataStructures: OrderedDict

@testset "bernoulli" begin

    # random
    x = bernoulli(0.5)

    # logpdf_grad
    f = (x::Bool, prob::Float64) -> logpdf(bernoulli, x, prob)
    args = (false, 0.3,)
    actual = logpdf_grad(bernoulli, args...)
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    args = (true, 0.3,)
    actual = logpdf_grad(bernoulli, args...)
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
end

@testset "beta" begin

    # random
    x = beta(0.5, 0.5)
    @test 0 < x < 1

    # out of support
    @test logpdf(beta, -1, 0.5, 0.5) == -Inf

    # avoid infinities, sanely
    @test logpdf(beta, eps(typeof(0.)), 0.5, 0.5) < logpdf(beta, 0., 0.5, 0.5) < Inf
    @test logpdf(beta, eps(typeof(0.)), 0.5, 1.5) < logpdf(beta, 0., 0.5, 1.5) < Inf
    @test logpdf(beta, 1-eps(typeof(0.)), 0.5, 0.5) < logpdf(beta, 1., 0.5, 0.5) < Inf
    @test logpdf(beta, 1-eps(typeof(0.)), 1.5, 0.5) < logpdf(beta, 1., 1.5, 0.5) < Inf

    # logpdf_grad
    f = (x, alpha, beta_param) -> logpdf(beta, x, alpha, beta_param)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(beta, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "categorical" begin

    # random
    x = categorical([0.2, 0.3, 0.5])
    @test 0 < x < 4

    # out of support
    @test logpdf(categorical, -1, [0.2, 0.3, 0.5]) == -Inf

    # integer-1 probability bug
    N = 20
    @test sum([rand(Gen.Distributions.Categorical([1, 0])) for i in 1:N]) == N

    # logpdf_grad
    f = (x, probs) -> logpdf(categorical, x, probs)
    args = (2, [0.2, 0.3, 0.5])
    actual = logpdf_grad(categorical, args...)
    @test actual[1] == nothing
    @test isapprox(actual[2][1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(actual[2][2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(actual[2][3], finite_diff_vec(f, args, 2, 3, dx))
end

@testset "gamma" begin

    # random
    x = gamma(1, 1)
    @test 0 < x

    # out of support
    @test logpdf(gamma, -1, 1, 1) == -Inf

    # logpdf_grad
    f = (x, shape, scale) -> logpdf(gamma, x, shape, scale)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(gamma, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "inv_gamma" begin

    # random
    x = inv_gamma(1, 1)
    @test 0 < x

    # out of support
    @test logpdf(inv_gamma, -1, 1, 1) == -Inf

    # logpdf_grad
    f = (x, shape, scale) -> logpdf(inv_gamma, x, shape, scale)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(inv_gamma, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "normal" begin

    # random
    x = normal(0, 1)

    # does not overflow
    @test logpdf(normal, 1e13, 0, 1) == -5e25
    @test logpdf_grad(normal, 1e13, 0, 1) == (-1e13, 1e13, 1e26)

    # logpdf_grad
    f = (x, mu, std) -> logpdf(normal, x, mu, std)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(normal, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "zero-dimensional broadcasted normal" begin

    # random
    x = broadcasted_normal(fill(0), fill(1))

    # logpdf_grad
    f = (x, mu, std) -> logpdf(broadcasted_normal, x, mu, std)
    args = (fill(0.4), fill(0.2), fill(0.3))
    actual = logpdf_grad(broadcasted_normal, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx; broadcast=true))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx; broadcast=true))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx; broadcast=true))
end

@testset "array normal (trivially broadcasted: all args have same shape)" begin

    mu =  [ 0.1   0.2   0.3  ;
            0.4   0.5   0.6  ]
    std = [ 0.01  0.02  0.03 ;
            0.04  0.05  0.06 ]
    x =   [ 1.    2.    3.   ;
            4.    5.    6.   ]

    # random
    broadcasted_normal(mu, std)

    # logpdf_grad
    f = (x, mu, std) -> logpdf(broadcasted_normal, x, mu, std)
    args = (mu, std, x)
    actual = logpdf_grad(broadcasted_normal, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx; broadcast=true))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx; broadcast=true))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx; broadcast=true))
end

@testset "broadcasted normal" begin

    ## Return shape of `broadcasted_normal`
    @test size(broadcasted_normal([0. 0. 0.], 1.)) == (1, 3)
    @test size(broadcasted_normal(zeros(1, 3, 4), ones(2, 1, 4))) == (2, 3, 4)
    @test_throws DimensionMismatch broadcasted_normal([0 0 0], [1 1])

    ## Return shape of `logpdf` and `logpdf_grad`
    @test size(logpdf(broadcasted_normal,
                      ones(2, 4), ones(2, 1), ones(1, 4))) == ()
    @test all(size(g) == ()
              for g in logpdf_grad(
                  broadcasted_normal, ones(2, 4), ones(2, 1), ones(1, 4)))
    # `x` has the wrong shape
    @test_throws DimensionMismatch logpdf(broadcasted_normal,
                                          ones(1, 2), ones(1,3), ones(2,1))
    @test_throws DimensionMismatch logpdf_grad(broadcasted_normal,
                                               ones(1, 2), ones(1,3), ones(2,1))
    # `x` has a shape that is broadcast-compatible with but not equal to the
    # right shape
    @test_throws DimensionMismatch logpdf(broadcasted_normal,
                                          ones(1, 3), ones(1,3), ones(2,1))
    @test_throws DimensionMismatch logpdf_grad(broadcasted_normal,
                                               ones(1, 3), ones(1,3), ones(2,1))
    # `mu` and `std` are broadcast-incompatible
    @test_throws DimensionMismatch logpdf(broadcasted_normal,
                                          ones(2, 1), ones(1,2), ones(1,3))
    @test_throws DimensionMismatch logpdf_grad(broadcasted_normal,
                                               ones(2, 1), ones(1,2), ones(1,3))

    ## Equivalence of broadcast to supplying bigger arrays for `mu` and `std`
    compact = OrderedDict(:x => reshape([ 0.2  0.3  0.4  0.5 ;
                                          0.5  0.4  0.3  0.2 ],
                                        (2, 4)),
                          :mu => reshape([0.7  0.7  0.8  0.6],
                                         (1, 4)),
                          :std => reshape([0.2, 0.1],
                                          (2, 1)))
    expanded = OrderedDict(:x => compact[:x],
                           :mu => repeat(compact[:mu], outer=(2, 1)),
                           :std => repeat(compact[:std], outer=(1, 4)))
    @test (logpdf(broadcasted_normal, values(compact)...) ==
           logpdf(broadcasted_normal, values(expanded)...))
    @test (logpdf_grad(broadcasted_normal, values(compact)...) ==
           logpdf_grad(broadcasted_normal, values(expanded)...))
end

@testset "multivariate normal" begin

    # random
    x = mvnormal([0.0, 0.0], [1.0 0.2; 0.2 1.4])
    @test length(x) == 2

    # logpdf_grad
    f = (x, mu, cov) -> logpdf(mvnormal, x, mu, cov)
    args = ([0.1, 0.2], [0.3, 0.4], [1.0 0.2; 0.2 1.4])
    actual = logpdf_grad(mvnormal, args...)
    @test isapprox(actual[1][1], finite_diff_vec(f, args, 1, 1, dx))
    @test isapprox(actual[1][2], finite_diff_vec(f, args, 1, 2, dx))
    @test isapprox(actual[2][1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(actual[2][2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(actual[3][1, 1], finite_diff_mat_sym(f, args, 3, 1, 1, dx))
    @test isapprox(actual[3][1, 2], finite_diff_mat_sym(f, args, 3, 1, 2, dx))
    @test isapprox(actual[3][2, 1], finite_diff_mat_sym(f, args, 3, 2, 1, dx))
    @test isapprox(actual[3][2, 2], finite_diff_mat_sym(f, args, 3, 2, 2, dx))
end

@testset "uniform" begin

    # random
    x = uniform(-0.5, 0.5)
    @test -0.5 < x < 0.5

    # out of support
    @test logpdf(uniform, -1, -0.5, 0.5) == -Inf

    # logpdf_grad
    f = (x, low, high) -> logpdf(uniform, x, low, high)
    args = (0.0, -0.5, 0.5)
    actual = logpdf_grad(uniform, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "uniform_discrete" begin

    # random
    x = uniform_discrete(10, 20)
    @test 10 <= x <= 20

    # out of support
    @test logpdf(uniform_discrete, -1, 10, 20) == -Inf

    # logpdf_grad
    args = (1, 1, 5)
    actual = logpdf_grad(uniform_discrete, args...)
    @test actual == (nothing, nothing, nothing)
end

@testset "piecewise_uniform" begin

    # random
    x = piecewise_uniform([-0.5, 0.5], [1.0])
    @test -0.5 < x < 0.5

    # out of support
    @test logpdf(piecewise_uniform, -1, [-0.5, 0.5], [1.0]) == -Inf

    # logpdf_grad
    f = (x, bounds, probs) -> logpdf(piecewise_uniform, x, bounds, probs)
    args = (0.5, [-1.0, 0.0, 1.0], [0.4, 0.6])
    actual = logpdf_grad(piecewise_uniform, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2][1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(actual[2][2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(actual[2][3], finite_diff_vec(f, args, 2, 3, dx))
    @test isapprox(actual[3][1], finite_diff_vec(f, args, 3, 1, dx))
    @test isapprox(actual[3][2], finite_diff_vec(f, args, 3, 2, dx))
end

@testset "beta uniform mixture" begin

    # random
    x = beta_uniform(0.5, 0.5, 0.5)
    @test 0 < x < 1

    # out of support
    @test logpdf(beta_uniform, -1, 0.5, 0.5, 0.5) == -Inf

    # logpdf_grad
    f = (x, theta, alpha, beta) -> logpdf(beta_uniform, x, theta, alpha, beta)
    args = (0.5, 0.4, 10., 2.)
    actual = logpdf_grad(beta_uniform, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
    @test isapprox(actual[4], finite_diff(f, args, 4, dx))
end

@testset "geometric" begin

    # random
    @test geometric(0.5) >= 0

    # out of support
    @test logpdf(geometric, -1, 0.5) == -Inf

    # logpdf_grad
    f = (x, p) -> logpdf(geometric, x, p)
    args = (4, 0.3)
    actual = logpdf_grad(geometric, args...)
    @test actual[1] == nothing
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
end

@testset "binom" begin

    # random
    @test binom(5, 0.3) >= 0

    # out of support
    @test logpdf(binom, -1, 5, 0.3) == -Inf

    # logpdf_grad
    f = (x, n, p) -> logpdf(binom, x, n, p)
    args = (2, 5, 0.3)
    actual = logpdf_grad(binom, args...)
    @test actual[1] == nothing
    @test actual[2] == nothing
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "neg_binom" begin

    # random
    @test neg_binom(5, 0.3) >= 0

    # out of support
    @test logpdf(neg_binom, -1, 5, 0.3) == -Inf

    # logpdf_grad
    f = (x, r, p) -> logpdf(neg_binom, x, r, p)
    args = (2, 5, 0.3)
    actual = logpdf_grad(neg_binom, args...)
    @test actual[1] == nothing
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "exponential" begin

    # random
    @test exponential(0.5) > 0

    # out of support
    @test logpdf(exponential, -1, 0.5) == -Inf

    # logpdf_grad
    f = (x, rate) -> logpdf(exponential, x, rate)
    args = (1.2, 0.5)
    actual = logpdf_grad(exponential, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
end

@testset "poisson" begin

    # random
    @test poisson(1.0) >= 0

    # out of support
    @test logpdf(poisson, -1, 1.0) == -Inf

    # logpdf_grad
    f = (x, lambda) -> logpdf(poisson, x, lambda)
    args = (4, 2.0)
    actual = logpdf_grad(poisson, args...)
    @test actual[1] == nothing
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
end

@testset "laplace" begin

    # random
    x = laplace(1, 1)

    # logpdf_grad
    f = (x, loc, scale) -> logpdf(laplace, x, loc, scale)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(laplace, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

@testset "cauchy" begin

    # random
    x = cauchy(0, 5)

    # logpdf_grad
    f = (x, x0, gamma) -> logpdf(cauchy, x, x0, gamma)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(cauchy, args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end


text:
using PH2514

# Question 1
# What is the Energy (in electron Volts - eV) associated with radiation at 13 micrometers? (Do not include units with your answer)
eV(Wavelength(13e-6))

# Question 2
# Transitions to the n = 5 state of the hydrogen atom are known as the
# Pfund

# Question 3
# What is the ground state energy for a Lithium ion, in eV? (Note that Z=3, the ground state energy for hydrogen is -13.58 eV). Enter the absolute value of your answer, i.e. don't enter the minus sign or list any units.
E_Zn(3,1)

# Question 4
# Calculate the energy for the hydrogen transition from N=5 to N=3, in electron Volts (eV). Do not enter units with your answer. Positive (absolute) value.
E_Zn(1,5)-E_Zn(1,3)

# Question 5
# Calculate the wavelength for the peak of blackbody radiation at T = 5000 K. (Do not include units with your answer).
Wavelength(WiensDisplacementLaw(5000))

# Question 6
# Use the Stefan Boltzman Law to calculate the radiated power, in Watts, of a gray body at 97 K (boiling point of liquid nitrogen), emissitivity = 0.9 and surface area = 2 m2. (Do not enter units with your answer). Note the difference between a gray body and a blackbody is the emissivity.
SteffannBoltzman(.9,97)* 2

# Question 7
# Kinetic Theory: For a helium gas, atoms with a temperature of 1100 K have a random (rms) velocity of:
Vrms(6.6464731e-27,1100)

# Question 8
# Plasma Physics: Gyroradius: For an electron with velocity of 527 km/s (5.27x105 m/s), Magnetic field of 150 x 10-9 Tesla (150 nanoTesla); calculate the gyroradius in meters. Do not enter the units with your answer.
GyroRadius(C[:me],5.27e5,C[:qe],150e-9)

# Question 9
# Plasma Physics: Gyrofrequency: For an electron with velocity of 527 km/s, Magnetic field of 150 x 10-9 Tesla (150 nT); calculate the gyrofrequency in Hz (not radians per second; f, not omega). Do not enter the units with your answer.
GyroFrequency(C[:qe],150e-9,C[:me]) / (2*π)

# Question 10
# Plasma Physics: Plasma Frequency:  Electrons with a density of 2x106 e/m3 are present.  Calculate the plasma frequency in KHz. (not radians per second). (Enter answer without units)
PlasmaFrequency(2e6)/2/π/1e3

# Question 11
# The direction of propagation of an electromagnetic wave is perpendicular to the magnetic field and parallel to the electric field.
False

# Question 12
# A magnetic field that is 'frozen' into a plasma will move with the plasma as it moves.
True

# Question 13
# A magnetic field can be considered as also exerting magnetic pressure.
True

# Question 14
# What is the drift velocity for a proton in a magnetic field of 35 micro Tesla with a perpendicular electric field of 1.5 mV/meter?
vD = 1.5e-3 / 35e-6

# Question 15
# If you double the magnetic field, what happens to the gyro radius?
# It halves


text:
using DeformableRegistration: ImageProcessing, Transformation, Visualization
 using PyPlot
 using Base.Test
#  using Logging

# setup logging

@testset "Visualization" begin
    testimage = dirname(Base.source_path()) * "/testdata/testimage.png"
    img = loadImage(testimage)
    figure();
    showImage(img,name="/testdata/testimage.png")
    @test typeof(gcf()) == Figure
end

close("all")
#
# # test showGrid
# gridSize = (256,256)
#  voxelsize = [0.1, 0.1]
#  shift = [5,5.0]
#  centeredGrid = getCellCenteredGrid(voxelsize, shift, gridSize)
#  affineParameters = [0.5,0.1,0,0,0.5,0.5]
#  deformationField = zeros(prod(gridSize)*2)
#  deformationField[prod(gridSize)+1:end] = 0.05*sin(5*centeredGrid[1:prod(gridSize)])
#  figure();
#  showGrid(spatialDomain,gridSize,centeredGrid,
#          affineParameters=affineParameters,
#          deformationField=deformationField)
# # check if figure was created
# @assert typeof(gcf()) == Figure
# Logging.info("Visualization: showGrid: Figure was created ✔")
# # check y and x limits
# @test ylim() == (2.1,-1.1)
# @test xlim() == (-1.0,2.21)
# Logging.info("Visualization: showGrid: x,y limits are alright ✔")


text:
export NoParams

struct NoParams{N,BT<:AbstractBlock{N}} <: TagBlock{BT,N}
    content::BT
end

YaoBlocks.parameters!(out, x::NoParams) = out
YaoBlocks.nparameters(x::NoParams) = 0
YaoBlocks.niparams(x::NoParams) = 0
YaoBlocks.mat(::Type{T}, x::NoParams) where {T} = mat(T, content(x))
YaoBlocks.PropertyTrait(x::NoParams) = PreserveAll()
YaoBlocks.dispatch!(f::Union{Function,Nothing}, x::NoParams, it::YaoBlocks.Dispatcher) = return x
Base.adjoint(x::NoParams) = NoParams(adjoint(content(x)))


text:
function eval_evenodd_trig(a, a², pos, neg, ::Val{even}) where even
	s = sign(a²)
	norm = sqrt(abs(a²))
	if even
		s > 0 ? pos(norm) : s < 0 ? neg(norm) : one(a)
	else
		s > 0 ? pos(norm)*a/norm : s < 0 ? neg(norm)*a/norm : a
	end
end

SERIES_ORDER = [150]

function eval_series(a, fn, order=first(SERIES_ORDER))
	T = eltype(a)
	series = fn(Taylor1(T, order))
	series(a)
end

isapproxzero(x) = abs(x) <= eps(float(one(x)))/2

crudenorm(a::Scalar) = a
crudenorm(a::Blade) = a.coeff
crudenorm(a::CompositeMultivector) = sum(abs.(a.components))

function exp_with_scalar_square(a, a²)
	isapproxzero(a²) && return one(a) + a

	norm = sqrt(abs(a²))
	if a² > 0
		cosh(norm) + sinh(norm)/norm*a
	else
		cos(norm) + sin(norm)/norm*a
	end
end

function exp_series(a, p=1)
	# TODO: figure out when exp(a) -> exp(a/p)^p is appropriate
	#       and how to find optimal p
	a /= p

	result = one(a)
	term = one(a)
	for i = 1:100
		term *= a/i
		if isapproxzero(crudenorm(term))
			break
		end
		result += term
	end

	result^p
end

function Base.exp(a::AbstractMultivector)
	# use fact that exp(scalar + a) = exp(scalar)exp(a)
	s = scalar(a)
	if !iszero(s)
		prefactor = exp(s)
		a -= s
	else
		prefactor = one(s)
	end

	a² = a*a
	if isscalar(a²)
		prefactor*exp_with_scalar_square(a, scalar(a²))
	else
		prefactor*exp_series(a - s)
	end
end

Base.exp(a::Blade) = exp_with_scalar_square(a, scalar(a*a))

for (fn,     pos,    neg,    even) ∈ [
	(:cosh,  :cosh,  :cos,   true),
	(:cos,   :cos,   :cosh,  true),
	(:sinh,  :sinh,  :sin,   false),
	(:sin,   :sin,   :sinh,  false),
	(:tanh,  :tanh,  :tan,   false),
	(:tan,   :tan,   :tanh,  false),
	(:asinh, :asinh, :asin,  false),
	(:asin,  :asin,  :asinh, false),
	(:atanh, :atanh, :atan,  false),
	(:atan,  :atan,  :atanh, false),
]
	@eval function Base.$fn(a::AbstractMultivector)
		a² = a*a
		if isscalar(a²)
			eval_evenodd_trig(a, scalar(a²), $pos, $neg, Val($even))
		else
			eval_series(a, $fn)
		end
	end
end

for (fn,    invfn) ∈ [
	(:sec,  :cos),
	(:sech, :cosh),
	(:csc,  :sin),
	(:csch, :sinh),
	(:cot,  :tan),
	(:coth, :tanh),
]
	@eval Base.$fn(a::AbstractMultivector) = inv(Base.$invfn(a))
end

# log(x) = 2atanh((x - 1)/(x + 1))
Base.log(a::AbstractMultivector) = 2atanh((a - 1)/(a + 1))


text:
# This example shows how to use your own field implementation for using them to compute on
# elliptic curves. For making things interesting we shall wrap GaloisFields implementeed binary extension field.
using Test
using CryptoGroups
import CryptoGroups: BinaryField, specialize, PB, bitlength
import GaloisFields: ExtensionField, GaloisField

struct GF₂{F <: ExtensionField} <: BinaryField
    x::F
    GF₂{F}(coeffs::NTuple{N, GaloisField(2)}) where {N, F <: ExtensionField} = new(F(coeffs))

    function GF₂{F}(xrev::BitVector) where F <: ExtensionField
        x = reverse(xrev)
        𝔽₂ = GaloisField(2)
        coeffs = Tuple(𝔽₂(i) for i in x) # May need a reverse
        return GF₂{F}(coeffs)
    end

    GF₂{F}(x::F) where F <: ExtensionField = new{F}(x)
    GF₂(x::F) where F <: ExtensionField = new{F}(x)

    function GF₂(f::BitVector)
        F, x = GaloisField(2, :x=>f)
        return GF₂{F}
    end

    GF₂(poly::Vector{Int}) = GF₂(BitVector(i in poly for i in 0:maximum(poly)))
end

bitlength(::Type{F}) where F <: ExtensionField = length(F.parameters[4]) - 1
bitlength(::Type{GF₂{F}}) where F <: ExtensionField = bitlength(F)

specialize(::Type{GF₂}, basis::PB) = GF₂(basis.f)

Base.convert(::Type{F}, a::BitVector) where F <: GF₂ = F(a)
Base.convert(::Type{BitVector}, a::GF₂) = reverse(BitVector(i.n for i in a.x.coeffs))

#tobits(a::GF₂) = reverse(BitVector(i.n for i in a.x.coeffs)) # For the sake of symmetry

Base.:+(a::F, b::F) where F <: GF₂ = F(a.x + b.x)
Base.:*(a::F, b::F) where F <: GF₂ = F(a.x * b.x)

Base.inv(a::F) where F <: GF₂ = F(inv(a.x))

Base.zero(::Type{GF₂{F}}) where F <: ExtensionField = GF₂(zero(F))
Base.one(::Type{GF₂{F}}) where F <: ExtensionField = GF₂(one(F))

########################## This one we can test easally as follows ##################

import CryptoGroups: BinaryCurve, generator, oncurve, order, @def, AffinePoint, <|, ECPoint, validate
import CryptoGroups: Specs
import .Specs: Curve_B_163_PB, Curve_K_163_PB

let
    B_163v3 = specialize(AffinePoint{BinaryCurve, GF₂}, Curve_B_163_PB)

    g = B_163v3 <| generator(Curve_B_163_PB)
    q = order(Curve_B_163_PB)

    @test oncurve(g)
    @test oncurve(g*3)
    @test g * q == zero(g)
end

@def K_163v3 ECPoint{AffinePoint{BinaryCurve, GF₂}} Curve_K_163_PB

let
    g = K_163v3 <| generator(Curve_K_163_PB)

    @test oncurve(g)
    @test oncurve(g*3)
    @test validate(g)
end


text:
using XGrad
# using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
using Test

function test_compare(f; inputs...)
    vals = ([val for (name, val) in inputs]...,)
    ctx = Dict()
    df = xdiff(f; ctx=ctx, inputs...)
    # compare_test runs in an older world age than the generated function
    # need to invoke the latest one in the test
    dvals = Base.invokelatest(df, vals...)
    dvals_a = [dvals...]

    # TODO: ReverseDiff doesn't support Julia 0.7 yet, use finite difference instead
    # f_tape = GradientTape(f, vals)
    # compiled_f_tape = compile(f_tape)
    # cfg = GradientConfig(vals)
    # results = map(similar, vals)
    # gradient!(results, compiled_f_tape, vals)
    # results_a = [results...]
    # @test isapprox(results_a, dvals_a[2:end]; atol=0.1)
end

logistic(x) = 1 ./ (1 + exp.(-x))
@diffrule logistic(x::Number) x (logistic(x) .* (1 .- logistic(x))) .* ds

include("tricky.jl")
include("aggr.jl")
include("linreg.jl")
include("ann.jl")
include("autoencoder.jl")
# include("vae.jl")  -- disabled until Distributions are updated to Julia 0.7
include("others.jl")
include("struct.jl")
include("xgrad.jl")


text:
#=

    This is a port of the R2 model TrueSkillSimple.cs

    With the datafiles
    - true_skill_simple_r2_players.csv
    - true_skill_simple_r2_games.csv

    Output from the R2 model (skills):
    ```
    [0] Mean: 107.265          skills[0]
    [0] Variance: 72.9736

    [1] Mean: 100.541          skills[1]
    [1] Variance: 86.9287

    [2] Mean: 95.7907          skills[2]
    [2] Variance: 84.9168
    ```

    Output from this Turing.jl model:

        parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
            Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

         skills[1]   105.8074    8.9089     0.0891    0.2314   1273.9109    1.0009      155.4119
         skills[2]    99.9198    8.8072     0.0881    0.2444   1261.0051    0.9999      153.8374
         skills[3]    94.0010    9.2404     0.0924    0.2506   1305.8709    0.9999      159.3108
  performance[1,1]   112.9550   14.6471     0.1465    0.4051   1246.4095    0.9999      152.0568
  performance[1,2]    92.8609   14.5794     0.1458    0.4111   1205.0434    1.0000      147.0103
  performance[2,1]   107.0177   14.7617     0.1476    0.3511   1777.9260    1.0001      216.8996
  performance[2,2]    86.8285   14.6496     0.1465    0.3437   1802.4896    1.0000      219.8963
  performance[3,1]   111.3551   15.1643     0.1516    0.2474   2981.5172    1.0002      363.7327
  performance[3,2]    88.0326   15.0679     0.1507    0.3018   2847.2484    1.0005      347.3525

    Cf true_skill.jl

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV, DataFrames
include("jl_utils.jl")

@model function true_skill_simple_r2(players,games)
    num_players = length(players)
    num_games, m = size(games)

    skills ~ filldist(Normal(100,10),num_players)

    performance = Matrix{Real}(undef, num_players,2)
    for i in 1:num_games
        performance[i,1] ~ Normal(skills[games[i,1]+1],15)
        performance[i,2] ~ Normal(skills[games[i,2]+1],15)

        if games[i,3] == 1
            true ~ Dirac(performance[i,1] > performance[i,2])
        elseif games[i,3] == 0
            true ~ Dirac(performance[i,1] == performance[i,2])
        else
            true ~ Dirac(performance[i,1] < performance[i,2])
        end
    end

end

# 0 1 2
players = parse.(Int64,split(readline("true_skill_simple_r2_players.csv"),","))

# 0  1  1
# 1  2  1
# 0  2  1
games = Matrix(CSV.read("true_skill_simple_r2_games.csv",header=false,DataFrame))
println("games:")
display(games)

model = true_skill_simple_r2(players,games)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

# chns = sample(model, NUTS(), 1_000)
# chns = sample(model, HMC(0.1,6), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function bench()
    print("( ")
    for p in PrimesPaged() p > 100 && break; print(p, " ") end
    println(")")
    countPrimesTo(Prime(100)) # warm up JIT
#=
    println(@time let count = 0
                      for p in PrimesPaged()
                          p > 1000000000 && break
                          count += 1
                      end; count end) # much slower counting over iteration
=#
    println(@time countPrimesTo(Prime(1000000000)))
end
bench()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
############################################################################################
### Types
############################################################################################

# Type to store Expectation Maximization result --------------------------------------------
mutable struct EmMVNModel{A, b, B}
    Σ::A
    μ::b
    fitted::B
end

"""
For observed data with missing values.

# Constructor

    SemObservedMissing(;
        specification,
        data,
        obs_colnames = nothing,
        kwargs...)

# Arguments
- `specification`: either a `RAMMatrices` or `ParameterTable` object (1)
- `data`: observed data
- `obs_colnames::Vector{Symbol}`: column names of the data (if the object passed as data does not have column names, i.e. is not a data frame)

# Extended help
## Interfaces
- `n_obs(::SemObservedMissing)` -> number of observed data points
- `n_man(::SemObservedMissing)` -> number of manifest variables

- `get_data(::SemObservedMissing)` -> observed data
- `data_rowwise(::SemObservedMissing)` -> observed data as vector per observation, with missing values deleted

- `patterns(::SemObservedMissing)` -> indices of non-missing variables per missing patterns
- `patterns_not(::SemObservedMissing)` -> indices of missing variables per missing pattern
- `rows(::SemObservedMissing)` -> row indices of observed data points that belong to each pattern
- `pattern_n_obs(::SemObservedMissing)` -> number of data points per pattern
- `pattern_nvar_obs(::SemObservedMissing)` -> number of non-missing observed variables per pattern
- `obs_mean(::SemObservedMissing)` -> observed mean per pattern
- `obs_cov(::SemObservedMissing)` -> observed covariance per pattern
- `em_model(::SemObservedMissing)` -> `EmMVNModel` that contains the covariance matrix and mean vector found via optimization maximization

## Implementation
Subtype of `SemObserved`

## Remarks
(1) the `specification` argument can also be `nothing`, but this turns of checking whether
the observed data/covariance columns are in the correct order! As a result, you should only
use this if you are shure your observed data is in the right format.

## Additional keyword arguments:
- `spec_colnames::Vector{Symbol} = nothing`: overwrites column names of the specification object
"""
mutable struct SemObservedMissing{
        A <: AbstractArray,
        D <: AbstractFloat,
        O <: AbstractFloat,
        P <: Vector,
        P2 <: Vector,
        R <: Vector,
        PD <: AbstractArray,
        PO <: AbstractArray,
        PVO <: AbstractArray,
        A2 <: AbstractArray,
        A3 <: AbstractArray,
        S <: EmMVNModel
        } <: SemObserved
    data::A
    n_man::D
    n_obs::O
    patterns::P # missing patterns
    patterns_not::P2
    rows::R # coresponding rows in data_rowwise
    data_rowwise::PD # list of data
    pattern_n_obs::PO # observed rows per pattern
    pattern_nvar_obs::PVO # number of non-missing variables per pattern
    obs_mean::A2
    obs_cov::A3
    em_model::S
end

############################################################################################
### Constructors
############################################################################################

function SemObservedMissing(;
        specification,
        data,

        obs_colnames = nothing,
        spec_colnames = nothing,

        kwargs...)

    if isnothing(spec_colnames) spec_colnames = get_colnames(specification) end

    if !isnothing(spec_colnames)
        if isnothing(obs_colnames)
            try
                data = data[:, spec_colnames]
            catch
                throw(ArgumentError(
                    "Your `data` can not be indexed by symbols. "*
                    "Maybe you forgot to provide column names via the `obs_colnames = ...` argument.")
                    )
            end
        else
            if data isa DataFrame
                throw(ArgumentError(
                    "You passed your data as a `DataFrame`, but also specified `obs_colnames`. "*
                    "Please make shure the column names of your data frame indicate the correct variables "*
                    "or pass your data in a different format.")
                    )
            end

            if !(eltype(obs_colnames) <: Symbol)
                throw(ArgumentError("please specify `obs_colnames` as a vector of Symbols"))
            end

            data = reorder_data(data, spec_colnames, obs_colnames)
        end
    end

    if data isa DataFrame
        data = Matrix(data)
    end

    # remove persons with only missings
    keep = Vector{Int64}()
    for i = 1:size(data, 1)
        if any(.!ismissing.(data[i, :]))
            push!(keep, i)
        end
    end
    data = data[keep, :]

    n_obs, n_man = size(data)

    # compute and store the different missing patterns with their rowindices
    missings = ismissing.(data)
    patterns = [missings[i, :] for i = 1:size(missings, 1)]

    patterns_cart = findall.(!, patterns)
    data_rowwise = [data[i, patterns_cart[i]] for i = 1:n_obs]
    data_rowwise = convert.(Array{Float64}, data_rowwise)

    remember = Vector{BitArray{1}}()
    rows = [Vector{Int64}(undef, 0) for i = 1:size(patterns, 1)]
    for i = 1:size(patterns, 1)
        unknown = true
        for j = 1:size(remember, 1)
            if patterns[i] == remember[j]
                push!(rows[j], i)
                unknown = false
            end
        end
        if unknown
            push!(remember, patterns[i])
            push!(rows[size(remember, 1)], i)
        end
    end
    rows = rows[1:length(remember)]
    n_patterns = size(rows, 1)

    # sort by number of missings
    sort_n_miss = sortperm(sum.(remember))
    remember = remember[sort_n_miss]
    remember_cart = findall.(!, remember)
    remember_cart_not = findall.(remember)
    rows = rows[sort_n_miss]

    pattern_n_obs = size.(rows, 1)
    pattern_nvar_obs = length.(remember_cart)

    cov_mean = [cov_and_mean(data_rowwise[rows]) for rows in rows]
    obs_cov = [cov_mean[1] for cov_mean in cov_mean]
    obs_mean = [cov_mean[2] for cov_mean in cov_mean]

    em_model = EmMVNModel(zeros(n_man, n_man), zeros(n_man), false)

    return SemObservedMissing(data, Float64(n_man), Float64(n_obs), remember_cart,
    remember_cart_not,
    rows, data_rowwise, Float64.(pattern_n_obs), Float64.(pattern_nvar_obs),
    obs_mean, obs_cov, em_model)
end

############################################################################################
### Recommended methods
############################################################################################

n_obs(observed::SemObservedMissing) = observed.n_obs
n_man(observed::SemObservedMissing) = observed.n_man

############################################################################################
### Additional methods
############################################################################################

get_data(observed::SemObservedMissing) = observed.data
patterns(observed::SemObservedMissing) = observed.patterns
patterns_not(observed::SemObservedMissing) = observed.patterns_not
rows(observed::SemObservedMissing) = observed.rows
data_rowwise(observed::SemObservedMissing) = observed.data_rowwise
pattern_n_obs(observed::SemObservedMissing) = observed.pattern_n_obs
pattern_nvar_obs(observed::SemObservedMissing) = observed.pattern_nvar_obs
obs_mean(observed::SemObservedMissing) = observed.obs_mean
obs_cov(observed::SemObservedMissing) = observed.obs_cov
em_model(observed::SemObservedMissing) = observed.em_model


text:
function read_to_check(x::String)
    y = read_file(x)
    d = Dict(k => sort(read_file(joinpath(SLIDE_DIR, y["Path"]..., ensurearray(v["name"])...)))
        for (k,v) in y["Input"])
    return Dict(Symbol(k) => edit_with(v, Rename.(propertynames(v), Symbol.(y["Input"][k]["col"])))
        for (k,v) in d)
end

# READ SETS
y = read_file(joinpath("data", "readfiles", "list_sets.yml"));
set = Dict((length(ensurearray(k)) == 1 ? Symbol(k) : Tuple(Symbol.(k))) =>
    sort(read_file(joinpath(y["Path"]..., ensurearray(v)...)))[:,1] for (k,v) in y["Input"])

# ******************************************************************************************
#   READ BLUENOTE OUTPUT RESULTS TO CHECK.
# ******************************************************************************************
println("  Reading i/o data...")
bio = read_to_check(joinpath(SLIDE_DIR,"calibrate","data_temp","check_io.yml"));

println("  Reading share data...")
bshr = read_to_check(joinpath(SLIDE_DIR,"calibrate","data_temp","check_share.yml"));
set[:notrd] = bshr[:notrd][:,:s];

println("  Reading calibration data...")
bcal = read_to_check(joinpath(SLIDE_DIR,"calibrate","data_temp","check_cal.yml"));
if :va in propertynames(bcal[:va0])
    valcols = Symbol.(unique(bcal[:va0][:,:va]))
    bcal[:va0] = unstack(bcal[:va0], :va, :value)
    bcal[:va0] = edit_with(bcal[:va0], Replace.(valcols,missing,0.0))
end

println("  Read disaggregation intermediary data...")
bdisagg_int = read_to_check(joinpath(SLIDE_DIR,"calibrate","data_temp","check_disagg_int.yml"));

println("  Read disaggregation output (this will take a minute)...")
bdisagg_out = read_to_check(joinpath(SLIDE_DIR,"calibrate","data_temp","check_disagg.yml"));
bdisagg = merge(bdisagg_int, bdisagg_out)

# disagg_check = Dict(k => missing for k in keys(bdisagg))

# k = :netval
# d = bshr[k]
# for (k,d) in bshr
#     if :gdpcat in propertynames(d)
#         global bshr[k] = unstack(d, :gdpcat, :value)
#         global bshr[k] = edit_with(bshr[k],
#             Replace.(setdiff(propertynames(bshr[k]),[:yr,:r,:g]), missing, 0.0))
#     end
# end

# # GSP -- LABOR
# bshr[:wg][!,:wg] .= true
# bshr[:hw][!,:hw] .= true

# bshr[:labor] = leftjoin(bshr[:labor], bshr[:hw], on = Pair.([:r,:g], [:r,:g]))
# bshr[:labor] = leftjoin(bshr[:labor], bshr[:wg], on = Pair.([:yr,:r,:g], [:yr,:r,:g]))
# bshr[:labor] = edit_with(bshr[:labor], Replace.([:hw,:wg], missing, false))

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Images, ImageFiltering, FileIO
Base.isless(a::RGB{T}, b::RGB{T}) where T =
    red(a) < red(b) || green(a) < green(b) || blue(a) < blue(b)
Base.middle(x::RGB) = x

img = load("data/lenna100.jpg")
mapwindow(median!, img, (3, 3))


text:
struct CuCodeGen
    eltyp::Type
end

const FT = Float32

const CUDA_NATIVE_RULES = [
    :(log.(x)) => :(CUDAnative.log.(x)),
    :(exp.(x)) => :(CUDAnative.exp.(x)),
    :(sqrt.(x)) => :(CUDAnative.sqrt.(x)),
    :(x .^ n) => :(CUDAnative.pow.(x, __FLOAT_TYPE__(n))),
    :(ones(n)) => :(CuArray(ones(__FLOAT_TYPE__, n))),
    # :(transpose(x)) => :(permutedims(x, (2,1))),  -- seems to cauase segfault in complex cases
]

function cuda_buffer_expr(var, eltyp, buffer_var, sz)
    if isempty(sz)
        return :($var = $eltyp(0.0))
    else
        return :($var = @get_or_create($buffer_var, $(Expr(:quote, var)),
                                       CuArray(zeros($eltyp, $sz))))
    end
end

function generate_code(codegen::CuCodeGen, g::ExGraph, nd::ExNode)
    # nd = cast_const_type(nd, codegen.eltyp)
    ex = to_buffered(g, nd)
    ex = rewrite_all(ex, CUDA_NATIVE_RULES; phs=[:x, :y, :z, :n])
    ex = subs(ex, Dict(:__FLOAT_TYPE__ => codegen.eltyp))
    return ex
end

function generate_code(codegen::CuCodeGen, g::ExGraph)
    g = eliminate_common(g)
    # g = cast_const_type(g, codegen.eltyp)
    ex = to_buffered(g)
    ex = rewrite_all(ex, CUDA_NATIVE_RULES; phs=[:x, :y, :z, :n])
    ex = subs(ex, Dict(:__FLOAT_TYPE__ => codegen.eltyp))
    init_exs = [cuda_buffer_expr(var, codegen.eltyp, :mem, sz) for (var, sz) in g.ctx[:rsizes]
                if haskey(g, var) && getcategory(g[var]) != :input]
    res = Expr(:block, init_exs..., ex.args...)
    return res
end

eval_codegen(codegen::CuCodeGen) = CuVecCodeGen(codegen.eltyp)


text:
#
# Parrot.jl/src/utils.jl
#
# By Elijah Mathews (me@elijahmathews.com)
# openpgp4fpr:ac1d3fb1e8a5eb7d14bd587b2932c725055a90d8
#

"""
    splitindices(frac::Real, N::Integer)

Generate two sets of unique indices of cumulative length `N`.
Returns a tuple:

    (indices1, indices2)

where `indices1` is an array of unique indices of length `frac × N` and
`indices2` is an array of unique indices of length `(1 - frac) × N`. Can
be used for generating indices to separate training and test datasets, for
example.
"""
function splitindices(frac, N::Integer)

    @assert 0 <= frac <= 1

    # Generate unique random indices.
    base = randperm(N)

    indices1 = base[1:ceil(Int, frac * N)]
    indices2 = base[(ceil(Int, frac * N) + 1):end]

    return indices1, indices2

end

"""
    fractionalerrorquantiles(fractionalerror; quants=[0.0005, 0.005, 0.025, 0.5, 0.975, 0.995, 0.9995])

Compute the 0.05%, 0.5%, 2.5%, 50% (median), 97.5%, 99.5%, and 99.95% quantiles
for a set of spectral energy distributions (SEDs) predicted by Parrot given an
array of the fractional errors `fractionalerror`.

The selected quantiles can be changed from their defaults by passing an array to
the `quants` keyword argument.
"""
function fractionalerrorquantiles(fractionalerror; quants=[0.0005, 0.005, 0.025, 0.5, 0.975, 0.995, 0.9995])

    returnval = []

    for i in 1:length(fractionalerror)

        result = zeros(length(quants), size(fractionalerror[i], 1))

        for j in 1:size(fractionalerror[i], 1)

            for k in 1:length(quants)

                result[k,j] = quantile(fractionalerror[i][j,:], quants[k])

            end

        end

        push!(returnval, result)

    end

    return returnval

end

"""
    convert(MultivariateStats.PCA{Float32}, pca::MultivariateStats.PCA{Float64})

Convert a `Float64` PCA model `pca` to a `Float32` PCA model.
"""
function Base.convert(::Type{MultivariateStats.PCA{Float32}}, pca::MultivariateStats.PCA{Float64})
    MultivariateStats.PCA(
        convert(Array{Float32,1}, pca.mean    ),
        convert(Array{Float32,2}, pca.proj    ),
        convert(Array{Float32,1}, pca.prinvars),
        convert(Float32,          pca.tprinvar),
        convert(Float32,          pca.tvar    ),
    )
end


text:
using DrWatson, PyCall
@quickactivate "GCN_HM_GRN-Integration"

pushfirst!(PyVector(pyimport("sys")."path"), "");
run_sim = pyimport("port_batched").run_sim

allparams = Dict(
    :layer => ["arma", "sage", "tag"],
    :lr => [0.0005, 0.0001],
    :cl => ["E116"]
)

dicts = dict_list(allparams)

function makesim(d::Dict)
    @unpack layer, lr, cl = d
    best_auc, losses_test, losses_train, auc_l = run_sim(cl, lr, layer)
    fulld = copy(d)
    fulld[:auc] = best_auc
    fulld[:losses_test] = losses_test
    fulld[:losses_train] = losses_train
    fulld[:auc_list] = auc_l
    return fulld
end

for (i, d) in enumerate(dicts)
    f = makesim(d)
    print("Finished $(i)")
    @tagsave(datadir("port_batched", savename(d, "bson"; digits=5)), f; safe = true)
end


text:
using ConcreteStructs
using ComponentArrays
using DifferentialEquations
using DiffEqFlux
using DiffEqSensitivity
using DiffEqUncertainty
using Flux
using Optim
using Setfield
using UnPack

include("general_utils.jl")
include("neural_utils.jl")
include("pendulum_model.jl")
include("controllers.jl")


text:
# multigrid.jl : define multigrid cycle operations

# Multigrid cycles
abstract type MultigridCycle end

struct V <: MultigridCycle
    ν₁::Int
    ν₂::Int
end
V() = V(2, 1)

struct W <: MultigridCycle
    ν₁::Int
    ν₂::Int
end
W() = W(2, 1)

struct F <: MultigridCycle
    ν₀::Int
    ν₁::Int
    ν₂::Int
end
F() = F(2, 2, 1)

# Multigrid struct
struct MultigridIterable{C<:MultigridCycle, Gs, S, V}
    grids::Gs
    max_iter::Int
    cycle_type::C
    smoother::S
    resnorm::V
end

show(io::IO, mg::MultigridIterable) = print(io, print_grid_sizes(mg.grids), "-grid Multigrid method")

print_grid_sizes(grids::Array{<:Grid}) = join(size(grids), "x")

"""
    MultigridMethod(A, sz, cycle_type)
    MultigridMethod(A, sz, cycle_type; kwargs...)
    MultigridMethod(f, sz, cycle_type)
    MultigridMethod(f, sz, cycle_type; kwargs...)

Geometric Multigrid method of type `cycle_type` for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh. A Galerkin approach is used to compose the coarse matrices, unless a function `f` is provided for direct discretization.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
f          : Function, function used for direct-discretization of the PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`
cycle_type : Multigrid cycle type, can be `V(ν₁,ν₂)` (V-cycle), `W(ν₁,ν₂)` (W-cycle) or `F(ν₀,ν₁,ν₂)` (FMG)

Options
=======
* max_iter : maximum number of iterations. Default value is `20` For FMG using the `F()`-cycle, this is set to 1.
* R_op     : restriction operator type, can be `Injection()` or `FullWeighting()` (default)
* P_op     : interpolation operator type, can be `FullWeighting()` (default), or `Cubic()`
* ngrids   : total number of grids to use, default is `min.(⌊log₂(sz)⌋)`
* smoother : smoother, can be `GaussSeidel()` of `Jacobi()`
"""
MultigridMethod(A::Union{AbstractMatrix,Function}, sz::NTuple, cycle_type::MultigridCycle; max_iter::Int=20, R_op::TransferKind=FullWeighting(), P_op::TransferKind=FullWeighting(), ngrids::Int=minimum(factor_twos.(sz)), smoother::Smoother=GaussSeidel()) = MultigridIterable(coarsen(A, sz, R_op, P_op, ngrids), max_iter, cycle_type, smoother, Vector{Float64}(undef, 0))

"""
V_cycle(A, sz)
V_cycle(A, sz; kwargs...)

Geometric Multigrid V(2,1)-cycle for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
V_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A ,sz, V(); kwargs...)

"""
W_cycle(A, sz)
W_cycle(A, sz; kwargs...)

Geometric Multigrid W(2,1)-cycle for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
W_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A, sz, W(); kwargs...)

"""
F_cycle(A, sz)
F_cycle(A, sz; kwargs...)

Geometric Full Multigrid for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
F_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A, sz, F(); max_iter=1, kwargs...)

# solver
"""
\\(M,b)

Solve the matrix problem `Mx=b` using Multigrid.

Inputs
======
M : MultigridIterable
b : right-hand side
"""
function \(mg::MultigridIterable, b::AbstractVector)
    length(b) == length(mg.grids[1].b) || throw(DimensionMismatch(string("Right-hand side b has length ", length(b), " but needs ", length(mg.grids[1].b)))) # check dimensions
    mg.grids[1].b .= b # copy rhs
    ϵ = 1/prod(mg.grids[1].sz)

    for (i, iter) in enumerate(Iterators.take(mg, mg.max_iter))
        mg.resnorm[end] < ϵ && break
    end

    mg.resnorm[end] ≥ ϵ && @warn @sprintf("maximum number of iterations reached, norm of residual is %6.3e > %6.3e", mg.resnorm[end], ϵ) # check convergence with max_iter iterations

    return mg.grids[1].x
end

# iterator commands
function iterate(iter::MultigridIterable, count=0)
    count == 0 && push!(iter.resnorm, norm_of_residu(iter.grids[1]))
    cycle!(iter)
    push!(iter.resnorm, norm_of_residu(iter.grids[1])) # log convergence history
    nothing, count+1
end

# multigrid cycles
cycle!(mg::MultigridIterable{C} where {C<:V}) = μ_cycle!(mg.grids, 1, mg.cycle_type.ν₁, mg.cycle_type.ν₂, 1, mg.smoother)
cycle!(mg::MultigridIterable{C} where {C<:W}) = μ_cycle!(mg.grids, 2, mg.cycle_type.ν₁, mg.cycle_type.ν₂, 1, mg.smoother)
cycle!(mg::MultigridIterable{C} where {C<:F}) = F_cycle!(mg.grids, mg.cycle_type.ν₀, mg.cycle_type.ν₁, mg.cycle_type.ν₂, 1, mg.smoother)

function μ_cycle!(grids::Vector{G} where {G<:Grid}, μ::Int, ν₁::Int, ν₂::Int, grid_ptr::Int, smoother::Smoother)
    smooth!(grids[grid_ptr], ν₁, smoother)
    if grid_ptr == length(grids)
        grids[grid_ptr].x .= grids[grid_ptr].A\grids[grid_ptr].b # exact solve
    else
        grids[grid_ptr+1].b .= grids[grid_ptr].R * residu(grids[grid_ptr])
        grids[grid_ptr+1].x .= zero(grids[grid_ptr+1].x)
        for i in 1:μ
            μ_cycle!(grids, μ, ν₁, ν₂, grid_ptr+1, smoother)
        end
        coarse_grid_correction!(grids, grid_ptr, grids[grid_ptr+1].P * grids[grid_ptr+1].x)
    end
    smooth!(grids[grid_ptr], ν₂, smoother)
end

# damping factor chosen to minimize energy norm
function coarse_grid_correction!(grids, grid_ptr, c)
    d = residu(grids[grid_ptr])
    α = c'*d/(c'*grids[grid_ptr].A*c)
    α = isnan(α) ? one(eltype(c)) : min(1.1, max(0.7, α))
    grids[grid_ptr].x .+= α * c
end

function F_cycle!(grids::Vector{G} where {G<:Grid}, ν₀::Int, ν₁::Int, ν₂::Int, grid_ptr::Int, smoother::Smoother)
    if grid_ptr == length(grids)
        grids[grid_ptr].x .= zero(grids[grid_ptr].x)
    else
        grids[grid_ptr+1].b .= grids[grid_ptr].R * grids[grid_ptr].b
        F_cycle!(grids, ν₀, ν₁, ν₂, grid_ptr+1, smoother)
        grids[grid_ptr].x .= P(Cubic(),grids[grid_ptr+1].sz...) * grids[grid_ptr+1].x # FMG with cubic interpolation
    end
    for i in 1:ν₀
        μ_cycle!(grids, 1, ν₁, ν₂, grid_ptr, smoother)
    end
end


text:
module Peter
  import Base.+

  struct PeterNumber <: Number
    x::Number
  end

  Base.:+(p::PeterNumber, y::Number) = PeterNumber(p.x + y)

  export PeterNumber
end

module Jarrett
  import Base.+

  struct JarrettNumber <: Number
    y::Number
  end

  Base.:+(x::Number, j::JarrettNumber) = JarrettNumber(x + j.y)

  export JarrettNumber
end

using .Peter
using .Jarrett

p = PeterNumber(1.0) + 3.0
j = 2.0 + JarrettNumber(2.0)
f = p + j

println("Peter and Jarrett are friends! $f")


text:
function ueberpruefeAbsorptionskoeffizienten2x2(Zeile1,Zeile2,Spalte1,Spalte2,f1,f2,f3,f4)

    Zeile1 = 20-Zeile1;
    Zeile2 = 20-Zeile2;
    Spalte1 = 20-Spalte1;
    Spalte2 = 20-Spalte2;

    if (Zeile1 == f1+f3) && (Zeile2 == f2+f4) && (Spalte1 == f1+f2) && (Spalte2 == f3+f4)
        println("Ihr habt eine mögliche Verteilung der Absorptionskoeffizienten gefunden!")
    else
        println("Überprüft noch einmal, ob wirklich alle Zeilen- und Spaltensummen erfüllt sind.")
    end

end


text:
function MinMaxModel(A,b,C,d,XLB,XUB,VariableDef,YR)

    m,n = size(A);
    nObj = size(C)[1];
    MyModel = Model();
    @variable(MyModel,XLB[i] <= X[i=1:n] <= XUB[i],category = VariableType(VariableDef,i))#
    @variable(MyModel,Z)

    @objective(MyModel,Min,Z)
    # @constraintref Ucons[1:m];
    # @constraintref LAMBDAcons[1:nObj];
    for i = 1:m
         @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    for j = 1:nObj
         @constraint(MyModel,sum(C[j,i]*X[i] for i = 1:n) + d[j] - Z <= YR[j] )
    end
    for j = 1:nObj
         @constraint(MyModel,sum(C[j,i]*X[i] for i = 1:n) + d[j] >= YR[j])
    end
    return MyModel;

end

function MinMaxModel(A,b,C,d,XLB,XUB,VariableDef,YR; Cuts = [])

    m,n = size(A);
    nObj = size(C)[1];
    MyModel = Model();
    @variable(MyModel,XLB[i] <= X[i=1:n] <= XUB[i],category = VariableType(VariableDef,i))#
    @variable(MyModel,Z)

    @objective(MyModel,Min,Z)
    # @constraintref Ucons[1:m];
    # @constraintref LAMBDAcons[1:nObj];
    for i = 1:m
         @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    for j = 1:nObj
         @constraint(MyModel,sum(C[j,i]*X[i] for i = 1:n) + d[j] - Z <= YR[j] )
    end
    for j = 1:nObj
         @constraint(MyModel,sum(C[j,i]*X[i] for i = 1:n) + d[j] >= YR[j])
    end
    if !isempty(Cuts)
        for i in Cuts
            ZeroVals = setdiff(VariableDef[1]+VariableDef[2]+1:n,i.OneVals)
            @constraint(MyModel, sum(1-X[Int(j+VariableDef[1]+VariableDef[2])] for j in i.OneVals)
            + sum(X[Int(j)] for j = ZeroVals) >= 1)
        end
    end
    return MyModel;

end


text:
# varinterm.jl - basic operations to compute price change arrays. Operaciones
# básicas para computar vectores o matrices de variaciones intermensuales de
# índices de precios

"""
    varinterm(idx::AbstractVecOrMat, base_index = 100)

Función para computar un vector o una matriz de variaciones intermensuales de los
índices de precios en `idx`, utilizando como índice base `base_index` en la
primera observación.

Ver también: [`varinterm!`](@ref)
"""
function varinterm end

"""
    varinterm!([v::AbstractVecOrMat, ] idx::AbstractVecOrMat, base_index = 100)

Función para computar un vector o una matriz de variaciones intermensuales de
los índices de precios en `idx`, utilizando como índice base `base_index` en la
primera observación. Si `idx` es una matriz, `v` es opcional y el cómputo se
realiza sobre la misma matriz. Si `idx` es un vector, es necesario proporcionar
`v` para realizar el cómputo.

Ver también: [`varinterm`](@ref).
"""
function varinterm! end

## Definición de métodos

function varinterm!(v::AbstractVector, idx::AbstractVector, base_index::Real = 100)
    l = length(v)
    for i in l:-1:2
        @inbounds v[i] = 100 * (idx[i] / idx[i-1] - 1)
    end
    v[1] = 100 * (idx[1] / base_index - 1)
end

function varinterm(idx::AbstractVector, base_index::Real = 100)
    v = similar(idx)
    varinterm!(v, idx, base_index)
    v
end

function varinterm(cpimat::AbstractMatrix, base_index::Real = 100)
    c = size(cpimat, 2)
    vmat = similar(cpimat)
    for j in 1:c
        vcol = @view vmat[:, j]
        idxcol = @view cpimat[:, j]
        varinterm!(vcol, idxcol, base_index)
    end
    vmat
end

function varinterm(cpimat::AbstractMatrix, base_index::AbstractVector)
    c = size(cpimat, 2)
    vmat = similar(cpimat)
    for j in 1:c
        vcol = @view vmat[:, j]
        idxcol = @view cpimat[:, j]
        varinterm!(vcol, idxcol, base_index[j])
    end
    vmat
end

function varinterm!(cpimat::AbstractMatrix, base_index::Real = 100)
    c = size(cpimat, 2)
    for j in 1:c
        idxcol = @view cpimat[:, j]
        varinterm!(idxcol, idxcol, base_index)
    end
end

"""
    varinterm(base::IndexCPIBase)

Devuelve una nueva copia de tipo `VarCPIBase` de un `IndexCPIBase`.
"""
function varinterm(base::IndexCPIBase)
    VarCPIBase(base)
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
using JuMP, Ipopt

# Use nonlinear optimization to compute the maximum likelihood estimate (MLE)
# of the parameters of a normal distribution
# aka the sample mean and variance

n = 1000
data = randn(n)

m = Model(solver=IpoptSolver(print_level=0))

@variable(m, μ, start = 0.0)
@variable(m, σ >= 0.0, start = 1.0)

@NLobjective(m, Max, (n/2)*log(1/(2π*σ^2))-sum((data[i]-μ)^2 for i=1:n)/(2σ^2))

solve(m)

println("μ = ", getvalue(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", getvalue(σ)^2)
println("var(data) = ", var(data))
println("MLE objective: ", getobjectivevalue(m))

# constrained MLE?
@NLconstraint(m, μ == σ^2)

solve(m)
println("\nWith constraint μ == σ^2:")
println("μ = ", getvalue(μ))
println("σ^2 = ", getvalue(σ)^2)

println("Constrained MLE objective: ", getobjectivevalue(m))

max_stars_repo_name: econtoolkit/tomlab

max_issues_repo_name: econtoolkit/tomlab

max_forks_repo_name: econtoolkit/tomlab


text:
using GLMakie
# using WGLMakie

fig = Figure()
ax = LScene(fig[1, 1], scenekw=(show_axis=false,))

x = LinRange(0, π, 50);
y = LinRange(0, 2π, 100);
z = LinRange(0, 3π, 150);
t = LinRange(0.1, 1., 50);

lsgrid = labelslidergrid!(
  fig,
  ["yz plane - x axis", "xz plane - y axis", "xy plane - z axis", "time t"],
  [1:length(x), 1:length(y), 1:length(z), 1:length(t)]
);
fig[2, 1] = lsgrid.layout;

vol4 = [cos(X)*sin(Y)*sin(Z)*T for X ∈ x, Y ∈ y, Z ∈ z, T ∈ t];
vol = @view vol4[:, :, :, 1];
plt = volumeslices!(ax, x, y, z, vol)

# connect sliders to `volumeslices` update methods
sl_yz, sl_xz, sl_xy, sl_t = lsgrid.sliders

on(sl_yz.value) do v; plt[:update_yz][](v) end
on(sl_xz.value) do v; plt[:update_xz][](v) end
on(sl_xy.value) do v; plt[:update_xy][](v) end
on(sl_t.value) do v;
  plt[:volume][] = @view vol4[:, :, :, v];
  plt[:update_xy][](sl_xy.value[]);
  plt[:update_xz][](sl_xz.value[]);
  plt[:update_yz][](sl_yz.value[]);
end

set_close_to!(sl_yz, .5length(x))
set_close_to!(sl_xz, .5length(y))
set_close_to!(sl_xy, .5length(z))

# cam3d!(ax.scene, projectiontype=Makie.Orthographic)

fig

try
  while true
    for t_ ∈ eachindex(t)
      sl_t.value[] = t_
      sleep(0.1)
    end
    for t_ ∈ reverse(eachindex(t))
      set_close_to!(sl_t, t_)
      sleep(0.1)
    end
  end
catch InterruptException
  nothing
end


text:
############################################################################################
"""
$(TYPEDEF)
Initial distribution for state and duration of semi-Markov kernel.

# Fields
$(TYPEDFIELDS)
"""
struct SemiMarkovInitiation{A,B}
    "Initial distribution of state variable."
    state::A
    "Initial distribution of duration variable."
    duration::B
    function SemiMarkovInitiation(state::A, duration::B) where {A,B}
        return new{A,B}(state, duration)
    end
end

############################################################################################
"""
$(TYPEDEF)
Transition distributions for state and duration of semi-Markov kernel.

# Fields
$(TYPEDFIELDS)
"""
struct SemiMarkovTransition{A,B}
    "Transition distribution of state variable."
    state::A
    "Transition distribution of duration variable."
    duration::B
    function SemiMarkovTransition(state::A, duration::B) where {A,B}
        return new{A,B}(state, duration)
    end
end

############################################################################################
"""
$(TYPEDEF)
Semi-Markov Kernel for particle propagation.

# Fields
$(TYPEDFIELDS)
"""
struct SemiMarkov{A<:SemiMarkovInitiation,B<:SemiMarkovTransition,C} <: ParticleKernel
    "Initial distribution, function of iter only."
    initial::A
    "Transition distribution, function of full particle trajectory and current iteration count."
    transition::B
    "Data distribution to weight particles. Function of full data, particle trajectory and current iteration count."
    evidence::C
    function SemiMarkov(
        initial::A, transition::B, evidence::C
    ) where {A<:SemiMarkovInitiation,B<:SemiMarkovTransition,C}
        return new{A,B,C}(initial, transition, evidence)
    end
end

############################################################################################
function initial(_rng::Random.AbstractRNG, kernel::SemiMarkov)
    s = rand(_rng, kernel.initial.state)
    d = rand(_rng, kernel.initial.duration(s))
    return (s, d)
end

############################################################################################
function transition(
    _rng::Random.AbstractRNG, kernel::SemiMarkov, val::AbstractArray{P}, iter::Integer
) where {P}
    if val[iter - 1][2] > 0
        return (val[iter - 1][1], val[iter - 1][2] - 1)
    else
        #!NOTE: Inconsistent that s depends on whole Vector, but d only on current s - but no other way as of now
        sₜ = rand(_rng, kernel.transition.state(val, iter))
        dₜ = rand(_rng, kernel.transition.duration(sₜ, iter))
        return (sₜ, dₜ)
    end
end
############################################################################################
function ℓtransition(
    valₜ::Union{P,AbstractArray{P}},
    kernel::SemiMarkov,
    val::AbstractArray{P},
    iter::Integer,
) where {P}
    sₜ₋₁, dₜ₋₁ = val[iter - 1]
    sₜ, dₜ = valₜ
    #!NOTE: If duration at t-1 is 0, and particles states are not the same from t-1 to t
    if dₜ₋₁ == 0 && sₜ₋₁ != sₜ
        ℓπ = logpdf(kernel.transition.state(val, iter), sₜ) #current state given past state
        ℓπ += logpdf(kernel.transition.duration(sₜ, iter), dₜ) #duration given current state
        return ℓπ
    elseif (dₜ₋₁ - 1) == (dₜ) && sₜ₋₁ == sₜ
        return 0.0 #log(1.0) = 0
    else
        return -Inf
    end
end

############################################################################################
export SemiMarkovInitiation, SemiMarkovTransition, SemiMarkov


text:
function internal_logpdf(d::CompoundDirichlet, b_lens::Array{T},
                         int_leave_map::Vector{Int64})::T where T<:Real

    blen_int = 0.0
    blen_leave = 0.0
    blen_int_log = 0.0
    blen_leave_log = 0.0
    nterm = 0.0

    @views for i in eachindex(int_leave_map)
        if int_leave_map[i] === 1
            @inbounds blen_int += b_lens[i]
            @inbounds blen_int_log += log(b_lens[i])
        else
            @inbounds blen_leave += b_lens[i]
            @inbounds blen_leave_log += log(b_lens[i])
            nterm += 1
        end
    end

    t_l = blen_int+blen_leave
    n_int = nterm-3.0

    first = (d.alpha*log(d.beta))-log(gamma(d.alpha)) - (t_l*d.beta)
    second = -log(gamma(d.a))-log(gamma(d.c))+log(gamma(d.a+d.c))
    third = blen_leave_log*(d.a-1.0) + blen_int_log*(d.a*d.c-1.0)
    fourth = (d.alpha-d.a*nterm-d.a*d.c*n_int)*log(t_l)

    r2 = first + second + third + fourth

    return r2

end

function gradlogpdf(d::CompoundDirichlet, x::GeneralNode)
    int_ext = internal_external(x)
    blv = get_branchlength_vector(x)
    # use let block for proper capturing of variables
    f = let d=d, int_ext=int_ext
        y -> internal_logpdf(d, y, int_ext)
    end
    r = Zygote.pullback(f, blv)
    return r[1],r[2](1.0)[1]
end

function gradlogpdf(d::exponentialBL, x::GeneralNode)
    bl = get_branchlength_vector(x)
    g(y) = sum(logpdf.(Exponential(d.scale), y))
    r = Zygote.pullback(g, bl)
    r[1],r[2](1.0)[1]
end

function gradlogpdf(t::Union{UniformConstrained, UniformTopology, UniformBranchLength},
                    x::GeneralNode)::Tuple{Float64, Vector{Float64}}

    blv = get_branchlength_vector(x)
    0.0, zeros(length(blv))
end

function logpdf(d::CompoundDirichlet, x::GeneralNode)
    internal_logpdf(d, get_branchlength_vector(x), internal_external(x))
end

function logpdf(t::Union{UniformConstrained, UniformTopology, UniformBranchLength}, x::GeneralNode)
    0.0
end

function logpdf(ex::exponentialBL, x::GeneralNode)
    _logpdf(ex, x)
end

function _logpdf(d::exponentialBL, x::GeneralNode)
    bl = get_branchlength_vector(x)
    sum(logpdf.(Exponential(d.scale), bl))
end

function logpdf_sub(d::CompoundDirichlet, x::GeneralNode, transform::Bool)
    insupport(LengthDistribution(d), x) ? logpdf(d, x) : -Inf
end

function insupport(l::LengthDistribution, x::GeneralNode)
    bl = get_branchlength_vector(x)
    all(isfinite.(bl)) && topo_placeholder(x, l) && !any(isnan.(bl)) && all(0.0 .< bl)
     #&&
end

function insupport(t::UniformConstrained, x::GeneralNode)::Bool
    topological.constraint_dict(t, x)
end

function insupport(t::UniformTopology, x::GeneralNode)::Bool
    true
end

function insupport(t::UniformTopology, x::AbstractArray)::Bool
    true
end

function topo_placeholder(x::GeneralNode , l::LengthDistribution)
    true
end

function relistlength(d::CompoundDirichlet, x::AbstractArray)
  n = length(x)
  (Array(x), n)
end

"""
Strict Molecular Clock - BirthDeath
Implemented following Yang & Rannala 1997
doi.org/10.1093/oxfordjournals.molbev.a025811
"""
mutable struct BirthDeath <: ContinuousMultivariateDistribution
    s::Int64
    rho::Float64
    mu::Float64
    lambd::Float64
end # mutable struct

length(d::BirthDeath) = d.s - 1

function insupport(d::BirthDeath, t::AbstractVector{T}) where {T<:Real}
    length(d) == length(t) && all(isfinite.(t)) && all(0 .< t)
end # function

function _logpdf(d::BirthDeath, t::AbstractVector{T}) where {T<:Real}
    numerator::Float64 = (d.rho*(d.lambd-d.mu))/(d.rho*d.lambd +
                         (d.lambd*(1.0-d.rho)-d.mu)*exp(d.mu-d.lambd))
    denum::Float64 = d.rho*(d.lambd-d.mu)
    vt1::Float64 = log(1.0-((denum*exp(d.mu-d.lambd))/(d.rho * numerator)))
    f::Float64 = log((2.0^(s-1.0))/(factorial(s)*(s-1.0)))
    for i in t
        f += ((d.lambd+(1.0-d.rho)+2.0*(denum-numerator))+(d.mu-d.lambd)*i)-vt1
    end # for
    return f
end # function

"""
Strict Molecular Clock - Simplified Birth Death
Implemented folloing Yang & Rannala 1996
doi.org/10.1007/BF02338839
"""
mutable struct BirthDeathSimplified <: ContinuousMultivariateDistribution
    s::Int64
    mu::Float64
    lambd::Float64
end # mutable struct

length(d::BirthDeathSimplified) = d.s - 1

function insupport(d::BirthDeathSimplified, t::AbstractVector{T}) where {T<:Real}
    length(d) == length(t) && all(isfinite.(t)) && all(0 .< t)
end # function

function _logpdf(d::BirthDeathSimplified, t::AbstractVector{T}) where {T<:Real}

    f::Float64 = log(2.0)*(d.s-1.0)+log(d.mu)*(d.s-2.0)
    p0::Float = (d.s-2.0)*log((d.mu*(1.0-exp(-(d.lambd-d.mu)))) /
                (d.lambd-d.mu*exp(-(d.lambd-d.mu))))
    p0 += log(factorial(d.s)*(d.s-1.0))
    f -= po

    con::Float64 = 2.0*log(d.lambd - d.mu)

    for i in t
        f += (con-(d.lambd-d.mu)*i) - log(d.lambd-d.mu*exp(-(d.lambd-d.mu)*i))*2.0
    end # for
    return f
end # function

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
export CDS_MVp

"""
compute single-thread matrix vector product with vector x, output is vector y: y=A*x
MVP is in the compressed diagonal format.
R is a tall matrix N by ndiagonals, corresponding to a square matrix A
offsets indicate offset of diagonal compared to the main diagonal in A (which is 0)
"""
function CDS_MVp(
                        N::Integer,
                        ndiags::Integer,
                        R::Array{TF,2},
                        offset::Vector{TI},
                        x::Vector{TF},
                        y::Vector{TF}) where {TF<:Real,TI<:Integer}

  for i = 1 : ndiags
      d = offset[i]
      r0 = max(1, 1-d)
      r1 = min(N, N-d)
      c0 = max(1, 1+d)
       for r = r0 : r1
         c = r - r0 + c0 #original
       @inbounds y[r] = y[r] + R[r,i] * x[c]#original
      end
  end
  return y
end


text:
include("../../testcases.jl")
using PyCall
using DFTK

atoms = [ElementPsp(:Fe, psp=load_psp(iron_bcc.psp)) => iron_bcc.positions]
model = model_DFT(iron_bcc.lattice, atoms, :lda_xc_teter93,
                  temperature=iron_bcc.temperature,
                  smearing=Smearing.FermiDirac(), spin_polarization=:collinear)

abinitpseudos = [joinpath(@__DIR__, "Fe-q8-lda.abinit.hgh")]
DFTK.run_abinit_scf(model, @__DIR__; abinitpseudos=abinitpseudos,
                    Ecut=20, kgrid=[4, 4, 4], nshiftk=1, shiftk=[0.5, 0.5, 0.5],
                    spinat=[0 0 4.0], n_bands=10, tol=1e-10, iscf=3)
# iscf == 3: Use Anderson mixing instead of minimization


text:
#    Problem :
#    *********
#
#    The Boundary Value problem.
#    This is the nonlinear least-squares version without fixed variables.
#
#
#    Source:  problem 28 in
#    J.J. More', B.S. Garbow and K.E. Hillstrom,
#    "Testing Unconstrained Optimization Software",
#    ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
#
#    Implementation translated from:
#    See also Buckley#17 (p. 94):
#
#
#    SIF input: Ph. Toint, Dec 1989 and Nick Gould, Oct 1992.
#
#    classification SUR2-MN-V-0
#    The number of variables is N.
#
# Daniel Henderson, 08/2021

# f = x -> begin

# end

# g! = (g, x) -> begin
# 	n = lastindex(x)
# end

# fg! = (g, x) -> begin

# end

# init = (n::Int=5000) -> begin

# end

# @warn "FRUEROTH needs to be implemented"

# TestSet["FREUROTH"] = UncProgram("FREUROTH",  f, g!, fg!, init)


text:
using ReverseDiff
using CatViews

x = [-1. -2; -1 -2]
y = [1.; 2]
w = [2.; 3]
b = [-3.]

type Vars
    lookup::Dict
    grad_lookup::Dict
    flat
    grad_flat

    Vars() = Vars(Dict(), Dict(), nothing, nothing)
end
getindex(vars::Vars, varname::String) = vars.lookup[varname]
setindex!(vars::Vars, varname::String, var) = vars.lookup[varname] = var
function compile!(vars::Vars)
    for (k, v) in vars.lookup
        vars.grad_lookup[k] = similar(v)
    end
    vars.flat = CatView(map(x -> view(x, :), values(vars.lookup)))
    vars.grad_flat = CatView(map(x -> view(x, :), values(vars.grad_lookup)))
end

vars = Vars()
vars[:w] = w
vars[:b] = b
compile!(vars)

ReverseDiff.@forward sigmoid(x) = 1. ./ (1. + exp(-x))
mse(a, y) = mean((a - y) .* (a - y))
#= mse(a, y) = mean((a - y) .^ 2) =#
f1(x, y, vars::Vars) = mse(sigmoid(x*vars[:w] .+ vars[:b]), y)

out = f1(x, y, vars)
println(out)

inputs = (x, y, vars)
result = map(similar, inputs)
result_view = CatView(map(x -> view(x, :), result))
f1! = ReverseDiff.compile_gradient(f1, result)

@time f1!(result, inputs)
@time f1!(result, inputs)

for r in result
    println(r)
end

println(result_view)

#= Expected output =#
#=  =#
#= [array([[ -3.34017280e-05,  -5.01025919e-05], =#
#=        [ -6.68040138e-05,  -1.00206021e-04]]), array([[ 0.9999833], =#
#=        [ 1.9999833]]), array([[  5.01028709e-05], =#
#=        [  1.00205742e-04]]), array([ -5.01028709e-05])] =#


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "QPALM"
version = v"0.3.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Benny44/QPALM_vLADEL.git", "e272603e833db9c6b6de0a7ca4630da2a7fd7268"),
    GitSource("https://github.com/Benny44/LADEL.git", "5af7a8d6a7ad76aeb00edf2f673fb3e6b703c5f2"),
    GitSource("https://github.com/kul-optec/QPALM.jl.git", "7b4424ce17ba232aa6e89648032c93572ca5a3e1"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd QPALM_vLADEL/
ln -s ${WORKSPACE}/srcdir/LADEL/* ${WORKSPACE}/srcdir/QPALM_vLADEL/LADEL/
ln -s ${WORKSPACE}/srcdir/QPALM.jl/* ${WORKSPACE}/srcdir/QPALM_vLADEL/interfaces/QPALM.jl/
mkdir -p build
cd build
mkdir lib debug
cd debug  # just following the instructions here

[[ $target == *-mingw* ]] && BLAS=libopenblas.${dlext}.a || BLAS=libopenblas.${dlext}
[[ $target == *-linux-* ]] && std="gnu99" || std="c99"

cmake ../.. \
      -DCMAKE_BUILD_TYPE=release \
      -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DUSE_LADEL=ON \
      -DUSE_CHOLMOD=ON \
      -DINTERFACES=OFF \
      -DPYTHON=OFF \
      -DCOVERAGE=OFF \
      -DUNITTESTS=OFF \
      -DJULIA=ON \
      -DLAPACKE=${prefix}/lib/${BLAS} \
      -DLAPACK=${prefix}/lib/${BLAS} \
      -DCMAKE_C_FLAGS="-I${includedir} -std=${std}"

make

# make install doesn't appear to work...
if [[ $target == *-mingw* ]]; then
  # headers don't make it to build/include
  cp ../../include/* $prefix/include
  cp bin/* $prefix/bin
else
  cp ../include/* $prefix/include
  cp ../lib/* $prefix/lib
fi
cd ../..
mv LICENSE LICENSE.QPALM
install_license LICENSE.QPALM
mv ${WORKSPACE}/srcdir/LADEL/LICENSE ./LICENSE.LADEL
install_license LICENSE.LADEL
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
                   LibraryProduct("libladel", :libladel),
                   LibraryProduct("libqpalm", :libqpalm),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="SuiteSparse_jll", uuid="bea87d4a-7f5b-5778-9afe-8cc45184846c"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.3")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil