text:
## make options for the stationary GP
fileprefix = "MT_only_$(string(F.stretch))_"
K = transD_GP.GP.OrstUhn()
nmin, nmax = 2, 40
fbounds = F.stretch ? [0 1.] : [-0.5 4]
sdev_pos = [0.05abs(diff([extrema(znall)...])[1])]
sdev_prop = 0.07*diff(fbounds, dims=2)[:]
demean = false
sdev_dc = 0.008*diff(fbounds, dims=2)[:]
sampledc = true
xall = permutedims(collect(znall))
xbounds = permutedims([extrema(znall)...])
λ, δ = [2], F.stretch ? 0.03 : 0.1
debug = false
## Initialize a stationary GP using these options
using Random
Random.seed!(12)
opt = transD_GP.OptionsStat(debug = debug,
                        nmin = nmin,
                        nmax = nmax,
                        xbounds = xbounds,
                        fbounds = fbounds,
                        fdataname = fileprefix,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        save_freq = 50,
                        demean = demean,
                        sampledc = sampledc,
                        sdev_dc = sdev_dc,
                        quasimultid = false,
                        K = K
                        )

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
mutable struct Figure
    # NOTE: couldn't use @kw_args here, clashed with the format of the base constructor.
    g ::GS     # description buffer
    id::String # id of the figure
    # ---
    axes     ::Vector{Axes}   # subplots
    size     ::T2F               # (width, heigth)
    textstyle::TextStyle         # parent font etc
    bgcolor  ::Option{Colorant}  # background col, nothing=transparent
    # ---
    texlabels   ::Option{Bool}   # true if has tex
    texscale    ::Option{String} # scale latex * hei (def=1)
    texpreamble ::Option{String} # latex preamble
    transparency::Option{Bool}   # if true, use cairo device
    # ---
    subroutines::Dict{String,String}
end
Figure(g::GS, id::String) =
    Figure(
        g, id,
        # --
        Vector{Axes}(),
        (12., 9.),
        TextStyle(),
        c"white",
        # --
        ∅,
        ∅,
        ∅,
        ∅,
        # --
        Dict{String,String}()
    )

"""
    reset!(f)

Internal function to completely refresh the figure `f` only keeping its id and size.
"""
function reset!(f::Figure)
    take!(f.g) # empty the buffer
    f.axes         = Vector{Axes}() # clean axes
    f.textstyle    = TextStyle(font="texcmss", hei=0.35) # default fontstyle
    f.bgcolor      = c"white" # default bg color
    f.texlabels    = ∅
    f.texscale     = ∅
    f.texpreamble  = ∅
    f.transparency = ∅
    f.subroutines  = Dict{String,String}()
    GP_ENV["CUR_FIG"]  = f
    GP_ENV["CUR_AXES"] = nothing
    return f
end

ext: jl


text:
"""
A structure which defines a hypercube in `N`-dimensional searchspace, with a `center` location
and conrners defined by `ranges`.
"""
struct Block{N,T<:Real}
    ranges::NTuple{N, Tuple{T,T}}
    center::NTuple{N,T}
    lowerbound::T
    upperbound::T
end
Block{N,T}() where {N,T} = Block{N,T}(ntuple(x->(zero(T),zero(T)),N), ntuple(x->(zero(T)),N), typemax(T), typemax(T))

numbertype(::Block{N,T}) where {N,T} = T
dims(::Block{N,T}) where {N,T} = N
Base.size(::Block{N,T}) where {N,T} = (N,)

# for speeding up hashing and performance of the priority queue in the branch and bound procedure
const hash_block_seed = UInt === UInt64 ? 0x03f7a7ad5ef46a89 : 0xa9bf8ce0
function hash(B::Block, h::UInt)
    h += hash_block_seed
    h = hash(B.ranges, h)
    return h
end

CoordinateTransformations.AffineMap(bl::Block) = AffineMap(bl.center)

#
function boxranges(center::Union{Tuple, NTuple, AbstractArray}, widths::Union{Tuple, NTuple, AbstractArray})
    t = eltype(center)
    return NTuple{length(center),Tuple{t,t}}([(center[i]-widths[i], center[i]+widths[i]) for i=1:length(center)])
end

"""
    sbrngs = subranges(ranges, nsplits)

Takes `ranges`, a nested tuple describing intervals for each dimension in rigid-rotation space
defining a hypercube, and splits the hypercube into `nsplits` even components along each dimension.
Since the space is 6-dimensional, the number of returned sub-cubes will be `nsplits^6`.
"""
function subranges(ranges, nsplits::Int)
    ndims = length(ranges)
    t = eltype(eltype(ranges))

    # calculate even splititng points for each dimension
    splitvals = [range(r[1], stop=r[2], length=nsplits+1) |> collect for r in ranges]
    splits = [[(splitvals[i][j], splitvals[i][j+1]) for j=1:nsplits] for i=1:ndims]
    f(x) = splits[x[1]][x[2]]
    children = fill(ranges, nsplits^ndims)
    for (i,I) in enumerate(CartesianIndices(NTuple{ndims,UnitRange{Int}}(fill(1:nsplits, ndims))))
        children[i] = NTuple{ndims,Tuple{t,t}}(map(x->f(x), enumerate(Tuple(I))))
    end
    return children
end

"""
    lim = translation_limit(gmmx, gmmy)

Computes the largest translation needed to ensure that the searchspace contains the best alignment transformation.
"""
function translation_limit(gmmx::AbstractSingleGMM, gmmy::AbstractSingleGMM)
    trlim = typemin(promote_type(numbertype(gmmx),numbertype(gmmy)))
    for gaussians in (gmmx.gaussians, gmmy.gaussians)
        if !isempty(gaussians)
            trlim = max(trlim, maximum(gaussians) do gauss
                    maximum(abs, gauss.μ) end)
        end
    end
    return trlim
end

function translation_limit(mgmmx::AbstractMultiGMM, mgmmy::AbstractMultiGMM)
    trlim = typemin(promote_type(numbertype(mgmmx),numbertype(mgmmy)))
    for key in keys(mgmmx.gmms) ∩ keys(mgmmy.gmms)
        trlim = max(trlim, translation_limit(mgmmx.gmms[key], mgmmy.gmms[key]))
    end
    return trlim
end

# Block for GOGMA procedure in full transformational space
function fullBlock(gmmx::AbstractGMM, gmmy::AbstractGMM, ranges=nothing, pσ=nothing, pϕ=nothing, rot=nothing, trl=nothing)
    # get center and uncertainty region
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    if isnothing(ranges)
        trlim = translation_limit(gmmx, gmmy)
        ranges = ((-t(π),t(π)), (-t(π),t(π)), (-t(π),t(π)), (-trlim,trlim), (-trlim,trlim), (-trlim,trlim))
    end
    center = NTuple{length(ranges),t}([sum(dim)/2 for dim in ranges])
    rwidth = ranges[1][2] - ranges[1][1]
    twidth = ranges[4][2] - ranges[4][1]

    # calculate objective function bounds for the block
    lb, ub = get_bounds(gmmx, gmmy, rwidth, twidth, center, pσ, pϕ)

    return Block(ranges, center, lb, ub)
end

# Block for rotational space only (i.e. the first stage of TIV-GOGMA)
function rotBlock(gmmx::AbstractGMM, gmmy::AbstractGMM, ranges=nothing, pσ=nothing, pϕ=nothing, rot=nothing, trl=nothing)
    # get center and uncertainty region
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    if isnothing(ranges)
        ranges = ((-t(π),t(π)), (-t(π),t(π)), (-t(π),t(π)))
    end
    center = NTuple{length(ranges),t}([sum(dim)/2 for dim in ranges])
    rwidth = ranges[1][2] - ranges[1][1]

    # calculate objective function bounds for the block
    if isnothing(trl)
        trl = (zero(t), zero(t), zero(t))
    end
    lb, ub = get_bounds(gmmx, gmmy, rwidth, zero(t), (center..., trl...), pσ, pϕ)

    return Block(ranges, center, lb, ub)
end

# Block for translational space only (i.e. the second stage of TIV-GOGMA)
function trlBlock(gmmx::AbstractGMM, gmmy::AbstractGMM, ranges=nothing, pσ=nothing, pϕ=nothing, rot=nothing, trl=nothing)
    # get center and uncertainty region
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    if isnothing(ranges)
        trlim = translation_limit(gmmx, gmmy)
        ranges = ((-trlim,trlim), (-trlim,trlim), (-trlim,trlim))
    end
    center = NTuple{length(ranges),t}([sum(dim)/2 for dim in ranges])
    twidth = ranges[1][2] - ranges[1][1]

    # calculate objective function bounds for the block
    if isnothing(rot)
        rot = (zero(t), zero(t), zero(t))
    end
    lb, ub = get_bounds(gmmx, gmmy, zero(t), twidth, (rot..., center...), pσ, pϕ)

    return Block(ranges, center, lb, ub)
end


text:
"
Creates the experiment directory form the parameter structure.
"
function create_experiment_directory(params)
    if !isnothing(params.data_params)
        data_tag = @sprintf("m%d-σ%1.3f-rs%d", params.data_params.data_num, params.data_params.noise_sigma, params.random_seed)
    else
        data_tag = "known-system"
    end
    for dim_key in keys(params.discretization_step)
        data_tag = @sprintf("%s-%0.3f", data_tag, params.discretization_step[dim_key])
    end
    exp_dir = @sprintf("%s/modes/%s/%s", params.experiment_directory, params.system_params.mode_tag, data_tag)
    !isdir(exp_dir) && mkpath(exp_dir)
    return exp_dir
end

function initialize_log(exp_dir; logging=Logging.Info)
    glogger = SimpleLogger(stdout, logging)
    logfile = open("$exp_dir/log.txt", "w+")
    text_logger = SimpleLogger(logfile, logging)
    demux_logger = TeeLogger(glogger, text_logger)
    global_logger(demux_logger)
    @info "Experiment directory: " exp_dir
    return logfile
end

function save_time_info(exp_dir, time_info)
    time_filename = "$exp_dir/time_info.bson"
    bson(time_filename, time_info)
end


text:
# Page not found

notfound(s = "Not found") = mux(status(404), respond(s))


text:
"""
    Coveralls(; config_file::Union{AbstractString, Nothing}=nothing) -> Coveralls

Add `Coveralls` to a template's plugins to optionally add a `.coveralls.yml` configuration
file to generated repositories, and an appropriate badge to the README. Also updates the
`.gitignore` accordingly.

# Keyword Arguments:
* `config_file::Union{AbstractString, Nothing}=nothing`: Path to a custom `.coveralls.yml`.
  If left unset, no file will be generated.
"""
struct Coveralls <: GenericPlugin
    gitignore::Vector{String}
    src::Union{String, Nothing}
    dest::String
    badges::Vector{Badge}
    view::Dict{String, Any}

    function Coveralls(; config_file::Union{AbstractString, Nothing}=nothing)
        if config_file != nothing
            config_file = if isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError("File $(abspath(config_file)) does not exist"))
            end
        end
        new(
            ["*.jl.cov", "*.jl.*.cov", "*.jl.mem"],
            config_file,
            ".coveralls.yml",
            [
                Badge(
                    "Coveralls",
                    "https://coveralls.io/repos/github/{{USER}}/{{PKGNAME}}.jl/badge.svg?branch=master",
                    "https://coveralls.io/github/{{USER}}/{{PKGNAME}}.jl?branch=master",
                ),
            ],
            Dict{String, Any}(),
        )
    end
end

interactive(::Type{Coveralls}) = interactive(Coveralls; file=nothing)


text:
"""Implements the CreditMetrics model for portfolio credit risk"""
module CreditMetrics

using StatsFuns

include("threshold_links.jl")
include("thresholds.jl")

end # module


text:
using MixedModels
using MixedModels: dataset

@isdefined(gfms) || const global gfms = Dict(
    :cbpp => [@formula((incid/hsz) ~ 1 + period + (1|herd))],
    :contra => [@formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist))],
    :grouseticks => [@formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location))],
    :verbagg => [@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item))],
)

@isdefined(fms) || const global fms = Dict(
    :oxide => [@formula(Thickness ~ 1 + (1|Lot/Wafer)),
               @formula(Thickness ~ 1 + Source + (1+Source|Lot) + (1+Source|Lot&Wafer))],
    :dyestuff => [@formula(yield ~ 1 + (1|batch))],
    :dyestuff2 => [@formula(yield ~ 1 + (1|batch))],
    :d3 => [@formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i))],
    :insteval => [
        @formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)),
        @formula(y ~ 1 + service*dept + (1|s) + (1|d)),
    ],
    :kb07 => [
        @formula(rt_trunc ~ 1+spkr+prec+load+(1|subj)+(1|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1|subj)+(1+prec|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1+spkr+prec+load|subj)+(1+spkr+prec+load|item)),
    ],
    :pastes => [
        @formula(strength ~ 1 + (1|batch&cask)),
        @formula(strength ~ 1 + (1|batch/cask)),
    ],
    :penicillin => [@formula(diameter ~ 1 + (1|plate) + (1|sample))],
    :sleepstudy => [
        @formula(reaction ~ 1 + days + (1|subj)),
        @formula(reaction ~ 1 + days + zerocorr(1+days|subj)),
        @formula(reaction ~ 1 + days + (1|subj) + (0+days|subj)),
        @formula(reaction ~ 1 + days + (1+days|subj)),
    ],
)

# for some reason it seems necessary to prime the pump in julia-1.6.0-DEV
@isdefined(fittedmodels) || const global fittedmodels = Dict{Symbol,Vector{LinearMixedModel}}(
    :dyestuff => [fit(MixedModel, only(fms[:dyestuff]), dataset(:dyestuff))]
);

function models(nm::Symbol)
    get!(fittedmodels, nm) do
        [fit(MixedModel, f, dataset(nm)) for f in fms[nm]]
    end
end


text:
module HOHQMesh

import HOHQMesh_jll

export generate_mesh

"""
    generate_mesh(control_file;
                  output_directory="out",
                  mesh_filename=nothing, plot_filename=nothing, stats_filename=nothing,
                  verbose=false)

Generate a mesh based on the `control_file` with the HOHQMesh mesh generator and store resulting
files in `output_directory`.

You can set the mesh filename, the plot filename, and the statistics filename using the keyword
arguments `mesh_filename`, `plot_filename`, and `stats_filename`, respectively. If set to `nothing`,
the filenames for the mesh file, plot file, and statistics file are generated automatically from the
control file name. For example, `path/to/ControlFile.control` will result in output files
`ControlFile.mesh`, `ControlFile.tec`, and `ControlFile.txt`.

You can activate verbose output from HOHQMesh that prints additional messages and debugging
mesh information with the keyword argument `verbose`.

This function returns the output to `stdout` of the HOHQMesh binary when generating the mesh.
"""
function generate_mesh(control_file;
                       output_directory="out",
                       mesh_filename=nothing, plot_filename=nothing, stats_filename=nothing,
                       verbose=false)
  @assert isfile(control_file) "'$control_file' is not a valid path to an existing file"

  # Determine output filenames
  filebase = splitext(basename(control_file))[1]
  if isnothing(mesh_filename)
    mesh_file_format = extract_mesh_file_format(control_file)
    if mesh_file_format == "ISM" || mesh_file_format == "ISM-V2" || mesh_file_format == "ISM-v2"
      mesh_filename = filebase * ".mesh"
    elseif mesh_file_format == "ABAQUS"
      mesh_filename = filebase * ".inp"
    else
      error("Unknown mesh file format: ", mesh_file_format, " (must be one of ISM, ISM-V2, or ABAQUS)")
    end
  end
  if isnothing(plot_filename)
    plot_filename = filebase * ".tec"
  end
  if isnothing(stats_filename)
    stats_filename = filebase * ".txt"
  end

  # Determine output filepaths
  mesh_filepath = joinpath(output_directory, mesh_filename)
  plot_filepath = joinpath(output_directory, plot_filename)
  stats_filepath = joinpath(output_directory, stats_filename)

  # Create output directory if it does not exist
  if !isdir(output_directory)
    mkdir(output_directory)
  end

  # Create temporary copy of mesh file to replace paths
  output = mktemp() do tmppath, tmpio
    # Update paths of output files
    lines = readlines(control_file, keep=true)
    for line in lines
      if occursin("mesh file name", line)
        write(tmpio, "      mesh file name   = " * mesh_filepath * "\n")
      elseif occursin("plot file name", line)
        write(tmpio, "      plot file name   = " * plot_filepath * "\n")
      elseif occursin("stats file name", line)
        write(tmpio, "      stats file name  = " * stats_filepath * "\n")
      else
        write(tmpio, line)
      end
    end
    flush(tmpio)

    # Run HOHQMesh and store output
    if verbose
      readchomp(`$(HOHQMesh_jll.HOHQMesh()) -verbose -f $tmppath`)
    else
      readchomp(`$(HOHQMesh_jll.HOHQMesh()) -f $tmppath`)
    end
  end

  String(output)
end

"""
    extract_mesh_file_format(control_file)

Return a string with the desired output format of the HOHQMesh generated mesh file.
This information is given within the `RUN_PARAMETERS` of the `CONTROL_INPUT` block
of the control file.
See the [`HOHQMesh` documentation](https://trixi-framework.github.io/HOHQMesh/) for details.
"""
function extract_mesh_file_format(control_file)
  # Find the file line that gives the mesh file format
  file_lines = readlines(open(control_file))
  line_index = findfirst(contains("mesh file format"), file_lines)
  # Extract the mesh file format keyword
  file_format = split(file_lines[line_index])[5]

  return file_format
end

"""
    examples_dir()

Return the path to the directory with some example mesh setups.
"""
examples_dir() = joinpath(pathof(HOHQMesh) |> dirname |> dirname, "examples")

end # module


text:
"""
Prototype for IPs based on "Machine-Learning", to be imported by
modules that either define basis sets or regression methods
"""
module MLIPs

using JuLIP:       AbstractCalculator, AbstractAtoms
using JuLIP.FIO:   decode_dict

import JuLIP:      energy, forces, virial, site_energy, site_energy_d
import Base:       Dict, convert, ==

export IPSuperBasis, IPCollection, combine

"""
`abstract type IPBasis` : A type derived from `IPBasis` defines
not an IP as such but a collection of IPs that can be understood as
a basis set. E.g., a 2B-potential could be obtained from
polynomials,
```
   ϕ(r) = ∑_k c_k exp(-k r)
```
then the `B::IPBasis` might specify the basis functions {r -> exp(-kr)} as
a collection of `PairPotential`s.

Then, `energy(at, B)` does not return a single energy but a vector
of energies.

## Developer Notes

* `length(basis::IPBasis)` must return the number of basis functions
"""
abstract type IPBasis end

# ========== wrap one more more calculators into a basis =================

struct IPCollection{T} <: IPBasis
   coll::Vector{T}
end

IPCollection(args...) = IPCollection([args...])
Base.length(coll::IPCollection) = length(coll.coll)
combine(basis::IPCollection, coeffs::AbstractVector{<:Number}) =
      SumIP([ c * b for (c, b) in zip(coeffs, basis.coll) ])

energy(coll::IPCollection, at::AbstractAtoms) =
         [ energy(B, at) for B in coll.coll ]
forces(coll::IPCollection, at::AbstractAtoms) =
         [ forces(B, at) for B in coll.coll ]
virial(coll::IPCollection, at::AbstractAtoms) =
         [ virial(B, at) for B in coll.coll ]
site_energy(coll::IPCollection, at::AbstractAtoms, i0::Integer) =
         [ site_energy(V, at, i0) for V in coll.coll ]
site_energy_d(coll::IPCollection, at::AbstractAtoms, i0::Integer) =
         [ site_energy_d(V, at, i0) for V in coll.coll ]

Dict(coll::IPCollection) = Dict(
      "__id__" => "JuLIP_IPCollection",
      "coll" => Dict.(coll.coll) )
IPCollection(D::Dict) = IPCollection( decode_dict.( D["coll"] ) )
convert(::Val{:JuLIP_IPCollection}, D::Dict) = IPCollection(D)
import Base.==
==(B1::IPCollection, B2::IPCollection) = all(B1.coll .== B2.coll)

# ========== IPSuperBasis : combine multiple sub-basis =================

"""
`struct IPSuperBasis:` a collection of IP basis sets, re-interpreted
as a large basis
"""
struct IPSuperBasis{TB <: IPBasis} <: IPBasis
   BB::Vector{TB}
end

_convert_basis_(b::IPBasis) = b
_convert_basis_(b::AbstractCalculator) = IPCollection(b)
_convert_basis_(b::AbstractVector) = IPCollection(b)

IPSuperBasis(args...) = IPSuperBasis([_convert_basis_(b) for b in args])

Base.length(super::IPSuperBasis) = sum(length.(super.BB))

function combine(basis::IPSuperBasis, coeffs::AbstractVector{<:Number})
   i0 = 0
   components = []
   for B in basis.BB
      lenB = length(B)
      push!(components, combine(B, coeffs[i0+1:i0+lenB]))
      i0 += lenB
   end
   return SumIP(components)
end

energy(superB::IPSuperBasis, at::AbstractAtoms) =
         vcat([ energy(B, at) for B in superB.BB ]...)
forces(superB::IPSuperBasis, at::AbstractAtoms) =
         vcat([ forces(B, at) for B in superB.BB ]...)
virial(superB::IPSuperBasis, at::AbstractAtoms) =
         vcat([ virial(B, at) for B in superB.BB ]...)
site_energy(superB::IPSuperBasis, at::AbstractAtoms, i0::Integer) =
         vcat([ site_energy(B, at, i0) for B in superB.BB ]...)
site_energy_d(superB::IPSuperBasis, at::AbstractAtoms, i0::Integer) =
         vcat([ site_energy_d(B, at, i0) for B in superB.BB ]...)

Dict(superB::IPSuperBasis) = Dict(
      "__id__" => "JuLIP_IPSuperBasis",
      "components" => Dict.(superB.BB) )
IPSuperBasis(D::Dict) = IPSuperBasis( decode_dict.( D["components"] ) )
convert(::Val{:JuLIP_IPSuperBasis}, D::Dict) = IPSuperBasis(D)
import Base.==
==(B1::IPSuperBasis, B2::IPSuperBasis) = all(B1.BB .== B2.BB)

# ========== SumIP =================
# a sum of several IPs.

struct SumIP{T} <: AbstractCalculator
   components::Vector{T}
end

SumIP(args...) = SumIP([args...])

energy(sumip::SumIP, at::AbstractAtoms; kwargs...) =
         sum(energy(calc, at; kwargs...) for calc in sumip.components)
forces(sumip::SumIP, at::AbstractAtoms; kwargs...) =
         sum(forces(calc, at; kwargs...) for calc in sumip.components)
virial(sumip::SumIP, at::AbstractAtoms; kwargs...) =
         sum(virial(calc, at; kwargs...) for calc in sumip.components)
site_energy(sumip::SumIP, at::AbstractAtoms, i0::Integer) =
         sum(site_energy(V, at, i0) for V in sumip.components)
site_energy_d(sumip::SumIP, at::AbstractAtoms, i0::Integer) =
         sum(site_energy_d(V, at, i0) for V in sumip.components)

Dict(sumip::SumIP) = Dict(
      "__id__" => "JuLIP_SumIP",
      "components" => Dict.(sumip.components) )
SumIP(D::Dict) = SumIP( decode_dict.( D["components"] ) )
convert(::Val{:JuLIP_SumIP}, D::Dict) = SumIP(D)

SumIP(V::AbstractCalculator, sumip::SumIP) =
   SumIP( [ [V]; sumip.components ]  )
SumIP(sumip::SumIP, V::AbstractCalculator) =
   SumIP( [ sumip.components; [V] ]  )
SumIP(sum1::SumIP, sum2::SumIP) =
   SumIP( [ sum1.components; sum2.components ]  )

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

mutable struct _SubexpressionStorage
    nodes::Vector{Nonlinear.Node}
    adj::SparseArrays.SparseMatrixCSC{Bool,Int}
    const_values::Vector{Float64}
    forward_storage::Vector{Float64}
    partials_storage::Vector{Float64}
    reverse_storage::Vector{Float64}
    forward_storage_ϵ::Vector{Float64}
    partials_storage_ϵ::Vector{Float64}
    reverse_storage_ϵ::Vector{Float64}
    linearity::Linearity

    function _SubexpressionStorage(
        expr::Nonlinear.Expression,
        subexpression_linearity,
        moi_index_to_consecutive_index,
    )
        nodes =
            _replace_moi_variables(expr.nodes, moi_index_to_consecutive_index)
        adj = Nonlinear.adjacency_matrix(nodes)
        N = length(nodes)
        linearity = _classify_linearity(nodes, adj, subexpression_linearity)
        return new(
            nodes,
            adj,
            expr.values,
            zeros(N),  # forward_storage,
            zeros(N),  # partials_storage,
            zeros(N),  # reverse_storage,
            Float64[],
            Float64[],
            Float64[],
            linearity[1],
        )
    end
end

mutable struct _FunctionStorage
    nodes::Vector{Nonlinear.Node}
    adj::SparseArrays.SparseMatrixCSC{Bool,Int}
    const_values::Vector{Float64}
    forward_storage::Vector{Float64}
    partials_storage::Vector{Float64}
    reverse_storage::Vector{Float64}
    grad_sparsity::Vector{Int}
    # Nonzero pattern of Hessian matrix
    hess_I::Vector{Int}
    hess_J::Vector{Int}
    rinfo::Coloring.RecoveryInfo # coloring info for hessians
    seed_matrix::Matrix{Float64}
    linearity::Linearity
    # subexpressions which this function depends on, ordered for forward pass.
    dependent_subexpressions::Vector{Int}

    function _FunctionStorage(
        nodes::Vector{Nonlinear.Node},
        const_values,
        num_variables,
        coloring_storage::Coloring.IndexedSet,
        want_hess::Bool,
        subexpressions::Vector{_SubexpressionStorage},
        dependent_subexpressions,
        subexpression_linearity,
        subexpression_edgelist,
        subexpression_variables,
        moi_index_to_consecutive_index,
    )
        nodes = _replace_moi_variables(nodes, moi_index_to_consecutive_index)
        adj = Nonlinear.adjacency_matrix(nodes)
        N = length(nodes)
        empty!(coloring_storage)
        _compute_gradient_sparsity!(coloring_storage, nodes)
        for k in dependent_subexpressions
            _compute_gradient_sparsity!(
                coloring_storage,
                subexpressions[k].nodes,
            )
        end
        grad_sparsity = sort!(collect(coloring_storage))
        empty!(coloring_storage)
        if want_hess
            linearity = _classify_linearity(nodes, adj, subexpression_linearity)
            edgelist = _compute_hessian_sparsity(
                nodes,
                adj,
                linearity,
                coloring_storage,
                subexpression_edgelist,
                subexpression_variables,
            )
            hess_I, hess_J, rinfo = Coloring.hessian_color_preprocess(
                edgelist,
                num_variables,
                coloring_storage,
            )
            seed_matrix = Coloring.seed_matrix(rinfo)
            return new(
                nodes,
                adj,
                const_values,
                zeros(N),  # forward_storage,
                zeros(N),  # partials_storage,
                zeros(N),  # reverse_storage,
                grad_sparsity,
                hess_I,
                hess_J,
                rinfo,
                seed_matrix,
                linearity[1],
                dependent_subexpressions,
            )
        else
            return new(
                nodes,
                adj,
                const_values,
                zeros(N),  # forward_storage,
                zeros(N),  # partials_storage,
                zeros(N),  # reverse_storage,
                grad_sparsity,
                Int[],
                Int[],
                Coloring.RecoveryInfo(),
                Array{Float64}(undef, 0, 0),
                NONLINEAR,
                dependent_subexpressions,
            )
        end
    end
end

"""
    NLPEvaluator(
        model::Nonlinear.Model,
        ordered_variables::Vector{MOI.VariableIndex},
    )

Return an `NLPEvaluator` object that implements the `MOI.AbstractNLPEvaluator`
interface.

!!! warning
    Before using, you must initialize the evaluator using `MOI.initialize`.
"""
mutable struct NLPEvaluator <: MOI.AbstractNLPEvaluator
    data::Nonlinear.Model
    ordered_variables::Vector{MOI.VariableIndex}

    objective::Union{Nothing,_FunctionStorage}
    constraints::Vector{_FunctionStorage}
    subexpressions::Vector{_SubexpressionStorage}
    subexpression_order::Vector{Int}
    # Storage for the subexpressions in reverse-mode automatic differentiation.
    subexpression_forward_values::Vector{Float64}
    subexpression_reverse_values::Vector{Float64}
    subexpression_linearity::Vector{Linearity}

    # A cache of the last x. This is used to guide whether we need to re-run
    # reverse-mode automatic differentiation.
    last_x::Vector{Float64}

    # Temporary storage for computing jacobians. This is also used as temporary
    # storage for the input of multivariate functions.
    jac_storage::Vector{Float64}
    # Temporary storage for the gradient of multivariate functions
    user_output_buffer::Vector{Float64}

    # storage for computing hessians
    # these Float64 vectors are reinterpreted to hold multiple epsilon components
    # so the length should be multiplied by the maximum number of epsilon components
    disable_2ndorder::Bool # don't offer Hess or HessVec
    want_hess::Bool
    forward_storage_ϵ::Vector{Float64} # (longest expression)
    partials_storage_ϵ::Vector{Float64} # (longest expression)
    reverse_storage_ϵ::Vector{Float64} # (longest expression)
    input_ϵ::Vector{Float64} # (number of variables)
    output_ϵ::Vector{Float64}# (number of variables)
    subexpression_forward_values_ϵ::Vector{Float64} # (number of subexpressions)
    subexpression_reverse_values_ϵ::Vector{Float64} # (number of subexpressions)
    hessian_sparsity::Vector{Tuple{Int64,Int64}}
    max_chunk::Int # chunk size for which we've allocated storage

    function NLPEvaluator(
        data::Nonlinear.Model,
        ordered_variables::Vector{MOI.VariableIndex},
    )
        return new(data, ordered_variables)
    end
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using CrystallographyBase: Lattice
using LinearAlgebra: I

export distortby, distort, strainstate

# See https://link.springer.com/content/pdf/10.1007%2F978-3-7091-0382-1_7.pdf and https://doi.org/10.2138/am-1997-1-207
distortby(lattice::Lattice, strain::TensorStrain) =
    Lattice((I + strain.data) * lattice.data)
distortby(lattice::Lattice, strain::EngineeringStrain) =
    distortby(lattice, TensorStrain(strain))
const distort = distortby  # For the sake of compatibility

strainstate(old::Lattice, new::Lattice) = TensorStrain(new.data / old.data - I)


text:
using Revise
using IMPES
using Plots

# using Plots, LinearAlgebra, IMPES
# using AnalyticalEOR

#* Define Relative Perms
kr_ww = SingleRelPerms(
        swr= 0.11,
        sor= 0.10,
        krw0= 0.2,
        kro0= 0.8,
        nw= 3.0,
        no= 2.0,
            )

kr_ow = SingleRelPerms(
        swr= 0.11,
        sor= 0.16,
        krw0= 0.4,
        kro0= 0.5,
        nw= 2.0,
        no= 3.0,
        )

#* Define reservoir
res = Reservoir(
    300.,
    200 * 0.3048,
    50 * 0.3048,
    100.0 * 9.869233e-16,
    0.2,
    )

#* Define fluids
fluid = Fluids(
            μw=1e-3,
            μo=5e-3,
)

#* Define grid
grid = Grid(
            res,
            100,
            )

grid = Grid(res,
                10,
                10)
#* Define BC

bc = BoundaryConditions(
                        [:neumann, :neumann, :neumann, :neumann],
                        [0., 0., 0., 0.]
                        )

sim = Simulation(zeros(grid.Nx, grid.Ny),
                0.2 * ones(grid.Nx, grid.Ny),
                zeros(grid.Nx, grid.Ny),
                100.,
                1.)

function Thalf(l::Int64, m::Int64, res::Reservoir,
            kr::RelPerms, fluid::Fluids, grid::Grid, sim::Simulation)
    Nx = grid.Nx
    Δx = grid.dx
    Δy = grid.dy
    k = res.k
    h = res.h
    μw = fluid.μw
    μo = fluid.μo
    s = sim.s
    p = sim.p
    θ = sim.θ

    # lj = ceil(Int64,l / Nx)
    # li = l - (lj-1) * Nx

    # mj = ceil(Int64, m / Nx)
    # mi = m - (mj-1) * Nx

    lj= l ÷ Nx + 1
    li= l - (lj - 1) * Nx

    mj= l ÷ Nx + 1
    mi= m - (mj - 1) * Nx

    if p[l] >= p[m]
        krw, kro = kr(s[l], θ)
    else
        krw, kro = kr(s[m], θ)
    end

    if lj==mj # Get x-direction transmissibilities
        Tw, To = (krw/μw, kro/μo ) .* (k * Δy * h / Δx)
    elseif li==mi # Get y-direction transmissibilities
        Tw, To = (krw/μw, kro/μo ) .* (k * Δx * h / Δy)
    elseif lj==mj && li==mi
        if l-1 % Nx == 0 || l % Nx == 0
            Tw, To = (krw/μw, kro/μo) .* (k * Δy * h / Δx)
        else
            Tw, To = (krw/μw, kro/μo) .* (k * Δx * h / Δy)
        end
    end
end

T(l,m) = Thalf(l,m, res, kr_ww, fluid, grid, sim)

Nx = grid.Nx
Ny = grid.Ny
N = grid.N

Tw = zeros(N, N)
To = zeros(N, N)

for l in 1:N

# l = 11

    if (l-1) % Nx != 0 #* Not in Left
        Tw[l, l-1], To[l, l-1] = .-T(l,l-1)
        Tw[l,l] -= Tw[l,l-1]
        To[l,l] -= To[l,l-1]
    end

    if l % Nx != 0 #* Not in Right
        Tw[l, l+1], To[l, l+1] = .-T(l,l+1)
        Tw[l,l] -= Tw[l,l+1]
        To[l,l] -= To[l,l+1]
    end

    if ceil(Int64, l / Nx) != 1 #* Not in Bottom
        Tw[l, l-Nx], To[l, l-Nx] = .-T(l,l-Nx)
        Tw[l,l] -= Tw[l,l-Nx]
        To[l,l] -= To[l,l-Nx]
    end

    if ceil(Int64, l / Nx) != Ny #* Not in Top
        Tw[l, l+Nx], To[l, l+Nx] = .-T(l,l+Nx)
        Tw[l,l] -= Tw[l,l+Nx]
        To[l,l] -= To[l,l+Nx]
    end
end

Tt = Tw .+ To
sparse(Tt)

Q = zeros(grid.N)
Q[1] = 426.5 * 0.3048 ^ 2 /3600 / 24
Q[end] = -426.5 * 0.3048 ^ 2 /3600 / 24

Q *= grid.dx * grid.dy

p = Tt \ Q
plot(reshape(p, grid.Nx, grid.Ny) / 6894, linetype=:contourf)

p

max_stars_repo_path: examples/Ex1 - WF2D.jl

max_issues_repo_path: examples/Ex1 - WF2D.jl

max_forks_repo_path: examples/Ex1 - WF2D.jl


text:
import ThreeJS
import Reactive: every

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"nested-props")

    map(every(1 / 30)) do t
        ThreeJS.outerdiv() <<
        (ThreeJS.initscene() <<
            [
                ThreeJS.mesh(0.0, 0.0, 0.0) <<
                [
                    ThreeJS.plane(1.0, 1.0),
                    ThreeJS.shadermaterial(open(readall, "assets/colorplane.vert", "r"), open(readall, "assets/colorplane.frag", "r");
                        :uniforms => Dict(:b => Dict(:type => "f", :value => sin(t)^2), :a => Dict(:type => "f", :value => cos(t)^2)))
                ],
                ThreeJS.camera(0.0, 0.0, 2.0)
            ]
        )
    end
end


text:
# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

using CloudMicrophysics, Documenter
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

pages = Any[
    "Home" => "index.md",
    "0-moment microphysics" => "Microphysics_0M.md",
    "1-moment microphysics" => "Microphysics_1M.md",
    "Aerosol activation" => "Aerosol_activation.md",
    "API" => "API.md",
    "References" => "References.md",
]

mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
)

makedocs(
    bib,
    sitename = "CloudMicrophysics.jl",
    strict = true,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [CloudMicrophysics],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/CloudMicrophysics.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
@doc """
       Newton method for local optimization (multivariable)

       Supported by:

       Ph D Kelvyn B. Sánchez
       e-mail: kelvyn.baruc@gmail.com
       GitHub: Kelvyn88
""" ->

function RandomSearch(f, vars, lb, up; maxiter = 100, iter = true)
      #try
            const version = Pkg.installed("OptimJulia")
            length(vars) == 0 ? len = 1 : len = length(vars)

            length(vars) == 0 ? sep_var = 1 : sep_var = length(vars)
            const separator = "="^(9 + (sep_var + 2)*16)

            # Initial values
            maxfn = 0
            varX = zeros(1,len)
            varXX = zeros(1,len)
            count = 1

            # Numeric iterations
            iter == true ? ti = time_ns() : nothing # To measure algorithm performance

            for i=1:maxiter
                  fn = f
                  # Numeric evaluations for simbolic variables
                  for j=1:len
                        varX[j] = lb[j] + (up[j] - lb[j])*rand()
                        fn = fn |> subs(vars[j], varX[1,j])
                  end

                  fn = convert(Float32, fn)
                  if fn > maxfn
                        maxfn = fn
                        varXX = varX
                  end

                  # Print iteration for the case of iter == true
                  if iter == true
                        if count == 1
                              # Header message and version
                              println(separator)
                              msg1 = @sprintf(" RandomSearch Solver for Optimization %s ", version)
                              msg2 = @sprintf(" Ph D Kelvyn B. Sánchez ")
                              header_left = ceil((length(separator) - length(msg1))/2)
                              header_right = (length(separator) - length(msg1)) - header_left
                              msg1_left = ":"^convert(Int32,header_left)
                              msg1_right = ":"^convert(Int32,header_right)
                              msg1 = string(msg1_left, msg1, msg1_right)
                              println(msg1)

                              header_left = ceil((length(separator) - length(msg2))/2)
                              header_right = (length(separator) - length(msg2)) - header_left
                              msg2_left = ":"^convert(Int32,header_left)
                              msg2_right = ":"^convert(Int32,header_right)
                              msg2 = string(msg2_left, msg2, msg2_right)
                              println(msg2)
                        end

                        if count == 1 || ceil((count/1)/10) == (count/1)/10
                              println(separator)

                              # Header
                              msg1 = "Iter\t "

                              for i=1:len
                                    msg1 = string(msg1, vars[i], "\t\t ")
                              end

                              msg1 = string(msg1, "f(x)\t\t")
                              println(msg1)
                              println(separator)
                        end

                        # Custom message for iterator
                        msg2 = @sprintf("%i\t ", count)
                        for i=1:len
                              msg2 = string(msg2, @sprintf("%4.6e\t ", varX[i]))
                        end

                        msg2 = string(msg2, @sprintf("%4.6e\t ", fn))
                        println(msg2)
                  end
                  count += 1

            end

            iter == true ? println(separator) : nothing
            iter == true ? tf = time_ns() : nothing # To measure algorithm performance
            iter == true ? println("elapse time: ", (tf - ti)/1.0e9, " seconds") : nothing

            return varXX, maxfn
      #catch Msg
      #      println(separator)
      #      println(Msg)
      #      return NaN, NaN, NaN
      #end
end


text:
using CFTime
using Dates
using Formatting
using JLD2
using Distributed
using SharedArrays
using MPI

if !(:LogSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "LogSystem.jl")))
end
using .LogSystem

if !(:ModelClockSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "ModelClockSystem.jl")))
end
using .ModelClockSystem

if !(:ConfigCheck in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "ConfigCheck.jl")))
end
using .ConfigCheck

if !(:appendLine in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "share", "AppendLine.jl")))
end

function runModel(
    OMMODULE      :: Any,
    coupler_funcs :: Any,
    t_start       :: AbstractCFDateTime,
    t_simulation  :: Any,
    read_restart  :: Bool,
    config        :: Dict,
)

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    ntask = MPI.Comm_size(comm)

    writeLog("===== [ Master Created ] =====")
    writeLog("Number of total tasks       : {:d}", ntask)
    writeLog("Number of total worker tasks: {:d}", ntask-1)

    MPI.Barrier(comm)

    is_master = rank == 0

    if is_master
        writeLog("Validate driver config.")
        config[:DRIVER] = validateConfigEntries(config[:DRIVER], getDriverConfigDescriptor())
    end

    writeLog("Setting working directory: {:s}", config[:DRIVER][:caserun])
    cd(config[:DRIVER][:caserun])

    if is_master

        global t_start, t_end

        if read_restart
            writeLog("read_restart is true.")
            restart_info = JLD2.load("model_restart.jld2")
            t_start = restart_info["timestamp"]
        end

        t_end = t_start + t_simulation

        writeLog("### Simulation time start: {:s}", Dates.format(t_start, "yyyy-mm-dd HH:MM:SS"))
        writeLog("### Simulation time end  : {:s}", Dates.format(t_end, "yyyy-mm-dd HH:MM:SS"))

    end

    writeLog("Passing t_start and config to slaves.")
    if ! is_master
        t_start = nothing
        config = nothing
    end
    t_start = MPI.bcast(t_start, 0, comm)
    config = MPI.bcast(config, 0, comm)

    # copy the start time
    beg_datetime = t_start + Dates.Second(0)

    writeLog("Read restart  : {}", read_restart)

    # Construct model clock
    clock = ModelClock("Model", beg_datetime)

    # initializing
    is_master && println("===== INITIALIZING MODEL: ", OMMODULE.name , " =====")

    OMDATA = OMMODULE.init(
        casename     = config[:DRIVER][:casename],
        clock        = clock,
        config      = config,
        read_restart = read_restart,
    )

    coupler_funcs.after_model_init!(OMMODULE, OMDATA)

    is_master && println("Ready to run the model.")
    step = 0
    while true

        step += 1

        writeLog("Current time: {:s}", clock2str(clock))

        stage, Δt, write_restart = coupler_funcs.before_model_run!(OMMODULE, OMDATA)

        if stage == :RUN

            cost = @elapsed let

                OMMODULE.run!(
                    OMDATA;
                    Δt = Δt,
                )

                MPI.Barrier(comm)

            end
            writeLog("Computation cost: {:f} secs.", cost)
            coupler_funcs.after_model_run!(OMMODULE, OMDATA)

            advanceClock!(clock, Δt)

        elseif stage == :END
            writeLog("stage == :END. Break loop now.")
            break
        end

    end

    OMMODULE.final(OMDATA)
    coupler_funcs.final(OMMODULE, OMDATA)

    if is_master

        writeLog("Writing restart file of driver")

        JLD2.save("model_restart.jld2", "timestamp", clock.time)
        archive_list_file = joinpath(
            config[:DRIVER][:caserun],
            config[:DRIVER][:archive_list],
        )

        timestamp_str = format(
            "{:s}-{:05d}",
            Dates.format(clock.time, "yyyy-mm-dd"),
            floor(Int64, Dates.hour(clock.time)*3600+Dates.minute(clock.time)*60+Dates.second(clock.time)),
        )

        appendLine(archive_list_file,
            format("cp,{},{},{}",
                "model_restart.jld2",
                config[:DRIVER][:caserun],
                joinpath(config[:DRIVER][:archive_root], "rest", timestamp_str)
            )
        )

    end

    if is_master
        archive(joinpath(
            config[:DRIVER][:caserun],
            config[:DRIVER][:archive_list],
        ))
    end

    is_master && println("Program Ends.")

end

function getDriverConfigDescriptor()

    return [
            ConfigEntry(
                :casename,
                :required,
                [String,],
            ),

            ConfigEntry(
                :caseroot,
                :required,
                [String,],
            ),

            ConfigEntry(
                :caserun,
                :required,
                [String,],
            ),

            ConfigEntry(
                :archive_root,
                :required,
                [String,],
            ),

            ConfigEntry(
                :archive_list,
                :optional,
                [String,],
                "archive_list.txt",
            ),
   ]
end

function archive(
    archive_list_file :: String,
)

    println("===== Archiving files BEGIN =====")

    for line in eachline(archive_list_file)

        args = split(line, ",")

        if length(args) == 0
            continue
        end

        action = args[1]
        args = args[2:end]

        if action in ["mv", "cp"]

            fname, src_dir, dst_dir = args

            if ! isdir(dst_dir)
                mkpath(dst_dir)
            end

            src_file = joinpath(src_dir, fname)
            dst_file = joinpath(dst_dir, fname)

            if isfile(src_file)

                if action == "mv"
                    mv(src_file, dst_file, force=true)
                    println(format("Moving file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                elseif action == "cp"
                    cp(src_file, dst_file, force=true)
                    println(format("Copying file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                end

            else
                println("File does not exist: ", src_file)
            end

        elseif action == "rm"
            fname, fdir = args
            rm(joinpath(fdir, fname), force=true)
            println(format("Removing file: {:s} in {:s}", fname, fdir))
        else
            throw(ErrorException(format("Unknown action in archive list: {:s}", action)))
        end

    end

    println("===== Archiving files END =====")

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
@testset "phyLiNC" begin

fastasimple = joinpath(@__DIR__, "..", "examples", "simple.aln")
fasta8sites = joinpath(@__DIR__, "..", "examples", "Ae_bicornis_8sites.aln") # 8 sites only
fastacontig = joinpath(@__DIR__, "..", "examples", "Ae_bicornis_Tr406_Contig10132.aln")
fasta1missing = joinpath(@__DIR__, "..", "examples", "simple_missingone.aln")
fastaindiv    = joinpath(@__DIR__, "..", "examples", "individuals.aln")
mappingfile   = joinpath(@__DIR__, "..", "examples", "mappingIndividuals.csv")
#= for local trials only, not for automatic testing
  pkgpath = dirname(Base.find_package("PhyloNetworks"))
  fastasimple   = abspath(joinpath(pkgpath, "..", "examples", "simple.aln"))
  fasta8sites   = abspath(joinpath(pkgpath, "..", "examples", "Ae_bicornis_8sites.aln"))
  fastacontig   = abspath(joinpath(pkgpath, "..", "examples", "Ae_bicornis_Tr406_Contig10132.aln"))
  fasta1missing = abspath(joinpath(pkgpath, "..", "examples", "simple_missingone.aln"))
  fastaindiv    = abspath(joinpath(pkgpath, "..", "examples", "individuals.aln"))
  mappingfile   = abspath(joinpath(pkgpath, "..", "examples", "mappingIndividuals.csv"))
=#
emptyconstraint = PhyloNetworks.TopologyConstraint[]

@testset "optimize local BL & gammas, simple example" begin
Random.seed!(99)
net_simple = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net_simple, fastasimple, :JC69)

## Local BL
obj.loglik = -Inf64
obj.net.edge[2].length= 0.0 # unzipping
e = PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[6],
        PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 10))
@test length(e) == 4     # not 5: edge below hybrid was excluded
@test obj.loglik > -40.0 # starting likelihood = -43.95468386633092

# ## Local Gamma
γcache = PhyloNetworks.CacheGammaLiNC(obj)
hybridmajorparent = PhyloNetworks.getMajorParentEdge(obj.net.hybrid[1])
@test_nowarn PhyloNetworks.optimizelocalgammas_LiNC!(obj, hybridmajorparent, 1e-6, γcache)
@test hybridmajorparent.gamma != 0.9
@test PhyloNetworks.getMinorParentEdge(obj.net.hybrid[1]).gamma != 0.1
end

@testset "optimize local BL & gammas, complex network and 8 sites" begin
Random.seed!(98)
dna_dat, dna_weights = readfastatodna(fasta8sites, true); # 22 species, 3 hybrid nodes, 103 edges
net = readTopology("((((((((((((((Ae_caudata_Tr275,Ae_caudata_Tr276),Ae_caudata_Tr139))#H1,#H2),(((Ae_umbellulata_Tr266,Ae_umbellulata_Tr257),Ae_umbellulata_Tr268),#H1)),((Ae_comosa_Tr271,Ae_comosa_Tr272),(((Ae_uniaristata_Tr403,Ae_uniaristata_Tr357),Ae_uniaristata_Tr402),Ae_uniaristata_Tr404))),(((Ae_tauschii_Tr352,Ae_tauschii_Tr351),(Ae_tauschii_Tr180,Ae_tauschii_Tr125)),(((((((Ae_longissima_Tr241,Ae_longissima_Tr242),Ae_longissima_Tr355),(Ae_sharonensis_Tr265,Ae_sharonensis_Tr264)),((Ae_bicornis_Tr408,Ae_bicornis_Tr407),Ae_bicornis_Tr406)),((Ae_searsii_Tr164,Ae_searsii_Tr165),Ae_searsii_Tr161)))#H2,#H4))),(((T_boeoticum_TS8,(T_boeoticum_TS10,T_boeoticum_TS3)),T_boeoticum_TS4),((T_urartu_Tr315,T_urartu_Tr232),(T_urartu_Tr317,T_urartu_Tr309)))),(((((Ae_speltoides_Tr320,Ae_speltoides_Tr323),Ae_speltoides_Tr223),Ae_speltoides_Tr251))H3,((((Ae_mutica_Tr237,Ae_mutica_Tr329),Ae_mutica_Tr244),Ae_mutica_Tr332))#H4))),Ta_caputMedusae_TB2),S_vavilovii_Tr279),Er_bonaepartis_TB1),H_vulgare_HVens23);");
PhyloNetworks.fuseedgesat!(93, net)
obj = (@test_logs (:warn, r"taxa with no data") PhyloNetworks.StatisticalSubstitutionModel(net, fasta8sites, :JC69))
@test length(obj.net.leaf) == 22
io = IOBuffer();
PhyloNetworks.showdata(io, obj)
@test String(take!(io)) == "data:\n  22 species\n  8 sites"
PhyloNetworks.showdata(io, obj, true)
@test String(take!(io)) ==
"data:
  22 species
  8 sites
  0 sites with no data (0.0%)
  2 invariant sites (25.0%)
  6 sites with 2 distinct states (75.0%)
  6 parsimony-informative sites (75.0%)
  6 sites with 1 or more missing values (75.0%)
  3.41% missing values overall"
close(io)

preorder!(obj.net)
PhyloNetworks.checknetwork_LiNC!(obj.net, 3, true, true, emptyconstraint)
# checknetwork removes degree-2 nodes (including root) and 2- and 3-cycles
# and requires that the network is preordered.
PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
PhyloNetworks.unzip_canonical!(obj.net)
## Local BL
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5)
obj.loglik = +Inf64
PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[27], lcache)
obj.loglik = -Inf64
@test_nowarn PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[27], lcache)
@test obj.net.edge[27].length ≈ 1e-8 atol=.0001
@test obj.net.edge[32].length ≈ 0.032 atol=.0005
@test obj.loglik > -80.4
# Local BL constrained edge
lengthe = obj.net.edge[44].length
@test_nowarn PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[44], lcache)
@test obj.net.edge[44].length == 0.0
@test obj.loglik > -78.95

# ## Local Gamma
# edge[4] = major parent edge of hybrid[1]
γcache = PhyloNetworks.CacheGammaLiNC(obj)
@test_nowarn PhyloNetworks.optimizelocalgammas_LiNC!(obj, obj.net.edge[4], 1e-6,γcache)
@test obj.net.edge[4].gamma == 0.0
@test PhyloNetworks.getMajorParentEdge(obj.net.hybrid[1]).gamma == 1.0

# gamma at a hybrid ladder: when some displayed trees don't have the focus edge
# 2 unzipped reticulations in a hybrid ladder, reasonable (small) branch lengths
net = readTopology("(#H2:0.0001::0.0244,((C:0.0262,((B:0.0)#H1:0.0::0.6)#H2:0.03::0.9756):0.4812,(#H1:0.0001::0.4,A:0.1274):0.0001):0.0274,D:0.151);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
obj.trait[1][3] = 3 # to create discordance across sites
obj.trait[4][4] = 2 #    for estimated γ to be within (0,1)
γcache = PhyloNetworks.CacheGammaLiNC(obj)
PhyloNetworks.discrete_corelikelihood!(obj) # -29.82754754416619
ll = PhyloNetworks.optimizegamma_LiNC!(obj, obj.net.edge[4], .001, γcache, 3)
@test ll ≈ -29.655467369763585
@test obj.net.edge[4].gamma ≈ 0.8090635871910823
@test γcache.hase == [true, true, false]
@test all(γcache.clikn .== 0.0)
ll = PhyloNetworks.optimizegamma_LiNC!(obj, obj.net.edge[1], .001, γcache, 3)
@test ll ≈ -29.468517891012983
@test obj.net.edge[1].gamma ≈ 0.19975558937688737
@test γcache.hase[[1,2]] == [false,true]
@test γcache.hase[3] === missing
@test γcache.clikn[[6,7]] ≈ [0.25, 0.0007] rtol=0.01
end #of local branch length and gamma optimization with localgamma! localBL! with 8 sites

@testset "global branch length and gamma optimization" begin
# to run locally on complex network:
# net = readTopology("(H_vulgare_HVens23:0.5,(((Ae_speltoides_Tr251:0.5):0.5,(Ae_mutica_Tr237:0.0)#H4:1.0::0.7):0.5,((((((Ae_caudata_Tr139:0.5,Ae_caudata_Tr275:0.5):0.0)#H1:1.0::0.7,#H2:1.0::0.3):0.5,#H1:1.0::0.3):0.5,((Ae_comosa_Tr271:0.5,Ae_comosa_Tr272:0.5):0.5,((Ae_uniaristata_Tr403:0.5,Ae_uniaristata_Tr357:0.5):0.5,Ae_uniaristata_Tr402:0.5):0.5):0.5):0.5,(((Ae_tauschii_Tr352:0.5,Ae_tauschii_Tr351:0.5):0.5,Ae_tauschii_Tr125:0.5):0.5,(((((((Ae_longissima_Tr241:0.5,Ae_longissima_Tr242:0.5):0.5,Ae_longissima_Tr355:0.5):0.5,Ae_sharonensis_Tr265:0.5):0.5,((Ae_bicornis_Tr408:0.5,Ae_bicornis_Tr407:0.5):0.5,Ae_bicornis_Tr406:0.5):0.5):0.5,(Ae_searsii_Tr164:0.5,Ae_searsii_Tr165:0.5):0.5):0.0)#H2:1.0::0.7,#H4:1.0::0.3):0.5):0.5):0.5):0.5);");
# obj = PhyloNetworks.StatisticalSubstitutionModel(net, fasta8sites, :JC69);
net = readTopology("(((A:0.5,(B:0.0)#H1:1.0::0.9):0.5,(C:0.5,#H1:1.0::0.1):0.5):0.5,D:0.5);")
# branch lengths set to 0.5, then unzipped -> some BL are 0, some 1, most 0.5
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69);

## optimizeBL
Random.seed!(5);
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5)
obj.loglik = -Inf64
@test_nowarn PhyloNetworks.optimizealllengths_LiNC!(obj, lcache);
@test all(e.length != 1.0 for e in obj.net.edge)
@test [e.length for e in obj.net.edge] ≈ [0.3727, 0.0, 1.0e-8, 1.0e-8, 1.0e-8, 1.0e-8, 1.0e-8, 0.5201, 1.0e-8] rtol=.001

## optimizegammas -- and delete hybrid edges with γ=0
γcache = PhyloNetworks.CacheGammaLiNC(obj)
@test_nowarn PhyloNetworks.optimizeallgammas_LiNC!(obj,1e-6,γcache,100)
@test obj.net.numHybrids == 0
end

@testset "data to SSM pruning: complex network" begin
net = readTopology("((((((((((((((Ae_caudata_Tr275,Ae_caudata_Tr276),Ae_caudata_Tr139))#H1,#H2),(((Ae_umbellulata_Tr266,Ae_umbellulata_Tr257),Ae_umbellulata_Tr268),#H1)),((Ae_comosa_Tr271,Ae_comosa_Tr272),(((Ae_uniaristata_Tr403,Ae_uniaristata_Tr357),Ae_uniaristata_Tr402),Ae_uniaristata_Tr404))),(((Ae_tauschii_Tr352,Ae_tauschii_Tr351),(Ae_tauschii_Tr180,Ae_tauschii_Tr125)),(((((((Ae_longissima_Tr241,Ae_longissima_Tr242),Ae_longissima_Tr355),(Ae_sharonensis_Tr265,Ae_sharonensis_Tr264)),((Ae_bicornis_Tr408,Ae_bicornis_Tr407),Ae_bicornis_Tr406)),((Ae_searsii_Tr164,Ae_searsii_Tr165),Ae_searsii_Tr161)))#H2,#H4))),(((T_boeoticum_TS8,(T_boeoticum_TS10,T_boeoticum_TS3)),T_boeoticum_TS4),((T_urartu_Tr315,T_urartu_Tr232),(T_urartu_Tr317,T_urartu_Tr309)))),(((((Ae_speltoides_Tr320,Ae_speltoides_Tr323),Ae_speltoides_Tr223),Ae_speltoides_Tr251))H3,((((Ae_mutica_Tr237,Ae_mutica_Tr329),Ae_mutica_Tr244),Ae_mutica_Tr332))#H4))),Ta_caputMedusae_TB2),S_vavilovii_Tr279),Er_bonaepartis_TB1),H_vulgare_HVens23);");
PhyloNetworks.fuseedgesat!(93, net)
for edge in net.edge # reset network
    setLength!(edge,1.0)
end
for h in net.hybrid
    setGamma!(PhyloNetworks.getMajorParentEdge(h),0.6)
end
obj = (@test_logs (:warn, r"pruned") PhyloNetworks.StatisticalSubstitutionModel(net, fasta8sites, :JC69))
@test length(obj.net.leaf) == 22
@test length(obj.net.edge) == 52
@test length(obj.net.hybrid) == 3
@test !PhyloNetworks.hashybridladder(obj.net)
end

@testset "checknetwork LiNC" begin
tree = readTopology("(A:3.0,(B:2.0,(C:1.0,D:1.0):1.0):1.0);");
@test any(length(n.edge) == 2 for n in tree.node) # one node of degree 2
preorder!(tree)
PhyloNetworks.checknetwork_LiNC!(tree, 1, true, true)
@test all(length(n.edge) != 2 for n in tree.node) # no nodes of degree 2
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
@test_throws ErrorException PhyloNetworks.checknetwork_LiNC!(net, 0, true, true)

netstr = "(#H2:0.1::0.2,((C:0.2,((B:0.3)#H1:0.4)#H2:0.5::0.8):0.6,(#H1:0.7,((A1:0.8)#H3:0.01,(A2:0.9,#H3:0.02):0.03):1.0):1.1):1.2,O:1.3);"
net = readTopology(netstr)
# 2 reticulation in a hybrid ladder, and another isolated reticulation
undoinfo = PhyloNetworks.unzip_canonical!(net)
@test all(PhyloNetworks.getChildEdge(h).length == 0.0 for h in net.hybrid) # unzipped
@test writeTopology(net, round=true) == "(#H2:0.8::0.2,((C:0.2,((B:0.0)#H1:0.0)#H2:1.2::0.8):0.6,(#H1:1.0,((A1:0.0)#H3:0.81,(A2:0.9,#H3:0.82):0.03):1.0):1.1):1.2,O:1.3);"
PhyloNetworks.rezip_canonical!(undoinfo...)
@test writeTopology(net, round=true) == netstr
end

@testset "update root in SSM displayed trees" begin
# W structure
net = readTopology("(C:0.0262,(B:0.0)#H2:0.03::0.9756,(((D:0.1,A:0.1274):0.0)#H1:0.0::0.6,(#H2:0.0001::0.0244,#H1:0.151::0.4):0.0274):0.4812);")
# "((C:0.0262,(B:0.0)#H2:0.03::0.9756):0.4812,((D:0.1,A:0.1274):0.0)#H1:0.0::0.6,(#H2:0.0001::0.0244,#H1:0.151::0.4):0.0274);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
@test [t.node[t.root].number for t in obj.displayedtree] == [5,5,5,5]
# obj.displayedtree[1]: (C:0.026,(B:0.0):0.03,(((D:0.1,A:0.127):0.0):0.0):0.481);
# move the root to place the W structure at the root:
# the network's root node will be missing from some displayed trees.
rootatnode!(obj.net, 7) # node 7 = tree node whose 2 children are both hybrids
# "(#H2:0.0001::0.0244,((C:0.0262,(B:0.0)#H2:0.03::0.9756):0.4812,((D:0.1,A:0.1274):0.0)#H1:0.0::0.6):0.0274,#H1:0.151);"
PhyloNetworks.updateSSM_root!(obj) # re-root displayed trees in the same way
@test [t.node[t.root].number for t in obj.displayedtree] == [6,7,7,7]
@test writeTopology(obj.displayedtree[1]) == "(((D:0.1,A:0.1274):0.0)H1:0.0,(C:0.0262,(B:0.0)H2:0.03):0.4812);"
end

@testset "skip γ and lengths optimization when needed" begin
# W structure, with middle γs = 0
net = readTopology("(C:0.0262,(B:0.0)#H2:0.03::1,(((D:0.1,A:0.1274):0.0)#H1:0.004::1,(#H2:0.0001,#H1:0.151):0.0274):0.4812);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
for i in [8,9] setGamma!(obj.net.edge[i], 0.0); end
PhyloNetworks.updateSSM!(obj)
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5);
e = PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[10], lcache)
@test length(e) == 3 # not 5: the 2 edges with γ = 0 were excluded
@test obj.net.edge[10].length == 0.0274
# hybrid ladder, with lower γ = 0
net = readTopology("(#H2:0.0001::0.0244,((C:0.0262,((B:0.0)#H1:0.0::0.6)#H2:0.03::0.9756):0.4812,(#H1:0.0001::0.4,A:0.1274):0.0001):0.0274,D:0.151);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
setGamma!(obj.net.edge[4], 0.0); PhyloNetworks.updateSSM!(obj)
e = PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[5], lcache)
@test length(e) == 4
@test [obj.net.edge[i].length for i in [1,5]] == [0.0001, 0.03]
γcache = PhyloNetworks.CacheGammaLiNC(obj);
ll = PhyloNetworks.optimizegamma_LiNC!(obj, obj.net.edge[5], .001, γcache, 3)
@test ll ≈ -31.124547305074074 # same as before optimization
end

@testset "optimizestructure with simple example" begin
#= network pre-processing: resulting net hard-coded below.
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69; maxhybrid=1)
PhyloNetworks.checknetwork_LiNC!(obj.net, 1, true, true)
PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
PhyloNetworks.unzip_canonical!(obj.net)
writeTopology(obj.net)
=#
net = readTopology("((A:0.3399824481995197,(B:0.0)#H1:0.08353360676474617::0.9):0.0001,(C:0.0001,#H1:0.048844990600034506::0.1):0.10871530327558311,D:0.33998306091744424);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69; maxhybrid=1)
PhyloNetworks.discrete_corelikelihood!(obj)
@test obj.loglik ≈ -29.7762035
maxmoves = 2
Random.seed!(96)
γcache = PhyloNetworks.CacheGammaLiNC(obj)
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5)
PhyloNetworks.optimizestructure!(obj, maxmoves, 1, true, true, 0,100,
                                emptyconstraint, 1e-6, γcache, lcache)
@test obj.loglik > -27.6 # -27.42 with RNG from julia 1.5

# allow hybrid ladders
Random.seed!(110)
PhyloNetworks.optimizestructure!(obj, maxmoves, 1, true, false, 0,100,
                                emptyconstraint, 1e-6, γcache, lcache)
@test obj.loglik > -27.5 # -27.42 with RNG from julia 1.5
end # of optimizestructure with simple example

@testset "phyLiNCone with simple net, no constraints" begin
no3cycle = true
net = readTopology("((A:0.3399824481995197,(B:0.0)#H1:0.08353360676474617::0.9):0.0001,(C:0.0001,#H1:0.048844990600034506::0.1):0.10871530327558311,D:0.33998306091744424);")
seed = 102
for nohybridladder in [true, false]
    #=
    net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
    obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
    PhyloNetworks.checknetwork_LiNC!(obj.net, 1, no3cycle, nohybridladder)
    PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
    PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
    PhyloNetworks.unzip_canonical!(obj.net)
    writeTopology(obj.net) # result hard-coded above. independent of nohybridladder
    =#
    obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
    obj.loglik = -Inf # missing otherwise, which would cause an error below
    γcache = PhyloNetworks.CacheGammaLiNC(obj)
    lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5)
    @test_nowarn PhyloNetworks.phyLiNCone!(obj, 1, no3cycle,
            nohybridladder, 3, 2, false, false,
            devnull, seed, 0.5, emptyconstraint,
            1e-2, 1e-2, 1e-2, 1e-2, 0.0, 25.0, 0.01,0.9,
            γcache, lcache)
    @test obj.loglik > -27.45
end
end

@testset "phyLiNC no constraints: HKY, rate variation" begin
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
obj = @test_nowarn PhyloNetworks.phyLiNC(net, fastasimple, :JC69, :G, 2; maxhybrid=2, # no missing BLs, so they're not re-estimated
                    no3cycle=true, nohybridladder=true, maxmoves=2,
                    nreject=1, nruns=1, filename="", verbose=false, seed=108)
@test obj.loglik > -27.5 # depends on RNG. -27.4 with julia 1.5, -27.27 earlier
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
obj = @test_nowarn PhyloNetworks.phyLiNC(net, fastasimple, :HKY85; maxhybrid=2,
                    no3cycle=true, nohybridladder=true, maxmoves=2, probST=1.0, # not enough moves to get back to a good topology
                    nreject=1, nruns=1, filename="phyLiNC2", verbose=false, seed=0)
@test obj.loglik > -24.21
@test read("phyLiNC2.err", String) == ""
@test startswith(read("phyLiNC2.log", String), "PhyLiNC network estimation starting")
rm("phyLiNC2.log")
rm("phyLiNC2.err")

net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
addprocs(1) # multiple cores
@everywhere using PhyloNetworks
#using Distributed; @everywhere begin; using Pkg; Pkg.activate("."); using PhyloNetworks; end
originalstdout = stdout  # verbose=true below
redirect_stdout(devnull)
obj = PhyloNetworks.phyLiNC(net, fastasimple, :JC69; maxhybrid=2, no3cycle=true,
                        nohybridladder=true, maxmoves=2, nreject=1, nruns=2,
                        filename="phyLiNCmult", verbose=true, seed=106)
redirect_stdout(originalstdout)
@test obj.loglik > -28.27
rmprocs(workers()) # remove extra processors
@test occursin("using 1 worker", read("phyLiNCmult.log", String))
rm("phyLiNCmult.log")

# phyLiNC w/ maxhybrid = 0
net_h0 = readTopology("(((A:2.0,B:1.0):1.5,C:0.6):0.5,D:2.0);");
obj = @test_nowarn PhyloNetworks.phyLiNC(net_h0, fastasimple, :JC69, :G, 2; maxhybrid=0,
                    no3cycle=true, nohybridladder=true, maxmoves=2,
                    nreject=1, nruns=1, filename="", verbose=false, seed=115)
@test obj.loglik > -27.4
end

@testset "phyLiNC with simple net and one constraint" begin
net_level1_s = readTopology("(((S8,S9),((((S1,S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));") # S1A S1B S1C go on leaf 1
# 3-cycle at degree-2 root -> 2-cycle after root deletion, removed within LiNC
# constraint
net_level1_i, c_species = PhyloNetworks.mapindividuals(net_level1_s, mappingfile)
PhyloNetworks.resetNodeNumbers!(net_level1_i)
net_level1_i.node[22].number = 100
PhyloNetworks.updateconstraints!(c_species, net_level1_i)
@test c_species[1].taxonnums == Set([8,9,100])
@test c_species[1].node.number == 21
@test PhyloNetworks.getParent(net_level1_i.node[22].edge[1]).number == 21

obj = PhyloNetworks.StatisticalSubstitutionModel(net_level1_i,fastaindiv,:JC69,:GI,2)
# obj.net = deepcopy of input net, so we need to rebuild the constraints
c_species[1] = PhyloNetworks.TopologyConstraint(0x01, c_species[1].taxonnames, obj.net)
# obj.net = deepcopy of input net, so we need to rebuild the constraints if done after
@test_logs (:warn, r"no 3-cycle") match_mode=:any PhyloNetworks.checknetwork_LiNC!(obj.net, 2,
                                                    true, true, c_species, true)
PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)

#= assign good starting branch lengths: find them outside of tests
for e in obj.net.edge e.length = 0.1; end # was -1.0 for missing
PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
PhyloNetworks.unzip_canonical!(obj.net)
print([e.length for e in obj.net.edge])
=#
startingbl = [0.07817851911808402, 0.16346448006948466, 0.007466318403444228, 0.34338840607739174, 0.17023758685593676, 0.03060241610507292, 0.0, 0.5775895221009606, 0.05156530617279971, 0.33292850647444944, 0.13056644275564092, 0.27314259657043044, 0.07441709795812174, 0.0001, 0.0909334086743112, 0.017615628089233237, 0.13268439288408584, 0.14375488372292924, 0.3199334263745607]
for (i,e) in enumerate(obj.net.edge)
  e.length = startingbl[i]
end
PhyloNetworks.setalpha!(obj.ratemodel, 0.48438)
obj.loglik = -Inf # actual likelihood -56.3068141288164. Need something non-missing
seed = 103
γcache = PhyloNetworks.CacheGammaLiNC(obj)
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-2,1e-2,1e-2,1e-2, 5)
@test_nowarn PhyloNetworks.phyLiNCone!(obj, 2, true, true,
        3, 2, false, false, devnull,
        seed, 0.5, c_species, 1e-2, 1e-2,
        1e-2, 1e-2, 0.0,50.0, 0.01,.9, γcache, lcache)
@test obj.loglik > -65.4 # -65.0 with RNG from julia 1.5

obj = PhyloNetworks.phyLiNC(net_level1_s, # missing BLs, so BLs are re-estimated before starting
            fastaindiv, :JC69, :Inv; maxhybrid=2, no3cycle=true, nohybridladder=true,
            verbose=false, filename="", speciesfile=mappingfile, seed=138, nruns=1,
            maxmoves=10, nreject=2)
@test obj.loglik > -67.7 # -69.83824 with :noRV
@test 0.19 < obj.ratemodel.pinv[1] < 0.3
# test that species stayed together after optimization, as the only polytomy
function polytomyS1(node)
    length(node.edge) > 3 || return false
    return Set(n.name for n in PhyloNetworks.getChildren(node)) == Set(["S1A", "S1B", "S1C"])
end
@test sum(polytomyS1(nod) for nod in obj.net.node) == 1
end

# hybrid flip tests
@testset "hybrid flip basics" begin
no3cycle = true
# nohybridladder = true w/ simple 1 hybrid starting network
Random.seed!(123)
# below: net already checked for PhyLiNC, startingBL and unzipped
net = readTopology("((A:0.3399824481995197,(B:0.0)#H1:0.08353360676474617::0.9):0.0001,(C:0.0001,#H1:0.048844990600034506::0.1):0.10871530327558311,D:0.33998306091744424);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
obj.loglik = -Inf # loglik missing otherwise, which would cause an error below
γcache = PhyloNetworks.CacheGammaLiNC(obj)
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5)
@test PhyloNetworks.fliphybridedgeLiNC!(obj, obj.loglik, true, emptyconstraint, 1e-6, γcache, lcache)
@test obj.loglik ≈ -27.7273003 atol=.0001
@test !PhyloNetworks.deletehybridedgeLiNC!(obj, obj.loglik, no3cycle, emptyconstraint, γcache, lcache)
@test obj.loglik ≈ -27.7273003 atol=.0001

# nohybridladder = false w/ hybrid ladder starting network
Random.seed!(2)
# below: net already checked for PhyLiNC and unzipped
net = readTopology("(#H2:0.02495259889870113::0.0244,((C:1e-4,((B:0.0)#H1:0.0::0.6)#H2:0.034190897863530335::0.9756):0.24434924848805456,(#H1:0.01539513240840275::0.4,A:0.2864250860992079):1.0e-8):1e-4,D:0.2716998373895161);")
obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69)
obj.loglik = -Inf # loglik missing otherwise, which would cause an error below
γcache = PhyloNetworks.CacheGammaLiNC(obj);
lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-6,1e-6,1e-2,1e-3, 5);
@test PhyloNetworks.fliphybridedgeLiNC!(obj, obj.loglik, false, emptyconstraint, 1e-6, γcache, lcache)
@test obj.loglik ≈ -29.36982 atol=.01
@test PhyloNetworks.deletehybridedgeLiNC!(obj, obj.loglik, no3cycle, emptyconstraint, γcache, lcache)
@test obj.loglik ≈ -28.30294 atol=.01
@test !PhyloNetworks.fliphybridedgeLiNC!(obj, obj.loglik, false, emptyconstraint, 1e-6, γcache, lcache)
end # hybrid flip basics
end # of overall phyLiNC test set


text:
abstract type AbstractAugmentedLattice{T,N_ARR,N_CDT,L} <: AbstractLattice{T,N_ARR,N_CDT} end
abstract type AbstractEmbeddedLattice{T,N_ARR,N_CDT,L} <: AbstractAugmentedLattice{T,N_ARR,N_CDT,L} end

function coordinates(lattice::AbstractEmbeddedLattice)
    lattice.coordinates
end
function size(lattice::AbstractEmbeddedLattice)
    size(lattice.lattice)
end

function difference(aug_lattice::AbstractEmbeddedLattice{T,N_ARR,N_CDT,L},
                    edge::Tuple{PT,PT}) where {T,N_ARR,N_CDT,L_N_CDT,
                                               L<:AbstractLattice{T,N_ARR,L_N_CDT},
                                               PT<:NTuple{N_CDT,T}
                                               }
    edge_first_dims = (edge[1][1:L_N_CDT], edge[2][1:L_N_CDT])
    edge_trailing_dims = (edge[1][L_N_CDT+1:end], edge[2][L_N_CDT+1:end])
    return (difference(aug_lattice.lattice, edge_first_dims)...,
        difference(aug_lattice.embedded_lattice, edge_trailing_dims)...)
end

function Base.step(aug_lattice::AbstractEmbeddedLattice)
    (step(aug_lattice.lattice)..., step(aug_lattice.embedded_lattice)...)
end

function unembed_values(lattice::AbstractEmbeddedLattice{T,N_ARR,N_CDT}, values::AbstractArray{T,N_ARR}) where {T,N_ARR,N_CDT}
    inner_coords = [coord[N_ARR+1:end] for coord in coordinates(lattice)]
    return [values[findall(map((x) -> all(isapprox.(embedded_coord, x)), inner_coords))] for embedded_coord in coordinates(lattice.embedded_lattice)]
end

linear_next(num::Int) = num + 1
#@recipe function f(lattice::AbstractEmbeddedLattice, values; layout=nothing, subplot=nothing)
#    @series begin
#        subplot := subplot
#        (lattice.lattice, values)
#    end
#    @series begin
#        subplot := linear_next(subplot)
#        (lattice.embedded_lattice, unembed_values(lattice, values))
#    end
#end

struct RandomlyEmbeddedLattice{T,N_ARR,N_CDT,L<:AbstractLattice{T,N_ARR},E<:AbstractSpace{T}} <: AbstractEmbeddedLattice{T,N_ARR,N_CDT,L}
    lattice::L
    embedded_lattice::E
    coordinates::Array{NTuple{N_CDT,T},N_ARR}
end
function RandomlyEmbeddedLattice(; lattice::L, embedded_lattice::E) where {T,N_ARR,L<:AbstractLattice{T,N_ARR},E<:AbstractSpace{T}}
    embedded_coordinates = embed_randomly(lattice, embedded_lattice)
    RandomlyEmbeddedLattice(lattice, embedded_lattice, embedded_coordinates)
end
function embed_randomly(lattice, embedded_lattice)
    [(lattice_coord..., sample(embedded_lattice)...) for lattice_coord in coordinates(lattice)]
end
function sample(lattice::AbstractLattice)
    (rand(length(extent(lattice))...) .* extent(lattice)) .- (extent(lattice) ./ 2)
end


text:
splat(f) = args -> f(args...)

struct Calibration{T}
    video::SystemPath
    extrinsic::Float64
    intrinsic::T
    checker_size::Float64
end

const Intrinsic = Pair{Float64, Float64}

struct ExtrinsicCalibration
    tform
    itform
    ϵ::Vector{Float64}
end

struct BothCalibration
    matlab
    ϵ::Vector{Float64}
end

extract(::Missing, _, path) = missing
function extract(extrinsic::Float64, video, path)
    to = joinpath(path, "extrinsic.png")
    ffmpeg_exe(`-loglevel 8 -ss $extrinsic -i $video -vf format=gray,yadif=1,scale=sar"*"iw:ih -pix_fmt gray -vframes 1 $to`)
    to
end
function extract(intrinsic::Intrinsic, video, path)
    ss, t2 = intrinsic
    t = t2 - ss
    r = 25/t
    files = joinpath(path, "intrinsic%03d.png")
    ffmpeg_exe(`-loglevel 8 -ss $ss -i $video -t $t -r $r -vf format=gray,yadif=1,scale=sar"*"iw:ih -pix_fmt gray $files`)
    readdir(path, join = true)
end

@memoize function build_calibration(c)
    # @debug "building calibration" c
    mktempdir() do path
        # path = mktempdir(cleanup = false)
        extrinsic = extract(c.extrinsic, c.video, path)
        intrinsic = extract(c.intrinsic, c.video, path)
        # spawnmatlab(c.checker_size, extrinsic, intrinsic)
        buildcalibration(c.checker_size, extrinsic, intrinsic)
    end
end

calibrate!(poi, c) = map!(xy -> calibrate(c, xy), space(poi), space(poi))

function createAffineMap(poic, expected)
    X = vcat((poic[k]' for k in keys(expected))...)
    X = hcat(X, ones(3))
    Y = hcat(values(expected)...)'
    c = (X \ Y)'
    A = c[:, 1:2]
    b = c[:, 3]
    AffineMap(SMatrix{2,2,Float64}(A), SVector{2, Float64}(b))
end

function build_extra_calibration(c, e)
    k = filter_collinearity(e)
    deleteat!(c, k)
    deleteat!(e, k)
    npoints = length(e)
    tform = if npoints < 2
        IdentityTransformation()
    elseif npoints < 3
        s = norm(e[2] - e[1])/norm(c[2] - c[1])
        # LinearMap(SDiagonal(s, s))
        LinearMap(UniformScaling(s))
    else
        createAffineMap(c[1:3], e[1:3])
    end
    itform = inv(tform)
    (; tform, itform)
end

function find_collinear(xy)
    inds = ((1, 2), (1, 3), (2, 3))
    Δ = [norm(xy[i1] - xy[i2]) for (i1, i2) in inds]
    M, i = findmax(Δ)
    j, l = setdiff(1:3, i)
    Δ[j] + Δ[l] ≈ M && return only(setdiff(1:3, inds[i]))
    nothing
end

updatek!(k, j, ::Nothing) = nothing
updatek!(k, j, i) = push!(k, j[i])

function filter_collinearity(xy)
    k = Int[]
    for j in combinations(1:length(xy), 3)
        if !any(∈(k), j)
            i = find_collinear(xy[j])
            updatek!(k, j, i)
        end
    end
    k
end


text:
import Pkg
Pkg.activate(@__DIR__)

using BenchmarkTools
using PlanarConvexHulls
using StaticArrays
using Random

const T = Float64

suite = BenchmarkGroup()

for n in [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000]
    rng = MersenneTwister(n)
    hull = ConvexHull{CCW, T}()
    sizehint!(hull, n)
    suite["jarvis_march! $n"] = @benchmarkable jarvis_march!($hull, points) setup = begin
        points = [rand($rng, SVector{2, T}) for _ = 1 : $n]
    end
end

overhead = BenchmarkTools.estimate_overhead()
results = run(suite, verbose=true, overhead=overhead, gctrial=false)
for result in results
    println("$(first(result)):")
    display(last(result))
    println()
end


text:
import Interpolations: AbstractExtrapolation

"""
A `TransformedArray` is an `AbstractArray` type with an affine
coordinate shift: `A[i,j]` evaluates the "parent" array `P` used to
construct `A` at a location `x,y` given by `[x, y] = tfm*[i,j]`.  It
therefore allows you to use lazy-evaluation to perform affine
coordinate transformations.

```
A = TransformedArray(etp, tfm)
```
where `etp` is an extrapolation object as defined by the
Interpolations package, and `tfm` is an `AffineTransform`.
"""
immutable TransformedArray{T,N,E<:AbstractExtrapolation,Tf<:AffineTransform} <: AbstractArray{T,N}
    data::E
    tform::Tf
end
TransformedArray{T,N}(etp::AbstractExtrapolation{T,N}, a::AffineTransform) =
    TransformedArray{T,N,typeof(etp),typeof(a)}(etp, a)

function TransformedArray(A::AbstractInterpolation, a::AffineTransform)
    etp = extrapolate(A, NaN)
    TransformedArray(etp, a)
end

function TransformedArray(A::AbstractArray, a::AffineTransform)
    itp = interpolate(A, BSpline(Linear()), OnGrid())
    TransformedArray(itp, a)
end

Base.size(A::TransformedArray) = size(A.data)

function Base.getindex{T}(A::TransformedArray{T,2}, i::Number, j::Number)
    x, y = tformfwd(A.tform, i, j)
    A.data[x, y]
end

function Base.getindex{T}(A::TransformedArray{T,3}, i::Number, j::Number, k::Number)
    x, y, z = tformfwd(A.tform, i, j, k)
    A.data[x, y, z]
end

Base.similar{T}(A::TransformedArray, ::Type{T}, dims::Dims) = Array{T}(dims)

"""
`transform(A, tfm; origin_dest=center(A), origin_src=center(A)`
computes the transformed `A` over its entire domain.  By default the
transformation is assumed to operate around the center of the input
array, and output coordinates are referenced relative to the center of
the output.

If `A` is a TransformedArray, then the syntax is just `transform(A;
origin_dest=center(A), origin_src=center(A))`. This is different from
the behavior of `A[:,:]`, which assumes the origin of coordinates to
be all-zeros.  To obtain behavior equivalent to `getindex`, supply
zero-vectors for both of them. Alternatively to make `getindex` behave
as `transform`, offset the origin of the transform used to construct
`A` by
```
origin_src - tform.scalefwd*origin_dest
```
"""
function transform{_,N}(A::TransformedArray{_,N}; kwargs...)
    y = A.tform * ones(Int, N)
    yt = (y...)::NTuple{N,eltype(y)}
    a = A.data[yt...]
    dest = Array{typeof(a)}(size(A))
    transform!(dest, A; kwargs...)
    dest
end

transform(A, a::AffineTransform; kwargs...) = transform(TransformedArray(A, a); kwargs...)

"""
`transform!(dest, src, tfm; origin_dest=center(dest),
origin_src=center(src))` is like `transform`, but using a
pre-allocated output array `dest`.

If `src` is already a TransformedArray, use `transform!(dest, src;
kwargs...)`.
"""
function transform!{S,T,N}(dest::AbstractArray{S,N},
                           src::TransformedArray{T,N};
                           origin_dest = center(dest),
                           origin_src = center(src))
    tform = src.tform
    if tform.scalefwd == eye(T, N) && tform.offset == zeros(N) && size(dest) == size(src) && origin_dest == origin_src
        copy!(dest, src)
        return dest
    end
    offset = tform.offset - tform.scalefwd*origin_dest + origin_src
    _transform!(dest, src, offset)
end

transform!(dest, src, a::AffineTransform; kwargs...) = transform!(dest, TransformedArray(src, a); kwargs...)

@require ImageMetadata begin
    transform(A::ImageMetadata.ImageMeta, a::AffineTransform; kwargs...) = ImageMetadata.copyproperties(A, transform(ImageMetadata.data(A), a; kwargs...))
    transform!(dest, A::ImageMetadata.ImageMeta, a::AffineTransform; kwargs...) = ImageMetadata.copyproperties(A, transform!(dest, ImageMetadata.data(A), a; kwargs...))
end

# For a FilledExtrapolation, this is designed to (usually) avoid evaluating
# the interpolation unless it is in-bounds.  This often improves performance.
@generated function _transform!{S,T,N,E<:Interpolations.FilledExtrapolation}(
                                dest::AbstractArray{S,N},
                                src::TransformedArray{T,N,E},
                                offset)
    # Initialize the final column of s "matrix," e.g., s_1_3 = A_1_3 + o_3
    # s stands for source-coordinates. The first "column," s_d_1, corresponds
    # to the actual interpolation position. The later columns simply cache
    # previous computations.
    sN = ntuple(i->Expr(:(=), Symbol(string("s_",i,"_$N")), Expr(:call, :+, Symbol(string("A_",i,"_$N")), Symbol(string("o_",i)))), N)
    quote
        tform = src.tform
        data = src.data
        @nexprs $N d->(o_d = offset[d])
        @nexprs $N j->(@nexprs $N i->(A_i_j = tform.scalefwd[i,j]))
        fill!(dest, data.fillvalue)
        $(sN...)
        @nloops($N,i,d->(d>1 ? (1:size(dest,d)) : (imin:imax)),
                # The pre-expression chooses the range within each column that will be in-bounds
                d->(d >  2 ? (@nexprs $N e->s_e_{d-1}=A_e_{d-1}+s_e_d) :
                    d == 2 ? begin
                        imin = 1
                        imax = size(dest,1)
                        @nexprs $N e->((imin, imax) = irange(imin, imax, A_e_1, s_e_2, size(data, e)))
                        @nexprs $N e->(s_e_1=A_e_1*imin+s_e_d)
                    end :
                    nothing), # pre
            d->(@nexprs $N e->(s_e_d += A_e_d)), # post
            # Perform the interpolation of the source data
            @inbounds (@nref $N dest i) = (@nref $N data d->s_d_1)
        )
        dest
    end
end

center(A::AbstractArray) = [(size(A,d)+1)/2 for d = 1:ndims(A)]

# Find i such that
#      imin <= i <= imax
#      1 <= coef*i+offset <= upper
function irange(imin::Int, imax::Int, coef, offset, upper)
    const thresh = 10^4/typemax(Int)  # needed to avoid InexactError for results with abs() bigger than typemax
    if coef > thresh
        return max(imin, floor(Int, (1-offset)/coef)), min(imax, ceil(Int, (upper-offset)/coef))
    elseif coef < -thresh
        return max(imin, floor(Int, (upper-offset)/coef)), min(imax, ceil(Int, (1-offset)/coef))
    else
        if 1 <= offset <= upper
            return imin, imax
        else
            return 1, 0   # empty range
        end
    end
end


text:
using ITensors
using OMEinsum

const Cup = [0.0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]
const Cdn = [0.0 0 1 0; 0 0 0 1; 0 0 0 0; 0 0 0 0]
const Nup = [0.0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]
const Ndn = [0.0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]
abstract type HamiltonianModel end

@doc raw"
    hamiltonian(model<:HamiltonianModel)

return the hamiltonian of the `model` as a two-site tensor operator.
"
function hamiltonian end

@doc raw"
    Hubbard(t::Real,U::Real,μ::Real)

return a struct representing Hubbard model
"
struct Hubbard{T<:Real} <: HamiltonianModel
    t::T
	U::T
    μ::T
end

"""
	hamiltonian(model::Hubbard)
"""
function hamiltonian(model::Hubbard)
	t = model.t
	U = model.U
    μ = model.μ
    ampo = AutoMPO()
    sites = siteinds("Electron",2)
    ampo .+= -t, "Cdagup",1,"Cup",2
    ampo .+= -t, "Cdagup",2,"Cup",1
    ampo .+= -t, "Cdagdn",1,"Cdn",2
    ampo .+= -t, "Cdagdn",2,"Cdn",1

    if U ≠ 0
        ampo .+= 1/4*U, "Nupdn", 1
        ampo .+= 1/4*U, "Nupdn", 2
    end
    if μ ≠ 0
        ampo .+= -1/4*μ, "Nup", 1
        ampo .+= -1/4*μ, "Ndn", 1
        ampo .+= -1/4*μ, "Nup", 2
        ampo .+= -1/4*μ, "Ndn", 2
    end
    H = MPO(ampo,sites)

    H1 = Array(H[1],inds(H[1])...)
    H2 = Array(H[2],inds(H[2])...)
    h = reshape(ein"aij,apq->ipjq"(H1,H2),16,16)

    return h
end

@doc raw"
    THubbard(t::Real,U::Real,μ::Real)

return a struct representing THubbard model
"
struct THubbard{T<:Real} <: HamiltonianModel
    t::T
	U::T
    μ::T
end

"""
	hamiltonian(model::Hubbard)
"""
function hamiltonian(model::THubbard)
	t = model.t
	U = model.U
    μ = model.μ
    ampo = AutoMPO()
    sites = siteinds("Electron",2)
    ampo .+= -t, "Cdagup",1,"Cdn",2
    ampo .+= -t, "Cdagdn",2,"Cup",1
    ampo .+= t, "Cdagdn",1,"Cup",2
    ampo .+= t, "Cdagup",2,"Cdn",1

    if U ≠ 0
        ampo .+= 1/4*U, "Nupdn", 1
        ampo .+= 1/4*U, "Nupdn", 2
    end
    if μ ≠ 0
        ampo .+= -1/4*μ, "Nup", 1
        ampo .+= -1/4*μ, "Ndn", 1
        ampo .+= -1/4*μ, "Nup", 2
        ampo .+= -1/4*μ, "Ndn", 2
    end
    H = MPO(ampo,sites)

    H1 = Array(H[1],inds(H[1])...)
    H2 = Array(H[2],inds(H[2])...)
    h = reshape(ein"aij,apq->ipjq"(H1,H2),16,16)

    return h
end

@doc raw"
    hop_pair(t::Real,γ::Real)

return a struct representing hop_pair model
"
struct hop_pair{T<:Real} <: HamiltonianModel
    t::T
	γ::T
end

"""
	hamiltonian(model::hop_pair)
"""
function hamiltonian(model::hop_pair)
	t = model.t
    γ = model.γ
    ampo = AutoMPO()
    sites = siteinds("Electron",2)
    ampo .+= -t, "Cdagup",1,"Cup",2
    ampo .+= -t, "Cdagup",2,"Cup",1
    ampo .+= -t, "Cdagdn",1,"Cdn",2
    ampo .+= -t, "Cdagdn",2,"Cdn",1

    ampo .+= γ, "Cdagup",1,"Cdagdn",2
    ampo .+= -γ, "Cdagdn",1,"Cdagup",2
    ampo .+= γ, "Cdn",2,"Cup",1
    ampo .+= -γ, "Cup",2,"Cdn",1
    H = MPO(ampo,sites)

    H1 = Array(H[1],inds(H[1])...)
    H2 = Array(H[2],inds(H[2])...)
    h = reshape(ein"aij,apq->ipjq"(H1,H2),16,16)

    return h
end

function Hubbard_hand(model::Hubbard) where T <: Real
    t = model.t
	U = model.U
    μ = model.μ
    H = zeros(4,4,4,4)
    H .+= -t * (ein"ab,cd -> acbd"(Cup',Cup) .+ ein"ab,cd -> acbd"(Cup,Cup'))
    H .+= -t * (ein"ab,cd -> acbd"(Cdn',Cdn) .+ ein"ab,cd -> acbd"(Cdn,Cdn'))
    H .+= U/4 * (ein"ab,cd -> acbd"(Nup*Ndn,I(4)) .+ ein"ab,cd -> acbd"(I(4),Nup*Ndn))
    H .+= -μ/4 * (ein"ab,cd -> acbd"(Nup.+Ndn,I(4)) .+ ein"ab,cd -> acbd"(I(4),Nup.+Ndn))
    return reshape(H,16,16)
end

function hop_pair_hand(t::T,γ::T) where T <: Real
    H = zeros(4,4,4,4)
    H .+= -t * (ein"ab,cd -> acbd"(Cup',Cup) .+ ein"ab,cd -> acbd"(Cup,Cup'))
    H .+= -t * (ein"ab,cd -> acbd"(Cdn',Cdn) .+ ein"ab,cd -> acbd"(Cdn,Cdn'))
    H .+= γ * (ein"ab,cd -> acbd"(Cup',Cdn') .- ein"ab,cd -> acbd"(Cdn',Cup'))
    H .+= γ * (ein"ab,cd -> acbd"(Cup,Cdn) .- ein"ab,cd -> acbd"(Cdn,Cup))
    return reshape(H,16,16)
end


text:
include("euler/euler.jl")
using .Grid: max_path_sum

compute(triangle::String)::Int = max_path_sum([parse.(Int, split(line, " ")) for line ∈ split(triangle, "\n")])

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Minimum budget.

"""
    struct MinimumBudget{CI <: CombinatorialInstance, T <: Real}

Implement a minimum-budget constraint, i.e. a minimum quantity of a resource
that must be used.

``\\sum_i x_i \\times \\mathtt{weights}_i \\geq \\mathtt{min\\_budget}``
"""
struct MinimumBudget{CI <: CombinatorialInstance, T <: Real} <:
       CombinatorialVariation
    instance::CI
    weights::Union{Vector{T}, Dict{<:Any, T}}
    min_budget::T
    compute_all_values::Bool # Solutions will be output for all budget values,
    # usually with the same time complexity as a single solve.

    function MinimumBudget(instance::CI, weights::Vector{T}, min_budget::T, compute_all_values::Bool) where {CI <: CombinatorialInstance, T <: Real}
        # There should be as many weights as the inner instance has rewards.
        if length(weights) != dimension(instance)
            error("The MinimumBudget constraint has $(length(weights)) weights, but the combinatorial instance has dimension $(dimension(instance)).")
        end

        # Ensure that important properties of the algorithms are respected.
        if any(weights .< zero(T))
            error("Some weights are negative, which is not allowed.")
        end
        if min_budget < zero(T)
            error("The minimum budget is negative, which is not allowed.")
        end

        return new{CI, T}(instance, weights, min_budget, compute_all_values)
    end

    function MinimumBudget(instance::CI, weights::Dict{<:Any, T}, min_budget::T, compute_all_values::Bool) where {CI <: CombinatorialInstance, T <: Real}
        # There should be as many weights as the inner instance has rewards.
        if length(weights) != dimension(instance)
            error("The MinimumBudget constraint has $(length(weights)) weights, but the combinatorial instance has dimension $(dimension(instance)).")
        end

        # Ensure that important properties of the algorithms are respected.
        if any(values(weights) .< zero(T))
            error("Some weights are negative, which is not allowed.")
        end
        if min_budget < zero(T)
            error("The minimum budget is negative, which is not allowed.")
        end

        return new{CI, T}(instance, weights, min_budget, compute_all_values)
    end
end

function MinimumBudget(
    i::CI,
    weights::Union{Vector{T}, Dict{<:Any, T}},
    min_budget::T=zero(T);
    compute_all_values::Bool=false,
) where {CI, T}
    return MinimumBudget(i, weights, min_budget, compute_all_values)
end

function copy(
    i::MinimumBudget;
    instance::CI=i.instance,
    weights::Union{Vector{T}, Dict{<:Any, T}}=i.weights,
    min_budget::T=i.min_budget,
    compute_all_values::Bool=i.compute_all_values,
) where {CI <: CombinatorialInstance, T <: Real}
    return MinimumBudget(instance, weights, min_budget, compute_all_values)
end

"""
    abstract type MinBudgetedSolution

A type of solution that is specifically tailored for the minimum-budget
variant.
"""
abstract type MinBudgetedSolution <: CombinatorialSolution end

"""
    abstract type SingleMinBudgetedSolution

A type of `MinBudgetedSolution` when only one solution is available, to the
single value of the budget that is defined in the related instance. Objects
of this type do not have to implement more fields than any
`CombinatorialSolution`.
"""
abstract type SingleMinBudgetedSolution <: MinBudgetedSolution end

"""
    abstract type MultipleMinBudgetedSolution

A type of `MinBudgetedSolution` when several solutions are available, with
the several values of the budget that are specified by the combinatorial
instance (typically, from zero to a maximum value).

Objects of this type have to implement the usual fields of
`CombinatorialSolution`, but also `solutions`, a mapping from a budget value
to the corresponding solution (`solutions`[i]` corresponds to a minimum budget
of `i`).
"""
abstract type MultipleMinBudgetedSolution <: MinBudgetedSolution end

function reward(i::MinimumBudget, var)
    return reward(i.instance, var)
end

# TODO: weight interface for a solution (like value).

# # Maximum budget.

# """
#     struct MaximumBudget{CI <: CombinatorialInstance, T <: Real}

# Implement a maximum-budget constraint, i.e. a maximum quantity of a resource
# that can be used.

# ``\\sum_i x_i \\times \\mathtt{weights}_i \\leq \\mathtt{max\\_budget}``
# """
# struct MaximumBudget{CI <: CombinatorialInstance, T <: Real} <:
#        CombinatorialVariation
#     instance::CI
#     weights::Union{Vector{T}, Dict{<:Any, T}}
#     max_budget::T
#     compute_all_values::Bool
# end

# function MaximumBudget(
#     i::CI,
#     weights::Union{Vector{T}, Dict{<:Any, T}},
#     max_budget::T;
#     compute_all_values::Bool=false,
# ) where {CI, T}
#     return MaximumBudget(i, weights, max_budget, compute_all_values)
# end

# function copy(
#     i::MaximumBudget;
#     instance::CI=i.instance,
#     weights::Union{Vector{T}, Dict{<:Any, T}}=i.weights,
#     max_budget::T=i.max_budget,
#     compute_all_values::Bool=i.compute_all_values,
# ) where {CI <: CombinatorialInstance, T <: Real}
#     return MaximumBudget(ci, weights, max_budget, compute_all_values)
# end

# """
#     abstract type MaxBudgetedSolution

# A type of solution that is specifically tailored for the maximum-budget
# variant.
# """
# abstract type MaxBudgetedSolution <: CombinatorialSolution end

# """
#     abstract type SingleMaxBudgetedSolution

# A type of `MaxBudgetedSolution` when only one solution is available, to the
# single value of the budget that is defined in the related instance. Objects
# of this type do not have to implement more fields than any
# `CombinatorialSolution`.
# """
# abstract type SingleMaxBudgetedSolution <: MaxBudgetedSolution end

# """
#     abstract type MultipleMaxBudgetedSolution

# A type of `MaxBudgetedSolution` when several solutions are available, with
# the several values of the budget that are specified by the combinatorial
# instance (typically, from zero to a maximum value).

# Objects of this type have to implement the usual fields of
# `CombinatorialSolution`, but also `solutions`, a mapping from a budget value
# to the corresponding solution (`solutions`[i]` corresponds to a maximum budget
# of `i`).
# """
# abstract type MultipleMaxBudgetedSolution <: MaxBudgetedSolution end


text:
######################################################################
# Scene graph UTs
# -----
# Licensed under Apache License 2.0
using Test
using ExtraFun
using SceneGraphs
using StaticArrays

abstract type AbstractEntity end

const TestScene = Scene3D{AbstractEntity, Float64}

EntityTransform = SpriteTransform{AbstractEntity, Float64}

struct Entity <: AbstractEntity
    name::Symbol
    transform::EntityTransform
end
Entity(name::Symbol) = Entity(name, EntityTransform())

function transform_matrix(translation, rotation, scale)
    mt = @SMatrix([1 0 0 translation[1];
                   0 1 0 translation[2];
                   0 0 1 translation[3];
                   0 0 0 1])

    rcos = cos(deg2rad(rotation))
    rsin = sin(deg2rad(rotation))
    mr = @SMatrix([rcos -rsin 0 0;
                   rsin  rcos 0 0;
                   0     0    1 0;
                   0     0    0 1])

    ms = @SMatrix([scale[1] 0 0 0;
                   0 scale[2] 0 0;
                   0 0 1 0;
                   0 0 0 1])

    return mt*mr*ms
end

macro test_nothrow(expr)
    esc(:(@test begin $(expr); true end))
end

@testset "graph" begin
    let scene = TestScene(), entity1 = Entity(:first), entity2 = Entity(:second), entity3 = Entity(:third), entity4 = Entity(:fourth)
        push!(scene, entity1)
        push!(scene, entity2)
        @test entity1 ∈ scene && entity2 ∈ scene

        @test sceneof(entity1) === scene && sceneof(entity2) === scene && sceneof(entity3) === nothing

        @test parent!(entity4, entity3) === entity4
        @test parentof(entity4) === entity3 && sceneof(entity4) === nothing && sceneof(entity3) === sceneof(entity4)

        @test parent!(entity3, entity1) === entity3
        @test parentof(entity3) === entity1 && parentof(entity4) === entity3 && sceneof(entity3) === scene && sceneof(entity1) === sceneof(entity3) && sceneof(entity3) === sceneof(entity4)
        @test entity3 ∈ childrenof(entity1) && entity3 ∉ childrenof(entity2)

        @test parent!(entity3, entity2) === entity3
        @test parentof(entity3) === entity2 && parentof(entity4) === entity3 && sceneof(entity3) === scene && sceneof(entity2) === sceneof(entity3) && sceneof(entity3) === sceneof(entity4)
        @test entity3 ∉ childrenof(entity1) && entity3 ∈ childrenof(entity2)

        @test deparent!(entity3) === entity3
        @test parentof(entity3) === nothing && parentof(entity4) === entity3
        @test sceneof(entity3) === nothing && sceneof(entity4) === nothing
        @test entity3 ∉ childrenof(entity1) && entity3 ∉ childrenof(entity2)
    end
end

@testset "transformation" begin
    @testset "translate" begin
        let scene = TestScene(), parent = Entity(:parent), interim = Entity(:interim), child = Entity(:child)
            push!(scene, parent)

            parent!(child, interim)
            parent!(interim, parent)

            @test_nothrow translate!(parent, 1, 2, 3)
            @test_nothrow translate!(interim, 2, 3, 4)
            update!(transformof(parent))

            @test worldlocationof(parent)  ≈ [1, 2, 3]
            @test worldlocationof(interim) ≈ [3, 5, 7]
            @test worldlocationof(child)   ≈ [3, 5, 7]
        end
    end

    @testset "rotate" begin
        let scene = TestScene(), parent = Entity(:parent), interim = Entity(:interim), child = Entity(:child)
            push!(scene, parent)

            parent!(child, interim)
            parent!(interim, parent)

            @test_nothrow rotate!(parent, deg2rad(45))
            @test_nothrow rotate!(interim, deg2rad(30))
            update!(transformof(parent))

            @test worldrotationof(parent)  ≈ deg2rad(45)
            @test worldrotationof(interim) ≈ deg2rad(75)
            @test worldrotationof(child)   ≈ deg2rad(75)
        end
    end

    @testset "scale" begin
        let scene = TestScene(), parent = Entity(:parent), interim = Entity(:interim), child = Entity(:child)
            push!(scene, parent)

            parent!(child, interim)
            parent!(interim, parent)

            @test_nothrow scale!(parent, 1, 2)
            @test_nothrow scale!(interim, 4, 2)
            update!(transformof(parent))

            @test worldscaleof(parent)  ≈ [1, 2]
            @test worldscaleof(interim) ≈ [4, 4]
            @test worldscaleof(child)   ≈ [4, 4]
        end
    end

    @testset "combined" begin
        let scene = TestScene(), parent = Entity(:parent), interim = Entity(:interim), child = Entity(:child)
            push!(scene, parent)

            parent!(child, interim)
            parent!(interim, parent)

            transform!(parent,  translate=(1, 2, 3), rotate=deg2rad(45), scale=(1, 2))
            transform!(interim, translate=(2, 3, 4), rotate=deg2rad(30), scale=(4, 2))
            update!(transformof(parent))

            r45 = deg2rad(45)
            r30 = deg2rad(30)
            m_p = transform_matrix((1, 2, 3), 45, (1, 2))
            m_i = transform_matrix((2, 3, 4), 30, (4, 2))
            m_c = transform_matrix((0, 0, 0),  0, (1, 1))

            @test obj2world(parent)  ≈ m_p
            @test obj2world(interim) ≈ m_p*m_i
            @test obj2world(child)   ≈ m_p*m_i*m_c
        end
    end
end