text:
#: Optimization functions and solvers

export OptFISTA, opt_fista

# Options type for FISTA

mutable struct OptFISTA{T<:Real}<:OptimOptions
    steplength::T
    niter::Int64
    tol_x::Union{Nothing,T}
    nesterov::Bool
end

opt_fista(steplength::T; niter::Int64=1000, tol_x::Union{Nothing,T}=nothing, nesterov::Bool=true) where {T<:Real} = OptFISTA{T}(steplength, niter, tol_x, nesterov)

# Generic FISTA solver

"""
Solver for the regularized problem via FISTA-like gradient-projections:
```math
min_x f(x)+g(x)
```
Reference: Beck, A., and Teboulle, M., 2009, A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems
https://www.ceremade.dauphine.fr/~carlier/FISTA
"""
function minimize_fista!(fun::DiffPlusProxFun{CT,N}, initial_estimate::AbstractArray{CT,N}, steplength::T, niter::Int64, nesterov::Bool, tol_x::Union{Nothing,T}, x0::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}

    # Initialization
    x0   .= initial_estimate
    xtmp  = similar(x0)
    x     = similar(x0)
    t0 = T(1)

    # Optimization loop
    for i = 1:niter

        grad!(fun.f, x0, xtmp; eval=false)         # Compute gradient
        xtmp .= x0-steplength*xtmp               # Gradient update
        # xtmp .*= -steplength; xtmp .+= x0          # Gradient update
        proxy!(xtmp, steplength, fun.g, x)         # Compute proxy
        if nesterov                                # > Nesterov two-step update:
            t = T(0.5)*(T(1)+sqrt(T(1)+T(4)*t0^2)) # - compute dynamic step size
            x .+= (t0-T(1))/t*(x-x0)               # - update momentum
            # x .*= (t+t0-T(1))/t; x .-= (t0-T(1))/t*x0 # - update momentum
            t0 = t                                 # - update step size
        end                                        # <

        # Update unknowns
        ~isnothing(tol_x) && norm(x-x0)/norm(x0) < tol_x && break
        x0 .= x

    end

    # Return output
    return x0

end

function minimize_fista_debug!(fun::DiffPlusProxFun{CT,N}, initial_estimate::AbstractArray{CT,N}, steplength::T, niter::Int64, nesterov::Bool, tol_x::Union{Nothing,T}, x0::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}

    # Initialization
    x0 .= initial_estimate
    xtmp = similar(x0)
    x = similar(x0)
    t0 = T(1)
    fval_hist = Array{T,1}(undef, niter)
    err_rel = Array{T,1}(undef, niter)

    # Optimization loop
    for i = 1:niter

        fval = grad!(fun.f, x0, xtmp; eval=true)   # Compute gradient
        xtmp .= x0-steplength*xtmp                 # Gradient update
        proxy!(xtmp, steplength, fun.g, x)         # Compute proxy
        if nesterov                                # > Nesterov two-step update:
            t = T(0.5)*(T(1)+sqrt(T(1)+T(4)*t0^2)) # - compute dynamic step size
            x .= x+(t0-T(1))/t*(x-x0)              # - update momentum
            t0 = t                                 # - update step size
        end                                        # <

        # Convergence check
        fval_hist[i] = fval
        err_rel[i] = norm(x-x0)/norm(x0)
        ~isnothing(tol_x) && err_rel[i] < tol_x && break

        # Update unknowns
        x0 .= x

    end

    # Return output
    return x0, fval_hist, err_rel

end

minimize!(fun::DiffPlusProxFun{CT,N}, x0::AbstractArray{CT,N}, opt::OptFISTA{T}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = minimize_fista!(fun, x0, opt.steplength, opt.niter, opt.nesterov, opt.tol_x, x)

minimize_debug!(fun::DiffPlusProxFun{CT,N}, x0::AbstractArray{CT,N}, opt::OptFISTA{T}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = minimize_fista_debug!(fun, x0, opt.steplength, opt.niter, opt.nesterov, opt.tol_x, x)

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
using Test
using TreeTools
using TreeKnit, TreeKnit.SplitGraph

println("##### main #####")

nwk1 = "((1:1,2:1):1,((3:1,4:1):1,5:1,6:1):1)"
nwk2 = "((6:1,2:1):1,((3:1,4:1):1,5:1,1:1):1)"

t1 = node2tree(TreeTools.parse_newick(nwk1))
t2 = node2tree(TreeTools.parse_newick(nwk2))

solutions = [
	[["1"], ["2"], ["3", "4", "5", "6"]],
	[["2"], ["6"], ["1", "3", "4", "5"]],
	[["1"], ["6"], ["2", "3", "4", "5"]],
	[["1"], ["2", "6"], ["3", "4", "5"]],
	[["6"], ["1", "2"], ["3", "4", "5"]],
]
@testset "Dyn-resolve" begin
	MCCs = TreeKnit.computeMCCs(
		t1, t2, OptArgs(likelihood_sort=false, γ=3);
	)
	@test in(MCCs, solutions)
end

ext: jl

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
function update_german_id()
    item("Update German ID")
    dic = Dict()
    for line in eachline(joinpath(work_dir, "data/pedigree/german.ref"))
        a, b = split(line)
        dic[a] = b
    end
    til = joinpath(work_dir, "data/genotypes/step-6.plk/german.fam")
    fra = joinpath(work_dir, "tmp/tt")
    mv(til, fra, force=true)
    open(til, "w") do io
        for line in eachline(fra)
            a, b, c, d, e, f = split(line)
            write(io, "$a ", dic[b], " $c $d $e $f\n")
        end
    end
end

"""
    merge_by_country(country, platform_list)
---
Merge the platforms by country.
"""
function merge_by_country(country::AbstractString, list)
    cd(work_dir)
    tmp = joinpath(work_dir, "tmp")
    empty_dir(tmp)
    fra = "data/genotypes/step-5.plk"
    dto = "data/genotypes/step-6.plk"
    ref = joinpath(work_dir, "data/maps/ref.allele")
    lst = joinpath(tmp, "merge.lst")
    isdir(dto) || mkdir(dto)

    item("Merge $country data")
    open("tmp/merge.lst", "w") do io
        for l in list
            write(io, "$fra/$country-$l\n")
        end
    end

    _ = read(`plink --cow
                    --merge-list "tmp/merge.lst"
                    --reference-allele $ref
                    --make-bed
                    --out $dto/$country`,
             String)
    #merge_beds("tmp/merge.lst", "$dto/$country")
    done()
end

"""
    imputation_with_beagle(countries)
---
1. Convert `bed` to `vcf`
2. Impute with `beagle.jar`
3. Convert back them to `bed`
"""
function imputation_with_beagle()
    cd(work_dir)
    fra = joinpath(work_dir, "data/genotypes/step-6.plk")
    dto = joinpath(work_dir, "data/genotypes/step-7.plk")
    tmp = joinpath(work_dir, "tmp")
    ref = joinpath(work_dir, "data/maps/ref.allele")
    isdir(dto) || mkdir(dto)
    empty_dir(tmp)
    countries = ["dutch", "german", "norge", "two"]

    title("Impute with beagle.jar")
    for country in countries
        item("Impute $country data")
        plink_2_vcf("$fra/$country", "tmp/miss")
        _ = read(`java -jar $beagle
                        gt=tmp/miss.vcf
                        ne=100
                        out=tmp/imp`,
                 String)
        _ = read(`plink --cow
			--vcf $tmp/imp.vcf.gz
                        --reference-allele $ref
			--const-fid
                        --make-bed
			--out $dto/$country`,
                 String)
        done()
    end
end

"""
    merge_dg()
---
Merge Dutch and German data together for imputation, as they are connected.
"""
function merge_dg()
    cd(work_dir)
    fra = joinpath(work_dir, "data/genotypes/step-6.plk")
    tmp = joinpath(work_dir, "tmp")
    ref = joinpath(work_dir, "data/maps/ref.allele")
    empty_dir(tmp)

    title("Merge dutch and german data")
    lst = joinpath(tmp, "list")
    write(lst,
          joinpath(fra, "dutch"), '\n',
          joinpath(fra, "german"), '\n')
    out = joinpath(fra, "two")
    _ = read(`plink --cow
                    --merge-list $lst
                    --reference-allele $ref
                    --make-bed
                    --out $out`,
             String)
    done()
end

"""
    rm_country_specific_snp()
---
The final step for `G` calculation. After imputation, remove country specific SNP
to calculate cross-country `G` matrix.
"""
function rm_country_specific_snp()
    tmp = joinpath(work_dir, "tmp")
    fra = joinpath(work_dir, "data/genotypes/step-7.plk")
    til = joinpath(work_dir, "data/genotypes/step-8.plk")
    ref = joinpath(work_dir, "data/maps/ref.allele")
    isdir(til) || mkdir(til)
    lst = joinpath(tmp, "snp.list")
    empty_dir(tmp)

    function snp_set(bim)
        snp = AbstractString[]
        for line in eachline(bim)
            push!(snp, split(line)[2])
        end
        return Set(snp)
    end

    title("Remove country specific SNP")
    item("Create shared SNP list")
    sd = snp_set(joinpath(fra, "dutch.bim"))
    sg = snp_set(joinpath(fra, "german.bim"))
    sn = snp_set(joinpath(fra, "norge.bim"))
    ss = intersect(sd, sg, sn)

    open(lst, "w") do io
        for snp in ss
            write(io, snp, '\n')
        end
    end
    for country in ["dutch", "german", "norge", "two"]
        # bed_snp_subset(joinpath(fra, country),
        #                lst,
        #                joinpath(til, country)
        #                )
        src = joinpath(fra, country)
        out = joinpath(til, country)
        _ = read(`plink --cow
                        --bfile $src
                        --extract $lst
                        --reference-allele $ref
                        --make-bed
                        --out $out`,
                 String)
    end
end


text:
include("BoundingBoxes.jl")
include("SynapseTables.jl")
include("VectorClouds.jl")
include("RangeIndexingArrays.jl")


text:
#!/usr/bin/env julia
using Luxor, Colors
fname = "/tmp/arc-twopoints.pdf"
Drawing(1200, 1200, fname)
origin()
background("white")
sethue("magenta")
setopacity(0.6)
setline(2)
tiles = Tiler(1200, 1200, 8, 8)
for (pos, n) in tiles
  _, pt2, pt3 = ngon(pos, rand(10:tiles.tilewidth/2), 3, rand(0:pi/12:2pi), vertices=true)
  sethue("black")
  if n % 4 == 0
    pt3 *= 1.1 # test for points not on arc
  end
  map(pt -> circle(pt, 4, :fill), [pos, pt3])
  sethue("red")
  circle(pt2, 3, :fill)
  randomhue()
  arc2r(pos, pt2, pt3, :stroke)
end
finish()
println("finished test: output in $(fname)")


text:
using JULES, Oceananigans
using Plots
using Printf

Nx = Ny = 1
Nz = 32
L = 10e3

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(2, 2, 2),
                            x=(0, L), y=(0, L), z=(0, L))

pₛ = 100000
Tₐ = 293.15
g  = 9.80665

buoyancy = IdealGas()
Rᵈ, cₚ = buoyancy.Rᵈ, buoyancy.cₚ

####
#### Isothermal atmosphere
####

H = Rᵈ * Tₐ / g    # Scale height [m]
ρₛ = pₛ / (Rᵈ*Tₐ)  # Surface density [kg/m³]

p₀(x, y, z) = pₛ * exp(-z/H)
ρ₀(x, y, z) = ρₛ * exp(-z/H)

θ₀(x, y, z) = Tₐ * exp(z/H * Rᵈ/cₚ)
Θ₀(x, y, z) = ρ₀(x, y, z) * θ₀(x, y, z)

####
#### Sponge layer forcing to damp vertically propagating acoustic waves at the top.
####

const τ⁻¹ = 1    # Damping/relaxation time scale [s⁻¹]. This is very strong damping.
const Δμ = 0.1L  # Sponge layer width [m] set to 10% of the domain height.
@inline μ(z, Lz) = τ⁻¹ * exp(-(Lz-z) / Δμ)

@inline Fw(i, j, k, grid, t, Ũ, C̃, p) = @inbounds -μ(grid.zF[k], grid.Lz) * Ũ.ρw[i, j, k]
forcing = ModelForcing(w=Fw)

####
#### Create model
####

model = CompressibleModel(
                      grid = grid,
                  buoyancy = buoyancy,
        reference_pressure = pₛ,
    thermodynamic_variable = ModifiedPotentialTemperature(),
                   tracers = (:Θᵐ,),
                   forcing = forcing
)

####
#### Set initial conditions
####

set!(model.density, ρ₀)
set!(model.tracers.Θᵐ, Θ₀)

# Initial profiles
Θ₀_prof = model.tracers.Θᵐ[1, 1, 1:Nz]
ρ₀_prof = model.density[1, 1, 1:Nz]

# Arrays we will use to store a time series of ρw(z = L/2).
times = [model.clock.time]
ρw_ts = [model.momenta.ρw[1, 1, Int(Nz/2)]]

####
#### Time step and keep plotting vertical profiles of ρθ′, ρw, and ρ′.
####

# @animate for i=1:100
while model.clock.time < 500
    time_step!(model; Δt=0.5, Nt=10)

    @show model.clock.time
    Θ_prof = model.tracers.Θᵐ[1, 1, 1:Nz] .- Θ₀_prof
    W_prof = model.momenta.ρw[1, 1, 1:Nz+1]
    ρ_prof = model.density[1, 1, 1:Nz] .- ρ₀_prof

    Θ_plot = plot(Θ_prof, grid.zC, xlim=(-2.5, 2.5), xlabel="Theta_prime", ylabel="z (m)", label="")
    W_plot = plot(W_prof, grid.zF, xlim=(-1.2, 1.2), xlabel="rho*w", label="")
    ρ_plot = plot(ρ_prof, grid.zC, xlim=(-0.01, 0.01), xlabel="rho", label="")

    t_str = @sprintf("t = %d s", model.clock.time)
    display(plot(Θ_plot, W_plot, ρ_plot, title=["" "$t_str" ""], layout=(1, 3), show=true))

    push!(times, model.clock.time)
    push!(ρw_ts, model.momenta.ρw[1, 1, Int(Nz/2)])
end

####
#### Plot the initial profile and the hydrostatically balanced profile.
####

ρ∞_prof = model.density[1, 1, 1:Nz]

θ₀_prof = Θ₀_prof ./ ρ₀_prof
θ∞_prof = model.tracers.Θᵐ[1, 1, 1:Nz] ./ ρ∞_prof

θ_plot = plot(θ₀_prof, grid.zC, xlabel="theta (K)", ylabel="z (m)", label="initial", legend=:topleft)
plot!(θ_plot, θ∞_prof, grid.zC, label="balanced")

ρ_plot = plot(ρ₀_prof, grid.zC, xlabel="rho (kg/m³)", ylabel="z (m)", label="initial")
plot!(ρ_plot, ρ∞_prof, grid.zC, label="balanced")

display(plot(θ_plot, ρ_plot, show=true))

####
#### Plot timeseries of maximum ρw
####

plot(times, ρw_ts, linewidth=2, xlabel="time (s)", ylabel="rho*w(z=$(Int(L/2)) m)", label="", show=true)


text:
# White Noise kernel

"""
    Noise <: Kernel

Noise kernel (covariance)
```math
k(x,x') = σ²δ(x-x'),
```
where ``δ`` is the Kronecker delta function and ``σ`` is the signal standard deviation.
"""
mutable struct Noise{T<:Real} <: Kernel
    "Signal variance"
    σ2::T
    "Priors for kernel parameters"
    priors::Array
end

"""
White Noise kernel

    Noise(lσ::Real)

# Arguments
  - `lσ::Real`: signal standard deviation (given on log scale)
"""
Noise(lσ::T) where T = Noise{T}(exp(2 * lσ), [])

cov(noise::Noise, sameloc::Bool) = sameloc ? noise.σ2 : zero(noise.σ2)
cov(noise::Noise, x::AbstractVector, y::AbstractVector) = cov(noise, x ≈ y)
@inline @inbounds function cov_ij(noise::Noise, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, dim::Int)
    @inbounds for z in 1:dim
        if !(X1[z,i] ≈ X2[z,j])
            return zero(noise.σ2)
        end
    end
    return noise.σ2
    # return cov(noise, all(X1[z,i] ≈ X2[z,j] for z in 1:dim)) # this is inefficient for some reason
end

get_params(noise::Noise{T}) where T = T[log(noise.σ2) / 2]
get_param_names(noise::Noise) = [:lσ]
num_params(noise::Noise) = 1

function set_params!(noise::Noise, hyp::AbstractVector)
    length(hyp) == 1 || throw(ArgumentError("Noise kernel has one parameter, received $(length(hyp))."))
    noise.σ2 = exp(2 * hyp[1])
end

@inline dKij_dθp(noise::Noise, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int, dim::Int) =
    # 2 * cov(noise, view(X, :, i), view(X, :, j))
    2 * cov_ij(noise, X1, X2, i, j, dim)


text:
"""1INFILTR case from Hydrus-1D, horizontal."""
module Example1INFILTR

using Fronts
using Fronts.PorousModels
using Plots

ϵ = 1e-6

Ks = 25  # cm/h
α = 0.01433  # 1/cm
n = 1.506

θs = 0.3308

model = VanGenuchten(n=n, α=α, Ks=Ks, θs=θs)

prob = DirichletProblem(model, i=0.1003, b=0.3308-ϵ)

θ = solve(prob)

r = range(0, 200, length=500)

plt = plot(r, r -> θ(r,1), label="t=1 h", xguide="r [cm]", yguide="θ")
plot!(r, r -> θ(r,4), label="t=4 h")
display(plt)

plt = plot(r, r -> flux(θ,r,1), label="t=1 h", xguide="r [cm]", yguide="U [m/s]")
plot!(r, r -> flux(θ,r,4), label="t=4 h")
display(plt)

end


text:
@testset "accumulation.jl" begin
    @testset "is_inplaceable_destination" begin
        is_inplaceable_destination = ChainRulesCore.is_inplaceable_destination

        @test is_inplaceable_destination([1, 2, 3, 4])
        @test !is_inplaceable_destination(1:4)

        @test is_inplaceable_destination(Diagonal([1, 2, 3, 4]))
        @test !is_inplaceable_destination(Diagonal(1:4))

        @test is_inplaceable_destination(view([1, 2, 3, 4], :, :))
        @test !is_inplaceable_destination(view(1:4, :, :))

        @test is_inplaceable_destination(falses(4))
        @test is_inplaceable_destination(spzeros(4))
        @test is_inplaceable_destination(spzeros(2, 2))

        @test !is_inplaceable_destination(1.3)
        @test !is_inplaceable_destination(@SVector [1, 2, 3])
        @test !is_inplaceable_destination(Hermitian([1 2; 2 4]))
        @test !is_inplaceable_destination(Symmetric([1 2; 2 4]))
    end

    @testset "add!!" begin
        @testset "scalar" begin
            @test 16 == add!!(12, 4)
        end

        @testset "misc AbstractTangent subtypes" begin
            @test 16 == add!!(12, @thunk(2 * 2))
            @test 16 == add!!(16, ZeroTangent())

            @test 16 == add!!(16, NoTangent())  # Should this be an error?
        end

        @testset "add!!(::AbstractArray, ::AbstractArray)" begin
            @testset "LHS Array (inplace)" begin
                @testset "RHS Array" begin
                    A = [1.0 2.0; 3.0 4.0]
                    accumuland = -1.0 * ones(2, 2)
                    ret = add!!(accumuland, A)
                    @test ret === accumuland  # must be same object
                    @test accumuland == [0.0 1.0; 2.0 3.0]
                end

                @testset "RHS StaticArray" begin
                    A = @SMatrix [1.0 2.0; 3.0 4.0]
                    accumuland = -1.0 * ones(2, 2)
                    ret = add!!(accumuland, A)
                    @test ret === accumuland  # must be same object
                    @test accumuland == [0.0 1.0; 2.0 3.0]
                end

                @testset "RHS Diagonal" begin
                    A = Diagonal([1.0, 2.0])
                    accumuland = -1.0 * ones(2, 2)
                    ret = add!!(accumuland, A)
                    @test ret === accumuland  # must be same object
                    @test accumuland == [0.0 -1.0; -1.0 1.0]
                end
            end

            @testset "add!!(::StaticArray, ::Array) (out of place)" begin
                A = [1.0 2.0; 3.0 4.0]
                accumuland = @SMatrix [-1.0 -1.0; -1.0 -1.0]
                ret = add!!(accumuland, A)
                @test ret == [0.0 1.0; 2.0 3.0]  # must return right answer
                @test ret !== accumuland  # must not be same object
                @test accumuland == [-1.0 -1.0; -1.0 -1.0]  # must not have changed
            end

            @testset "add!!(::Diagonal{<:Vector}, ::Diagonal{<:Vector}) (inplace)" begin
                A = Diagonal([1.0, 2.0])
                accumuland = Diagonal([-2.0, -2.0])
                ret = add!!(accumuland, A)
                @test ret === accumuland  # must be same object
                @test accumuland == Diagonal([-1.0, 0.0])
            end

            @testset "Unhappy Path" begin
                # wrong length
                @test_throws DimensionMismatch add!!(ones(4, 4), ones(2, 2))
                # wrong shape
                @test_throws DimensionMismatch add!!(ones(4, 4), ones(16))
                # wrong type (adding scalar to array)
                @test_throws MethodError add!!(ones(4), 21.0)
            end
        end

        @testset "AbstractThunk $(typeof(thunk))" for thunk in (
            @thunk(-1.0 * ones(2, 2)),
            InplaceableThunk(x -> x .-= ones(2, 2), @thunk(-1.0 * ones(2, 2))),
        )
            @testset "in place" begin
                accumuland = [1.0 2.0; 3.0 4.0]
                ret = add!!(accumuland, thunk)
                @test ret == [0.0 1.0; 2.0 3.0]  # must return right answer
                @test ret === accumuland  # must be same object
            end

            @testset "out of place" begin
                accumuland = @SMatrix [1.0 2.0; 3.0 4.0]

                ret = add!!(accumuland, thunk)
                @test ret == [0.0 1.0; 2.0 3.0]  # must return right answer
                @test ret !== accumuland  # must not be same object
                @test accumuland == [1.0 2.0; 3.0 4.0]  # must not have mutated
            end
        end

        @testset "not actually inplace but said it was" begin
            # thunk should never be used in this test
            ithunk = InplaceableThunk(@thunk(@assert false)) do x
                77 * ones(2, 2)  # not actually inplace (also wrong)
            end
            accumuland = ones(2, 2)
            @assert ChainRulesCore.debug_mode() == false
            # without debug being enabled should return the result, not error
            @test 77 * ones(2, 2) == add!!(accumuland, ithunk)

            ChainRulesCore.debug_mode() = true  # enable debug mode
            # with debug being enabled should error
            @test_throws ChainRulesCore.BadInplaceException add!!(accumuland, ithunk)
            ChainRulesCore.debug_mode() = false  # disable it again
        end
    end

    @testset "showerror BadInplaceException" begin
        BadInplaceException = ChainRulesCore.BadInplaceException
        ithunk = InplaceableThunk(x̄ -> nothing, @thunk(@assert false))
        msg = sprint(showerror, BadInplaceException(ithunk, [22], [23]))
        @test occursin("22", msg)

        msg_equal = sprint(showerror, BadInplaceException(ithunk, [22], [22]))
        @test occursin("equal", msg_equal)
    end
end


text:
@testset "Distance" begin
    a = randn(Float,10)
    b = randn(Float,10)
    c1 = [a b a b]
    c2 = [a b]
    c3 = [a b b b]

    @test VectorDataUtils.sql2 <: VectorDataUtils.Distance
    @test VectorDataUtils.l2 <: VectorDataUtils.Distance
    @test VectorDataUtils.l1 <: VectorDataUtils.Distance

    @test VectorDataUtils.evaluate(VectorDataUtils.sql2(),a,b) ≈ VectorDataUtils.sql2dist(a,b)
    @test VectorDataUtils.evaluate(VectorDataUtils.l2(),a,b) ≈ VectorDataUtils.l2dist(a,b)
    @test VectorDataUtils.evaluate(VectorDataUtils.l1(),a,b) ≈ VectorDataUtils.l1dist(a,b)

    distances = (VectorDataUtils.sql2(), VectorDataUtils.l2(), VectorDataUtils.l1())
    for dm ∈ distances

        d = VectorDataUtils.colwise(dm,c1,c3)
        @test length(d) == size(c1,2)
        @test d[1] == 0.0
        @test d[2] == 0.0
        @test d[3] == VectorDataUtils.evaluate(dm,a,b)

        d = VectorDataUtils.colwise(dm,a,c1)
        @test length(d) == size(c1,2)
        @test d[1] == 0.0
        @test d[2] == VectorDataUtils.evaluate(dm,a,b)
        @test d[3] == 0.0
        @test d[4] == d[2]

        d = VectorDataUtils.pairwise(dm,a,c1)
        @test length(d) == size(c1,2)
        @test d[1] == 0.0
        @test d[2] == VectorDataUtils.evaluate(dm, a, b)

        d = VectorDataUtils.pairwise(dm,c1,a)
        @test length(d) == size(c1,2)
        @test d[3] == 0.0
        @test d[4] == VectorDataUtils.evaluate(dm, a, b)

        dd = VectorDataUtils.pairwise(dm, c2, c1)
        @test size(dd,1) == size(c2,2)
        @test size(dd,2) == size(c1,2)
        @test dd[1,1] == 0.0
        @test dd[2,2] == 0.0
        @test dd[1,2] == VectorDataUtils.evaluate(dm, a, b)
        @test dd[2,1] == VectorDataUtils.evaluate(dm, b, a)
    end
end


text:
# %%
using Tmatrix
using Profile

rx, rz = 1e-7, 1.3e-7
n_max = 8
k1 = Complex(6.2e6)
k2 = 7e6 + 1e3 * im
rotationally_symmetric = true

@profile calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2; rotationally_symmetric=rotationally_symmetric)
Profile.print()

"""
using VectorSphericalWaves
n_θ_points, n_ϕ_points = 10000,20000
θ_array, ϕ_array = Tmatrix.meshgrid_θ_ϕ(n_θ_points, n_ϕ_points; min_θ=1e-16, min_ϕ=1e-16, rotationally_symmetric=rotationally_symmetric)
r_array, n̂_array = Tmatrix.ellipsoid(rx, rz, θ_array);
VectorSphericalWaves.N_mn_wave_SVector.(0,1,k2*r_array, θ_array,ϕ_array; kind= "regular")
VectorSphericalWaves.N_mn_wave_SVector(0,1,k2*r_array[1], θ_array[1],ϕ_array[1]; kind= "regular")
@time for i = 1:1e3; VectorSphericalWaves.N_mn_wave_SVector.(0,1,k2*r_array, θ_array,ϕ_array; kind= "regular"); end
@time for i = 1:1e3; VectorSphericalWaves.N_mn_wave.(0,1,k2*r_array, θ_array,ϕ_array; kind= "regular"); end
"""


text:
"""
The LittleManComputer module used to assemble, disassamble and run code
written for the imaginary little man computer. This is a computer used
for educational purposed, to teach beginner assembly code programming.

# Example
    code = assemble("counter.lmc")
    inputs = [4, 3]
    outputs = simulate!(code, inputs)

Runs code stored in "counter.lmc" file with input `[4, 3]`.
"""
module LittleManComputer

const opcodes = Dict(
    "ADD" => 100,
    "SUB" => 200,
    "STA" => 300,
    "LDA" => 500,
    "BRA" => 600,
    "BRZ" => 700,
    "BRP" => 800,
    "INP" => 901,
    "OUT" => 902,
    "HLT" => 000
)

const mnemonics = Dict(value => key for (key, value) in opcodes)

include("assembler.jl")
include("disassembler.jl")
include("simulator.jl")

end # module


text:
# """
#     refine_scalar_potential(p::ScalarPotential{T}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T};
#         only2d::Val{only_2d} = Val(size(p.data, 2)==1)) where {T, only_2d}
#
# Refine any scalar potential `p`.
#
# 1. Extent the grid to be a closed grid in all dimensions.
# 2. Refine the axis of the grid based on `max_diffs` and `minimum_applied_potential`:
#    Insert N new ticks between to existing ticks such that the potential difference between each tick becomes
#    smaller than `max_diff[i]` (i -> dimension) but that the distances between the ticks stays larger than `minimum_distances[i]`.
# 3. Create the new data array for the refined grid and fill it by interpolation of the the initial (coarse) grid.
# """
function refine_scalar_potential(p::ScalarPotential{T}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T};
        only2d::Val{only_2d} = Val(size(p.data, 2)==1)) where {T, only_2d}
    closed_potential = _get_closed_potential(p)
    new_grid = _create_refined_grid(closed_potential, T.(max_diffs), T.(minimum_distances))
    new_data = Array{T, 3}(undef, size(new_grid))
    if only_2d
        int = interpolate_closed_potential(closed_potential, only2d)
        for i3 in axes(new_data, 3)
            x3 = new_grid.axes[3].ticks[i3]
            for i1 in axes(new_data, 1)
                x1 = new_grid.axes[1].ticks[i1]
                new_data[i1, 1, i3] = int(x1, x3)
            end
        end
    else
        int = interpolate_closed_potential(closed_potential, only2d)
        for i3 in axes(new_data, 3)
            x3 = new_grid.axes[3].ticks[i3]
            for i2 in axes(new_data, 2)
                x2 = new_grid.axes[2].ticks[i2]
                for i1 in axes(new_data, 1)
                    x1 = new_grid.axes[1].ticks[i1]
                    new_data[i1, i2, i3] = int(x1, x2, x3)
                end
            end
        end
    end
    return _convert_to_original_potential(p, new_data, new_grid)
end

function interpolate_closed_potential(p::ScalarPotential, ::Val{true}) where {T}
    interpolate((p.grid.axes[1], p.grid.axes[3]), p.data[:,1,:], Gridded(Linear()))
end
function interpolate_closed_potential(p::ScalarPotential, ::Val{false}) where {T}
    interpolate(p.grid.axes, p.data, Gridded(Linear()))
end

_get_closed_ticks(ticks::Vector{T}, int::ClosedInterval{T}) where {T} = ticks
_get_closed_ticks(ticks::Vector{T}, int::Interval{:closed, :open, T}) where {T} = vcat(ticks, int.right)
_get_closed_ticks(ticks::Vector{T}, int::Interval{:open, :closed, T}) where {T} = vcat(int.left, ticks)

_get_closed_values(values::Array{T,3}, dim::Int, int::ClosedInterval{T}) where {T} = values

function _get_closed_values(values::Array{T,3}, dim::Int, int::Interval{:closed, :open, T})::Array{T,3} where {T}
    if dim == 1
        cat(values, values[1:1,:,:], dims=dim)
    elseif dim == 2
        cat(values, values[:,1:1,:], dims=dim)
    else # dim == 3
        cat(values, values[:,:,1:1], dims=dim)
    end
end
function _get_closed_values(values::Array{T,3}, dim::Int, int::Interval{:open, :closed, T})::Array{T,3} where {T}
    if dim == 1
        cat(values[end:end,:,:], values, dims=dim)
    elseif dim == 2
        cat(values[:,end:end,:], values, dims=dim)
    else # dim == 3
        cat(values[:,:,end:end], values, dims=dim)
    end
end

function _get_closed_axis(ax::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    ticks = _get_closed_ticks(ax.ticks, ax.interval)
    return DiscreteAxis{T, BL, BR, ClosedInterval{T}}(ClosedInterval(ax.interval.left, ax.interval.right), ticks)
end

# """
#     _get_closed_potential(p::ScalarPotential{T,3,CS}) where {T, CS}
#
# Returns an closed Grid & Potential:
# E.g. if one of the axis is {:closed,:open} it will turn this into {:closed,:closed}
# and also extend the `data` field of the potential in the respective dimension and fill
# it with the respective values.
# """
function _get_closed_potential(p::ScalarPotential{T,3,CS}) where {T, CS}
    grid = p.grid
    closed_values = _get_closed_values(p.data, 1, grid.axes[1].interval)
    closed_values = _get_closed_values(closed_values, 2, grid.axes[2].interval)
    closed_values = _get_closed_values(closed_values, 3, grid.axes[3].interval)
    closed_axes = broadcast(i -> _get_closed_axis(grid.axes[i]), (1, 2, 3))
    AT = typeof(closed_axes)
    closed_grid = Grid{T, 3, CS, AT}(closed_axes)
    ScalarPotential(p, closed_values, closed_grid)
end

# """
#     _convert_to_original_potential(::Type{P}, data, grid) where {P, T, CS}
#
# Basically the counterpart to `_get_closed_potential`.
# """
function _convert_to_original_potential(p::ScalarPotential{T,3,CS,ATO}, data::Array{T, 3}, grid::Grid{T, 3, CS, AT}) where {T, CS, AT, ATO}
    ATs = get_axes_type(ATO)
    axs = broadcast(i -> _convert_closed_axis(ATs[i], grid.axes[i]), (1, 2, 3))
    new_data = _reduce_closed_values(ATs, data)
    new_grid = Grid{T,3,CS,ATO}(axs)
    ScalarPotential(p, new_data, new_grid)
end

_convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, I}}, ax::DiscreteAxis{T, BL, BR, I}) where {T, BL, BR, I} = ax

function _convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}}, ax::DiscreteAxis{T, BL, BR, ClosedInterval{T}}) where {T, BL, BR}
    DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}(Interval{:closed, :open, T}(ax.interval.left, ax.interval.right), ax.ticks[1:end-1])
end
function _convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}}, ax::DiscreteAxis{T, BL, BR, ClosedInterval{T}}) where {T, BL, BR}
    DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}(Interval{:open, :closed, T}(ax.interval.left, ax.interval.right), ax.ticks[2:end])
end

function _reduce_closed_values(axTypes::Tuple, a::Array{T, 3}) where {T}
   inds = broadcast(i -> _get_ind_range_of_closed_axis(size(a, i), axTypes[i]), (1, 2, 3))
   a[inds...]
end

_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, ClosedInterval{T}}}) where {T, BL, BR} = 1:l
_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}}) where {T, BL, BR} = 1:(l-1)
_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}}) where {T, BL, BR} = 2:l

function _create_refined_grid(p::ScalarPotential{T,3}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T}) where {T}
    max_diffs = broadcast(md -> iszero(md) ? Inf : md, max_diffs)
    n_1 = floor.(Int, [maximum(abs.(p.data[i+1,:,:] .- p.data[i,:,:])) for i in 1:size(p.data, 1)-1] ./ max_diffs[1])
    n_2 = floor.(Int, [maximum(abs.(p.data[:,i+1,:] .- p.data[:,i,:])) for i in 1:size(p.data, 2)-1] ./ max_diffs[2])
    n_3 = floor.(Int, [maximum(abs.(p.data[:,:,i+1] .- p.data[:,:,i])) for i in 1:size(p.data, 3)-1] ./ max_diffs[3])
    ns = (n_1, n_2, n_3)
    widths = diff.((p.grid.axes[1].ticks, p.grid.axes[2].ticks, p.grid.axes[3].ticks))
    sub_widths = broadcast(ia -> [widths[ia][i] / (ns[ia][i]+1) for i in eachindex(ns[ia])], (1,2,3))
    for ia in 1:3
        for i in eachindex(ns[ia])
            while sub_widths[ia][i] < minimum_distances[ia] && ns[ia][i] > 0
                ns[ia][i] -= 1
                sub_widths[ia][i] = widths[ia][i] / (ns[ia][i]+1)
            end
        end
    end
    for i in 1:3 # always add an even number of ticks
        if isodd(sum(ns[i]))
            i_max_width = findmax(sub_widths[i])[2]
            ns[i][i_max_width] += 1
            sub_widths[i][i_max_width] = widths[i][i_max_width] / (ns[i][i_max_width]+1)
        end
    end
    new_axes = broadcast(i -> _refine_axis(p.grid.axes[i], ns[i], sub_widths[i]), (1, 2, 3))
    return typeof(p.grid)(new_axes)
end

function _refine_axis(ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T}}, ns::Vector{Int}, sub_widths::Vector{T}) where {T, I}
    @assert length(ns) == length(ax.ticks)-1 # for ClosedInterval axis
    ticks = Vector{T}(undef, length(ax.ticks) + sum(ns))
    i = 1
    for j in eachindex(ns)
        ticks[i] = ax.ticks[j]
        for k in 1:ns[j]
            i += 1
            ticks[i] = ax.ticks[j] + k * sub_widths[j]
        end
        i += 1
    end
    ticks[end] = ax.ticks[end]
    typeof(ax)(ax.interval, ticks)
end

_extend_refinement_limits(rl::Real) = (rl, rl, rl )
_extend_refinement_limits(rl::Tuple{<:Real,<:Real,<:Real}) = rl


text:
export iauApio
"""
For a terrestrial observer, prepare star-independent astrometry
parameters for transformations between CIRS and observed
coordinates.  The caller supplies the Earth orientation information
and the refraction constants as well as the site coordinates.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   sp     double      the TIO locator s' (radians, Note 1)
   theta  double      Earth rotation angle (radians)
   elong  double      longitude (radians, east +ve, Note 2)
   phi    double      geodetic latitude (radians, Note 2)
   hm     double      height above ellipsoid (m, geodetic Note 2)
   xp,yp  double      polar motion coordinates (radians, Note 3)
   refa   double      refraction constant A (radians, Note 4)
   refb   double      refraction constant B (radians, Note 4)

Returned:
   astrom iauASTROM*  star-independent astrometry parameters:
      pmt    double       unchanged
      eb     double[3]    unchanged
      eh     double[3]    unchanged
      em     double       unchanged
      v      double[3]    unchanged
      bm1    double       unchanged
      bpn    double[3][3] unchanged
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Notes:

   1. sp, the TIO locator s', is a tiny quantity needed only by the
   most precise applications.  It can either be set to zero or
   predicted using the SOFA function iauSp00.

   2. The geographical coordinates are with respect to the WGS84
   reference ellipsoid.  TAKE CARE WITH THE LONGITUDE SIGN:  the
   longitude required by the present function is east-positive
   (i.e. right-handed), in accordance with geographical convention.

   3. The polar motion xp,yp can be obtained from IERS bulletins.  The
   values are the coordinates (in radians) of the Celestial
   Intermediate Pole with respect to the International Terrestrial
   Reference System (see IERS Conventions 2003), measured along the
   meridians 0 and 90 deg west respectively.  For many applications,
   xp and yp can be set to zero.

   Internally, the polar motion is stored in a form rotated onto the
   local meridian.

   4. The refraction constants refa and refb are for use in a
   dZ = A*tan(Z)+B*tan^3(Z) model, where Z is the observed
   (i.e. refracted) zenith distance and dZ is the amount of
   refraction.

   5. It is advisable to take great care with units, as even unlikely
   values of the input parameters are accepted and processed in
   accordance with the models used.

   6. In cases where the caller does not wish to provide the Earth
   rotation information and refraction constants, the function
   iauApio13 can be used instead of the present function.  This
   starts from UTC and weather readings etc. and computes suitable
   values using other SOFA functions.

   7. This is one of several functions that inserts into the astrom
   structure star-independent parameters needed for the chain of
   astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

   The various functions support different classes of observer and
   portions of the transformation chain:

         functions         observer        transformation

      iauApcg iauApcg13    geocentric      ICRS <-> GCRS
      iauApci iauApci13    terrestrial     ICRS <-> CIRS
      iauApco iauApco13    terrestrial     ICRS <-> observed
      iauApcs iauApcs13    space           ICRS <-> GCRS
      iauAper iauAper13    terrestrial     update Earth rotation
      iauApio iauApio13    terrestrial     CIRS <-> observed

   Those with names ending in "13" use contemporary SOFA models to
   compute the various ephemerides.  The others accept ephemerides
   supplied by the caller.

   The transformation from ICRS to GCRS covers space motion,
   parallax, light deflection, and aberration.  From GCRS to CIRS
   comprises frame bias and precession-nutation.  From CIRS to
   observed takes account of Earth rotation, polar motion, diurnal
   aberration and parallax (unless subsumed into the ICRS <-> GCRS
   transformation), and atmospheric refraction.

   8. The context structure astrom produced by this function is used by
   iauAtioq and iauAtoiq.

Called:
   iauPvtob     position/velocity of terrestrial station
   iauAper      astrometry parameters: update ERA

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauApio(sp::Real, theta::Real,
               elong::Real, phi::Real, hm::Real,
               xp::Real, yp::Real,
               refa::Real, refb::Real)
   # Allocate return value
   ref_astrom = Ref{iauASTROM}(iauASTROM())

   ccall((:iauApio, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble,
         Cdouble, Cdouble, Cdouble,
         Cdouble, Cdouble, Cdouble,
         Ref{iauASTROM}),
         convert(Float64, sp),
         convert(Float64, theta),
         convert(Float64, elong),
         convert(Float64, phi),
         convert(Float64, hm),
         convert(Float64, xp),
         convert(Float64, yp),
         convert(Float64, refa),
         convert(Float64, refb),
         ref_astrom)

    return ref_astrom[]
end


text:
###include("Crystal.jl")
###include("DFToutMod.jl")
#using XMLDict

"""
    module Force_Stress

Module for calculating force and stress
"""
module Force_Stress
"""
Scripts to calculate force and stress
"""

#calc for testing non-autodiff forces only
#using Calculus

using Base.Threads
using FFTW
using LinearAlgebra
using ForwardDiff
#using ReverseDiff
###using Optim
using LineSearches
using ..CrystalMod:crystal
using ..CrystalMod:makecrys

using ..CalcTB:calc_tb_fast
using ..CalcTB:distances_etc_3bdy_parallel
using ..TB:calc_energy_charge_fft
using ..TB:tb_crys
using ..TB:types_energy
using ..TB:make_kgrid
using ..TB:get_dq
using ..TB:get_h1
using ..TB:ewald_energy
using ..CalcTB:calc_frontier

using ..Ewald:electrostatics_getgamma
using ..Ewald:estimate_best_kappa
using ..SCF:scf_energy

using ..BandTools:gaussian
using ..CrystalMod:get_grid

export get_energy_force_stress
#export relax_structure

"""
    function get_energy_force_stress(crys::crystal, database; smearing = 0.01, grid = missing)

Get force and stress, non-fft algorithm. Generally use the fft algorithm.

`return energy_tot,  f_cart, stress`

Returns Ryd units. Generally users should use the `scf_energy_force_stress` function

Uses automatic differentation for gradient.
"""
function get_energy_force_stress(crys::crystal, database; smearing = 0.01, grid = missing)

    println("crys")
#    tbc = []
    tbc = calc_tb_fast(crys, database)
    return get_energy_force_stress_fft(tbc, database, do_scf=tbc.scf, grid = grid, smearing=smearing)
end

"""
    function get_energy_force_stress(crys::crystal, database; smearing = 0.01, grid = missing)

Get force and stress, non-fft algorithm
"""
function get_energy_force_stress(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing, cs = 4)

    if ismissing(grid)
        grid = get_grid(tbc.crys)
    end

    kgrid, kweights = make_kgrid(grid)
    nk = size(kgrid)[1]

    tooshort, energy_tot = safe_mode_energy(tbc.crys, database)

    if !(tooshort)
        if !ismissing(vv)
            VECTS, VALS, efermi = vv
            energy_tot = 0.0
        else

            #prepare eigenvectors / values
            error_flag = false
            if do_scf
                energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=e_den0, conv_thr = 1e-9)
            else
                energy_tot, efermi, e_den, VECTS, VALS, error_flag =  calc_energy_charge_fft(tbc, grid=grid, smearing=smearing)
            end
            if error_flag
                println("warning, trouble with eigenvectors/vals in initial step get_energy_force_stress")
            end

        end

        h1, dq = get_h1(tbc)

#        println("energy_tot $energy_tot")

        OCCS = gaussian.(VALS.-efermi, smearing)

    end

    ct = deepcopy(tbc.crys)
#=
    #reshape function
    function reshape_vec(x, nat)
        T=typeof(x[1])

        print("size x ", size(x))
        x_r = zeros(T, ct.nat, 3)
        for n = 1:ct.nat
            for j = 1:3
                x_r[n,j] = x[3*(n-1) + j]
            end
        end

        x_r_strain = zeros(T, 3,3)

        x_r_strain[1,1] = x[3*nat+1]
        x_r_strain[2,2] = x[3*nat+2]
        x_r_strain[3,3] = x[3*nat+3]

        x_r_strain[2,3] = x[3*nat+4]
        x_r_strain[3,2] = x[3*nat+4]

        x_r_strain[1,3] = x[3*nat+5]
        x_r_strain[3,1] = x[3*nat+5]

        x_r_strain[1,2] = x[3*nat+6]
        x_r_strain[2,1] = x[3*nat+6]

        return x_r, x_r_strain
    end

    #reshape function
    function inv_reshape_vec(x, strain, nat)
        T=typeof(x[1])
        x_r = zeros(T, nat*3 + 6)
        for n = 1:nat
            for j = 1:3
                x_r[3*(n-1) + j] = x[n,j]
            end
        end

        x_r[3*ct.nat + 1] = strain[1,1]
        x_r[3*ct.nat + 2] = strain[2,2]
        x_r[3*ct.nat + 3] = strain[3,3]
        x_r[3*ct.nat + 4] = strain[2,3]+strain[3,2]
        x_r[3*ct.nat + 5] = strain[1,3]+strain[3,1]
        x_r[3*ct.nat + 6] = strain[1,2]+strain[2,1]

        return x_r
    end
=#

    function f(x::Vector)
        T=typeof(x[1])

        x_r, x_r_strain = reshape_vec(x, ct.nat, strain_mode=true)

        A = ct.A * (I(3) + x_r_strain)
        #A = deepcopy(ct.A)

        crys_dual = makecrys( A , ct.coords + x_r, ct.types, units="Bohr")

        #this deals with cases where the distances between atoms become very short, which can happen during relaxations
        #currently we just have an artificial repulsive force in this case
        if tooshort
            tooshort, energy_short = safe_mode_energy(crys_dual, database, var_type=T)
            return energy_short
        end

        if database["scf"] == true
            scf = true
            kappa = estimate_best_kappa(ct.A)
            gamma_dual = electrostatics_getgamma(crys_dual, kappa=kappa)
        else
            scf = false
            gamma_dual=zeros(T, ct.nat,ct.nat)
        end

        tbc_dual = calc_tb_fast(crys_dual, database; verbose=false, var_type=T, use_threebody=true, use_threebody_onsite=true, gamma=gamma_dual, check_frontier=false)

        nwan = tbc.tb.nwan

        hk = zeros(Complex{T}, nwan, nwan)
        hk0 = zeros(Complex{T}, nwan, nwan)
        hka = zeros(Complex{T}, nwan, nwan)
        sk = zeros(Complex{T}, nwan, nwan)

        twopi_i = -1.0im*2.0*pi

        VALS0 = zeros(T, nk, nwan)

        #analytic fourier transform
        for k = 1:nk
            #vect, vals, hk, sk, vals0 = Hk(hktemp, sktemp, tbc.tb, grid[k,:])

            hk0[:,:] .= 0.0
            sk[:,:] .= 0.0

            kmat = repeat(kgrid[k,:]', tbc.tb.nr,1)
            exp_ikr = exp.(twopi_i * sum(kmat .* tbc.tb.ind_arr,dims=2))

            for m in 1:tbc.tb.nwan
                 for n in 1:tbc.tb.nwan
                     hk0[m,n] = tbc_dual.tb.H[m,n,:]'*exp_ikr[:]
                     sk[m,n]  = tbc_dual.tb.S[m,n,:]'*exp_ikr[:]
                 end
            end
            hk0 = 0.5*(hk0 + hk0')
            sk = 0.5*(sk + sk')

            if scf
                hk0 = hk0 + h1 .* sk
            end

            for a = 1:nwan

                hka[:,:] = hk0 - ( VALS[k,a] )  * sk
                #hka[:,:] = hk0 - ( VALS[k,a] )  * sk

                VALS0[k,a] += real.(VECTS[k,:,a]'*hka*VECTS[k,:,a])
            end
        end

        energy0 = sum(OCCS .* VALS0) / nk * 2.0

        if scf
            eewald, pot = ewald_energy(crys_dual, gamma_dual, dq)
        else
            eewald = 0.0
        end

        etypes = types_energy(tbc.crys)

        return energy0 + etypes + eewald

    end

#    x0 = inv_reshape_vec(ct.coords, ct.nat)

    chunksize=min(cs, 3*ct.nat + 6)
    cfg = ForwardDiff.GradientConfig(f, zeros(3*ct.nat + 6), ForwardDiff.Chunk{chunksize}())
    g = ForwardDiff.gradient(f, zeros(3*ct.nat + 6)  )

    x, stress = reshape_vec(g, ct.nat)

    f_cart = -1.0 * x
    f_cart = f_cart * inv(ct.A)'

    stress = -stress / abs(det(ct.A))

    for i = 1:3
        for j = 1:3
            if abs(stress[i,j]) < 1e-12
                stress[i,j] = 0.0
            end
        end
    end

    return energy_tot,  f_cart, stress

end

#primarily for testing
"""
    function finite_diff(crys::crystal, database, ind1, ind2; stress_mode=false, step = 0.0002, smearing = 0.01, grid = missing)

Finite differences force/stress, for testing.
# Arguments
- `crys::crystal` Crystal structure
- `database` Database of fitting coefficents.
- `ind1`  atom index for first stress index
- `ind2` cartesian index or second stress index
- `stress_mode=false` true for stress, otherwise force.
- `step = 0.0002` step_size for finite steps.
- `smearing = 0.01` smearing energy
- `grid = missing` kpoint grid
"""
function finite_diff(crys::crystal, database, ind1, ind2; stress_mode=false, step = 0.0002, smearing = 0.01, grid = missing)
    if ismissing(grid)
        grid = get_grid(crys)
    end

    tbc0 = calc_tb_fast(crys, database, verbose=false, check_frontier=false)

    energy_tot0, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc0, smearing=smearing, grid=grid, conv_thr = 1e-9)

    if stress_mode == false

        println("force mode")

        crys1 = deepcopy(crys)
        cart1 = crys1.coords * crys1.A
        cart1[ind1, ind2] += step

        crys1.coords = cart1 * inv(crys1.A)

        tbc1 = calc_tb_fast(crys1, database, verbose=false, check_frontier=false)

        energy_tot1, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc1, smearing=smearing, grid=grid, conv_thr = 1e-9)

        crys2 = deepcopy(crys)
        cart2 = crys2.coords * crys2.A
        cart2[ind1, ind2] -= step

        crys2.coords = cart2 * inv(crys2.A)

        tbc2 = calc_tb_fast(crys2, database, verbose=false, check_frontier=false)

        energy_tot2, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc2, smearing=smearing, grid=grid, conv_thr = 1e-9)

        force = - (energy_tot1 - energy_tot2) / (2 * step)

        return energy_tot0, force

    else
        println("stress mode")

        crys1 = deepcopy(crys)

        strain = zeros(3,3)
        strain[ind1,ind2] = step
        strain[ind2,ind1] = step

        crys1.A = crys1.A *(I(3) + strain)

        tbc1 = calc_tb_fast(crys1, database, verbose=false, check_frontier=false)

        energy_tot1, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc1, smearing=smearing, grid=grid, conv_thr = 1e-7)

        crys2 = deepcopy(crys)
        strain = zeros(3,3)
        strain[ind1,ind2] = -step
        strain[ind2,ind1] = -step

        crys2.A = crys2.A *(I(3) + strain)

        tbc2 = calc_tb_fast(crys2, database, verbose=false, check_frontier=false)

        energy_tot2, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc2, smearing=smearing, grid=grid, conv_thr = 1e-7)

        stress = -1.0* (energy_tot1 - energy_tot2) / (2 * step) / abs(det(crys.A))

        if ind1 != ind2
            stress = stress / 2.0
        end

        return energy_tot0, stress
    end
end

function cell_force(A, stress)

    stress = (stress + stress')/2.0

#    println("cf stress ", stress)
#    println("cf A      ", A)

    Ainv = inv(A)
    cf = zeros(3,3)
    for j= 1:3
        for i = 1:3
            cf[i,j] = sum(Ainv[j,:].*stress[i,:])
        end
    end

    return cf * abs(det(A))

end

######################################################
#=
"""
    function relax_structure(crys::crystal, database; smearing = 0.01, grid = missing, mode="vc-relax", nsteps=100, update_grid=true, conv_thr=2e-4)

Relax structure. Primary user function is relax_structure in ThreeBodyTB.jl, which calls this one.
"""
function relax_structure(crys::crystal, database; smearing = 0.01, grid = missing, mode="vc-relax", nsteps=100, update_grid=true, conv_thr = 2e-4)

    if update_grid==false
        grid = get_grid(crys)
    end

    eden = missing #starts off missing

    #do this ahead of first iteration, to get memory in correct place
    tbc = calc_tb_fast(deepcopy(crys), database)
    energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, conv_thr=1e-7)

    if error_flag
        println("warning error computing scf in relax_structure, zeroth iteration")
    end

    eden = deepcopy(tbc.eden)

    x0 = inv_reshape_vec(crys.coords, crys.A, crys.nat, strain_mode=false)
    crys_working = deepcopy(crys)

    fcall = 0
    firstiter = true

    nat = crys.nat

    function fn(x)
        println("CALL FN", x)
        coords, A = reshape_vec(x, nat)

#        println("coords ", coords)
#        println("A ", A)

        crys_working.coords = coords
        if  mode == "vc-relax"
            crys_working.A = A
        end
        println("FN crys")
        println(crys_working)

        tooshort, energy_short = safe_mode_energy(crys_working, database)
#        println("fn too short ", tooshort)
        if tooshort
            return energy_short
        end

        if firstiter
            firstiter = false
        else
            if crys_working != tbc.crys
#                println("yes calc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
#                println(crys_working)
                tbc = calc_tb_fast(deepcopy(crys_working), database, verbose=false)
            else
#                println("nocalc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
#                println(crys_working)
#                println(tbc.crys)

            end
#            println(crys_working)
            energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, verbose=false, conv_thr=1e-6)
            eden = deepcopy(tbcx.eden)
        end
#        println("fn $energy_tot fnffnfnfnffnfnfnfnfnfnfnfnfnfnfffffff")
        return energy_tot

    end

    f_cart_global = []
    stress_global = []

    function grad(storage, x)
        println("CALL GRAD", x)

#        println("typeof x ", typeof(x))
        fcall += 1

        coords, A = reshape_vec(x, nat)
        crys_working.coords = coords

        if  mode == "vc-relax"
            crys_working.A = A
        end

        println("GRAD crys")
        println(crys_working)

#        println("crys_working $fcall")
#        println(crys_working)

        tooshort, energy_short = safe_mode_energy(crys_working, database)

#        println("too short ", tooshort)

        if crys_working != tbc.crys && !tooshort
#            println("yes calc forces -----------------------------------------------------------------------------------------")
#            println(crys_working)
            tbc = calc_tb_fast(deepcopy(crys_working), database, verbose=false)
            energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, verbose=false, conv_thr=1e-6)

            eden = deepcopy(tbcx.eden)

            #        else
#            println("no calc forces ------------------------------------------------------------------------------------------")
#            println(crys_working)
#            println(tbc.crys)
#            println("__")
        end

        energy_tmp,  f_cart, stress =  get_energy_force_stress_fft(tbc, database; do_scf = false, smearing = smearing, grid = grid, vv=[VECTS, VALS, efermi] )
        #energy_tmp,  f_cart, stress =  get_energy_force_stress_fft(tbc, database; do_scf = true, smearing = smearing, grid = grid )

        f_cart_global = f_cart
        stress_global = stress

        fsum = sum((f_cart).^2)^0.5
        ssum= sum((stress).^2)^0.5

#        println("f_cart")
#        println(f_cart)

        println("FCALL $fcall en:  $energy_tot (Ryd)  fsum:  $fsum  ssum:  $ssum    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")

#        println("grad stress: ", stress)

        stress_units = cell_force( crys_working.A, stress)
        #        stress_units = (stress * abs(det(crys_working.A)) *  inv(crys_working.A))'
        #stress_units = (stress * abs(det(crys_working.A)) )'

#        println("grad stress units: ", -stress_units)

        #        f_crys = f_cart * inv(crys_working.A)
        f_crys = f_cart * crys_working.A

        g = inv_reshape_vec(-f_crys, -stress_units, crys_working.nat, strain_mode=false)

#        println("g ", g)

        if mode != "vc-relax" #don't need stress
            g[3*nat+1:end] = zeros(9)
        end
        storage[:] = g
#        println("ret storage $fcall")
#        println(storage)
    end

#=
    function fn_fake(x)
        println("fn_fake ", x[1]*2)
        return sum(x.^2)
    end

    function grad_fake(storage, x)
        storage[:] =  2.0 * x

        println("grad_fake  ", x[1]*2)
        println("grad_fakex ", x)
        println("grad_fakes ", storage)
        return storage

    end
=#

#    opts = Optim.Options(g_tol = 1e-3,f_tol = 1e-3, x_tol = 1e-3,
#                         iterations = 10,
#                         store_trace = true,
#                         show_trace = false)

#    println("X000000000000000000, ", x0)
#    println()
#    println(fn(x0))
#    println("X000000000000000001, ", x0)
#    storage = zeros(size(x0))
#    println(grad(x0, storage))

#    println("storage")
#    println(storage)

    println("starting vec ")
    println(x0)
    println()

#    opts = Optim.Options(g_tol = 7e-4,f_tol = 7e-4, x_tol = 7e-4,
#                         iterations = nsteps,
#                         store_trace = true,
#                         show_trace = false)

    opts = Optim.Options(g_tol = conv_thr,f_tol = conv_thr, x_tol = conv_thr,
                             iterations = nsteps,
                             store_trace = true,
                             show_trace = false)

    #res = optimize(fn,grad, x0, ConjugateGradient(), opts)

    #preconditioner, guess for inv Hess is based on the  metric. see qe bfgs_module.f90
    function init(x)
        num = 3*nat + 9

#        factor=0.1
        factor = 1.0

        P = zeros(eltype(x), num, num)

        A = crys.A
        g = A'*A
        ginv = inv(g)

        for a = 1:nat
            aa = (a-1)*3
            for i = 1:3
                for j = 1:3
                    P[aa+i ,aa+j] = g[i,j] * factor
                end
            end
        end

        vol = abs(det(A))

#        if mode == "vc-relax"
        for b = 1:3
            bb = 3 * nat + (b - 1)*3
            for i = 1:3
                for j = 1:3
                    P[bb+i,bb+j] = 0.04 * vol * ginv[i,j] * factor
                end
            end
        end
#    end

        return inv(P)
    end

#    P = init(x0)
#    for i = 1:size(P)[1]
#        for j = 1:size(P)[2]
#            println("P $i $j ", P[i,j])
#        end
#    end
#    println("P")
#    println(P)

    #x -> Matrix{eltype(x)}(I, length(x), length(x))

#    res = optimize(fn,x0, NelderMead(), opts)

    #    res = optimize(fn,grad, x0, ConjugateGradient())

#    res = optimize(fn,grad, x0, BFGS( ), opts)

    #res = optimize(fn,grad, x0, BFGS(  initial_invH = init ), opts)

    res = missing

    try
        res = optimize(fn,grad, x0, BFGS(  initial_invH = init, linesearch=LineSearches.MoreThuente() ), opts)

        #res = optimize(fn,grad, x0, BFGS(initial_invH = init ), opts)

    catch e
        if e isa InterruptException
            println("user interrupt")

            return tbc.crys
        else
            println("unknown error")
            return tbc.crys

        end
    end

    #res = optimize(fn,grad, x0, BFGS(  initial_invH = init, linesearch=LineSearches.HagerZhang() ), opts)

    energy = res.minimum

    # bad : res = optimize(fn,grad, x0, ConjugateGradient(linesearch=LineSearches.BackTracking(order=3) ), opts)
    #linesearch=LineSearches.BackTracking(order=2),
    #LineSearches.MoreThuente()
    #LineSearches.HagerZhang()
    #linesearch=LineSearches.BackTracking(order=2)
    #linesearch=LineSearches.StrongWolfe()

#for testing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#    storage = zeros(size(x0))
#    g  = grad(storage, x0)

#    println()

#    res = Calculus.gradient(fn, x0)

#    println("grad x0")
#    println(storage)
#    println("calc")
#    println(res)

#    println([storage res])
#
#    return [storage res]
#end for testing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#    res = optimize(fn, x0, LBFGS())
#
#    res = optimize(fn,grad, x0, LBFGS())

#    res = optimize(fn, x0, ConjugateGradient())

#    res = optimize(fn,grad, x0, LBFGS(), opts)

#    res = optimize(fn, x0, LBFGS(), opts)

#    f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
#    function g!(storage, x)
 #       println("g $x");
 #       storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
 #       storage[2] = 200.0 * (x[2] - x[1]^2)
 #       println("g storage", storage)
 #   end

#    res = optimize(f, g!,  [0.0, 0.0], LBFGS())

    println()
    println("res")
    println(res)

    minvec = Optim.minimizer(res)
#    println("minvec")
#    println(minvec)

    coords, A = reshape_vec(minvec, nat)
    cfinal = makecrys(A, coords, crys.types, units="Bohr")
    return cfinal, tbc, energy, f_cart_global, stress_global

#    return res

end
=#

"""
    function reshape_vec(x, nat; strain_mode=false)

The force and relax algorithms from outside codes take in vectors, not `crystal` 's. So we have to reshape
vectors to and from `crystals`
"""
function reshape_vec(x, nat; strain_mode=false)
#    println("RESHAPEVEC RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRr ", strain_mode)

    T=typeof(x[1])

#    println("size x ", size(x))
    x_r = zeros(T, nat, 3)
    for n = 1:nat
        for j = 1:3
            x_r[n,j] = x[3*(n-1) + j]
        end
    end

    x_r_strain = zeros(T, 3,3)

    if length(x) == 3*nat+6 && strain_mode
        x_r_strain[1,1] = x[3*nat+1]
        x_r_strain[2,2] = x[3*nat+2]
        x_r_strain[3,3] = x[3*nat+3]

        x_r_strain[2,3] = 0.5*x[3*nat+4]
        x_r_strain[3,2] = 0.5*x[3*nat+4]

        x_r_strain[1,3] = 0.5*x[3*nat+5]
        x_r_strain[3,1] = 0.5*x[3*nat+5]

        x_r_strain[1,2] = 0.5*x[3*nat+6]
        x_r_strain[2,1] = 0.5*x[3*nat+6]
    elseif length(x) == 3*nat+6
        x_r_strain[1,1] = x[3*nat+1]
        x_r_strain[2,2] = x[3*nat+2]
        x_r_strain[3,3] = x[3*nat+3]

        x_r_strain[2,3] = x[3*nat+4]
        x_r_strain[3,2] = x[3*nat+4]

        x_r_strain[1,3] = x[3*nat+5]
        x_r_strain[3,1] = x[3*nat+5]

        x_r_strain[1,2] = x[3*nat+6]
        x_r_strain[2,1] = x[3*nat+6]

    elseif length(x) == 3*nat+9
        x_r_strain[1,1] = x[3*nat+1]
        x_r_strain[1,2] = x[3*nat+2]
        x_r_strain[1,3] = x[3*nat+3]
        x_r_strain[2,1] = x[3*nat+4]
        x_r_strain[2,2] = x[3*nat+5]
        x_r_strain[2,3] = x[3*nat+6]
        x_r_strain[3,1] = x[3*nat+7]
        x_r_strain[3,2] = x[3*nat+8]
        x_r_strain[3,3] = x[3*nat+9]
    else
        println("I'm confusing about the length reshape_vec $nat ", length(x) )
    end

    return x_r, x_r_strain
end

"""
    function inv_reshape_vec(x, strain, nat; strain_mode=true)

The force and relax algorithms from outside codes take in vectors, not `crystal` 's. So we have to reshape
vectors to and from `crystals`
"""
function inv_reshape_vec(x, strain, nat; strain_mode=true)
    T=typeof(x[1])
    if strain_mode
        x_r = zeros(T, nat*3 + 6)
    else
        x_r = zeros(T, nat*3 + 9)
    end

    for n = 1:nat
        for j = 1:3
            x_r[3*(n-1) + j] = x[n,j]
        end
    end
    if strain_mode
        x_r[3*nat + 1] = strain[1,1]
        x_r[3*nat + 2] = strain[2,2]
        x_r[3*nat + 3] = strain[3,3]
        x_r[3*nat + 4] = (strain[2,3]+strain[3,2])
        x_r[3*nat + 5] = (strain[1,3]+strain[3,1])
        x_r[3*nat + 6] = (strain[1,2]+strain[2,1])
    else
        x_r[3*nat+1] = strain[1,1]
        x_r[3*nat+2] = strain[1,2]
        x_r[3*nat+3] = strain[1,3]
        x_r[3*nat+4] = strain[2,1]
        x_r[3*nat+5] = strain[2,2]
        x_r[3*nat+6] = strain[2,3]
        x_r[3*nat+7] = strain[3,1]
        x_r[3*nat+8] = strain[3,2]
        x_r[3*nat+9] = strain[3,3]
    end

    return x_r
end

"""
    function safe_mode_energy(crys::crystal, database; var_type=Float64)

Relaxation can accidently lead to very small atom-atom distances during the relaxation precedure if too large of step is taken. This function is a repulsive energy
function at short range to make sure the relaxtion doesn't get stuck at very short distances
where the fitting doesn't apply.
"""
function safe_mode_energy(crys::crystal, database; var_type=Float64, check=true)

    diststuff = distances_etc_3bdy_parallel(crys,10.0, 0.0, var_type=var_type)
    R_keep, R_keep_ab, array_ind3, array_floats3, dist_arr, c_zero, dmin_types, dmin_types3 = diststuff
    nkeep = size(R_keep_ab)[1]

    energy = 1.0
    tooshort = false

    warned = zeros(Bool, crys.nat)

    for a1 = 1:crys.nat
        t1 = crys.stypes[a1]
        for a2 = 1:crys.nat
            t2 = crys.stypes[a2]
            dmin = database[(t1,t2)].min_dist * 0.979
            for c = 1:nkeep
                cind = R_keep_ab[c,1]
                if dist_arr[a1,a2,cind,1] < dmin*2 && dist_arr[a1,a2,cind,1] > 1e-7
                    energy += 1e-3/dist_arr[a1,a2,cind,1]
                end

                if dist_arr[a1,a2,cind,1] < dmin*1.01999 && dist_arr[a1,a2,cind,1] > 1e-7
                    tooshort = true
                    energy += 0.02 * (dist_arr[a1,a2,cind,1] - dmin)^2 + 0.3 * abs(dist_arr[a1,a2,cind,1] - dmin)
                    if var_type == Float64 && warned[a1] == false
                        println("WARNING, SAFE MODE $a1 $t1 $a2 $t2 $c ", dist_arr[a1,a2,cind,1])
                        warned[a1] = true
                    end
                end
            end
        end
    end
    if tooshort
        println("WARNING, safe mode activated, minimum distances < fitting data * 0.98")
        return tooshort, energy
    elseif check==true
        #violation_list, vio_bool = calc_frontier(crys, database, test_frontier=true, diststuff=diststuff, verbose=false)
        violation_list, vio_bool = calc_frontier(crys, database, test_frontier=true, verbose=false)
#        println("vio_vool $vio_bool")
        if vio_bool == false
            println("ACTIVATE SAFE MODE")
            return true, energy
        else
            return false, energy
        end
    else
        return tooshort, energy
    end

end

##############################################################################################################
"""
    function get_energy_force_stress_fft(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing)

Calculate energy/force/stress using fft algorithm. Users should use `scf_energy_force_stress`, which calls this. Uses automatic differentation for jacobian.
"""
function get_energy_force_stress_fft(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing)

    println("get_energy_force_stress_fft")

    if ismissing(grid)
        grid = get_grid(tbc.crys)
    end

    kgrid, kweights = make_kgrid(grid)
    nk = size(kgrid)[1]

    ct = deepcopy(tbc.crys)

#    println("test safe get_energy_force_stress_fft")
    tooshort, energy_tot = safe_mode_energy(tbc.crys, database)

    if tooshort ##########################
        println("safemode")
        function f(x :: Vector)
            T=typeof(x[1])
            x_r, x_r_strain = reshape_vec(x, ct.nat, strain_mode=true)
            A = ct.A * (I(3) + x_r_strain)
            crys_dual = makecrys( A , ct.coords + x_r, ct.types, units="Bohr")

            tooshort, energy_short = safe_mode_energy(crys_dual, database, var_type=T, check=false)
            return energy_short
        end
        garr = ForwardDiff.gradient(f, zeros(3*ct.nat + 6) )

    else #not too short ##################

        scf = database["scf"]
        #    println("scf ", scf)
        if !(tooshort)
            if !ismissing(vv)
                VECTS, VALS, efermi = vv
                energy_tot = 0.0
            else
                #prepare eigenvectors / values
                error_flag = false
                if do_scf
                    println("doing scf aaaaaaaaaaa")
                    energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=e_den0, conv_thr = 1e-8)
                else
                    energy_tot, efermi, e_den, VECTS, VALS, error_flag =  calc_energy_charge_fft(tbc, grid=grid, smearing=smearing)
                end
                if error_flag
                    println("warning, trouble with eigenvectors/vals in initial step get_energy_force_stress")
                end
            end
            h1, dq = get_h1(tbc)
            OCCS = gaussian.(VALS.-efermi, smearing)
        end

        size_ret = tbc.tb.nwan * tbc.tb.nwan * tbc.tb.nr

        function ham(x :: Vector)
            T=typeof(x[1])

            x_r, x_r_strain = reshape_vec(x, ct.nat, strain_mode=true)

            #x_r, x_r_strain = reshape_vec(x, 0, strain_mode=true)

            A = ct.A * (I(3) + x_r_strain)
            crys_dual = makecrys( A , ct.coords + x_r, ct.types, units="Bohr")

            #crys_dual = makecrys( A , ct.coords , ct.types)

            #        gamma_dual=zeros(T, ct.nat,ct.nat)
            #println("gamma")
            if database["scf"] == true
                scf = true
                kappa = estimate_best_kappa(ct.A)
                gamma_dual = electrostatics_getgamma(crys_dual, kappa=kappa)
            else
                scf = false
                gamma_dual=zeros(T, ct.nat,ct.nat)
            end

            tbc_dual = calc_tb_fast(crys_dual, database; verbose=false, var_type=T, use_threebody=true, use_threebody_onsite=true, gamma=gamma_dual , check_frontier=false)
            ret = zeros(T, size_ret * 2 + 1)

            ret[1:size_ret] = real.(tbc_dual.tb.H[:])
            ret[size_ret+1:size_ret*2] = real.(tbc_dual.tb.S[:])

            #we stick eewald in the last entry of ret
            if scf
                eewald, pot = ewald_energy(crys_dual, gamma_dual, dq)
                ret[end] = eewald
            end

            return ret
        end

#        println("jac")
        @time begin

            ret = ham(zeros(3*ct.nat + 6))

            #ret = ham(zeros( 6))

            #chunksize=min(12, 3*ct.nat + 6)
            chunksize=min(9, 3*ct.nat + 6)
            #chunksize=min(6, 3*ct.nat + 6)

            cfg = ForwardDiff.JacobianConfig(ham, zeros(3*ct.nat + 6), ForwardDiff.Chunk{chunksize}())
            #cfg = ForwardDiff.JacobianConfig(ham, zeros( 6), ForwardDiff.Chunk{chunksize}())
            g = ForwardDiff.jacobian(ham, zeros(3*ct.nat + 6) , cfg ) ::  Array{Float64,2}
#            g = ForwardDiff.jacobian(ham, zeros( 6) , cfg ) ::  Array{Float64,2}

        end
#        println("end jac")
        #    function f_es(x::Vector)
        #        x_r, x_r_strain = reshape_vec(x, ct.nat, strain_mode=true)
        #        A = ct.A * (I(3) + x_r_strain)
        #
        #        crys_dual = makecrys( A , ct.coords + x_r, ct.types)
        #
        #
        #
        #        if scf
        #            kappa = estimate_best_kappa(ct.A)
        #            gamma_dual = electrostatics_getgamma(crys_dual, kappa=kappa)
        #            eewald, pot = ewald_energy(crys_dual, gamma_dual, dq)
        #        else
        #            eewald = 0.0
        #        end
        #        return eewald
        #    end

        #    cfg_grad = ForwardDiff.GradientConfig(f_es, zeros(3*ct.nat + 6), ForwardDiff.Chunk{3*ct.nat + 6}())
        #    g_es = ForwardDiff.gradient(f_es, zeros(3*ct.nat + 6) , cfg_grad )

        #    @time gr = ReverseDiff.jacobian(ham, zeros(3*ct.nat + 6)  )

        #    return 0,0,0

#        ham_orig = ham(zeros(3*ct.nat + 6))

        #    println("size of g ", size(g))
        #    println("eltype ", eltype(g))

        hr_g = zeros(Complex{eltype(g)},  tbc.tb.nwan, tbc.tb.nwan,  grid[1], grid[2], grid[3],3*ct.nat + 6 )
        sr_g = zeros(Complex{eltype(g)},  tbc.tb.nwan, tbc.tb.nwan,  grid[1], grid[2], grid[3],3*ct.nat + 6 )

        #    hr = zeros(Complex{eltype(g)}, tbc.tb.nwan, tbc.tb.nwan,  grid[1], grid[2], grid[3])

        ind = zeros(Int64, 3)
        new_ind = zeros(Int64, 3)

        for c in 1:size(tbc.tb.ind_arr)[1]

            ind[:] = tbc.tb.ind_arr[c,:]

            new_ind[1] = mod(ind[1], grid[1])+1
            new_ind[2] = mod(ind[2], grid[2])+1
            new_ind[3] = mod(ind[3], grid[3])+1

            for na = 1:tbc.tb.nwan
                for nb = 1:tbc.tb.nwan
                    hr_g[na,nb,new_ind[1], new_ind[2], new_ind[3], :] += g[na + (nb-1) * tbc.tb.nwan + tbc.tb.nwan^2 * (c-1) , :]

                    #                hr[na,nb,new_ind[1], new_ind[2], new_ind[3]] += ham_orig[na + (nb-1) * tbc.tb.nwan + tbc.tb.nwan^2 * (c-1) ]

                    sr_g[na,nb,new_ind[1], new_ind[2], new_ind[3], :] += g[size_ret + na + (nb-1) * tbc.tb.nwan + tbc.tb.nwan^2 * (c-1),: ]

                end
            end
        end

        begin
            hk_g = fft(hr_g, [3,4,5])
            sk_g = fft(sr_g, [3,4,5])

            VALS0 = zeros(Float64,  prod(grid), tbc.tb.nwan,3*ct.nat+6)

            #        hka = zeros(Complex{Float64}, tbc.tb.nwan, tbc.tb.nwan)
            #        hka = zeros(Complex{Float64}, tbc.tb.nwan, tbc.tb.nwan, 3*ct.nat+6)

            # hk = fft(hr, [3,4,5])
            #    VALS00 = zeros(Float64, prod(grid), tbc.tb.nwan)
        end

        #    if false
        #
        #        pVECTS = permutedims(VECTS, [2,3,1])
        #
        #
        #        hk_g2 = zeros(Complex{Float64}, prod(grid), tbc.tb.nwan, tbc.tb.nwan,  3*ct.nat+6)
        #        sk_g2 = zeros(Complex{Float64}, prod(grid), tbc.tb.nwan, tbc.tb.nwan,  3*ct.nat+6)
        #
        #
        #        for gc = 1:3*ct.nat+6
        #            c=0
        #            for k1 = 1:grid[1]
        #                for k2 = 1:grid[2]
        #                    for k3 = 1:grid[3]
        #                        c += 1
        #                        hk_g2[c, :,:,gc] = hk_g[:,:,k1,k2,k3,gc]
        #                        sk_g2[c, :,:,gc] = sk_g[:,:,k1,k2,k3,gc]
        #                    end
        #                end
        #            end
        #        end
        #
        #        hka = zeros(Complex{Float64}, prod(grid), tbc.tb.nwan, tbc.tb.nwan,  tbc.tb.nwan,3*ct.nat+6)
        #        for gc = 1:3*ct.nat+6
        #            for a = 1:tbc.tb.nwan
        #                for a2 = 1:tbc.tb.nwan
        #                    for a1 = 1:tbc.tb.nwan
        #                        hka[:,a1,a2,a, gc] = hk_g2[:,a1,a2,gc] - VALS[:,a] .*  sk_g2[:,a1,a2,gc]
        #                    end
        #                end
        #            end
        #        end
        #    end
        #    hka[:,:] =  hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]

        #    if scf
        #        hk0 = hk0 + h1 .* sk
        #    end

#        println("old vals")
        if false
            pVECTS = permutedims(VECTS, [2,3,1])

            if scf
                @threads for a1 = 1:tbc.tb.nwan
                    for a2 = 1:tbc.tb.nwan
                        hk_g[a2,a1,:,:,:,:] += h1[a2,a1]*sk_g[a2,a1,:,:,:,:]
                    end
                end
            end
            @threads for gc = 1:3*ct.nat+6
                c=0
                for k1 = 1:grid[1]
                    for k2 = 1:grid[2]
                        for k3 = 1:grid[3]
                            c += 1
                            for a =1:tbc.tb.nwan
                                #                        hka[:,:] =  hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]

                                VALS0[c, a,gc] = real(pVECTS[:,a,c]' * ( hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]   )  * pVECTS[:,a,c])                  #+ h1 .*  sk_g[:,:,k1,k2,k3, gc]
                            end
                        end
                    end
                end
            end
        end

#        println("new vals")

        psi_gradH_psi(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, tbc.tb.nwan, ct.nat, grid, OCCS)

        #        @time psi_gradH_psi(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, tbc.tb.nwan, ct.nat, grid)

        #    @time for gc = 1:3*ct.nat+6
        #        c=0
        #        for k1 = 1:grid[1]
        #            for k2 = 1:grid[2]
        #                for k3 = 1:grid[3]
        #                    c += 1
        #                    for a =1:tbc.tb.nwan
        #                        hka[:,:] =  hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]
        #                        VALS0[c, a,gc] = real(pVECTS[:,a,c]' * hka[c,:,:,a, gc]  * pVECTS[:,a,c])
        #                    end
        #                end
        #            end
        #        end
        #    end

        #end
        #    end

        garr =zeros(Float64, 3*ct.nat+6)
        for gc = 1:3*ct.nat+6
            garr[gc] = sum(OCCS .* VALS0[:,:,gc]) / nk * 2.0
        end

        if scf
            garr += g[end,:][:]
        end

    end

    x, stress = reshape_vec(garr, ct.nat)
    f_cart = -1.0 * x
    f_cart = f_cart * inv(ct.A)'
    stress = -stress / abs(det(ct.A))

    for i = 1:3
        for j = 1:3
            if abs(stress[i,j]) < 1e-12
                stress[i,j] = 0.0
            end
        end
    end
    for i = ct.nat
        for j = 1:3
            if abs(f_cart[i,j]) < 1e-12
                f_cart[i,j] = 0.0
            end
        end
    end

    return energy_tot,  f_cart, stress

end

#############################################################################################################

"""
    function psi_gradH_psi(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, nwan, nat, grid, OCCS)

Helper function for <psi | grad_Ham | psi >
"""
function psi_gradH_psi(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, nwan, nat, grid, OCCS)

    pVECTS = permutedims(VECTS, [2,3,1])

    if scf
        @threads for a1 = 1:nwan
            for a2 = 1:nwan
                @inbounds hk_g[a2,a1,:,:,:,:] +=  h1[a2,a1]*(@view sk_g[a2,a1,:,:,:,:])
            end
        end
    end
    occs_max = ones(Int64, prod(grid))*nwan

    @threads for c = 1:prod(grid)
        for a = 1:nwan
            if OCCS[c,a] < 1e-4
                occs_max[c] = max(a-1, 1)
                break
            end
        end
    end

    @threads for gc = 1:3*nat+6
        c=0
        for k1 = 1:grid[1]
            for k2 = 1:grid[2]
                for k3 = 1:grid[3]
                    c += 1
                    for a =1:occs_max[c]
                        #                        hka[:,:] =  hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]

                        #VALS0[c, a,gc] = real(pVECTS[:,a,c]' * ( hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]   )  * pVECTS[:,a,c])                  #+ h1 .*  sk_g[:,:,k1,k2,k3, gc]
                        @inbounds VALS0[c, a,gc] = real( (@view pVECTS[:,a,c])' * ( (@view hk_g[:,:,k1,k2,k3,gc]) - ( VALS[c,a] )  *   (@view sk_g[:,:,k1,k2,k3, gc])   )  * (@view pVECTS[:,a,c]))                  #+ h1 .*  sk_g[:,:,k1,k2,k3, gc]
                    end
                end
            end
        end
    end

end

"""
    function psi_gradH_psi(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, nwan, nat, grid, OCCS)

Helper function for <psi | grad_Ham | psi >

This version isn't used.
"""
function psi_gradH_psi2(VALS0, VECTS, hk_g, sk_g, h1, VALS, scf, nwan, nat, grid, OCCS)

    pVECTS = permutedims(VECTS, [2,3,1])

    if scf
        @threads for a1 = 1:nwan
            for a2 = 1:nwan
                hk_g[a2,a1,:,:,:,:] += h1[a2,a1]*(@view sk_g[a2,a1,:,:,:,:])
            end
        end
    end
    occs_max = ones(Int64, prod(grid))*nwan

    @threads for c = 1:prod(grid)
        for a = 1:nwan
            if OCCS[c,a] < 1e-4
                occs_max[c] = max(a-1, 1)
                break
            end
        end
    end

    @threads for gc = 1:3*nat+6
        c=0
        for k1 = 1:grid[1]
            for k2 = 1:grid[2]
                for k3 = 1:grid[3]
                    c += 1
                    for a =1:occs_max[c]
                        #                        hka[:,:] =  hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]

                        #                        VALS0[c, a,gc] = real(pVECTS[:,a,c]' * ( hk_g[:,:,k1,k2,k3,gc] - ( VALS[c,a] )  *   sk_g[:,:,k1,k2,k3, gc]   )  * pVECTS[:,a,c])                  #+ h1 .*  sk_g[:,:,k1,k2,k3, gc]

                        VALS0[c, a,gc] = real( (@view pVECTS[:,a,c])' * ( (@view hk_g[:,:,k1,k2,k3,gc]) - ( VALS[c,a] )  *   (@view sk_g[:,:,k1,k2,k3, gc])   )  * (@view pVECTS[:,a,c]))                  #+ h1 .*  sk_g[:,:,k1,k2,k3, gc]
                    end
                end
            end
        end
    end

end

end #end module


text:
export FunctionActor
export subscribe!

"""
    FunctionActor{D} <: Actor{D}

FunctionActor provides a simple interface to use a single function as a `next!` callback, `error!` callback throws an `ErrorException` and `complete!` callback does nothing.
Should not be used explicitly because it will be created automatically when passing a `Function` object as an actor in `subscribe!` function.

# Examples
```jldoctest
using Rocket

source = from(1:5)
subscribe!(source, (t) -> println(t))
;

# output
1
2
3
4
5
```

See also: [`Actor`](@ref), [`subscribe!`](@ref)
"""
struct FunctionActor{D, F} <: Actor{D}
    on_next :: F
end

on_next!(actor::FunctionActor{D}, data::D) where D = actor.on_next(data)
on_error!(actor::FunctionActor, err)               = error(err)
on_complete!(actor::FunctionActor)                 = begin end

struct FunctionActorFactory{F} <: AbstractActorFactory
    on_next :: F
end

create_actor(::Type{L}, factory::FunctionActorFactory{F}) where { L, F } = FunctionActor{L, F}(factory.on_next)

subscribe!(source, fn::F) where { F <: Function } = subscribe!(source, FunctionActorFactory{F}(fn))


text:
using PoincareInvariants
using GeometricIntegrators
using GeometricIntegrators.Utils
using Test

include("poincare_invariant_1st_tests.jl")
include("poincare_invariant_2nd_unittests.jl")
include("poincare_invariant_2nd_tests.jl")


text:
#=
Solving the optimal growth problem via value function iteration.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-07-05

References
----------

Simple port of the file quantecon.models.optgrowth

http://quant-econ.net/jl/dp_intro.html
=#

#=
    This type defines the primitives representing the growth model. The
    default values are

        f(k) = k**alpha, i.e, Cobb-Douglas production function
        u(c) = ln(c), i.e, log utility

    See the constructor below for details
=#
"""
Neoclassical growth model

##### Fields

- `f::Function` : Production function
- `bet::Real` : Discount factor in (0, 1)
- `u::Function` : Utility function
- `grid_max::Int` : Maximum for grid over savings values
- `grid_size::Int` : Number of points in grid for savings values
- `grid::FloatRange` : The grid for savings values

"""
type GrowthModel
    f::Function
    bet::Real
    u::Function
    grid_max::Int
    grid_size::Int
    grid::FloatRange
end

default_f(k) = k^0.65
default_u(c) = log(c)

"""
Constructor of `GrowthModel`

##### Arguments

- `f::Function(k->k^0.65)` : Production function
- `bet::Real(0.95)` : Discount factor in (0, 1)
- `u::Function(log)` : Utility function
- `grid_max::Int(2)` : Maximum for grid over savings values
- `grid_size::Int(150)` : Number of points in grid for savings values

"""
function GrowthModel(f=default_f, bet=0.95, u=default_u, grid_max=2,
                     grid_size=150)
    grid = linspace_range(1e-6, grid_max, grid_size)
    return GrowthModel(f, bet, u, grid_max, grid_size, grid)
end

"""
$(____bellman_main_docstring).

##### Arguments

- `g::GrowthModel` : Instance of `GrowthModel`
- `w::Vector`: Current guess for the value function
- `out::Vector` : Storage for output.
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function bellman_operator!(g::GrowthModel, w::Vector, out::Vector;
                           ret_policy::Bool=false)
    # Apply linear interpolation to w
    Aw = CoordInterpGrid(g.grid, w, BCnan, InterpLinear)

    for (i, k) in enumerate(g.grid)
        objective(c) = - g.u(c) - g.bet * Aw[g.f(k) - c]
        res = optimize(objective, 1e-6, g.f(k))
        c_star = res.minimum

        if ret_policy
            # set the policy equal to the optimal c
            out[i] = c_star
        else
            # set Tw[i] equal to max_c { u(c) + beta w(f(k_i) - c)}
            out[i] = - objective(c_star)
        end
    end

    return out
end

function bellman_operator(g::GrowthModel, w::Vector;
                          ret_policy::Bool=false)
    out = similar(w)
    bellman_operator!(g, w, out, ret_policy=ret_policy)
end

"""
$(____greedy_main_docstring).

##### Arguments

- `g::GrowthModel` : Instance of `GrowthModel`
- `w::Vector`: Current guess for the value function
- `out::Vector` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
function get_greedy!(g::GrowthModel, w::Vector, out::Vector)
    bellman_operator!(g, w, out, ret_policy=true)
end

get_greedy(g::GrowthModel, w::Vector) = bellman_operator(g, w, ret_policy=true)

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
module StaticBumperFrostHelper

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/StaticBumper" StaticBumper(x::Integer, y::Integer, respawnTime::Number=0.6, moveTime::Number=1.81818187, sprite::String="bumper", wobble::Bool=false, notCoreMode::Bool=false, easing::String="CubeInOut")

const easings = String[
    "BackIn",
    "BackInOut",
    "BackOut",
    "BounceIn",
    "BounceInOut",
    "BounceOut",
	"CubeIn",
    "CubeInOut",
	"CubeOut",
	"ElasticIn",
    "ElasticInOut",
	"ElasticOut",
	"ExpoIn",
    "ExpoInOut",
	"ExpoOut",
    "Linear",
	"QuadIn",
    "QuadInOut",
    "QuadOut",
    "SineIn",
    "SineInOut",
    "SineOut"
]

const placements = Ahorn.PlacementDict(
    "Custom Bumper (Frost Helper)" => Ahorn.EntityPlacement(
        StaticBumper
    )
)

Ahorn.editingOptions(entity::StaticBumper) = Dict{String, Any}(
    "easing" => easings
)

Ahorn.nodeLimits(entity::StaticBumper) = 0, 1

sprite = "objects/Bumper/Idle22.png"

function Ahorn.selection(entity::StaticBumper)
    x, y = Ahorn.position(entity)
	nodes = get(entity.data, "nodes", ())

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        return [Ahorn.getSpriteRectangle(sprite, x, y), Ahorn.getSpriteRectangle(sprite, nx, ny)]
    end

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::StaticBumper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        theta = atan(y - ny, x - nx)
        Ahorn.drawArrow(ctx, x, y, nx + cos(theta) * 8, ny + sin(theta) * 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::StaticBumper, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
# Autogenerated wrapper script for LIBMINCd_jll for x86_64-linux-gnu
export libminc2, libminc2_simple

using HDF5_jll
using NetCDF_jll
JLLWrappers.@generate_wrapper_header("LIBMINCd")
JLLWrappers.@declare_library_product(libminc2, "libminc2.so.5.3.0")
JLLWrappers.@declare_library_product(libminc2_simple, "libminc2-simple.so")
function __init__()
    JLLWrappers.@generate_init_header(HDF5_jll, NetCDF_jll)
    JLLWrappers.@init_library_product(
        libminc2,
        "lib/libminc2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libminc2_simple,
        "lib/libminc2-simple.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# =========================================================================== #
# Compliant julia 1.x

# Using the following packages
using JuMP, GLPK
using LinearAlgebra

include("loadSPP.jl")
include("setSPP.jl")
include("getfname.jl")

# =========================================================================== #

# Setting the data
fname = "Desktop/solveSPP/Data/didactic.dat"  # path for a standard config on macOS
cost, matrix = loadSPP(fname)

# Proceeding to the optimization
solverSelected = GLPK.Optimizer
ip, ip_x = setSPP(solverSelected, cost, matrix)
println("Solving..."); optimize!(ip)

# Displaying the results
println("z  = ", objective_value(ip))
print("x  = "); println(value.(ip_x))

# =========================================================================== #

# Collecting the names of instances to solve
target = "Desktop/solveSPP/Data"            # path for a standard config on macOS
fnames = getfname(target)

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
# This file was generated, do not modify it. # hide
X = source()
y = source(kind=:target)

model1 = linear
model2 = knn

m1 = machine(model1, X, y)
y1 = predict(m1, X)

m2 = machine(model2, X, y)
y2 = predict(m2, X)

yhat = 0.5*y1 + 0.5*y2

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
function embed(data::DataFrame;path::AbstractString="",xdim::Int64=20,ydim::Int64=20,perplexity=10,maxIter=10000,eta=200.0)
	# Kratochvíl M, Koladiya A and Vondrášek J. Generalized EmbedSOM on quadtree-structured self-organizing maps. F1000Research 2020, 8:2120

	imputed = Impute.srs(data)
	disallowmissing!(imputed)

	if isfile(path)
		@info """Using embedding saved at $path"""
		som = deserialize(path)

	else
		@info """Calculating embedding..."""
		som = initGigaSOM(imputed,xdim,ydim)
		som = trainGigaSOM(som,imputed)
		~isempty(path) && serialize(path,som)
	end

	######################## extract clusters and embedding
	try
		clusters = mapToGigaSOM(som,imputed)
		embedding = embedGigaSOM(som,imputed)

		rename!(clusters,"index"=>"cluster")
		return clusters, Matrix{Float64}(embedding')
	catch

		@info """Re-calculating embedding..."""
		som = initGigaSOM(imputed,xdim,ydim)
		som = trainGigaSOM(som,imputed)
		~isempty(path) && serialize(path,som)

		clusters = mapToGigaSOM(som,imputed)
		embedding = embedGigaSOM(som,imputed)

		rename!(clusters,"index"=>"cluster")
		return clusters, Matrix{Float64}(embedding')
	end
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
using DataFrames
using EMIRT

"""
transform error curve to dataframe
"""
function ec2df(ec::ScoreCurve, ecname::Symbol=:ecname)
    df = DataFrame()
    for (k,v) in ec
        df[k] = v
    end
    df[:tag] = ecname
    return df
end

"""
transfer error curves to dataframe
"""
function ecs2df(ecs::ScoreCurves)
    df = DataFrame()
    for (ecname,ec) in ecs
        if isempty(df)
            df = ec2df(ec, ecname)
        else
          metrics = collect(keys(first(values(ecs))))
          @show metrics
          df = join(df, ec2df(ec, ecname), on=metrics, kind=:outer)
        end
    end
    return df
end


text:
"Represents a point in physical space"
type Node{dim} <: DofObject
    "A unique ID for the Node"
    id::Int64

    "x,y,z coordinates"
    coords::Vec{2, Float64}

    "Degrees of freedom assigned to this Node"
    dofs::Array{Int64}

    "ProcessorID this Node is assigned to"
    processor_id::Int64
end

import Base.show

function show(io::IO, node::Node)
    println("Node: ", node.id)
    println("  Coords: ", node.coords)
    println("  Dofs: ", node.dofs)
    println("  pid: ", node.processor_id)
end


text:
function PCA(X)
    mean_vec = mean(X, dims=2)
    X = X .- mean_vec
    M = X'*X
    e, EV = eigen(M)
    tmp = (X*EV)'
    V = tmp[end:-1:1,:]
    S = reverse(sqrt.(e))
    for i in 1:size(V, 2)
        V[:,i] ./= S
    end
    V, S, mean_vec
end

function get_projection(image_matrix, d)
    V, S, mean_vec = PCA(image_matrix)
    V'[:,1:d], mean_vec
end

function image_to_eigenface_rep(image, proj, mean_vec)
    proj'*(reshape(convert(Array{Float64}, img), :, 1) - mean_vec)
end

function eigenface_rep_to_image(eigenface_rep, proj, mean_vec, dims)
    image = proj*eigenface_rep + mean_vec
    Gray.(N0f8.(clamp01.(reshape(image, dims[1], dims[2]))))
end

function get_eigenfaces(proj, d, dims)
    eig(n) = rescale(proj[:,n])
    eigenfaces = reduce(hcat,[eig(i) for i in 1:d])
    image_matrix_to_images(eigenfaces, dims)
end


text:
using MLDatasets
using Flux, Statistics
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser
using Flux: onehotbatch, onecold, @epochs
using Flux.Losses: logitcrossentropy
using Base: @kwdef
using CUDA
using MLDatasets
using Images
using ProgressBars

function get_data(args)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    # load train and test dataset
    x_train, y_train = CIFAR10.traindata(Float32)
    x_test,  y_test  = CIFAR10.testdata(Float32)

    # reshape
    x_train = reshape(x_train, 32, 32, 3, :)
    x_test = reshape(x_test, 32, 32, 3, :)

    # resize
    x_train = imresize(x_train, (64, 64, 3))
    x_test = imresize(x_test, (64, 64, 3))

    # one-hot-encode the labels
    y_train = onehotbatch(y_train, 0:9)
    y_test = onehotbatch(y_test, 0:9)

    train_loader = DataLoader((x_train, y_train), batchsize=args.batchsize, partial=false, shuffle=true)
    test_loader = DataLoader((x_test, y_test), batchsize=args.batchsize, partial=false)

    return train_loader, test_loader

end

function loss_function(ŷ, y)
    logitcrossentropy(ŷ, y)
end

function evaluation_loss_accuracy(loader, model)

    loss, accuracy, counter = 0f0, 0f0, 0

    for (x,y) in loader
        ŷ = model(x)
        loss += loss_function(ŷ,y)
        accuracy += sum(onecold(ŷ) .== onecold(y))
        counter +=  1
    end

    return loss / counter, accuracy / counter
end

# AdaptiveMeanPool cannot be used with 64x64 because it will throw a
# DimensionMismatch. You can use it with 256x256 or bigger images
# For more information refer to: https://github.com/FluxML/model-zoo/issues/334
function set_model(; imgsize=(64, 64,3), num_classes=10)
    return Chain(
                Conv((11, 11), imgsize[end]=>64, stride=(4,4), relu, pad=(2,2)),
                MaxPool((3, 3), stride=(2,2)),
                Conv((5, 5), 64=>192, relu, pad=(2,2)),
                MaxPool((3, 3), stride=(2,2)),
                Conv((3, 3), 192=>384, relu, pad=(1,1)),
                Conv((3, 3), 384=>256, relu, pad=(1,1)),
                Conv((3, 3), 256=>256, relu, pad=(1,1)),
                MaxPool((3, 3), stride=(2,2)),
                # AdaptiveMeanPool((6, 6)),
                flatten,
                Dropout(0.5),
                Dense(256*1*1, 4096, relu), # With AdaptiveMeanPool((6, 6)) set 256 * 6 * 6
                Dropout(0.5),
                Dense(4096, 4096, relu),
                Dense(4096, num_classes))
end

Base.@kwdef mutable struct Args
    η = 1e-4             # learning rate
    batchsize = 128     # batch size
    epochs = 1          # number of epochs
end

function train(; kws...)
    args = Args(; kws...) # collect options in a struct for convenience

    @info "Getting data..."

    train_loader, test_loader = get_data(args)

    model = set_model() # my implmentation
    # model = alexnet() # Metalhead's implmentation
    # model = LeNet5() # test pipeline

    ps = Flux.params(model)
    opt = ADAM(args.η)

    @info "Start training..."

    for epoch in 1:args.epochs

        for (x,y) in ProgressBar(train_loader)
            @assert size(x) == (64, 64, 3, args.batchsize)
            @assert size(y) == (10, args.batchsize)

            gs = gradient(() -> loss_function(model(x), y), ps)
            Flux.Optimise.update!(opt, ps, gs)
        end

        train = evaluation_loss_accuracy(train_loader, model)
        test = evaluation_loss_accuracy(test_loader, model)

        println("Epoch $(epoch-1)")
        println("\t Train => loss = $(train[1]) \t acc = $(train[2])")
        println("\t Test => loss = $(test[1]) \t acc = $(test[2])")

    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    train()
end


text:
using Printf

struct MyTimeOld
    hour :: Int64
    minute :: Int64
    second :: Int64
end

# This is a copy constructor, it creates a copy of the original object.
function MyTimeOld(time::MyTimeOld)
    MyTimeOld(time.hour, time.minute, time.second)
end

println(MyTimeOld(1, 2, 3))

# Inner constructors. Note that there will be no default constructor.
struct MyTime
    hour :: Int64
    minute :: Int64
    second :: Int64
    function MyTime(hour::Int64=0, minute::Int64=0, second::Int64=0)
        @assert(0 ≤ minute < 60, "Minute is not between 0 and 60.")
        @assert(0 ≤ second < 60, "Second is not between 0 and 60.")
        new(hour, minute, second)
    end
end

a = MyTime(11)
b = MyTime(11, 23)
c = MyTime(11, 23, 45)

println(MyTime())
println(a)
println(b)
println(c)

mutable struct MMyTime
    hour :: Int
    minute :: Int
    second :: Int
    function MMyTime(hour::Int64=0, minute::Int64=0, second::Int64=0)
        @assert(0 ≤ minute < 60, "Minute is between 0 and 60.")
        @assert(0 ≤ second < 60, "Second is between 0 and 60.")
        time = new()
        time.hour = hour
        time.minute = minute
        time.second = second
        time
    end
end

function Base.show(io::IO, time::MyTime)
    @printf(io, "%02d:%02d:%02d", time.hour, time.minute, time.second)
end

time = MyTime(9, 45)
print(time)

# Overloading
function timetoint(time::MyTime)
    minutes = time.hour * 60 + time.minute
    seconds = minutes * 60 + time.second
end

function inttotime(seconds::Int64)
    (minutes, second) = divrem(seconds, 60)
    hour, minute = divrem(minutes, 60)
    MyTime(hour, minute, second)
end

import Base.+

function +(t1::MyTime, t2::MyTime)
    seconds = timetoint(t1) + timetoint(t2)
    inttotime(seconds)
end

start = MyTime(9, 45)
duration = MyTime(1, 35, 0)
println(start + duration)

function increment(time, seconds)
    return inttotime(timetoint(time) + seconds)
end

function +(time::MyTime, seconds::Int64)
    increment(time, seconds)
end

function +(seconds::Int64, time::MyTime)
    time + seconds
end

println(start + 5)
println(10 + start)

function histogram(s)
    d = Dict()
    for c in s
        if c ∉ keys(d)
            d[c] = 1
        else
            d[c] += 1
        end
    end
    d
end

println(histogram(("spam", "egg", "spam", "spam", "bacon", "spam")))

println(sum((a, b, c)))