text:
MaybeFunction = Nullable{Function}
const nullFunction = MaybeFunction()
typealias AkoFunction Union{Function, MaybeFunction}

Base.isnull(x::MaybeFunction) = x.isnull
Base.isnull(x::Function)      = false

abstract  Multivalued

abstract  Surface
abstract  SurfaceAction        <: Surface
abstract  Boundary
abstract  BoundaryMaterial     <: Boundary

abstract  MateralBounding      <: Boundary
abstract  Wood                 <: MaterialBounding
abstract  Clay                 <: MaterialBounding

abstract  BoundingContext      <: Boundary
abstract  Imprecise            <: BoundingContext
abstract  Precise              <: Imprecise
abstract  Inaccurate           <: BoundingContext
abstract  Accurate             <: Inaccurate

abstract  BoundingCoverness    <: Boundary
abstract  Spread               <: BoundingCoverness
abstract  Shared               <: Spread
abstract  Throughout           <: BoundingCoverness
abstract  Approached           <: Throughout

abstract  AnyBoundaries        <: Boundary

abstract  OneBoundary          <: AnyBoundaries
abstract  TwoBoundaries        <: AnyBoundaries
abstract  ThreeBoundaries      <: AnyBoundaries

abstract
AbsInterval{          Surface,       Tension                  } <: Number

                #   type covering   configures     perceptual
                #   value-context   boundaries     values map

abstract
AbsFlexibleInterval{  Surface,       Tension,        Action   } <: AbsInterval{ Surface, Tension }

immutable
FlexibleInterval{     Surface <: Real,
                      Tension <: AnyBoundaries,
                      Action  <: FunctionMaybe  } <: AbsFlexibleInterval{ Surface, Tension, Action }
    lo::Surface
    hi::Surface

    FlexibleInterval{S}(lo::S,hi::S) = (lo <= hi) ? new(lo,hi) : new(hi,lo)
end

FlexibleInterval{S<:Real,T<:AnyBoundaries,A<:FunctionMaybe}(boundaries::T, action::A, lo::S,hi::S) =
    FlexibleInterval{S,T,A}(lo,hi)

FlexibleInterval{S,T<:AnyBoundaries,A<:FunctionMaybe}(::Type{T}, ::Type{A}, lo::S, hi::S) = FlexibleInterval(lo,hi)
FlexibleInterval{S,T<:AnyBoundaries}(::Type{T}, lo::S, hi::S) = FlexibleInterval(T, nullFunction, lo, hi)
FlexibleInterval{S,A<:FunctionMaybe}(::Type{A}, lo::S, hi::S) = FlexibleInterval(ClosedCloseD, A, lo, hi)
FlexibleInterval{S}(lo::S, hi::S)   = FlexibleInterval(ClosedCloseD, nullFunction, lo, hi)


text:
## SymRange, an AbstractUnitRange that's symmetric around 0
# These are used as axes for CenterIndexedArrays
struct SymRange <: AbstractUnitRange{Int}
    n::Int  # goes from -n:n
end

function SymRange(r::AbstractUnitRange)
    first(r) == -last(r) || error("cannot convert $r to a SymRange")
    SymRange(last(r))
end

Base.first(r::SymRange) = -r.n
Base.last(r::SymRange) = r.n
Base.axes(r::SymRange) = (r,)

@inline Base.unsafe_indices(r::SymRange) = (r,)

function iterate(r::SymRange)
    r.n == 0 && return nothing
    first(r), first(r)
end

function iterate(r::SymRange, s)
    s == last(r) && return nothing
    copy(s+1), s+1
end

@inline function Base.getindex(v::SymRange, i::Int)
    @boundscheck abs(i) <= v.n || Base.throw_boundserror(v, i)
    return i
end

Base.intersect(r::SymRange, s::SymRange) = SymRange(min(last(r), last(s)))

@inline function Base.getindex(r::SymRange, s::SymRange)
    @boundscheck checkbounds(r, s)
    s
end

@inline function Base.getindex(r::SymRange, s::AbstractUnitRange{<:Integer})
    @boundscheck checkbounds(r, s)
    return s
end

# TODO: should we be worried about the mismatch in axes?
# And should `convert(SymRange, r)` fail if axes(r) isn't the same as the result?
Base.promote_rule(::Type{SymRange}, ::Type{UR}) where {UR<:AbstractUnitRange} =
    UR
Base.promote_rule(::Type{UnitRange{T2}}, ::Type{SymRange}) where {T2} =
    UnitRange{promote_type(T2, Int)}

Base.show(io::IO, r::SymRange) = print(io, "SymRange(", repr(last(r)), ')')

if isdefined(Base, :reduced_index)
    Base.reduced_index(r::SymRange) = SymRange(0)
end

ext: jl


text:
using CSV, DataFrames, Plots, Statistics
df = DataFrame(CSV.File("data/btc.csv"))
x = log.(df.cl[2:2001]) .- log.(df.cl[1:2000])
x .= x .- mean(x)

m1 = Std()
m2 = GJR()

res1 = fit!(m1, x)
res2 = fit!(m2, x)

v1 = predict(m1, 200, 240)
v2 = predict(m2, 200, 240)


text:
using MechanicalSketch
import MechanicalSketch: sethue, background, O, HE, color_with_lumin
import MechanicalSketch: finish, set_scale_sketch, arrow, °, m, kN, empty_figure
import MechanicalSketch: PALETTE

let
empty_figure(filename = joinpath(@__DIR__, "test_16.png"))
background(color_with_lumin(PALETTE[6], 80))
sethue(PALETTE[4])

function scaleindicator_force(α)
    p1 = O + 2m .* (cos(α), sin(α))
    arrow(p1, 1kN, α = α, labellength = true)
end

for i = 1:10
    sethue(PALETTE[i])
    set_scale_sketch(20kN,  HE * i)
    scaleindicator_force(36° * (i - 1))
end

finish()
# reset default:
set_scale_sketch()
nothing
end


text:
module alum_cyl_mode_examples
using FinEtools
using FinEtools.MeshExportModule
using Arpack

# Mesh alum_cyl.inp
# Abaqus with the Standard isoparametric C3D4 tetrahedron:
C3D4 = [1 0
2 0
3 6.63586E-005
4 0.000171053
5 0.000211299
6 0.000244378
7 2564.63
8 2568.09
9 2597.26
10 4094.38
11 4714.36
12 4717.19
13 5181.98
14 6865.13
15 6868.17
16 6962.86
17 6965.67
18 7024.97
19 7029.44
20 7108.54]
# Abaqus with the standard quadratic tetrahedron:
C3D10 = [1 0
2 0
3 0
4 0.000139365
5 0.000221551
6 0.000291805
7 2546.81
8 2546.81
9 2560.69
10 4100
11 4693.55
12 4693.56
13 5121.57
14 6841.21
15 6841.24
16 6914.22
17 6914.23
18 6950.64
19 6950.66
20 7000.64]

function alum_cyl_mode_nice_t4()

    E = 70000*phun("MPa");
    nu = 0.33;
    rho = 2700*phun("KG/M^3");
    radius = 0.5*phun("ft");
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (10.0*2*pi)^2;
    stabfact = 0.0062

    MR = DeforModelRed3D
    output = import_ABAQUS("alum_cyl.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm") # The input is provided in SI(mm) units
    fens, fes = T10toT4(fens, fes)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material, stabfact)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(Symmetric(K+OmegaShift*M), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    true

end # alum_cyl_modes

function alum_cyl_mode_esnice_t4()

    E = 70000*phun("MPa");
    nu = 0.33;
    rho = 2700*phun("KG/M^3");
    radius = 0.5*phun("ft");
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (10.0*2*pi)^2;

    MR = DeforModelRed3D
    output = import_ABAQUS("alum_cyl.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm") # The input is provided in SI(mm) units
    fens, fes = T10toT4(fens, fes)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(Symmetric(K+OmegaShift*M), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    true

end # alum_cyl_modes

function alum_cyl_modes_algo()
    println("""
    % Vibration modes of unit cube  of almost incompressible material.
    %
    % Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    % tetrahedral. International Journal for Numerical Methods in
    % Engineering 67: 841-867.""")
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho= 1*phun("KG/M^3");
    a=1*phun("M"); b=a; h= a;
    n1=2;# How many element edges per side?
    na= n1; nb= n1; nh =n1;
    neigvs=20                   # how many eigenvalues
    omega_shift=(0.1*2*pi)^2;

    fens,fes =H20block(a,b,h, na,nb,nh)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)),
    material), "femm_mass"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)),
    material))

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1],
    "omega_shift"=>omega_shift, "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtools.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    modeldata["postprocessing"] = FDataDict("file"=>"alum_cyl_mode",
    "mode"=>10)
    modeldata=FinEtools.AlgoDeforLinearModule.exportmode(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

    true

end # alum_cyl_modes_algo

function alum_cyl_modes_export()
    println("""
    Vibration modes of unit cube  of almost incompressible material.

    This example EXPORTS the model to Abaqus.

    Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    tetrahedral. International Journal for Numerical Methods in
    Engineering 67: 841-867.
    """)
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 5;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;

    MR = DeforModelRed3D
    fens,fes  = H20block(a,b,h, na,nb,nh)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material)

    K =stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M =mass(femm, geom, u)
    d,v,nev,nconv = eigs(Symmetric(K+OmegaShift*M), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    d = d - OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    mode = 7
    scattersysvec!(u, v[:,mode])
    File =  "alum_cyl_modes.vtk"
    vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
    @async run(`"paraview.exe" $File`)

    AE = AbaqusExporter("alum_cyl_modes_h20");
    # AE.ios = STDOUT;
    HEADING(AE, "Vibration modes of unit cube  of almost incompressible material.");
    COMMENT(AE, "The  first six frequencies are rigid body modes.");
    COMMENT(AE, "The  first nonzero frequency (7) should be around 0.26 Hz");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    COMMENT(AE, "The hybrid form of the serendipity hexahedron is chosen because");
    COMMENT(AE, "the material is  nearly incompressible.");
    ELEMENT(AE, "c3d20rh", "AllElements", 1, connasarray(fes))
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    DENSITY(AE, rho)
    STEP_FREQUENCY(AE, neigvs)
    END_STEP(AE)
    close(AE)

    true

end # alum_cyl_modes_export

function alum_cyl_modes_msh8_algo()
    println("""
    % Vibration modes of unit cube  of almost incompressible material.
    % Mean-strain hexahedron.
    % Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    % tetrahedral. International Journal for Numerical Methods in
    % Engineering 67: 841-867.""")
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho= 1*phun("KG/M^3");
    a=1*phun("M"); b=a; h= a;
    n1=8 # How many element edges per side?
    na= n1; nb= n1; nh =n1;
    neigvs=20                   # how many eigenvalues
    omega_shift=(0.1*2*pi)^2;

    fens,fes = H8block(a,b,h, na,nb,nh)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)),
    material), "femm_mass"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,3)),
    material))

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1],
    "omega_shift"=>omega_shift, "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtools.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    modeldata["postprocessing"] = FDataDict("file"=>"alum_cyl_mode",
    "mode"=>10)
    modeldata=FinEtools.AlgoDeforLinearModule.exportmode(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

    true

end # alum_cyl_modes_msh8_algo

function allrun()
    println("#####################################################")
    println("# alum_cyl_mode_nice_t4 ")
    alum_cyl_mode_nice_t4()
    println("#####################################################")
    println("# alum_cyl_mode_esnice_t4 ")
    alum_cyl_mode_esnice_t4()
    return true
end # function allrun

end # module alum_cyl_mode_examples


text:
# Constants
# =========

@enum htsLogLevel begin
    HTS_LOG_OFF
    HTS_LOG_ERROR
    HTS_LOG_WARNING = 3
    HTS_LOG_INFO
    HTS_LOG_DEBUG
    HTS_LOG_TRACE
end

const HTS_FMT_CSI  = Cint(0)
const HTS_FMT_BAI  = Cint(1)
const HTS_FMT_TBI  = Cint(2)
const HTS_FMT_CRAI = Cint(3)

# CIGAR operations
const BAM_CMATCH     = Cint(0)
const BAM_CINS       = Cint(1)
const BAM_CDEL       = Cint(2)
const BAM_CREF_SKIP  = Cint(3)
const BAM_CSOFT_CLIP = Cint(4)
const BAM_CHARD_CLIP = Cint(5)
const BAM_CPAD       = Cint(6)
const BAM_CEQUAL     = Cint(7)
const BAM_CDIFF      = Cint(8)
const BAM_CBACK      = Cint(9)

const BAM_CIGAR_STR   = b"MIDNSHP=XB"
const _BAM_CIGAR_STR_PADDED = b"MIDNSHP=XB??????"  # padded for safe lookup
const BAM_CIGAR_SHIFT = Cint(4)
const BAM_CIGAR_MASK  = 0xf
const BAM_CIGAR_TYPE  = 0x3C1A7

const bam_cigar_table = [unsafe_load(cglobal((:bam_cigar_table, libhts), UInt8), i) for i in 1:256]

# Record flags
"The read is paired in sequencing, no matter whether it is mapped in a pair"
const BAM_FPAIRED        = UInt16(1 << 0)
"The read is mapped in a proper pair"
const BAM_FPROPER_PAIR   = UInt16(1 << 1)
"The read itself is unmapped; conflictive with `BAM_FPROPER_PAIR`"
const BAM_FUNMAP         = UInt16(1 << 2)
"The mate is unmapped"
const BAM_FMUNMAP        = UInt16(1 << 3)
"The read is mapped to the reverse strand"
const BAM_FREVERSE       = UInt16(1 << 4)
"The mate is mapped to the reverse strand"
const BAM_FMREVERSE      = UInt16(1 << 5)
"This is read1"
const BAM_FREAD1         = UInt16(1 << 6)
"This is read2"
const BAM_FREAD2         = UInt16(1 << 7)
"Not primary alignment"
const BAM_FSECONDARY     = UInt16(1 << 8)
"QC failure"
const BAM_FQCFAIL        = UInt16(1 << 9)
"Optical or PCR duplicate"
const BAM_FDUP           = UInt16(1 << 10)
"Supplementary alignment"
const BAM_FSUPPLEMENTARY = UInt16(1 << 11)


text:
# Test Linear System Optimizaiton withh one switching time
function linsystest(solver=IpoptSolver(tol=1e-03); objtol = 1e-3, primaltol = 1e-3)
  println("Testing Switched Linear Systems Optimization ", string(typeof(solver)))

  # Size of the state space
  nx = 2;

  # Cost function Matrix
  Q = eye(2)

  # Initial State
  x0 = [1.0; 1.0]

  # Define initial and final time
  t0 = 0.0
  tf = 1.0

  ### Define System Dynamics
  A = zeros(nx, nx, 2)
  A[:, :, 1] = [-1 0;
                1  2]
  A[:, :, 2] = [1 1;
                1 -2]

  m = stoproblem(x0, A, solver=solver)
  solve!(m)
  @test_approx_eq_eps gettau(m) 0.26486646235103123 primaltol
  @test_approx_eq_eps getobjval(m) 5.2545429449272145 objtol

  println("Passed")
end


text:
# for Q dynamics set α = r-λ*(exp(μ+0.5*δ^2)-1)
@sde_model Merton begin
  dS = S*α*dt + S*σ*dW + S*(ξ-1.0)*dN
  J ~ LogNormal(μ, δ)
end α σ λ μ δ

function sample(model::Merton, scheme::Exact, t0, x0)
  @unpack α, σ, λ, μ, δ = model
  N = rand(Poisson(model.λ*scheme.Δt))
  rand(LogNormal(log(x0)+(α-0.5*σ^2)*scheme.Δt+N*μ, sqrt(N*δ^2+scheme.Δt*σ^2)))
end


text:
module SpecialMatrices

using LinearAlgebra
using FFTW
using Polynomials

import LinearAlgebra: Matrix, inv, det, ishermitian, isposdef, mul!

import Base: getindex, isassigned, size, *
import Base.\

include("cauchy.jl") #Cauchy matrix
include("companion.jl") #Companion matrix
include("frobenius.jl") #Frobenius matrix
include("hankel.jl") #Hankel matrix
include("hilbert.jl") #Hilbert matrix
include("kahan.jl") #Kahan matrix
include("riemann.jl") #Riemann matrix
include("strang.jl") #Strang matrix
include("toeplitz.jl") #Toeplitz matrix, Circulant matrix
include("vandermonde.jl") #Vandermonde matrix

end # module


text:
using Pardiso

export newton, L₁

L₁(V::Array) = sum(abs.(V)) / length(V)
L₂(V::Array, xs::Array) = trapz(xs, V.^2)

function newton(V₀::Array, rhs::Function, jac::Function, p::Dict, N::Int, Δx::Real)
    max_iter = 100
    tol = 1e-8

    V = copy(V₀)

    ps = MKLPardisoSolver()

    for i in 1:max_iter
        Y = rhs(V, p, N, Δx)
        error = L₁(Y)

        println("Running iteration n $i, error is $error")
        if error < tol
            break
        end

        J = jac(V, p, N, Δx)
        Nx, Ny = size(J)
        Nv = length(V)
        V += solve(ps, J, -Y)
    end
    V
end


text:
using ApproxFun
using Plots
using BenchmarkTools

function sir_eqn(S,I,u0,p)
  (β,γ) = p
  (S0,I0) = u0
   return [S(0)-S0,
           I(0)-I0,
           S' + β*S*I,
           I' - β*S*I + γ*I]
end;

tmax = 40.0
t=Fun(identity, 0..tmax);

u0 = [990.0,10.0]; # S,I

p = [0.0005,0.25]; # β,γ

S,I = newton((S,I)->sir_eqn(S,I,u0,p), u0 .* one(t); maxiterations=50);

plot(S,label="S",xlabel="Time",ylabel="Number")
plot!(I,label="I")

@benchmark newton((S,I)->sir_eqn(S,I,u0,p), u0 .* one(t); maxiterations=50)


text:
struct IdentityHandler end

is_hardware_primitive(::Type{typeof(identity)}, x) = true
is_hardware_primitive(::Type{typeof(Base.broadcasted)}, ::Type{typeof(identity)}, x) = true
gethandler(::Bool, ::Type{typeof(identity)}, x) = IdentityHandler()
init_state(::IdentityHandler) = nothing

function (handler::IdentityHandler)(buffer, netlist, state, inputs, outputs)
    if issigned(inputs[1])
        setsigned!(netlist, outputs[1], true)
        write(buffer, """
            $stdcomment
            // BEGIN identity
            assign $(name(outputs[1]))_p = $(name(inputs[1]))_p;
            assign $(name(outputs[1]))_m = $(name(inputs[1]))_m;
            // END identity
            \n""")
    else
        write(buffer, """
            $stdcomment
            // BEGIN identity
            assign $(name(outputs[1])) = $(name(inputs[1]));
            // END identity
            \n""")
    end

    return buffer, state
end


text:
function spabac!(a::Vector{Float64}, b::Vector{Float64}, kdiag::Vector{Int64})
  n = size(kdiag, 1)
  b[1] /= a[1]
  for i in 2:n
    ki = kdiag[i] - i
    l = kdiag[i-1] - ki + 1
    x = b[i]
    if l != i
      m = i - 1
      for j in l:m
        x -= a[ki+j] * b[j]
      end
    end
    b[i] = x / a[ki+i]
  end
  for it in 2:n
    i = n + 2 - it
    ki = kdiag[i] - i
    x = b[i] / a[ki+i]
    b[i] = x
    l = kdiag[i-1] - ki + 1
    if l != i
      m = i - 1
      for k in l:m
        b[k] -= x * a[ki+k]
      end
    end
  end
  b[1] /= a[1]
  nothing
end

function spabac(a::Vector{Float64}, b_in::Vector{Float64}, kdiag::Vector{Int64})
  b = deepcopy(b_in)
  n = size(kdiag, 1)
  b[1] /= a[1]
  for i in 2:n
    ki = kdiag[i] - i
    l = kdiag[i-1] - ki + 1
    x = b[i]
    if l != i
      m = i - 1
      for j in l:m
        x -= a[ki+j] * b[j]
      end
    end
    b[i] = x / a[ki+i]
  end
  for it in 2:n
    i = n + 2 - it
    ki = kdiag[i] - i
    x = b[i] / a[ki+i]
    b[i] = x
    l = kdiag[i-1] - ki + 1
    if l != i
      m = i - 1
      for k in l:m
        b[k] -= x * a[ki+k]
      end
    end
  end
  b[1] /= a[1]
  b
end


text:
"""
`FourierLayer(in, out, grid, modes, σ=identity, init=glorot_uniform)`
`FourierLayer(Wf::AbstractArray, Wl::AbstractArray, [bias_f, bias_l, σ])`

Create a Layer of the Fourier Neural Operator as proposed by Li et al.
arXiv: 2010.08895

The layer does a fourier transform on the grid dimension of the input array,
filters higher modes out by the weight matrix and transforms it to the
specified output dimension such that In x M x N -> Out x M x N.
The output though only contains the relevant Fourier modes with the rest padded to zero
in the last axis as a result of the filtering.

The input `x` should be a rank 3 tensor of shape
(num parameters (`in`) x num grid points (`grid`) x batch size (`batch`))
The output `y` will be a rank 3 tensor of shape
(`out` x num grid points (`grid`) x batch size (`batch`))

You can specify biases for the paths as you like, though the convolutional path is
originally not intended to perform an affine transformation.

# Examples
Say you're considering a 1D diffusion problem on a 64 point grid. The input is comprised
of the grid points as well as the IC at this point.
The data consists of 200 instances of the solution.
Beforehand we convert the two input channels into a higher-dimensional latent space with 128 nodes by using a regular `Dense` layer.
So the input takes the dimension `128 x 64 x 200`.
The output would be the diffused variable at a later time, which initially makes the output of the form `128 x 64 x 200` as well. Finally, we have to squeeze this high-dimensional ouptut into the one quantity of interest again by using a `Dense` layer.

We wish to only keep the first 16 modes of the input and work with the classic sigmoid function as activation.

So we would have:

```julia
model = FourierLayer(128, 128, 100, 16, σ)
```
"""
struct FourierLayer{F,Tc<:Complex{<:AbstractFloat},Tr<:AbstractFloat,Bf,Bl}
    # F: Activation, Tc/Tr: Complex/Real eltype
    Wf::AbstractArray{Tc,3}
    Wl::AbstractMatrix{Tr}
    grid::Int
    σ::F
    λ::Int
    bf::Bf
    bl::Bl
    # Constructor for the entire fourier layer
    function FourierLayer(
        Wf::AbstractArray{Tc,3}, Wl::AbstractMatrix{Tr},
        grid::Int, σ::F = identity,
        λ::Int = 12, bf = true, bl = true) where
        {F,Tc<:Complex{<:AbstractFloat},Tr<:AbstractFloat}

        # create the biases with one singleton dimension
        bf = Flux.create_bias(Wf, bf, 1, size(Wf,2), size(Wf,3))
        bl = Flux.create_bias(Wl, bl, 1, size(Wl,1), grid)
        new{F,Tc,Tr,typeof(bf),typeof(bl)}(Wf, Wl, grid, σ, λ, bf, bl)
    end
end

# Declare the function that assigns Weights and biases to the layer
# `in` and `out` refer to the dimensionality of the number of parameters
# `modes` specifies the number of modes not to be filtered out
# `grid` specifies the number of grid points in the data
function FourierLayer(in::Integer, out::Integer, grid::Integer, modes = 12,
                        σ = identity; initf = cglorot_uniform, initl = Flux.glorot_uniform,
                        bias_fourier=true, bias_linear=true)

    # Initialize Fourier weight matrix (only with relevant modes)
    Wf = initf(in, out, modes)
    # Make sure filtering works
    @assert modes <= floor(Int, grid/2 + 1) "Specified modes exceed allowed maximum.
    The number of modes to filter must be smaller than N/2 + 1"
    # Pad the fourier weight matrix with additional zeros
    Wf = pad_zeros(Wf, (0, floor(Int, grid/2 + 1) - modes), dims=3)

    # Initialize Linear weight matrix
    Wl = initl(out, in)

    # Pass the bias bools
    bf = bias_fourier
    bl = bias_linear

    # Pass the modes for output
    λ = modes

    return FourierLayer(Wf, Wl, grid, σ, λ, bf, bl)
end

# Only train the weight array with non-zero modes
Flux.@functor FourierLayer
Flux.trainable(a::FourierLayer) = (a.Wf[:,:,1:a.λ], a.Wl,
                                typeof(a.bf) != Flux.Zeros ? a.bf[:,:,1:a.λ] : nothing,
                                typeof(a.bl) != Flux.Zeros ? a.bl : nothing)

# The actual layer that does stuff
function (a::FourierLayer)(x::AbstractArray)
    # Assign the parameters
    Wf, Wl, bf, bl, σ, = a.Wf, a.Wl, a.bf, a.bl, NNlib.fast_act(a.σ, x)

    # Do a permutation: DataLoader requires batch to be the last dim
    # for the rest, it's more convenient to have it in the first one
    xp = permutedims(x, [3,1,2])

    # The linear path
    # x -> Wl
    @ein linear[batch, out, grid] := Wl[out, in] * xp[batch, in, grid]
    linear .+ bl

    # The convolution path
    # x -> 𝔉 -> Wf -> i𝔉
    # Do the Fourier transform (FFT) along the grid dimension of the input and
    # Multiply the weight matrix with the input using batched multiplication
    # We need to permute the input to (channel,batch,grid), otherwise batching won't work
    @ein 𝔉[batch, out, grid] := Wf[in, out, grid] * rfft(xp, 3)[batch, in, grid]
    𝔉 .+ bf

    # Do the inverse transform
    # We need to permute back to match the shape of the linear path
    i𝔉 = irfft(𝔉, size(xp,3),3)

    # Return the activated sum
    return permutedims(σ.(linear + i𝔉), [2,3,1])
end

# Print nicely
function Base.show(io::IO, l::FourierLayer)
    print(io, "FourierLayer with\nConvolution path: (", size(l.Wf, 2), ", ",
            size(l.Wf, 1), ", ", size(l.Wf, 3))
    print(io, ")\n")
    print(io, "Linear path: (", size(l.Wl, 2), ", ", size(l.Wl, 1))
    print(io, ")\n")
    print(io, "Fourier modes: ", l.λ)
    print(io, "\n")
    l.σ == identity || print(io, "Activation: ", l.σ)
end


text:
# TODO the two functions elemental and binary need to be put into one
# and comments need to be added.

function elemental(i, j, d, D, atoms_numbers)

    n_atoms = size(atoms_numbers, 1)

    i = i .+ 1
    j = j .+ 1

    n_i  = zeros(Int, n_atoms)
    for ii=1:size(i, 1)
        n_i[i[ii]] += 1
    end

    i_grouped = [Array{Int}(undef, n) for n in n_i]
    counter_group = ones(Int, n_atoms)

    for ii=1:size(i, 1)
        i_grouped[i[ii]][counter_group[i[ii]], 1] = ii
        counter_group[i[ii]] += 1
    end

    n_triplets::Int = 0
    for ii=1:n_atoms
        n_triplets += (counter_group[ii]-1) * (counter_group[ii]-2) / 2
    end

    D = copy(transpose(D))
    d_out = Array{Float64}(undef, 3, n_triplets)
    counter = 1
    for i_atom=1:n_atoms
      @inbounds for iii=1:n_i[i_atom] for jjj=iii+1:n_i[i_atom]
            @inbounds ii = i_grouped[i_atom][iii]
            @inbounds jj = i_grouped[i_atom][jjj]

            b = 0.
            for s=1:3
                @inbounds b += (D[s, ii] - D[s, jj])^2
            end
            b = sqrt(b)

            @inbounds d_out[1, counter] = d[ii]
            @inbounds d_out[2, counter] = d[jj]
            @inbounds d_out[3, counter] = b
            counter += 1
        end
        end
    end
    return d_out
end

function binary(i, j, d, D, atoms_numbers)

    atoms_tuples = sort(collect(Set(atoms_numbers)))
    A, B = atoms_tuples
    n_atoms = size(atoms_numbers, 1)
    s = size(atoms_tuples, 1)
    a= [(atoms_tuples[i], atoms_tuples[j], atoms_tuples[k]) for i=1:s for j=1:s for k=j:s]
    sizes_dict = Dict( aa => 0 for aa=a)

    i = i .+ 1
    j = j .+ 1

    bins = zeros(Int, 2, n_atoms)
    n_i  = zeros(Int, n_atoms)
    for ii=1:size(i, 1)
        n_i[i[ii]] += 1
        if atoms_numbers[j[ii]] == atoms_tuples[1]
            bins[1, i[ii]] += 1
        else
            bins[2, i[ii]] += 1
        end
    end

    i_grouped = [Array{Int}(undef, n) for n in n_i]
    counter_group = ones(Int, n_atoms)

    for ii=1:size(i, 1)
        i_grouped[i[ii]][counter_group[i[ii]], 1] = ii
        counter_group[i[ii]] += 1
    end

    for ii=1:n_atoms
        if bins[1, ii] > 1
            sizes_dict[(atoms_numbers[ii], atoms_tuples[1], atoms_tuples[1])] += (bins[1, ii] - 1)bins[1, ii] / 2
        end

        if (bins[1, ii] > 0) & (bins[2, ii] > 0)
            sizes_dict[(atoms_numbers[ii], atoms_tuples[1], atoms_tuples[2])] += bins[1, ii] * bins[2, ii]
        end

        if bins[2, ii] > 1
            sizes_dict[(atoms_numbers[ii], atoms_tuples[2], atoms_tuples[2])] += (bins[2, ii] - 1)bins[2, ii] / 2
        end
    end

    D = copy(transpose(D))
    ds = [Array{Float64}(undef, 3, sizes_dict[aa]) for aa=a]
    counter = [1 for aa=a]
    for i_atom=1:n_atoms
      @inbounds for iii=1:n_i[i_atom] for jjj=iii+1:n_i[i_atom]
            @inbounds ii = i_grouped[i_atom][iii]
            @inbounds jj = i_grouped[i_atom][jjj]

            @inbounds if atoms_numbers[j[ii]] < atoms_numbers[j[jj]]
                @inbounds tup = (atoms_numbers[i_atom], atoms_numbers[j[ii]], atoms_numbers[j[jj]])
            else
                @inbounds tup = (atoms_numbers[i_atom], atoms_numbers[j[jj]], atoms_numbers[j[ii]])
            end

            b = 0.
            for s=1:3
               @inbounds b += (D[s, ii] - D[s, jj])^2
            end
            b = sqrt(b)

            if     (tup[1] == A) & (tup[2] == A) & (tup[3] == A)
                idx = 1
            elseif (tup[1] == A) & (tup[2] == A) & (tup[3] == B)
                idx = 2
            elseif (tup[1] == A) & (tup[2] == B) & (tup[3] == B)
                idx = 3
            elseif (tup[1] == B) & (tup[2] == A) & (tup[3] == A)
                idx = 4
            elseif (tup[1] == B) & (tup[2] == A) & (tup[3] == B)
                idx = 5
            else
                idx = 6
            end
            c = counter[idx]
            @inbounds ds[idx][1, c] = d[ii]
            @inbounds ds[idx][2, c] = d[jj]
            @inbounds ds[idx][3, c] = b
            counter[idx] += 1
        end
        end
    end
    return ds
end

function get_3b_from_2b_desc(i, j, d, D, atoms_numbers)
    if length(Set(atoms_numbers)) == 1
        elemental(i, j, d, D, atoms_numbers)
    else
        binary(i, j, d, D, atoms_numbers)
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Serd_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Serd")
JLLWrappers.@generate_main_file("Serd", UUID("0d4e85e1-a563-5338-9f64-3e4dc8b67009"))
end  # module Serd_jll


text:
using FastGaussQuadrature
using LinearAlgebra

struct PBSpline{T}
    p :: Int       # spline degree
    nₕ :: Int      # number of splines
    L :: T         # domain length
    h :: T         # element width

    function PBSpline(p::Int, nₕ::Int, L::T, h::T) where {T}
        @assert nₕ > p
        new{T}(p, nₕ, L, h)
    end
end

# convenience constructor
PBSpline(p::Int, nₕ::Int, L) = PBSpline(p, nₕ, L, L/nₕ)

Base.size(S::PBSpline) = (1,)
Base.length(S::PBSpline) = 1
Base.iterate(S::PBSpline, state=1) = state > 1 ? nothing : (S, state+1)

@inline function eval_unitPBSpline(p::Int, x::T) where {T}
    if p == 0
        return zero(T) <= x < one(T) ? one(T) : zero(T)
    elseif x > p+1
        return zero(T)
    else
        return x/p * eval_unitPBSpline(p-1,x) + (p+1-x)/p * eval_unitPBSpline(p-1,x-1)
    end
end

@inline function eval_PBSpline(S::PBSpline{T}, i::Int, x::T) where {T}
    x /= S.h; x -= (i-1)
    # periodicity
    while x < zero(T)  x += S.nₕ end
    while x > S.nₕ     x -= S.nₕ end
    return eval_unitPBSpline(S.p, x)
end

function eval_PBSBasis(c::Array{T}, S::PBSpline{T}, x::T) where {T}
    @assert length(c) == S.nₕ
    r::T = 0
    for j in 1:S.nₕ
        r += c[j] * eval_PBSpline(S, j, x)
    end
    return r
end

function eval_PBSBasis(c::Array{T}, S::PBSpline{T}, x::Array{T}, r = zero(x)) where {T}
    @assert length(c) == S.nₕ
    for i in eachindex(r,x)
        r[i] = eval_PBSBasis(c, S, x[i])
    end
    return r
end

# evaluation point(s) x, degree p
function eval_deriv_PBSpline(S::PBSpline{T}, i::Int, x::T) where {T}
    x /= S.h; x -= i-1
    # periodicity
    while x < zero(T)  x += S.nₕ end
    while x > S.nₕ     x -= S.nₕ end
    1/S.h * (eval_unitPBSpline(S.p-1, x) - eval_unitPBSpline(S.p-1, x-1 < 0 ? x-1+S.nₕ : x-1) )
end

function eval_deriv_PBSBasis(c::Array{T}, S::PBSpline{T}, x::T) where {T}
    @assert length(c) == S.nₕ
    r::T = 0
    for j in 1:S.nₕ
        r += c[j] * eval_deriv_PBSpline(S, j, x)
    end
    return r
end

function eval_deriv_PBSBasis(c::Array{T}, S::PBSpline{T}, x::Array{T}, r = zero(x)) where {T}
    for i in eachindex(x)
        r[i] = eval_deriv_PBSBasis(c, S, x[i])
    end
    return r
end

function integrate_gausslegendre(f, a::T, b::T, nq::Int) where {T}
    nodes, weights = gausslegendre(nq)
    return 0.5*(b-a)*dot(weights, f( 0.5 .* (b-a) .* nodes .+ 0.5 .* (a+b) ) )
    # 0.5 * (b-a) * mapreduce(q -> q[2] * f( 0.5 * (b-a) * q[1] + 0.5 * (a+b) ), +, zip(gausslegendre(nq)...))
    # 0.5 * (b-a) * mapreduce((c,w) -> w * f( 0.5 * (b-a) * c + 0.5 * (a+b) ), +, zip(gausslegendre(nq)...))
end

function rhs_PBSBasis(f, S::PBSpline{T}, q = S.p+1, rhs = zeros(T, S.nₕ)) where {T}
    for i in 1:S.nₕ
        for k in 0:S.p
            if i+k <= S.nₕ
                rhs[i] += integrate_gausslegendre(x -> (eval_PBSpline(S, i, x) .* f(x)), (i-1+k)*S.h, (i+k)*S.h, q)
            else
                rhs[i] += integrate_gausslegendre(x -> (eval_PBSpline(S, i, x) .* f(x)), (i-1+k-S.nₕ)*S.h, (i+k-S.nₕ)*S.h, q)
            end
        end
    end
    return rhs
end

function rhs_particles_PBSBasis(x::Vector{T}, w::Vector{T}, S::PBSpline{T}, rhs = zeros(T, S.nₕ)) where {T}
    for i in eachindex(rhs)
        rhs[i] = 0
        for j in eachindex(w,x)
            rhs[i] += w[j] * eval_PBSpline(S, i, x[j])
        end
    end
    return rhs
end


text:
using Test
using VertexFinder: find

#    5--7--6
#    |  | /
# 4--8--2
# |  |  |
# 9--1--3

G = Dict(
    1 => (3, 8, 9),
    2 => (3, 8, 7),
    3 => (1, 2),
    4 => (8, 9),
    5 => (7, 8),
    6 => (2, 7),
    7 => (5, 2, 6),
    8 => (1, 2, 4, 5),
    9 => (1, 4))

w1 = find(G, 1)
w2 = find(G, w1)
w3 = find(G, w2)
@test w1 == 6
@test w2 == 9
@test w3 == 6


text:
include("./optimizers.jl")
include("./grad.jl")

function update!(opt, xs::Params, g::Grads)
    for x in xs
        x == nothing && continue
        update!(opt, x, g[x])
    end
end

function update!(opt, x, g)
    x .-= apply!(opt, x, g)
end

function train!(m::NetWork, loss, data, opt)
    ps = Params(m)
    #=
    for d in data
        try
            g = grad(m, d, loss, ps)
            update!(opt, ps, g)
        catch
            @warn "you can't train by one data because something wrong with the data"
        end
    end
    =#
    i = 1
    for d in data
        i += 1
        g = grad(m, d, loss, ps)
        update!(opt, ps, g)
    end
end


text:
using Test, QuHamiltonian

@testset "neighbor" begin
    @neighbor 2 i,j
    @test i == QuHamiltonian.Neighbor{2}(:i)
    @test j == QuHamiltonian.Neighbor{2}(:j)
end

@testset "nearest" begin
    @nearest i, j
    @test i == QuHamiltonian.Neighbor{1}(:i)
    @test j == QuHamiltonian.Neighbor{1}(:j)
end

@testset "vertex" begin
    @vertex i, j
    @test i == QuHamiltonian.Vertex(:i)
    @test j == QuHamiltonian.Vertex(:j)
end

@testset "surround" begin
    @surround i, j
    @test i == QuHamiltonian.Surround(:i)
    @test j == QuHamiltonian.Surround(:j)

    # NOTE: surround a vertex place holder is more readable?
    # @vertex k
    # @surround k
end

@testset "face" begin
    @face i, j
    @test i == QuHamiltonian.Face(:i)
    @test j == QuHamiltonian.Face(:j)
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module DataFramesTimings

using Random
using DataFrames
using DataFramesMeta

Random.seed!(1)
const n = 5_000_000
da = rand(n)
db = rand(n)
df = DataFrame(a = da, b = db)

function dot3(df::DataFrame)
    da, db = df.a, df.b
    T = eltype(da)
    x = 0.0
    for i in 1:length(da)
        x += df.a[i] * df.a[i]
    end
    return x
end

function dot4(df::DataFrame)
    da, db = df.a, df.b
    T = eltype(da)
    x = 0.0
    for i in 1:length(da)
        x += da[i]::T * db[i]::T
    end
    return x
end

@show t3 = @elapsed dot3(df)
@show t4 = @elapsed dot4(df)

end # module


text:
module FixStrings
export ASCIIStr, UTF8Str, ByteStr

# Deal with mess from #16058
if VERSION >= v"0.5.0-dev+3876"
    const ASCIIStr = String
    const UTF8Str = String
    const ByteStr = String
else
    const ASCIIStr = ASCIIString
    const UTF8Str = UTF8String
    const ByteStr = ByteString
end
end


text:
function StatisticalSummaryValues(arg0::jdouble, arg1::jdouble, arg2::jlong, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return StatisticalSummaryValues((jdouble, jdouble, jlong, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function equals(obj::StatisticalSummaryValues, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_max(obj::StatisticalSummaryValues)
    return jcall(obj, "getMax", jdouble, ())
end

function get_mean(obj::StatisticalSummaryValues)
    return jcall(obj, "getMean", jdouble, ())
end

function get_min(obj::StatisticalSummaryValues)
    return jcall(obj, "getMin", jdouble, ())
end

function get_n(obj::StatisticalSummaryValues)
    return jcall(obj, "getN", jlong, ())
end

function get_standard_deviation(obj::StatisticalSummaryValues)
    return jcall(obj, "getStandardDeviation", jdouble, ())
end

function get_sum(obj::StatisticalSummaryValues)
    return jcall(obj, "getSum", jdouble, ())
end

function get_variance(obj::StatisticalSummaryValues)
    return jcall(obj, "getVariance", jdouble, ())
end

function hash_code(obj::StatisticalSummaryValues)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::StatisticalSummaryValues)
    return jcall(obj, "toString", JString, ())
end


text:
#Julia v0.6.2
#Julia で Pkg.add("PyCall"), python3 でMDAnalysis MDAnalysisTests numpy が必要

using PyCall
@pyimport MDAnalysis as mda
@pyimport MDAnalysis.topology.guessers as mda_tg

#考慮対象外の水とイオン
#他に考慮しない残基名等あればここに追記してください
#最初に回転したりなんかしたりするなら不要？
const wation = ["HOH","H2O","WAT","Na+","Cl-"]

const boltzmann =  1.380649e-23 / 4.184
const boltzmannpermol = 1.380649 * 6.02214076 / 4.184

function dihedral_divide(p1, p2, p3, p4)
    v1 = p1-p2
    v2 = p3-p2
    v3 = p4-p3

    normalize!(v2)
    # v1/3 minus component that aligns with v2
    v1 -= (v1⋅v2)v2
    v3 -= (v3⋅v2)v2

    x = v1⋅v3
    y = (v2×v1)⋅v3

    return atan2(y,x) #(-π,π]
end
dihedral(p::Array{Float64,2}) = dihedral_divide(p[1,:], p[2,:], p[3,:], p[4,:])
dihedral(p::Array{Float64,2}, i1::Int, i2::Int, i3::Int, i4::Int) = dihedral_divide(p[i1,:], p[i2,:], p[i3,:], p[i4,:])
dihedral(p::Array{Float64,2}, iv::Vector{Int}) = dihedral(p, iv[1], iv[2], iv[3], iv[4])

function angle_divide(p1, p2, p3)
    v1 = p2-p1
    v2 = p2-p3
    normalize!(v1)
    normalize!(v2)
    return acos(v1⋅v2) #[0,π]
end
angle(p::Array{Float64,2}) = angle_divide(p[1,:], p[2,:], p[3,:])
angle(p::Array{Float64,2}, i1::Int, i2::Int, i3::Int) = angle_divide(p[i1,:], p[i2,:], p[i3,:])
angle(p::Array{Float64,2}, iv::Vector{Int}) = angle(p, iv[1], iv[2], iv[3])

function improper_divide(p1, p2, p3, p4)
    # p1 = root, p2 = hub, p3 & p4 = leaf
    # dihedral of plane(p1,p2,p3) and plane(p1,p2,p4)

    v1 = p2-p1
    v2 = p3-p2
    v3 = p4-p2

    normalize!(v1)
    v2 -= (v1⋅v2)v1
    v3 -= (v1⋅v3)v1

    x = v2⋅v3
    y = (v1×v2)⋅v3

    return atan2(y,x)
end
improper(p::Array{Float64,2}) = improper_divide(p[1,:], p[2,:], p[3,:], p[4,:])
improper(p::Array{Float64,2}, i1::Int, i2::Int, i3::Int, i4::Int) = improper_divide(p[i1,:], p[i2,:], p[i3,:], p[i4,:])
improper(p::Array{Float64,2}, iv::Vector{Int}) = improper(p, iv[1], iv[2], iv[3], iv[4])

bondlen(p::Array{Float64,2}, i1::Int, i2::Int) = vecnorm(p[i1,:] - p[i2,:])
bondlen(p::Array{Float64,2}, iv::Vector{Int}) = bondlen(p, iv[1], iv[2])

function connected_hydrogen(namelist::Vector{String}, connect, parent)
    ret = Int[]
    for i in connect[parent]
        if namelist[i][1] == 'H'
            push!(ret,i)
        end
    end
    return ret
end

const nm_conv = log(10.0^-9)
function jacobian(tm_val, angls_idx)
    tmpjac = 0.0
    #tmpjac += 2 * (log(tm_val[1]) )
    for i in 1:angls_idx
        #tmpjac += 2 * (log(tm_val[i+1]))
        tmpjac += log( sin(tm_val[i+angls_idx+1]) )
    end
    return tmpjac
end

function load_coord!(vals, io::IO, n_trajectory::Int, head::Int)
    seekstart(io)
    skip(io, 8 * (head-1) * n_trajectory )
    read!(io, vals)
end

function shannon(vals::Vector{Float64}, bins::Int, pitch::T where T<:AbstractFloat)
    nval = length(vals)
    ret = 0.0
    ctr = zeros(Int,bins+1)
    for i in vals
        ctr[floor(Int,i/pitch) + 1] += 1
    end
    for i in ctr
        if i == 0; continue; end
        p_i = Float64(i) / nval
        ret += p_i * log(p_i)
    end
    return -ret
end

function shannon_bond(vals::Vector{Float64}, bins::Int)
    vals .-= minimum(vals)
    pitch = maximum(vals) / bins
    return shannon(vals, bins, pitch)
end

function shannon_angle(vals::Vector{Float64}, bins::Int)
    vals .-= minimum(vals)
    pitch = maximum(vals) / bins
    return shannon(vals, bins, pitch)
end

function shannon_dihedral(vals::Vector{Float64}, bins::Int)
    vals .+= π
    pitch = 2π / bins
    return shannon(vals, bins, pitch)
end

function fast_cor(n_trajectory::Int, idx::Int, idy::Int)

    sx = unsafe_load(pmean,idx)
    sy = unsafe_load(pmean,idy)

    xy = 0.0

    @simd for i = 1:n_trajectory
        xi = unsafe_load(px,i)
        yi = unsafe_load(py,i)
        xi -= sx
        yi -= sy
        #xx += xi * xi
        #yy += yi * yi
        xy += xi * yi
    end

    xx = unsafe_load(psp,idx)
    yy = unsafe_load(psp,idy)

    return (xy / (sqrt(xx*yy)))
end

ff14types = Dict(
       "CA" => 1.5,
       "N3" => 1.1,
       "CN" => 1.5,
       "CT" => 1.5,
       "CO" => 1.5,
       "CB" => 1.5,
       "HA" => 0.4,
       "HC" => 0.4,
       "H4" => 0.4,
       "HO" => 0.4,
       "OH" => 1.05,
       "CX" => 1.5,
       "O2" => 1.05,
       "NA" => 1.1,
       "H1" => 0.4,
       "HP" => 0.4,
       "3C" => 1.5,
       "C8" => 1.5,
       "C*" => 1.5,
       "2C" => 1.5,
       "CW" => 1.5,
       "N2" => 1.1)

#テストのお約束
#PSF = mdad.PSF #data/adk.psf
#DCD = mdad.DCD #data/adk_dims.dcd
TPR = ARGS[1]
TRR = ARGS[2]
const test = mda.Universe(TPR,TRR)

const atomname = convert.(String, test[:atoms][:names]) #原子名
const resname = convert.(String, test[:atoms][:resnames]) #原子ごとの残基名
const n_atoms = length(atomname)

const allbonds = ARGS[4] == "0" ? test[:bonds] : mda_tg.guess_bonds(test[:atoms],test[:atoms][:positions],vdwradii=ff14types)

const n_bonds = length(allbonds)
@show const n_trajectory = length(test[:trajectory])
const temperature = float(ARGS[3])

#結合情報の格納　グラフの次数が高々4なのでリストで持つ
connect = [Int[] for i=1:n_atoms]
@inbounds for bond in allbonds
    bondno = bond[:atoms][:ix]
    push!(connect[bondno[1] + 1], bondno[2] + 1) #python : 0-index julia : 1-index
    push!(connect[bondno[2] + 1], bondno[1] + 1)
end

#参照先の格納
bonds = zeros(Int,(n_atoms-1,2))
angls = zeros(Int,(n_atoms-2,3))
dihds = zeros(Int,(n_atoms-3,4))
bonds_idx = angls_idx = dihds_idx = 1

#木のrootについて ad hoc な取扱い。
#グラフの次数で見て、0ならばイオンか何かなので飛ばす
#1ならばそれをrootに、4ならば多分N末N原子なので隣にある水素をどれか取って無理やりrootに
#タンパク質が存在しないと困るのでグラフの次数2-3はエラーを吐く　←改善の余地あり
rootidx = 1
seeking = true
while seeking
    if length(connect[rootidx]) == 4
        for i in connect[rootidx]
            if length(connect[i]) == 1
                rootidx = i
                seeking = false
                break
            end
        end
    elseif length(connect[rootidx]) == 1
        seeking = false
    elseif 2 <= length(connect[rootidx]) <= 3
        println( STDERR, "invalid connection order $(connect[rootidx]) in initializing the internal coordinate tree")
        return 1
    end
    rootidx += 1
end
rootidx -= 1

#rootからBFSして内部座標の構築
#rootを次数0として、次数1のノードは1つになるように取った
#グラフの次数2のノードが複数あるとdihedralのとり方が非自明になるので、どこか一個固定してimproperとしてみる
two_top = 0
#親
parent = zeros(Int,n_atoms)
parent[rootidx] = -1
#juliaだと特別dequeとかあるわけではなくvectorで全部ごり押すらしい
queue = [rootidx]
#念の為スコープの外で宣言
top = rootidx

finished = false
@inbounds while !finished
    if length(queue) != 0
        while length(queue) > 0
            top = shift!(queue)

            if parent[top] != -1
                bonds[bonds_idx,1] = top
                bonds[bonds_idx,2] = parent[top]
                bonds_idx += 1
            end
            if bonds_idx == 3
                two_top = top
            end

            if bonds_idx >= 3
                angls[angls_idx,1] = top
                angls[angls_idx,2] = parent[top]
                angls[angls_idx,3] = parent[parent[top]]
                angls_idx += 1
            end

            if bonds_idx >= 4

                #メチル基判定で使う
                #親ノードの接続している原子のうち原子名の1文字目が'H'であるインデックスのリスト
                #note:空になりうるので長さ判定を先に入れること
                cnh = connected_hydrogen(atomname, connect, parent[top])
                #次数2判定
                if parent[parent[parent[top]]] == -1
                    dihds[dihds_idx,1] = parent[parent[top]]
                    dihds[dihds_idx,2] = parent[top]
                    dihds[dihds_idx,3] = two_top
                    dihds[dihds_idx,4] = top
                #メチル基Hの判定　improperになるのはメチル基Hの2番手以降
                elseif atomname[top][1] == 'H' && atomname[parent[top]][1] == 'C' && length(cnh) == 3 && cnh[1] != top
                    dihds[dihds_idx,1] = parent[parent[top]]
                    dihds[dihds_idx,2] = parent[top]
                    dihds[dihds_idx,3] = cnh[1]
                    dihds[dihds_idx,4] = top
                #それ以外は通常のtorsion
                else
                    dihds[dihds_idx,1] = parent[parent[parent[top]]]
                    dihds[dihds_idx,2] = parent[parent[top]]
                    dihds[dihds_idx,3] = parent[top]
                    dihds[dihds_idx,4] = top
                end

                dihds_idx += 1
            end #add dihedral

            for i in connect[top]
                if parent[i] == 0
                    parent[i] = top
                    push!(queue,i)
                end
            end

        end #BFS while
    else # = queue is empty

        #チェインの切れ目とかが存在した場合、直前のところに無理やりつなぐ
        for i in 1:n_atoms
            if parent[i] == 0 && !(resname[i] in wation)
                parent[i] = top
                finished = !finished
                break
            end
        end
        #最後まで探索してしまった場合1回しかflipされないのでそのまま最外ループを抜ける
        finished = !finished

    end # if queue is not empty
end #while !finished

bonds_idx -= 1
angls_idx -= 1
dihds_idx -= 1
const n_coords = bonds_idx + angls_idx + dihds_idx
@show const n_nowation = (n_coords+6)/3

#S_{ind} の項
#各座標の時間発展を全部取ると確実にメモリからはみ出るので
#一度中間ファイルに全部書いてから一般化座標の一個ずつ足す
#水はないから多少はファイルサイズ減るはず
#-kBは括りだして後で掛ける

#S_{jac} の項 : jacobianがかかってくる項
#b2^2 ∏_{i=3}^{N} b_i^2 sin θ_i がjacobianで、それにlnがかかっているのでlnとって総和
#あとでpdf考慮するためにn_trajectoryで割ってkB掛ける

S_ind = 0.0
S_jac = 0.0

tempfile = open("mdentropy_tmp","w+") #w,r,create,truncate
unsafe_write(tempfile,Ref(UInt8(0)), 8 * n_coords * n_trajectory) #確保

tm_val = zeros(Float64, n_coords) #一般化座標の格納先
@inbounds for tm = 1:n_trajectory
    dump = test[:trajectory][tm] #フレーム移動のお約束

    coords = Float64.(test[:atoms][:positions])

    for i in 1:bonds_idx
        tm_val[i] = bondlen(coords,bonds[i,:])
    end
    for i in 1:angls_idx
        tm_val[i + bonds_idx] = angle(coords,angls[i,:])
    end
    for i in 1:dihds_idx
        i_dihds = dihds[i,:]
        #二面角がproper ⇔ 二面角計算の3,4番が結合している は逆が一般には偽だが
        #BFSで構築した都合上3員環が存在しない仮定を置くと成立する
        if parent[i_dihds[4]] == i_dihds[3]
            tm_val[i + bonds_idx + angls_idx] = dihedral(coords,dihds[i,:])
        else
            tm_val[i + bonds_idx + angls_idx] = improper(coords,dihds[i,:])
        end
    end

    #coord:元 tm_val:内部

    #情報落ち防止のためJ(b,θ)は一旦計算してからS_jacに足す
    S_jac += jacobian(tm_val,angls_idx)

    #各座標ごとにバラして格納
    #magic number の 8 は Float64 = (cのfloat) のbyte数
    seekstart(tempfile)
    skip(tempfile, 8*(tm-1) )
    p = pointer(tm_val)
    for i in 1:n_coords
        unsafe_write(tempfile,p,8)
        skip(tempfile,8*(n_trajectory-1))
        p += 8
    end
end

S_jac /= n_trajectory
S_jac *= boltzmannpermol
println()
println(STDOUT,"S_jac :", S_jac)
#seekstart(tempfile)

#15 = bin数
#あとで変数に変える
vals = zeros(Float64,n_trajectory)
for i = 1:bonds_idx
    load_coord!(vals, tempfile, n_trajectory, i)
    S_ind += shannon_bond(vals,15)
end
for i = 1:angls_idx
    load_coord!(vals, tempfile, n_trajectory, i + bonds_idx)
    S_ind += shannon_angle(vals,15)
end
for i = 1:dihds_idx
    load_coord!(vals, tempfile, n_trajectory, i + bonds_idx + angls_idx)
    S_ind += shannon_dihedral(vals,15)
end
S_ind *= boltzmannpermol
println()
println(STDOUT,"S_ind :", S_ind)

#ここから重いので対策
#繰り返し利用する平均と残差平方和を別に求めておき、相関行列を求めるときに都度持ってくる。

corr_matrix = zeros(Float64,(n_coords,n_coords))
pcorr = pointer(corr_matrix)

#作業スペース
x = zeros(Float64,n_trajectory)
y = zeros(Float64,n_trajectory)
const px = pointer(x)
const py = pointer(y)

const coef = 1.0

#平均
means = zeros(Float64,n_coords)
const pmean = pointer(means)
#平方和
sumpow = zeros(Float64,n_coords)
const psp = pointer(sumpow)

#xを使って上2つの作成
seekstart(tempfile)
for i=1:n_coords
    mx = 0.0
    read!(tempfile, x)
    @simd for j = 1:n_trajectory
        mx += unsafe_load(px,j)
    end
    mx /= n_trajectory
    unsafe_store!(pmean,mx,i)

    psx = 0.0
    @simd for j = 1:n_trajectory
        psx += (unsafe_load(px,j) - mx)^2
    end
    unsafe_store!(psp,psx,i)
end

S_indQH = 0.0
for i=1:n_coords
    @inbounds S_indQH += log(2*pi*e*sumpow[i])
end
S_indQH *= boltzmannpermol/2
println("S_indQH: ",S_indQH)
println()

for i = 1:n_coords
    #corr_matrix[i,i] = coef
    unsafe_store!(pcorr, coef, (i-1)*n_coords+i)

    load_coord!(x, tempfile, n_trajectory, i)

    for j=i+1:n_coords
        #一般化座標の連番に並んでいるので、yの大きさだけ読めばj番目の座標データが入る
        read!(tempfile, y)

        #ポインタを明示的に使って高速化
        #tmp = coef * cor(x,y)
        tmp = coef * fast_cor(n_trajectory, i, j)

        #debug
        #@assert abs( tmp - cor(x,y) * coef) < 1e-9 * coef

        #corr_matrix[j,i] = tmp
        unsafe_store!(pcorr, tmp, (i-1)*n_coords+j)
    end
end

corr_matrix = Symmetric(corr_matrix, :L)
#対称行列なのでBunch-Kauchman分解ができる
logdetcoefcorr, hoge = logabsdet(bkfact(corr_matrix))

logdetcorr = logdetcoefcorr - log(coef)*n_coords

S_corr = boltzmannpermol * logdetcorr / 2
println("S_corr: ",S_corr)

S_P = ((3n_nowation - 6.0) / 2 * log(temperature) + 1.5 * n_nowation ) * boltzmannpermol

println("S_P: ",S_P)

println("---------------------------------")
println("S_Conf: ", S_P+S_ind+S_corr+S_jac)

max_stars_repo_name: BILAB/scripts

max_issues_repo_name: BILAB/scripts

max_forks_repo_name: BILAB/scripts


text:
#####
##### Session management
#####
using Setfield

const DEFAULT_SESSIONS_DIR = "sessions"

"""
    SessionReport

The full collection of statistics and benchmark results
collected during a training session.

# Fields
- `iterations`: vector of ``n`` iteration reports with type
    [`Report.Iteration`](@ref)
- `benchmark`: vector of ``n+1`` benchmark reports with type
    [`Report.Benchmark`](@ref)
"""
struct SessionReport
  iterations :: Vector{Report.Iteration}
  benchmark :: Vector{Report.Benchmark}
  SessionReport() = new([], [])
end

function valid_session_report(r::SessionReport)
  length(r.benchmark) == length(r.iterations) == 0 && return true
  if length(r.benchmark) == length(r.iterations) + 1
    nduels = length(r.benchmark[1])
    return all(length(b) == nduels for b in r.benchmark)
  else
    return false
  end
end

"""
    Session{Env}

A wrapper on an AlphaZero environment that adds features such as:
- Logging and plotting
- Loading and saving of environments
In particular, it implements the [`Handlers`](@ref) interface.

# Public fields
- `env::Env` is the environment wrapped by the session
- `report` is the current session report, with type [`SessionReport`](@ref)
"""
mutable struct Session{Env}
  env :: Env
  dir :: String
  logger :: Logger
  autosave :: Bool
  save_intermediate :: Bool
  benchmark :: Vector{Benchmark.Evaluation}
  # Temporary state for logging
  progress :: Union{Progress, Nothing}
  report :: SessionReport

  function Session(env, dir, logger, autosave, save_intermediate, benchmark)
    return new{typeof(env)}(
      env, dir, logger, autosave, save_intermediate,
      benchmark, nothing, SessionReport())
  end
end

#####
##### Save and load environments
#####

const GSPEC_FILE       =  "gspec.data"
const BESTNN_FILE      =  "bestnn.data"
const CURNN_FILE       =  "curnn.data"
const MEM_FILE         =  "mem.data"
const ITC_FILE         =  "iter.txt"
const REPORT_FILE      =  "report.json"
const PARAMS_FILE      =  "params.data"
const PARAMS_JSON_FILE =  "params.json"      # not used when loading envs
const NET_PARAMS_FILE  =  "netparams.json"
const BENCHMARK_FILE   =  "benchmark.json"
const LOG_FILE         =  "log.txt"
const PLOTS_DIR        =  "plots"
const ITERS_DIR        =  "iterations"

iterdir(dir, i) = joinpath(dir, ITERS_DIR, "$i")

function valid_session_dir(dir)
  isfile(joinpath(dir, PARAMS_FILE)) &&
  isfile(joinpath(dir, BESTNN_FILE)) &&
  isfile(joinpath(dir, CURNN_FILE)) &&
  isfile(joinpath(dir, MEM_FILE)) &&
  isfile(joinpath(dir, ITC_FILE))
end

function save_env(env::Env, dir)
  isdir(dir) || mkpath(dir)
  serialize(joinpath(dir, GSPEC_FILE), env.gspec)
  serialize(joinpath(dir, PARAMS_FILE), env.params)
  open(joinpath(dir, PARAMS_JSON_FILE), "w") do io
    JSON3.pretty(io, JSON3.write(env.params))
  end
  open(joinpath(dir, NET_PARAMS_FILE), "w") do io
    JSON3.pretty(io, JSON3.write(Network.hyperparams(env.bestnn)))
  end
  serialize(joinpath(dir, BESTNN_FILE), Network.params(env.bestnn))
  serialize(joinpath(dir, CURNN_FILE), Network.params(env.curnn))
  serialize(joinpath(dir, MEM_FILE), get_experience(env))
  open(joinpath(dir, ITC_FILE), "w") do io
    JSON3.write(io, env.itc)
  end
end

function load_nn(location, e::Experiment)
  nn = e.mknet(e.gspec, e.netparams)
  params = deserialize(location)
  Network.set_params!(nn, params)
  return nn
end

function load_env(dir , e::Experiment)
  gspec = deserialize(joinpath(dir, GSPEC_FILE))
  params = deserialize(joinpath(dir, PARAMS_FILE))

  curnn = load_nn(joinpath(dir, CURNN_FILE), e)
  bestnn = load_nn(joinpath(dir, BESTNN_FILE), e)
  experience = deserialize(joinpath(dir, MEM_FILE))
  itc = open(JSON3.read, joinpath(dir, ITC_FILE), "r")
  return Env(gspec, params, curnn, bestnn, experience, itc)
end

#####
##### Save and load session reports
#####

function load_session_report(dir::String, niters)
  rep = SessionReport()
  for itc in 0:niters
    idir = iterdir(dir, itc)
    ifile = joinpath(idir, REPORT_FILE)
    bfile = joinpath(idir, BENCHMARK_FILE)
    # Load the benchmark report
    open(bfile, "r") do io
      push!(rep.benchmark, JSON3.read(io, Report.Benchmark))
    end
    # Load the iteration report
    if itc > 0
      open(ifile, "r") do io
        push!(rep.iterations, JSON3.read(io, Report.Iteration))
      end
    end
  end
  @assert valid_session_report(rep)
  return rep
end

function save_report_increment(session, bench, itrep, idir)
  open(joinpath(idir, BENCHMARK_FILE), "w") do io
    JSON3.pretty(io, JSON3.write(bench))
  end
  if session.env.itc > 0
    @assert !isnothing(itrep)
    open(joinpath(idir, REPORT_FILE), "w") do io
      JSON3.pretty(io, JSON3.write(itrep))
    end
  end
end

#####
##### Incremental saving
#####

autosave_enabled(s::Session) = s.autosave && !isnothing(s.dir)

function save_increment!(session::Session, bench, itrep=nothing)
  push!(session.report.benchmark, bench)
  isnothing(itrep) || push!(session.report.iterations, itrep)
  itc = session.env.itc
  if autosave_enabled(session)
    idir = iterdir(session.dir, itc)
    isdir(idir) || mkpath(idir)
    # Save the environment state,
    # both at the root and in the last iteration folder
    save(session, session.dir)
    session.save_intermediate && save(session, idir)
    # Save the collected statistics
    save_report_increment(session, bench, itrep, idir)
    # Do the plotting
    params = session.env.params
    plotsdir = joinpath(session.dir, PLOTS_DIR)
    isnothing(itrep) || plot_iteration(itrep, params, plotsdir, itc)
    plot_training(params, session.report.iterations, plotsdir)
    plot_benchmark(params, session.report.benchmark, plotsdir)
  end
end

#####
##### Run benchmarks
#####

function show_space_after_progress_bar(logger)
  Log.console_only(logger) do
    Log.sep(logger, force=true)
  end
end

function run_duel(env::Env, duel; logger)
  if isa(duel, Benchmark.Duel)
    player_name = Benchmark.name(duel.player)
    baseline_name = Benchmark.name(duel.baseline)
    legend = "$player_name against $baseline_name"
  else
    @assert isa(duel, Benchmark.Single)
    legend = Benchmark.name(duel.player)
  end
  Log.section(logger, 2, "Running benchmark: $legend")
  progress = Log.Progress(logger, duel.sim.num_games[env.itc])
  report = Benchmark.run(env, duel, progress)
  show_space_after_progress_bar(logger)
  print_report(
    logger, report,
    ternary_rewards=env.params.ternary_rewards)
  return report
end

function run_benchmark(session::Session)
  report = Report.Benchmark()
  for duel in session.benchmark
    outcome = run_duel(session.env, duel, logger=session.logger)
    push!(report, outcome)
  end
  return report
end

# return whether or not critical problems were found
function zeroth_iteration!(session::Session)
  session.env.itc = 0
  Log.section(session.logger, 2, "Initial report")
  report = initial_report(session.env)
  print_report(session.logger, report)
  isempty(report.errors) || return false
  bench = run_benchmark(session)
  save_increment!(session, bench)
  return true
end

function missing_zeroth_iteration(session::Session)
  return isempty(session.report.iterations) && isempty(session.report.benchmark)
end

#####
##### Session constructors
#####

default_session_dir(e::Experiment) = default_session_dir(e.name)
default_session_dir(exp_name::String) = joinpath(DEFAULT_SESSIONS_DIR, exp_name)

function session_logger(dir, nostdout, autosave)
  if autosave
    isdir(dir) || mkpath(dir)
    logfile = open(joinpath(dir, LOG_FILE), "a")
  else
    logfile = devnull
  end
  out = nostdout ? devnull : stdout
  return Logger(out, logfile=logfile)
end

"""
    Session(::Experiment; <optional kwargs>)

Create a new session from an experiment.

# Optional keyword arguments
- `dir="sessions/<experiment-name>"`: session directory in which all files and reports
    are saved.
- `autosave=true`: if set to `false`, the session won't be saved automatically nor
    any file will be generated
- `nostdout=false`: disables logging on the standard output when set to `true`
- `save_intermediate=false`: if set to true (along with `autosave`), all
    intermediate training environments are saved on disk so that
    the whole training process can be analyzed later. This can
    consume a lot of disk space.
"""
function Session(
    e::Experiment;
    dir=nothing,
    autosave=true,
    nostdout=false,
    save_intermediate=false,
    sim_num_games = e.params.self_play.sim.num_games,
    sim_num_workers = e.params.self_play.sim.num_workers,
    sim_batch_size = e.params.self_play.sim.batch_size,
    use_gpu = e.params.self_play.sim.use_gpu,
    mcts_iters = e.params.self_play.mcts.num_iters_per_turn,
    mcts_cpuct = e.params.self_play.mcts.cpuct,
    bench_num_games = e.benchmark[1].sim.num_games,
    bench_num_workers = e.benchmark[1].sim.num_workers,
    bench_batch_size = e.benchmark[1].sim.batch_size)

  isnothing(dir) && (dir = default_session_dir(e))
  logger = session_logger(dir, nostdout, autosave)
  if valid_session_dir(dir)
    Log.section(logger, 1, "Loading environment from: $dir")
    env = load_env(dir, e)
    # The parameters must be unchanged
    same_json(x, y) = JSON3.write(x) == JSON3.write(y)
    same_json(env.params, e.params) || @info "Using modified parameters"
    @assert same_json(Network.hyperparams(env.bestnn), e.netparams)
    session = Session(env, dir, logger, autosave, save_intermediate, e.benchmark)
    session.report = load_session_report(dir, env.itc)
  else
    e = @set e.params.self_play.sim.num_games = sim_num_games
    e = @set e.params.self_play.sim.num_workers = sim_num_workers
    e = @set e.params.self_play.sim.batch_size = sim_batch_size
    e = @set e.params.self_play.sim.use_gpu = use_gpu
    e = @set e.params.learning.use_gpu = use_gpu
    e = @set e.params.self_play.mcts.num_iters_per_turn = mcts_iters
    e = @set e.params.self_play.mcts.cpuct = mcts_cpuct
    benchs =  Vector{Benchmark.Evaluation}()
    for bench in e.benchmark
      bench = @set bench.sim.num_games = bench_num_games
      bench = @set bench.sim.num_workers = bench_num_workers
      bench = @set bench.sim.batch_size = bench_batch_size
      push!(benchs, bench)
    end
    network = e.mknet(e.gspec, e.netparams)
    env = Env(e.gspec, e.params, network)
    session = Session(env, dir, logger, autosave, save_intermediate, e.benchmark)
    Log.section(session.logger, 1, "Initializing a new AlphaZero environment")
    save(session, session.dir)
    session.save_intermediate && save(session, idir)
  end
  return session
end

#####
##### Public interface
#####

"""
    resume!(session::Session)

Resume a previously created or loaded session. The user can interrupt training
by sending a SIGKILL signal.
"""
function resume!(session::Session)
  try
    if missing_zeroth_iteration(session)
      success = zeroth_iteration!(session)
      success || return
    end
    train!(session.env, session)
  catch e
    isa(e, InterruptException) || rethrow(e)
    Log.section(session.logger, 1, "Interrupted by the user")
  end
end

"""
    save(session::Session)

Save a session on disk.

This function is called automatically by [`resume!`](@ref) after each
training iteration if the session was created with `autosave=true`.
"""
function save(session::Session, dir=session.dir)
  save_env(session.env, dir)
end

"""
    explore(session::Session, [mcts_params, use_gpu])

Start an explorer session for the current environment.
"""
function explore(session::Session; args...)
  Log.section(session.logger, 1, "Starting interactive exploration")
  explore(AlphaZeroPlayer(session.env; args...), session.env.gspec)
end

AlphaZero.AlphaZeroPlayer(s::Session; args...) = AlphaZero.AlphaZeroPlayer(s.env; args...)

#####
##### Utilities for printing reports
#####

const NUM_COL = Log.ColType(7, x -> fmt(".4f", x))
const BIGINT_COL = Log.ColType(10, n -> format(ceil(Int, n), commas=true))

const LEARNING_STATUS_TABLE = Log.Table([
  ("Loss",   NUM_COL,     s -> s.loss.L),
  ("Lv",     NUM_COL,     s -> s.loss.Lv),
  ("Lp",     NUM_COL,     s -> s.loss.Lp),
  ("Lreg",   NUM_COL,     s -> s.loss.Lreg),
  ("Linv",   NUM_COL,     s -> s.loss.Linv),
  ("Hp",     NUM_COL,     s -> s.Hp),
  ("Hpnet",  NUM_COL,     s -> s.Hpnet)])

const SAMPLES_STATS_TABLE = Log.Table([
  ("Loss",   NUM_COL,     s -> s.status.loss.L),
  ("Lv",     NUM_COL,     s -> s.status.loss.Lv),
  ("Lp",     NUM_COL,     s -> s.status.loss.Lp),
  ("Lreg",   NUM_COL,     s -> s.status.loss.Lreg),
  ("Linv",   NUM_COL,     s -> s.status.loss.Linv),
  ("Hpnet",  NUM_COL,     s -> s.status.Hpnet),
  ("Hp",     NUM_COL,     s -> s.status.Hp),
  ("Wtot",   BIGINT_COL,  s -> s.Wtot),
  ("Nb",     BIGINT_COL,  s -> s.num_boards),
  ("Ns",     BIGINT_COL,  s -> s.num_samples)])

function print_report(logger::Logger, status::Report.LearningStatus; kw...)
  Log.table_row(logger, LEARNING_STATUS_TABLE, status; kw...)
end

function print_report(logger::Logger, stats::Report.Memory)
  content, styles, comments = [], [], []
  # All samples
  push!(content, stats.all_samples)
  push!(styles, Log.BOLD)
  push!(comments, ["all samples"])
  # Latest batch
  push!(content, stats.latest_batch)
  push!(styles, Log.BOLD)
  push!(comments, ["latest batch"])
  # Per game stage
  for stage in stats.per_game_stage
    minrem = stage.min_remaining_length
    maxrem = stage.max_remaining_length
    push!(content, stage.samples_stats)
    push!(styles, Log.NO_STYLE)
    push!(comments, ["$minrem to $maxrem turns left"])
  end
  Log.table(
    logger, SAMPLES_STATS_TABLE, content, styles=styles, comments=comments)
end

function print_report(logger::Logger, report::Report.SelfPlay)
  sspeed = format(round(Int, report.samples_gen_speed), commas=true)
  Log.print(logger, "Generating $(sspeed) samples per second on average")
  avgdepth = fmt(".1f", report.average_exploration_depth)
  Log.print(logger, "Average exploration depth: $avgdepth")
  memf = format(report.mcts_memory_footprint, autoscale=:metric, precision=2)
  Log.print(logger, "MCTS memory footprint per worker: $(memf)B")
  mems = format(report.memory_size, commas=true)
  memd = format(report.memory_num_distinct_boards, commas=true)
  Log.print(logger, "Experience buffer size: $(mems) ($(memd) distinct boards)")
end

function print_report(logger::Logger, report::Report.Initial)
  nnparams = format(report.num_network_parameters, commas=true)
  Log.print(logger, "Number of network parameters: $nnparams")
  nnregparams = format(report.num_network_regularized_parameters, commas=true)
  Log.print(logger, "Number of regularized network parameters: $nnregparams")
  mfpn = report.mcts_footprint_per_node
  Log.print(logger, "Memory footprint per MCTS node: $(mfpn) bytes")
  for w in report.warnings
    Log.print(logger, crayon"yellow", "Warning: ", w)
  end
  for e in report.errors
    Log.print(logger, crayon"red", "Error: ", e)
  end
end

percentage(x, total) = round(Int, 100 * (x / total))

function print_report(
    logger::Logger,
    report::Report.Evaluation;
    nn_replaced=false,
    ternary_rewards=false)

  r = fmt("+.2f", report.avgr)
  if ternary_rewards
    n = length(report.rewards)
    stats = Benchmark.TernaryOutcomeStatistics(report)
    pwon = percentage(stats.num_won, n)
    pdraw = percentage(stats.num_draw, n)
    plost = percentage(stats.num_lost, n)
    details = ["$pwon% won, $pdraw% draw, $plost% lost"]
  else
    wr = round(Int, 100 * (report.avgr + 1) / 2)
    details = []
  end
  if nn_replaced
     push!(details, "network replaced")
  end
  details = isempty(details) ? "" : " (" * join(details, ", ") * ")"
  red = fmt(".1f", 100 * report.redundancy)
  msg = "Average reward: $r$details, redundancy: $red%"
  Log.print(logger, msg)
end

function print_report(logger::Logger, report::Report.Checkpoint; ternary_rewards=false)
  print_report(
    logger, report.evaluation,
    nn_replaced=report.nn_replaced,
    ternary_rewards=ternary_rewards)
end

#####
##### Event handlers
#####

function Handlers.iteration_started(session::Session)
  i = session.env.itc + 1
  Log.section(session.logger, 1, "Starting iteration $i")
end

function Handlers.self_play_started(session::Session)
  ngames = session.env.params.self_play.sim.num_games[session.env.itc]
  Log.section(session.logger, 2, "Starting self-play")
  session.progress = Log.Progress(session.logger, ngames)
end

function Handlers.game_played(session::Session)
  next!(session.progress)
end

function Handlers.self_play_finished(session::Session, report)
  show_space_after_progress_bar(session.logger)
  print_report(session.logger, report)
  session.progress = nothing
end

function Handlers.memory_analyzed(session::Session, report)
  Log.section(session.logger, 2, "Memory Analysis")
  print_report(session.logger, report)
end

function Handlers.learning_started(session::Session)
  Log.section(session.logger, 2, "Starting learning")
end

function Handlers.updates_started(session::Session, status)
  Log.section(session.logger, 3, "Optimizing the loss")
  print_report(session.logger, status, style=Log.BOLD)
end

function Handlers.updates_finished(session::Session, report)
  print_report(session.logger, report)
end

function Handlers.checkpoint_started(session::Session)
  Log.section(session.logger, 3, "Launching a checkpoint evaluation")
  num_games = session.env.params.arena.sim.num_games
  # In single player games, each game has to be played twice (with both networks)
  n = GI.two_players(session.env.gspec) ? num_games : 2 * num_games
  session.progress = Log.Progress(session.logger, n)
end

function Handlers.checkpoint_game_played(session::Session)
  next!(session.progress)
end

function Handlers.checkpoint_finished(session::Session, report)
  show_space_after_progress_bar(session.logger)
  ternary_rewards = session.env.params.ternary_rewards
  print_report(session.logger, report, ternary_rewards=ternary_rewards)
end

function Handlers.learning_finished(session::Session, report)
  return
end

function Handlers.iteration_finished(session::Session, report)
  bench = run_benchmark(session)
  save_increment!(session, bench, report)
  flush(Log.logfile(session.logger))
end

function Handlers.training_finished(session::Session)
  Log.section(session.logger, 1, "Training completed")
  close(Log.logfile(session.logger))
end

#####
##### Launch new experiments on an existing session
#####

function run_duel(session_dir::String, duel; logger)
  env = load_env(session_dir)
  run_duel(env, duel; logger)
end

function run_new_benchmark(session_dir, name, benchmark; logger, itcmax=nothing)
  outdir = joinpath(session_dir, name)
  isdir(outdir) || mkpath(outdir)
  logger = Logger()
  Log.section(logger, 1, "Computing benchmark")
  itc = 0
  reports = Report.Benchmark[]
  while valid_session_dir(iterdir(session_dir, itc))
    !isnothing(itcmax) && itc > itcmax && break
    Log.section(logger, 1, "Iteration: $itc")
    itdir = iterdir(session_dir, itc)
    env = load_env(itdir)
    report = [run_duel(env, duel; logger) for duel in benchmark]
    push!(reports, report)
    # Save the intermediate reports
    open(joinpath(outdir, BENCHMARK_FILE), "w") do io
      JSON3.pretty(io, JSON3.write(reports))
    end
    plot_benchmark(env.params, reports, outdir)
    itc += 1
  end
end

function regenerate_plots(session::Session)
  plotsdir = joinpath(session.dir, PLOTS_DIR)
  plot_training(session.env.params, session.report.iterations, plotsdir)
  plot_benchmark(session.env.params, session.report.benchmark, plotsdir)
  for (itc, itrep) in enumerate(session.report.iterations)
    plot_iteration(itrep, session.env.params, plotsdir, itc)
  end
end


text:
using Distributions
using Klara

CURRENTDIR, CURRENTFILE = splitdir(@__FILE__)
ROOTDIR = splitdir(splitdir(CURRENTDIR)[1])[1]
OUTDIR = joinpath(ROOTDIR, "output")

# OUTDIR = "../output"

SUBOUTDIR = "SMMALA"

nchains = 10
nmcmc = 110000
nburnin = 10000

function C(n::Int, c::Float64)
  X = eye(n)
  [(j <= n-i) ? X[i+j, i] = X[i, i+j] = c^j : nothing for i = 1:(n-1), j = 1:(n-1)]
  X
end

n = 20
μ = zeros(n)
Σ = C(n, 0.9)
ν = 30.

Σt = (ν-2)*Σ/ν
Σtinv = inv(Σt)

plogtarget(p::Vector, v::Vector) = logpdf(MvTDist(ν, zeros(n), (ν-2)*Σ/ν), p)

p = BasicContMuvParameter(:p, logtarget=plogtarget, nkeys=1, diffopts=DiffOptions(mode=:forward, order=2))

model = likelihood_model([p], isindexed=false)

sampler = SMMALA(0.02, H -> softabs(H, 1000.))

mcrange = BasicMCRange(nsteps=nmcmc, burnin=nburnin)

mctuner = AcceptanceRateMCTuner(0.7, score=x -> logistic_rate_score(x, 3.), verbose=false)

outopts = Dict{Symbol, Any}(:monitor=>[:value], :diagnostics=>[:accept])

times = Array{Float64}(nchains)
stepsizes = Array{Float64}(nchains)
i = 1

while i <= nchains
  v0 = Dict(:p=>rand(Normal(0, 2), n))

  job = BasicMCJob(model, sampler, mcrange, v0, tuner=mctuner, outopts=outopts)

  tic()
  run(job)
  runtime = toc()

  chain = output(job)
  ratio = acceptance(chain)

  if 0.55 < ratio < 0.85
    writedlm(joinpath(OUTDIR, SUBOUTDIR, "chain"*lpad(string(i), 2, 0)*".csv"), chain.value, ',')
    writedlm(joinpath(OUTDIR, SUBOUTDIR, "diagnostics"*lpad(string(i), 2, 0)*".csv"), vec(chain.diagnosticvalues), ',')

    times[i] = runtime
    stepsizes[i] = job.sstate.tune.step

    println("Iteration ", i, " of ", nchains, " completed with acceptance ratio ", ratio)
    i += 1
  end
end

writedlm(joinpath(OUTDIR, SUBOUTDIR, "times.csv"), times, ',')
writedlm(joinpath(OUTDIR, SUBOUTDIR, "stepsizes.csv"), stepsizes, ',')


text:
### Compare LombScargle.jl with Astropy

using LombScargle, Base.Test, PyCall
PyCall.@pyimport astropy.stats as ast

ntimes = 401
t = linspace(0.01, 10pi, ntimes)
t += step(t)*rand(ntimes)
for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
    s = f(t)
    # "psd" normalization in LombScargle slightly differ from that of
    # Astropy and the test would fail if we includ it.
    for fitmean in (true, false), nrm in ("standard", "model", "log")
        f_jl, p_jl = freqpower(lombscargle(t, s, fit_mean = fitmean,
                                           normalization=nrm,
                                           maximum_frequency=20))
        f_py, p_py =
            ast.LombScargle(t, s,
                            fit_mean = fitmean)[:autopower](method="cython",
                                                            normalization=nrm,
                                                            maximum_frequency=20)
        @test_approx_eq f_jl f_py
        @test_approx_eq p_jl p_py
    end
end

# Test the fast method with evenly spaced data.
t = linspace(0.01, 10pi, ntimes)
for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
    s = f(t)
    # "psd" normalization in LombScargle slightly differ from that of
    # Astropy and the test would fail if we includ it.
    for fitmean in (true, false), nrm in ("standard", "model", "log")
        f_jl, p_jl = freqpower(lombscargle(t, s, fit_mean = fitmean,
                                           normalization=nrm,
                                           maximum_frequency=10,
                                           samples_per_peak=10))
        f_py, p_py =
            ast.LombScargle(t, s,
                            fit_mean = fitmean)[:autopower](method="fast",
                                                            normalization=nrm,
                                                            maximum_frequency=10,
                                                            samples_per_peak=10)
        @test_approx_eq f_jl f_py
        @test_approx_eq p_jl p_py
    end
end

# Test the model functions
for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
    s = f(t)
    for fm in (true, false), cd in (true, false)
        m_jl = LombScargle.model(t, s, 1/2pi, fit_mean=fm, center_data=cd)
        m_py = ast.LombScargle(t, s, center_data=cd, fit_mean=fm)[:model](t, 1/2pi)
        @test_approx_eq m_jl m_py
    end
end

info("LombScargle.jl is consistent with Astropy")


text:
push!(LOAD_PATH, dirname(dirname(@__DIR__)))
using Documenter, MyMCMCSampler
makedocs(sitename="MyMCMCSampler Documentation")


text:
using MPI
using PencilArrays
using LinearAlgebra: transpose!
using CUDA

MPI.Init()
comm = MPI.COMM_WORLD       # we assume MPI.Comm_size(comm) == 12
rank = MPI.Comm_rank(comm)  # rank of local process, in 0:11

# Define MPI Cartesian topology: distribute processes on a 3×4 grid.
topology = MPITopology(comm, (3, 4))

# Let's decompose 3D arrays along dimensions (2, 3).
# This corresponds to the "x-pencil" configuration in the figure.
# This configuration is described by a Pencil object.
dims_global = (42, 31, 29)  # global dimensions of the array
decomp_dims = (2, 3)
pen_x = Pencil(topology, dims_global, decomp_dims)
s_local = size_local(pen_x)
# We can now allocate distributed arrays in the x-pencil configuration.
Ax = PencilArray{Float64}(undef, pen_x)
AxC = PencilArray(pen_x, CuArray{Float64}(undef, s_local ))
fill!(Ax,  rank * π)  # each process locally fills its part of the array
fill!(AxC, rank * π)  # each process locally fills its part of the array

MPI.Barrier(comm)

if rank == 0
 println( "Ax:  ", size(Ax)          )    # size of local part
 println( "Ax:  ", size_global(Ax)   )    # total size of the array = (42, 31, 29)
 println( "AxC: ", size(AxC)         )    # size of local part
 println( "AxC: ", size_global(AxC)  )    # total size of the array = (42, 31, 29)
end

# Create another pencil configuration, decomposing along dimensions (1, 3).
# We could use the same constructor as before, but it's recommended to reuse the
# previous Pencil instead.
pen_y = Pencil(pen_x, decomp_dims=(1, 3), permute = Permutation(2, 1, 3))
# pen_y = Pencil(pen_x, decomp_dims=(1, 3) )

# Now transpose from the x-pencil to the y-pencil configuration, redistributing
# the data initially in Ax.
Ay = PencilArray{Float64}(undef, pen_y)
s_local = size_local(pen_y, MemoryOrder() )
AyC = PencilArray(pen_y, CuArray{Float64}(undef, s_local ))

MPI.Barrier(comm)

if rank == 0
 println( "Ay:  ", size(Ay)          )    # size of local part
 println( "Ay:  ", size_global(Ay)   )    # total size of the array = (42, 31, 29)
 println( "AyC: ", size(AyC)         )    # size of local part
 println( "AyC: ", size_global(AyC)  )    # total size of the array = (42, 31, 29)
end

# Create another pencil configuration, decomposing along dimensions (1, 3).

transpose!(Ay, Ax)
transpose!(AyC, AxC)

localmaxAx=maximum(@. abs(AxC - Ax) )
localmaxAy=maximum(@. abs(AyC - Ay) )

maxAx=MPI.Allreduce( localmaxAx, +, comm )
maxAy=MPI.Allreduce( localmaxAy, +, comm )

MPI.Barrier(comm)

if rank == 0
 println("maxAx: ", maxAx)
 println("maxAy: ", maxAy)
end

# Test gather on transposed arrays
AyG=gather(Ay)
Ay.=AyC
if rank == 0
 println("typeof(Ay): ",typeof(Ay))
 println("typeof(AyC): ",typeof(AyC))
end
AyCG=gather(AyC)

if rank == 0
 println("typeof(AyC): ",typeof(AyC))
 println("typeof(AyCG): ",typeof(AyCG))
 println("size(AyC): ",size(AyC))
 println("size(AyCG): ",size(AyCG))
 println("maximum(AyCG): ",maximum(AyCG))
end

## if rank == 0
##  println( maximum(@. abs(AyG - AyCG) ) )
## end


text:
using NCBITaxonomy
using Test

my_tests = [
   "taxon.jl",
   "exceptions.jl",
   "lineages.jl",
   "rank.jl",
   "namefinders.jl",
   "divisions.jl"
]

global test_n
global anyerrors

test_n = 1
anyerrors = false

for my_test in my_tests
  try
    include(my_test)
    println("[TEST $(lpad(test_n,2))] \033[1m\033[32mPASS\033[0m $(my_test)")
  catch e
    global anyerrors = true
    println("[TEST $(lpad(test_n,2))] \033[1m\033[31mFAIL\033[0m $(my_test)")
    showerror(stdout, e, backtrace())
    println()
    throw("TEST FAILED")
  end
  global test_n += 1
end

if anyerrors
  throw("Tests failed")
end


text:
module TestDoctests

using DataUtils
using Documenter
using Test

@testset "docstest" begin
    doctest(DataUtils; manual=false)
end

end  # module


text:
"""
  hexlattice(width, height, r)

Build a 2D simplex from a hexagonal grid.

The grid will cover the given width and height. Each edge in the grid
will be of length `r`. This returns a tuple with vertices as a Float64
array of `(2, n)` and triangles as an Int64 array of `(3, m)`, where
triangle indices point into the vertices.
"""
function hexlattice(width, height, r)
    # Us a Bravais lattice.
    # Two vectors: r1, r2, of the same length.
    # The angle between them is 120 degrees.
    #
    # a1 = (r, 0)
    # a2 = (-r/sqrt(3), 2r/sqrt(3))
    # x,y = i a1 + j a2
    # x = ri - rj/2
    # y = 2jr/√3
    # m and n are the row count (x) and column count (y).
    m, n = Int64.(ceil.([width / r, height * sqrt(3) / 2]))
    # This is all about making a square region from a Bravais
    # lattice that is tilted. We do that by coming up with a rule
    # for which entries are saved.
    isvertex((i, j)) = (0 <= (i-1) - (j-1) ÷ 2 < m) && (0 <= j-1 < n)

    v = Dict{Tuple{Int64,Int64},Int64}()
    vidx = 1
    for j in 1:n
        ioffset = (j - 1) ÷ 2
        for i in (1 + ioffset):(m + ioffset)
            @assert isvertex((i, j))
            v[(i, j)] = vidx
            vidx += 1
        end
    end
    vertex = zeros(Float64, 2, length(v))
    for ((i, j), idx) in v
        vertex[:, idx] = [r * (i-1) - r * (j-1) / 2, 2 * r * (j-1) / 3]
    end

    # This is too many but we can trim at the end, faster than adding to array.
    triangles_pad = zeros(Int64, 3, length(v))
    triangle_cnt = 0
    for (i, j) in keys(v)
        # Each unit of the Bravais lattice is a parallelogram, so two triangles.
        triangle = [
                [(i, j), (i + 1, j), (i + 1, j + 1)]
                [(i, j), (i + 1, j + 1), (i, j + 1)]
                ]
        for t in 1:2
            if all(isvertex.(triangle[5]))
                triangle_cnt += 1
                triangles_pad[:, triangle_cnt] = [v[t] for t in triangle[t]]
            end
        end
    end
    triangles = triangles_pad[:, 1:triangle_cnt]
    vertex, triangles
end


text:
include("../src/processing/instance_generators.jl")

args = Args.parse_commandline()
Args.set(:level, 1)

# inst = Instance(Args.get(:infile))
print("Génération d'une instance... ")

inst = instance_build_mini10()

if lg2()
    println("\n", to_s_alpx(inst))
    print("Génération d'une instance... ")
end
print("=> fait\n")

print("Statistiques sur l'instance $(inst.name)... ")
@test inst.nb_planes == length(inst.planes)
@test inst.nb_planes == 10
@test inst.nb_kinds == 2
@test inst.name == "mini10"
print("=> fait\n")


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

__precompile__()
module StochasticPrograms

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using JuMP
using Compat
using OrderedCollections
using MathOptInterface
import Distributions: mean, std, quantile, ProbabilityWeights, TDist, Normal
using MacroTools
using MacroTools: @q, postwalk, prewalk
using Reexport
using ProgressMeter

import Base: getindex, length, in, issubset, show
import JuMP: optimize!, termination_status, index, value

const DenseAxisArray = JuMP.Containers.DenseAxisArray
const SparseAxisArray = JuMP.Containers.SparseAxisArray
const VectorizedProductIterator = JuMP.Containers.VectorizedProductIterator

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex
const AcceptableTermination = [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED]

import MutableArithmetics
const MA = MutableArithmetics

export
    StochasticModel,
    StochasticProgram,
    Probability,
    AbstractScenario,
    Scenario,
    StageParameters,
    parameters,
    AbstractSampler,
    Sampler,
    SampleSize,
    Confidence,
    ConfidenceInterval,
    StochasticSolution,
    AbstractStructuredOptimizerAttribute,
    AbstractSampledOptimizerAttribute,
    AbsoluteTolerance,
    RelativeTolerance,
    MasterOptimizer,
    MasterOptimizerAttribute,
    RawMasterOptimizerParameter,
    SubProblemOptimizer,
    SubProblemOptimizerAttribute,
    RawSubProblemOptimizerParameter,
    InstanceOptimizer,
    InstanceCrash,
    InstanceOptimizerAttribute,
    RawInstanceOptimizerParameter,
    NumSamples,
    NumEvalSamples,
    NumEWSSamples,
    NumEEVSamples,
    NumUpperTrials,
    NumLowerTrials,
    Execution,
    ExecutionParameter,
    RawExecutionParameter,
    set_execution_attribute,
    set_execution_attributes,
    Serial,
    Synchronous,
    Asynchronous,
    Quadratic,
    Linearized,
    InfNorm,
    ManhattanNorm,
    UnsupportedStructure,
    UnloadedStructure,
    UnloadableStructure,
    AbstractStructuredOptimizer,
    AbstractSampledOptimizer,
    optimizer_name,
    optimizer,
    optimizer_constructor,
    cache_solution!,
    AbstractStochasticStructure,
    StochasticInstantiation,
    UnspecifiedInstantiation,
    Deterministic,
    DeterministicEquivalent,
    StageDecomposition,
    Vertical,
    StageDecompositionStructure,
    ScenarioDecomposition,
    Horizontal,
    ScenarioDecompositionStructure,
    DistributedStageDecomposition,
    DistributedVertical,
    DistributedScenarioDecomposition,
    DistributedHorizontal,
    SampleAverageApproximation,
    Crash,
    distributed,
    initialize!,
    initialized,
    deferred,
    internal_model,
    add_scenario!,
    add_scenarios!,
    scenarioproblems,
    Decision,
    Decisions,
    SingleDecision,
    VectorOfDecisions,
    AffineDecisionFunction,
    QuadraticDecisionFunction,
    VectorAffineDecisionFunction,
    DecisionVariable,
    SPConstraintRef,
    ScenarioDependentModelAttribute,
    ScenarioDependentVariableAttribute,
    ScenarioDependentConstraintAttribute,
    DecisionRef,
    DecisionAffExpr,
    DecisionQuadExpr,
    decision,
    state,
    decisions,
    num_decisions,
    decision_names,
    all_decision_variables,
    all_known_decision_variables,
    all_auxiliary_variables,
    decision_by_name,
    recourse_length,
    first_stage,
    num_stages,
    scenario,
    scenarios,
    scenario_type,
    scenario_types,
    scenariotext,
    stage_parameters,
    stage_probability,
    probability,
    ExpectedScenario,
    expected,
    has_generator,
    generator,
    subproblem,
    subproblems,
    num_subproblems,
    masterterms,
    transfer_model!,
    num_scenarios,
    AbstractSampler,
    sample,
    sampler,
    sample!,
    generate!,
    stage_one_model,
    stage_two_model,
    outcome_model,
    evaluate_decision,
    recourse_decision,
    lower_bound,
    lower,
    upper,
    confidence,
    confidence_interval,
    lower_confidence_interval,
    upper_confidence_interval,
    gap,
    instantiate,
    set_optimizer,
    master_optimizer,
    set_masteroptimizer_attribute,
    set_masteroptimizer_attributes,
    subproblem_optimizer,
    set_suboptimizer_attribute,
    set_suboptimizer_attributes,
    set_instanceoptimizer_attribute,
    set_instanceoptimizer_attributes,
    optimize!,
    num_iterations,
    structure,
    default_structure,
    supports_structure,
    check_loadable,
    load_structure!,
    restore_structure!,
    load_model!,
    optimal_instance,
    decision,
    optimal_decision,
    optimal_recourse_decision,
    @define_scenario,
    @sampler,
    @first_stage,
    @second_stage,
    @stage,
    @stochastic_model,
    @zero,
    @expectation,
    @sample,
    @decision,
    @recourse,
    @known,
    @parameters,
    @uncertain,
    @scenario,
    WS,
    wait_and_see_decision,
    EWS,
    DEP,
    VRP,
    EVPI,
    EVP,
    expected_value_decision,
    EV,
    EEV,
    VSS,
    LShaped

macro exportJuMP()
    Expr(:export, names(JuMP)...)
end
@exportJuMP

include("types/types.jl")
include("methods/methods.jl")
include("macros/macros.jl")
include("optimizer_interface.jl")
include("crash.jl")
include("solvers/solvers.jl")
include("io/io.jl")

end # module


text:
using ConstrainedDynamics
using ConstrainedDynamics: vrotate, Vmat, GenericJoint, AbstractBody, orthogonalrows, deepcopy
using LinearAlgebra

function Fixed_abstract(body1::AbstractBody{T}, body2; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T
    vertices = (p1, p2)
    function gfixed(joint::GenericJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
       G= [
           vrotate(xb + vrotate(vertices[2], qb) - (xa + vrotate(vertices[1], qa)), inv(qa));
           Vmat(qa \ qb / qoffset)
       ]

       return G
   end
   function gfixed(joint::GenericJoint, xb::AbstractVector, qb::UnitQuaternion)
       G= [
           xb + vrotate(vertices[2], qb) - vertices[1];
           Vmat(qb / qoffset)
       ]

       return G
   end

   return GenericJoint{6}(body1, body2, gfixed)
end
function Revolute_abstract(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T
    vertices = (p1, p2)
    V1, V2, V3 = orthogonalrows(axis)
    V12 = [V1;V2]

    function grevolute(joint::GenericJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
       G= [
           vrotate(xb + vrotate(vertices[2], qb) - (xa + vrotate(vertices[1], qa)), inv(qa));
           V12*Vmat(qa \ qb / qoffset)
       ]

       return G
   end
   function grevolute(joint::GenericJoint, xb::AbstractVector, qb::UnitQuaternion)
       G= [
           xb + vrotate(vertices[2], qb) - vertices[1];
           V12*Vmat(qb / qoffset)
       ]

       return G
   end

   return GenericJoint{5}(body1, body2,grevolute)
end

joint_axis = [1.0;0.0;0.0]
l = 1.0
w, d = .1, .1

vert11 = [0.;0.;l / 2]
vert12 = -vert11
vert21 = [0.;0.;l / 2]

origin = Origin{Float64}()
origin_abstract = deepcopy(origin)
link1 = Box(w, d, l, l)
link1_abstract = deepcopy(link1)
link2 = deepcopy(link1)
link2_abstract = deepcopy(link1)

joint0to1 = EqualityConstraint(Fixed(origin, link1; p2=vert11))
joint1to2 = EqualityConstraint(Revolute(link1, link2,joint_axis; p1=vert12, p2=vert21))
joint0to1_abstract = EqualityConstraint(Fixed_abstract(origin_abstract, link1_abstract; p2=vert11))
joint1to2_abstract = EqualityConstraint(Revolute_abstract(link1_abstract, link2_abstract,joint_axis; p1=vert12, p2=vert21))

links = [link1;link2]
links_abstract = [link1_abstract;link2_abstract]
constraints = [joint0to1;joint1to2]
constraints_abstract = [joint0to1_abstract;joint1to2_abstract]

mech = Mechanism(origin, links, constraints)
mech_abstract= Mechanism(origin_abstract, links_abstract, constraints_abstract)

setPosition!(origin,link1,p2 = vert11)
setPosition!(origin_abstract,link1_abstract,p2 = vert11)

function compare(s1,s2)
    @test all(isapprox.(norm.(s1.x[1]-s2.x[1]), 0.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.q[1]./s2.q[1]), 1.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.x[2]-s2.x[2]), 0.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.q[2]./s2.q[2]), 1.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.ω[1]-s2.ω[1]), 0.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.v[1]-s2.v[1]), 0.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.ω[2]-s2.ω[2]), 0.0; atol = 1e-8)) &&
        all(isapprox.(norm.(s1.v[2]-s2.v[2]), 0.0; atol = 1e-8))
end

for i=1:10
    randang = pi/4*rand()
    setPosition!(link1,link2,p1 = vert12,p2 = vert21,Δq = UnitQuaternion(RotX(randang)))
    setPosition!(link1_abstract,link2_abstract,p1 = vert12,p2 = vert21,Δq = UnitQuaternion(RotX(randang)))
    storage = simulate!(mech, 1, record = true)
    storage_abstract = simulate!(mech_abstract, 1, record = true)

    compare(storage,storage_abstract)
end


text:
using .lib

module docs
	export help,version
	import ...clean

	function help()
		print(clean("""

			 使い方:
			  measure-jl [options] [command] [arg1] [arg2]…
			  measure-jl -multiple [options] "[command1]" "[command2]"…

			  [command] を実行し,最後にその所要時間を表示します

			  measure --version : バージョン表示
			  measure --help    : ヘルプを表示

			  オプション

			   -o,-out,-stdout [string]
			   -e,-err,-stderr [string]
			    標準出力,標準エラー出力の出力先を指定します
			    指定しなければ inherit になります
			    • inherit
			     stdoutはstdoutに,stderrはstderrにそれぞれ出力します
			    • discard
			     出力しません
			    • [file path]
			     指定したファイルに書き出します (追記)

			   -r,-result [string]
			    実行結果の出力先を指定します
			    指定しなければ stderr になります
			    • stdout,stderr
			    • [file path]
			     指定したファイルに書き出します (追記)

			   -m,-multiple [string?]
			    複数のコマンドを実行します。通常はシェル経由で実行されます。
			    例えば measure echo 1 と指定していたのを

			     measure -multiple "echo 1" "echo 2"

			    などと1つ1つのコマンドを1つの文字列として渡して実行します
				引数に次のいずれかの値を指定することができます (指定しなければserial)
				• none
				 単一のコマンドとして実行します (-mを指定しない場合と同じ)
				• serial
				 指定した複数のコマンドをその順に実行していきます
				• spawn,parallel
				 シェルの同時起動により並列実行します
				• thread
				 スレッドを利用して並列実行します

		"""))
		exit(0)
	end

	function version()
		print(clean("""

			 measure v2.4
			 Julia エディション (measure-jl)

		"""))
		exit(0)
	end

end

max_stars_repo_path: Time Measurement/Julia/docs.jl

max_issues_repo_path: Time Measurement/Julia/docs.jl

max_forks_repo_path: Time Measurement/Julia/docs.jl