text:
using Completion: critical_pairs, add_rule!, size, orient!, choose!, complete
using Rewrite: PatternRule

@testset "critical_pairs" begin
    @test critical_pairs(
        PatternRule{Term}(@term(f(f(x, y), z)), @term(f(x, f(y, z)))),
        PatternRule{Term}(@term(f(f(x, y), z)), @term(f(x, f(y, z)))),
    ) == [
        (@term(f(x, f(y, z))), @term(f(x, f(y, z)))),
        (@term(f(f(x, y), f(z, z₁))), @term(f(f(x, f(y, z)), z₁))),
    ]

    @test critical_pairs(
        PatternRule{Term}(@term(f(f(x, y), z)), @term(f(x, f(y, z)))),
        PatternRule{Term}(@term(f(g(x), x)), @term(h())),
    ) == [
        (@term(f(g(x), f(x, z₁))), @term(f(h(), z₁)))
    ]

    @test critical_pairs(
        PatternRule{Term}(@term(f(x, g(x))), @term(h())),
        PatternRule{Term}(@term(f(h(), x)), @term(x)),
    ) == [
        (@term(h()), @term(g(h())))
    ]
end

@testset "add_rule!" begin
    es, ss, rs = [], TRS(), TRS()
    @test add_rule!(PatternRule{Term}(@term(f()), @term(g())), es, ss, rs) ==
        ([], TRS(@term(f()) => @term(g())), TRS())
    @test isempty(es)
    @test !isempty(ss) && ss == TRS(@term(f()) => @term(g()))
    @test isempty(rs)

    @test add_rule!(PatternRule{Term}(@term(f(x, y)), @term(g(x, y))), [], [PatternRule{Term}(@term(f(x, x)), @term(x))], []) ==
        ([(@term(g(x, x)), @term(x))], [PatternRule{Term}(@term(f(x, y)), @term(g(x, y)))], [])
    @test add_rule!(PatternRule{Term}(@term(f(x, y)), @term(g(x, y))), [], [], Pair[@term(f(x, x)) => @term(x)]) ==
        ([(@term(g(x, x)), @term(x))], [PatternRule{Term}(@term(f(x, y)), @term(g(x, y)))], [])
end

@testset "orient!" begin
    a >ᵣ b = size(a) > size(b)

    @test orient!(>ᵣ, [(@term(f(x)), @term(f(x)))], TRS(), TRS(@term(g(x)) => @term(x))) ==
        (TRS(), TRS(@term(g(x)) => @term(x)))
    @test orient!(>ᵣ, [(@term(f(x)), @term(x))], TRS(), TRS(@term(f(x)) => @term(x))) ==
        (TRS(), TRS(@term(f(x)) => @term(x)))
    @test orient!(>ᵣ, [(@term(f(x, x)), @term(x))], TRS(), TRS(@term(f(x, y)) => @term(g(x)))) ==
        (TRS(@term(g(x)) => @term(x)), TRS(@term(f(x, y)) => @term(x)))
    @test orient!(>ᵣ, [(@term(x), @term(f(x, x)))], TRS(), TRS(@term(f(x, y)) => @term(g(x)))) ==
        (TRS(@term(g(x)) => @term(x)), TRS(@term(f(x, y)) => @term(x)))
    @test (@test_logs (:warn, "Unable to determine preferred form") orient!(
        >ᵣ, [(@term(f(x, x)), @term(h(x, x)))], TRS(), TRS()
    )) === nothing
end

@testset "choose!" begin
    rs = TRS(@term(f(f(f(x)))) => @term(x), @term(g(x)) => @term(x))
    @test choose!(rs) == PatternRule{Term}(@term(g(x)), @term(x))
    @test rs == TRS(@term(f(f(f(x)))) => @term(x))

    rs = TRS(@term(f(x)) => @term(x), @term(g(x)) => @term(x))
    @test choose!(rs) == PatternRule{Term}(@term(f(x)), @term(x))

    rs = TRS(@term(g(x)) => @term(x), @term(f(x)) => @term(x))
    @test choose!(rs) == PatternRule{Term}(@term(g(x)), @term(x))

    @test_throws ArgumentError choose!([])
end

@testset "complete" begin
    >ᵣ = LPO((:i, 1), (:*, 2), (:e, 0))
    with_context(AlgebraContext()) do
        @test Set(complete(>ᵣ, @term AXIOMS [
            ((x * y) * z  , x * (y * z)  )
            (i(x) * x     , e()          )
            (e() * x      , x            )
        ])) == Set(@term RULES [
                (x * y) * z  =>  x * (y * z)
                    e() * x  =>  x
                    x * e()  =>  x
                   x * i(x)  =>  e()
                   i(x) * x  =>  e()
                     i(e())  =>  e()
                    i(i(x))  =>  x
                 i(x₁ * y₁)  =>  i(y₁) * i(x₁)
            x * (i(x) * z₃)  =>  z₃
            i(x) * (x * z₁)  =>  z₁
        ])
    end
end


text:
module SimulatedNeuralMoments

# the type that holds the model specifics
struct SNMmodel
    modelname::String # name of model
    lb # vector of lower bounds. Can be -Inf, if desired
    ub # vector of upper bounds. Can be inf, if desired
    insupport::Function # function that checks if a draw is valid
    priordraw::Function # function that returns a draw from prior
    auxstat::Function # function that returns an array of draws of statistic
end

function TransformStats(data, info)
    q01,q50,q99,iqr = info
    data .= max.(data, q01)
    data .= min.(data, q99)
    data .= (data .- q50) ./ iqr
end

# neural moments given statistic
function NeuralMoments(z, nnmodel, nninfo)
    nnmodel(TransformStats(z, nninfo))
end

# mean and covariance of output of NN, given draw from prior (not transformed)
# remember that the NN outputs estimates of transformed parameters
using LinearAlgebra
function mΣ(θ, reps, model::SNMmodel, nnmodel, nninfo)
    z = model.auxstat(θ, reps)
    Zs = [NeuralMoments(z[i], nnmodel, nninfo) for i = 1:reps]
    mean(Zs), Symmetric(cov(Zs))
end

using PrettyTables, Printf
function prettyprint(a, cnames="", rnames="")
if rnames !=""
    rnames = rnames[:]
    a = [rnames a]
    if cnames != ""
        cnames = cnames[:]
        cnames = vcat("", cnames)
    end
end
if cnames !=""
    pretty_table(a; header=(cnames), formatters=ft_printf("%12.5f"))
else
    pretty_table(a; formatters=ft_printf("%12.5f"))
end
end

function dstats(x, rnames="";short=false, silent=false)
    k = size(x,2)
    if rnames==""
        rnames = 1:k
        rnames = rnames'
    end
    m = mean(x,dims=1)
    mm = median(x,dims=1)
    s = std(x,dims=1)
    sk = m-m
    kt = m-m
    mn = minimum(x,dims=1)
    mx = maximum(x,dims=1)
    q05 = fill(0.0,k)
    q25 = fill(0.0,k)
    q75 = fill(0.0,k)
    q95 = fill(0.0,k)
    if short == false
        for i = 1:size(x,2) q05[i], q25[i], q75[i],q95[i] = quantile(x[:,i], [0.05,0.25,0.75,0.95]) end
        cnames = ["  mean", " median","  std", "IQR", "min", "max", "q05", "q95"]
        stats = [m' mm' s' (q75-q25) mn' mx' q05 q95]
        if !silent prettyprint(stats, cnames, rnames) end
    else
        cnames = ["  mean", " median", "  std", "min", "max"]
        stats = [m' mm' s' mn' mx']
        if !silent prettyprint(stats, cnames, rnames) end
    end
    return stats
end

# simulates data from prior, trains and tests the net, and returns
# the trained net and the information for transforming the inputs
using Statistics, Flux
using Base.Iterators
function MakeNeuralMoments(model::SNMmodel, transf; TrainTestSize=1, Epochs=1000)
    data = 0.0
    datadesign = 0.0
    nParams = size(model.lb,1)
    # training and testing
    if (TrainTestSize == 1) TrainTestSize = Int64(2*nParams*1e4); end # use a default size if none provided
    params = zeros(TrainTestSize,nParams)
    statistics = zeros(TrainTestSize,size(model.auxstat(model.lb,1)[1],1))
    Threads.@threads for s = 1:TrainTestSize
        ok = false
        θ = model.priordraw()
        W = (model.auxstat(θ,1))[1]
        # repeat draw if necessary
        while any(isnan.(W))
            θ = model.priordraw()
            W = model.auxstat(θ,1)[1]
        end
        params[s,:] = transf(θ)
        statistics[s,:] = W
    end
    # transform stats to robustify against outliers
    q50 = zeros(size(statistics,2))
    q01 = similar(q50)
    q99 = similar(q50)
    iqr = similar(q50)
    for i = 1:size(statistics,2)
        q = quantile(statistics[:,i],[0.01, 0.25, 0.5, 0.75, 0.99])
        q01[i] = q[1]
        q50[i] = q[3]
        q99[i] = q[5]
        iqr[i] = q[4] - q[2]
    end
    nninfo = (q01, q50, q99, iqr)
    for i = 1: size(statistics,1)
        statistics[i,:] .= TransformStats(statistics[i,:], nninfo)
    end
    # train net
    TrainingProportion = 0.5 # size of training/testing
    params = Float32.(params)
    s = Float32.(std(params, dims=1)')
    statistics = Float32.(statistics)
    trainsize = Int(TrainingProportion*TrainTestSize)
    yin = params[1:trainsize, :]'
    yout = params[trainsize+1:end, :]'
    xin = statistics[1:trainsize, :]'
    xout = statistics[trainsize+1:end, :]'
    # define the neural net
    nStats = size(xin,1)
    NNmodel = Chain(
        Dense(nStats, 10*nParams, tanh),
        Dense(10*nParams, 3*nParams, tanh),
        Dense(3*nParams, nParams)
    )
    loss(x,y) = Flux.huber_loss(NNmodel(x)./s, y./s; δ=0.1) # Define the loss function
    # monitor training
    function monitor(e)
        println("epoch $(lpad(e, 4)): (training) loss = $(round(loss(xin,yin); digits=4)) (testing) loss = $(round(loss(xout,yout); digits=4))| ")
    end
    # do the training
    bestsofar = 1.0e10
    pred = 0.0 # define it here to have it outside the for loop
    batches = [(xin[:,ind],yin[:,ind])  for ind in partition(1:size(yin,2), 50)]
    bestmodel = 0.0
    for i = 1:Epochs
        if i < 20
            opt = Momentum() # the optimizer
        else
            opt = ADAMW() # the optimizer
        end
        Flux.train!(loss, Flux.params(NNmodel), batches, opt)
        current = loss(xout,yout)
        if current < bestsofar
            bestsofar = current
            bestmodel = NNmodel
            xx = xout
            yy = yout
            println("________________________________________________________________________________________________")
            monitor(i)
            pred = NNmodel(xx)
            error = yy .- pred
            results = [pred;error]
            rmse = sqrt.(mean(error.^Float32(2.0),dims=2))
            println(" ")
            println("RMSE for model parameters ")
            prettyprint(reshape(round.(rmse,digits=3),1,nParams))
            println(" ")
            println("dstats prediction of parameters:")
            dstats(pred')
            println(" ")
            println("dstats prediction error:")
            dstats(error')
        end
    end
    bestmodel, nninfo
end

export SNMmodel, MakeNeuralMoments, NeuralMoments, mΣ
end

ext: jl

max_stars_repo_name: mcreel/SimulatedNeuralMoments

max_issues_repo_name: mcreel/SimulatedNeuralMoments

max_forks_repo_name: mcreel/SimulatedNeuralMoments


text:
module SeqUMAP

using Distances
using FASTX
import MultivariateStats: fit, PCA, transform
import UMAP: umap

include("encoding.jl")
include("embedding.jl")
include("projection.jl")
include("io.jl")

#encoding.jl...
export AA_DICT,
NT_DICT,
IUPACbool,
resolve_base,
resolve_seq,
string2encoding,

#embedding.jl...
collect_kmers,
kmer_count!,
get_kmer_index,
kmer_embed,

#projection.jl...
CorrectedKmer,
sequmap,
seqpca,

#io.jl
read_fasta,
read_fastq,
clean_nt,
clean_aa

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthenticationV1TokenReviewStatus <: SwaggerModel
    audiences::Any # spec type: Union{ Nothing, Vector{String} } # spec name: audiences
    authenticated::Any # spec type: Union{ Nothing, Bool } # spec name: authenticated
    error::Any # spec type: Union{ Nothing, String } # spec name: error
    user::Any # spec type: Union{ Nothing, IoK8sApiAuthenticationV1UserInfo } # spec name: user

    function IoK8sApiAuthenticationV1TokenReviewStatus(;audiences=nothing, authenticated=nothing, error=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("audiences"), audiences)
        setfield!(o, Symbol("audiences"), audiences)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("authenticated"), authenticated)
        setfield!(o, Symbol("authenticated"), authenticated)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiAuthenticationV1TokenReviewStatus

const _property_map_IoK8sApiAuthenticationV1TokenReviewStatus = Dict{Symbol,Symbol}(Symbol("audiences")=>Symbol("audiences"), Symbol("authenticated")=>Symbol("authenticated"), Symbol("error")=>Symbol("error"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiAuthenticationV1TokenReviewStatus = Dict{Symbol,String}(Symbol("audiences")=>"Vector{String}", Symbol("authenticated")=>"Bool", Symbol("error")=>"String", Symbol("user")=>"IoK8sApiAuthenticationV1UserInfo")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthenticationV1TokenReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthenticationV1TokenReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1TokenReviewStatus[property_name]

function check_required(o::IoK8sApiAuthenticationV1TokenReviewStatus)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, name::Symbol, val)
end


text:
function rawmode(f, hide_cursor=true)
    rawenabled = enableRawMode()
    rawenabled && hide_cursor && cursor_hide(terminal.out_stream)
    try
        f()
    finally
        rawenabled && disableRawMode(); cursor_show(terminal.out_stream)
    end
end

function enableRawMode()
    try
        REPL.Terminals.raw!(terminal, true)
        return true
    catch err
        warn("TerminalMenus: Unable to enter raw mode: $err")
    end
    return false
end

function disableRawMode()
    try
        REPL.Terminals.raw!(terminal, false)
        return true
    catch err
        warn("TerminalMenus: Unable to disable raw mode: $err")
    end
    return false
end


text:
#Version Julia "1.1.1"
using DataFrames
using MathProgBase
using CSV
using LinearAlgebra
################################################################################
## Setting-up directory
tempdir1=@__DIR__
repdir=tempdir1[1:findfirst("ReplicationAK",tempdir1)[end]]
diroutput=repdir*"/Output_all"
dirdata=repdir*"/Data_all"

## Data
# Sample size
const n=154
# Number of time periods
const T=50
# Number of goods
const K=3
# Read csv files
dum0=CSV.read(dirdata*"/rationalitydata3goods.csv",DataFrame)
# break the dataset into subdatasets by individual
splitdum0=groupby(dum0,:id)
@eval  splitp=$splitdum0
#Initialize array of effective prices
rho=zeros(n,T,K)
#Initialize array of consumption
cve=zeros(n,T,K)
# Fill the arrays
# Columns 10:12 correspond to prices
# Columns 4:6 correspond to consumption bundles
for i=1:n
    dum0=convert(Array,splitp[i])
    rho[i,:,:]=dum0[1:T,10:12]
    cve[i,:,:]=dum0[1:T,4:6]
end

################################################################################
## Deterministic
ind=1
p=rho[ind,:,:]
q=cve[ind,:,:]

function rdmatrix(;p=p,q=q)
    rdmatrixfill=zeros(T,T)
    for t=1:T, s=1:T
      if (p[t,:]'*q[t,:])[1]>=(p[t,:]'*q[s,:])[1]
        rdmatrixfill[t,s]=1
      end
    end
    rdmatrixfill
end

function tclosure(;m=m)
    n0=size(m)[2]
    for k=1:n0, i=1:n0, j=1:n0
      if  m[i,k]==1 && m[k,j]==1
        m[i,j]=1
      end
    end
    m
end

function garpok(;p=p,q=q)
    garpmatrixfill=zeros(T,T)
    m=tclosure(m=rdmatrix(p=p,q=q))
    for t=1:T
        for s=1:T
            if m[t,s]==1 && ((p[s,:]'*q[s,:])[1]>(p[s,:]'*q[t,:])[1])
                garpmatrixfill[t,s]=1
            end
        end
    end

    maximum(garpmatrixfill)
    ## output of zero is a pass, no violation of GARP
end

function testgarp(;rho=rho,cve=cve)
    garpresults=zeros(n)
    for ind=1:n
      garpresults[ind]=garpok(p=rho[ind,:,:],q=cve[ind,:,:])
    end
    garpresults
end
garpresults=testgarp(rho=rho,cve=cve)
rate=1-sum(garpresults)/n

## Pass Rate
ratev=["pass-garp" rate]
DFsolv=convert(DataFrame,ratev)
CSV.write(diroutput*"/2App_dt_rr.csv",DFsolv)


text:
import Composites
using Test

@testset "Textbook ABD" begin
    e1 = fill(148e9, 2)
    e2 = fill(9.65e9, 2)
    g12 = fill(4.55e9, 2)
    nu12 = fill(0.3, 2)
    theta = [0.0, 45.0]

    # Single value getQ
    qtest = [148.87e9  2.91e9  0.0;
               2.91e9  9.71e9  0.0;
               0.0     0.0     4.55e9]
    q = Composites.getQ(e1[1], e2[1], g12[1], nu12[1], theta[1])
    for i = 1:size(qtest, 1)
        for j = 1:size(qtest, 2)
            @test round(q[i,j]/1e9, digits=2) == qtest[i,j]/1e9
        end
    end

    # Single value getQ with angle
    qbartest = [45.65e9 36.55e9 34.79e9;
                36.55e9 45.65e9 34.79e9;
                34.79e9 34.79e9 38.19e9]
    qbar = Composites.getQ(e1[2], e2[2], g12[2], nu12[2], theta[2])
    for i = 1:size(qbartest, 1)
        for j = 1:size(qbartest, 2)
            @test round(qbar[i,j]/1e9, digits=2) == qbartest[i,j]/1e9
        end
    end

    # getABD test
    nply = [10, 10]
    tply = [0.0001, 0.0001]
    matid = [1, 1]
    theta = [0.0, 45.0]
    Atest = [194.52 39.46 34.79;
              39.46 55.36 34.79;
              34.79 34.79 42.74]
    Btest = [-51.61 16.82 17.40;
              16.82 17.97 17.40;
              17.40 17.40 16.82]
    Dtest = [64.84 13.15 11.60;
             13.15 18.45 11.60;
             11.60 11.60 14.25]

    q = Composites.getQ.(e1, e2, g12, nu12)
    A, B, D = Composites.getABD(matid, nply, tply, theta, q)
    for i = 1:size(qbartest, 1)
        for j = 1:size(qbartest, 2)
            @test round(A[i,j]/1e6, digits=2) == Atest[i,j]
            @test round(B[i,j]/1e3, digits=2) == Btest[i,j]
            @test round(D[i,j], digits=2) == Dtest[i,j]
        end
    end

    alfatest = [13.44  -4.85  -7.14;
                -4.85  41.81 -21.23;
                -7.14 -21.23  64.95]
    betatest = [17.07  -6.01 -11.06;
                -6.01  -5.04 -11.06;
               -11.06 -11.06 -24.05]
    deltatest = [40.32 -14.56 -21.41;
                -14.56 125.42 -63.68;
                -21.41 -63.68 194.86]

    alfa, beta, delta = Composites.compliance(A, B, D)
    for i = 1:size(qbartest, 1)
        for j = 1:size(qbartest, 2)
            @test round(alfa[i,j]*1e9, digits=2) == alfatest[i,j]
            @test round(beta[i,j]*1e6, digits=2) == betatest[i,j]
            @test round(delta[i,j]*1e3, digits=2) == deltatest[i,j]
        end
    end
end

### Classical Laminate Theory ###

# materials
e1 = [181e9, 203e9, 38.6e9, 76e9]
e2 = [10.3e9, 11.2e9, 8.27e9, 5.5e9]
g12 = [7.17e9, 8.4e9, 4.14e9, 2.3e9]
nu12 = [0.28, 0.32, 0.26, 0.34]
rho = zeros(Float64, 4)
xt = [1500.0, 3500.0, 1062.0, 1400.0]*1e6
xc = [1500.0, 1540.0, 610.0, 235.0]*1e6
yt = [40.0, 56.0, 31.0, 12.0]*1e6
yc = [246.0, 150.0, 118.0, 53.0]*1e6
s = [68.0, 98.0, 72.0, 34.0]*1e6
t = fill(0.0001, 4)
mat = Composites.Material.(e1, e2, g12, nu12, rho, xt, xc, yt, yc, s, t)

# Laminate
matid = [4, 4, 4, 4, 4, 4, 4]
nply = [5, 5, 5, 10, 5, 5, 5]
tply = [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001]
theta = [0.0, 45.0, -45.0, 90.0, -45.0, 45.0, 0.0]
lam = Composites.Laminate{Float64}(matid, nply, tply, theta)

# getz(tply::Array{Float64,1},nply::Array{Int64,1))
z = Composites.getz(tply, nply)
@testset "getz" begin
    ztest = [-2.0, -1.5, -1.0, -0.5, 0.5, 1.0, 1.5, 2.0]*1e-3
    for i = 1:length(z)
        @test isapprox(z[i], ztest[i])
    end
    # getz(lam::Laminate)
    @test Composites.getz(lam) == z
end

# getQ(e1::Array{Float64,1},e2::Array{Float64,1},g12::Array{Float64,1},
#   nu12::Array{Float64,1},theta::Array{Float64,1}=zeros(Float64,length(e1)))
q = Composites.getQ.(e1, e2, g12, nu12)
@testset "getQ" begin
    qtest1 = [181.811 2.897 0.0; 2.897 10.346 0.0; 0.0 0.0 7.17] #GPa
    for i = 1:3
        for j = 1:3
            @test round.(q[1][i,j]/1e9, digits=3) == qtest1[i,j]
        end
    end
    qtest3 = [39.167 2.182 0.0; 2.182  8.392 0.0; 0.0 0.0 4.14] #GPa
    for i = 1:3
        for j = 1:3
            @test round.(q[3][i,j]/1e9, digits=3) == qtest3[i,j]
        end
    end
    for i = 1:length(mat)
        # getQ(e1::Float64,e2::Float64,g12::Float64,nu12::Float64,theta::Float64=0.0)
        @test Composites.getQ(e1[i], e2[i], g12[i], nu12[i]) == q[i][:,:]
        # getQ(mat::material,theta::Float64)
        @test Composites.getQ(mat[i]) == q[i][:,:]
    end
    # getQ(mat::Array{material,1},theta::Array{Float64,1}
    # getQ(mat::Array{material,1},lam::Laminate)
    @test Composites.getQ.(mat) == q
end

# rotQ(q::Array{Float64,2},theta::Float64)
qbar = [Composites.rotQ(q[matid[i]], theta[i]) for i in 1:length(theta)]
@testset "rotQ" begin
    qbartest1 = [76.641 1.886 0.0; 1.886 5.546 0.0; 0.0 0.0 2.3] #GPa
    qbartest2 = [23.79 19.19 17.774; 19.19 23.79 17.774; 17.774 17.774 19.604] #GPa
    qbartest3 = [23.79 19.19 -17.774; 19.19 23.79 -17.774; -17.774 -17.774 19.604] #GPa
    qbartest4 = [5.546 1.886 0.0; 1.886 76.641 0.0; 0.0 0.0 2.3] #GPa
    for i = 1:3
        for j = 1:3
            @test round.(qbar[1][i,j]/1e9, digits=3) == qbartest1[i,j]
            @test round.(qbar[2][i,j]/1e9, digits=3) == qbartest2[i,j]
            @test round.(qbar[3][i,j]/1e9, digits=3) == qbartest3[i,j]
            @test round.(qbar[4][i,j]/1e9, digits=3) == qbartest4[i,j]
            @test round.(qbar[5][i,j]/1e9, digits=3) == qbartest3[i,j]
            @test round.(qbar[6][i,j]/1e9, digits=3) == qbartest2[i,j]
            @test round.(qbar[7][i,j]/1e9, digits=3) == qbartest1[i,j]
        end
    end
    # rotQ(q11::Float64,q12::Float64,q22::Float64,q66::Float64,theta::Float64)
    for i = 1:size(qbar, 3)
        @test Composites.rotQ(q[matid[i]][1,1],q[matid[i]][1,2],q[matid[i]][2,2],
            q[matid[i]][3,3], theta[i]) == qbar[i][:,:]
    end
end

# getABD(matid::Array{Int64,1},nply::Array{Int64,1},tply::Array{Float64,1},
#   theta::Array{Float64,1},q::Array{Float64,3})
A,B,D = Composites.getABD(matid, nply, tply, theta, q)
@testset "getABD" begin
    Atest = [1.298e8 4.215e7 0.0; 4.215e7 1.298e8 0.0; 0.0 0.0 4.381e7] #N
    Btest = zeros(3, 3)
    Dtest = [288.317 47.549 17.774; 47.549 75.033 17.774; 17.774 17.774 49.759] #N*m
    for i = 1:3
        for j = 1:3
            @test isapprox(round(A[i,j], sigdigits=4),Atest[i,j])
            @test round(B[i,j], digits=3) == Btest[i,j]
            @test round(D[i,j], digits=3) == Dtest[i,j]
        end
    end
    # getABD(lam::Laminate,q::Array{Float64,3})
    @test Composites.getABD(lam, q) == (A,B,D)
end

# compliance(A::Array{Float64,2},B::Array{Float64,2},D::Array{Float64,2})
a, b, d = Composites.compliance(A, B, D)
@testset "compliance" begin
    atest = [8.615e-9 -2.798e-9 0.0; -2.798e-9  8.615e-9 0.0;0.0 0.0 2.283e-8] #N
    btest = zeros(3,3)
    dtest = [3.887e-3 -2.332e-3 -5.556e-4; -2.332e-3  1.596e-2 -4.867e-3;
        -5.556e-4 -4.867e-3  2.203e-2] #N*m
    for i = 1:3
        for j = 1:3
            @test isapprox(round(a[i,j],sigdigits=4), atest[i,j], atol = eps(Float64))
            @test round(b[i,j], digits=3) == btest[i,j]
            @test isapprox(round(d[i,j],sigdigits=4), dtest[i,j])
        end
    end
end

@testset "rotstress" begin
    # rotstress(stress::Array{Float64,1},theta::Float64)
    stress = rand(3)
    tsigmatest = [0.5 0.5 1.0;0.5 0.5 -1.0;-0.5 0.5 0.0]*stress
    tsigma = Composites.rotstress(stress, 45.0)
    for i = 1:size(tsigma, 1)
        for j = 1:size(tsigma, 2)
            @test isapprox(tsigma[i,j], tsigmatest[i,j])
        end
    end
end
@testset "rotstrain" begin
    # rotstrain(strain::Array{Float64,1},theta::Float64)
    strain = rand(3)
    tepstest = [0.5 0.5 0.5;0.5 0.5 -0.5;-1.0 1.0 0.0]*strain
    teps = Composites.rotstrain(strain, 45.0)
    for i = 1:size(teps, 1)
        for j = 1:size(teps, 2)
            @test isapprox(teps[i,j], tepstest[i,j])
        end
    end
end

force = [25000.0, 0.0, 0.0]
moment = [1000.0, 0.0, 0.0]
resultantstress = vcat(force, moment)
S = vcat(hcat(A, B), hcat(B', D))
resultantstrain = S\resultantstress

lowerplystrain,upperplystrain = Composites.getplystrain(lam, resultantstrain::Array{Float64,1})
lowerplystress = Composites.getplystress(lowerplystrain, q, matid)
upperplystress = Composites.getplystress(upperplystrain, q, matid)
@testset "getplystrain/getplystress" begin
    # getplystrain(lam::Laminate,resultantstrain::Array{Float64,1})
    # getplystrain(nply::Array{Int64,1},tply::Array{Float64,1},theta::Array{Float64,1},
    #   resultantstrain::Array{Float64,1})
    @test (lowerplystrain, upperplystrain) == Composites.getplystrain(nply, tply,
        theta, resultantstrain)
    # getplystress(plystrain::Array{Float64,2},q::Array{Float64,3},matid::Array{Int64,1})
    # getplystress(plystrain::Array{Float64,2},q::Array{Float64,3},lam::Laminate) = getplystress(plystrain,q,matid)
    @test Composites.getplystress(lowerplystrain, q, lam) == lowerplystress
    lowerplystresstest = [-570.672 -54.747 -76.129 80.735 46.629 46.014 456.66;
                          11.223  -9.655  -4.223 -7.519  2.922  7.337 -8.386;
                          2.556   20.799 -13.647 -0.639 7.808 -14.96  -1.917] #MPa
    for i = 1:3
        for j = 1:length(lowerplystress)
            @test round.(lowerplystress[j][i]/1e6, digits = 3) == lowerplystresstest[i,j]
        end
    end
    # getplystress(plystrain::Array{Float64,2},q::Array{Float64,3},lam::Laminate) = getplystress(plystrain,q,matid)
    @test Composites.getplystress(upperplystrain, q, lam) == upperplystress
    upperplystresstest = [-423.91 -34.595 -35.21  -90.646 87.548  66.167 603.422;
                           8.422  -6.257  -1.841  9.644  5.304  10.736 -11.187;
                           1.917  13.647  -6.496  0.639  14.96 -22.112 -2.556] #MPa
    for i = 1:3
        for j = 1:length(upperplystress)
            @test round(upperplystress[j][i]/1e6, digits=3) == upperplystresstest[i,j]
        end
    end
end

@testset "maxstress" begin
    failmaxstresstest = [
        -0.302793   -0.0247108  -0.0251497  -0.0647469   0.0625344   0.0472619   0.431016;
        1.80387     0.147213    0.149828    0.385726   -0.372545   -0.28156    -2.56775  ;
        0.701835   -0.521387   -0.153435    0.803671    0.441966    0.894661   -0.932254 ;
        -0.158906    0.11805     0.0347399  -0.181963   -0.100068   -0.202565    0.211076 ;
        0.0563767   0.401395   -0.191046    0.0187922   0.439999   -0.650348   -0.0751689;
        -0.0563767  -0.401395    0.191046   -0.0187922  -0.439999    0.650348    0.0751689]
    sfmaxstresstest = [
        -3.30258   -40.4682   -39.7618   -15.4447   15.9912   21.1587     2.3201;
        0.554362    6.79287    6.67431    2.59251  -2.68424  -3.55164   -0.389446;
        1.42484    -1.91796   -6.51743    1.24429   2.26262   1.11774   -1.07267;
        -6.29303     8.47099   28.7853    -5.49561  -9.99323  -4.93669    4.73762;
        17.7378      2.49131   -5.23433   53.2135    2.27273  -1.53764  -13.3034;
        -17.7378     -2.49131    5.23433  -53.2135   -2.27273   1.53764   13.3034]
    failmaxstress,sfmaxstress = Composites.getmatfail(upperplystress, mat,
        lam, "maxstress")
    for i = 1:6
        for j = 1:7
            @test round(failmaxstress[j][i], sigdigits=6) == failmaxstresstest[i,j]
            @test round(sfmaxstress[j][i], sigdigits=6) == sfmaxstresstest[i,j]
        end
    end
end

@testset "tsaiwu" begin
    failtsaiwutest = [2.95173,-0.069493,0.0470993,1.17469,0.260911,1.02618,-1.08208]
    sftsaiwutest = [0.413356,2.93886, 4.85947, 0.873163, 2.02113, 0.983487, 1.90494]
    failtsaiwu, sftsaiwu = Composites.getmatfail(upperplystress, mat, lam, "tsaiwu")
    @test isapprox(failtsaiwu, failtsaiwutest, atol=1e-5)
    @test isapprox(sftsaiwu, sftsaiwutest, atol=1e-5)
end

@testset "hashinrotem" begin
    failhashinrotemtest = [
    -0.407623  -0.0391053  -0.0543778   0.0576675   0.0333063   0.0328674   0.326186;
    2.42839    0.232968    0.323953   -0.343551   -0.198421   -0.195806   -1.94323;
    0.880392  -1.02161    -0.284952   -0.392936    0.112031    0.567472   -0.491518;
    -0.050493   0.407417    0.167466    0.0204784  -0.0557789  -0.212765    0.0282125]
    sfhashinrotemtest = [
    -2.45325   -25.572     -18.3899   17.3408   30.0243   30.4253    3.06574 ;
    0.411795    4.29244     3.08687  -2.91077  -5.0398   -5.10711  -0.514606;
    1.06577     0.989369    1.87333   1.59529   2.98766   1.32748   1.42636 ;
    -4.45025     1.56668     2.44364   6.98798  -4.23414  -2.16795   5.9536  ]
    failhashinrotem, sfhashinrotem = Composites.getmatfail(lowerplystress, mat,
        lam, "hashinrotem")
    for i = 1:4
        for j = 1:7
            @test round(failhashinrotem[j][i], sigdigits=6) == failhashinrotemtest[i,j]
            @test round(sfhashinrotem[j][i], sigdigits=6) == sfhashinrotemtest[i,j]
        end
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2018/21")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# Algorithm 6 from [Joldes, Muller, Popescu 2017]
# relative error < 3u², 20 FP Ops
# rel accuracy [102.4, 104.4] bits
# "The largest relative error we have obtain through many tests is around 2.25 × 2^−2p = 2.25 u^2.
#  For Float64, this is 2.25*2^(-106) [rel accuracy 104.8 bits]." pg 11
function Base.:(+)(x::Double64, y::Double64)
    hi, lo   = two_sum(Hi(x), Hi(y))
    thi, tlo = two_sum(Lo(x), Lo(y))
    c = lo + thi
    hi, lo = two_hilo_sum(hi, c)
    c = tlo + lo
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

biterror(::typeof(+)) = 2.25

function Base.:(+)(x::Double64, y::Float64)
    hi, lo = two_sum(Hi(x), y)
    c = lo + Lo(x)
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

function Base.:(+)(x::Float64, y::Double64)
    hi, lo = two_sum(x, Hi(y))
    c = lo + Lo(y)
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

# relative error < 3u², 20 FP Ops, 102.4 bits (relative)
# Algorithm 6 from [Joldes, Muller, Popescu 2017]
# reworked for subtraction
function Base.:(-)(x::Double64, y::Double64)
    hi, lo   = two_diff(Hi(x), Hi(y))
    thi, tlo = two_diff(Lo(x), Lo(y))
    c = lo + thi
    hi, lo = two_hilo_sum(hi, c)
    c = tlo + lo
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

biterror(::typeof(-)) = 2.25

function Base.:(-)(x::Double64, y::Float64)
    hi, lo = two_diff(Hi(x), y)
    c = lo + Lo(x)
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

function Base.:(-)(x::Float64, y::Double64)
    hi, lo = two_diff(x, Hi(y))
    c = lo + Lo(y)
    hi, lo = two_hilo_sum(hi, c)
    return Double64((hi, lo))
end

# relative error < 5u², 9 FP Ops, 101.6 bits (relative)
# "The largest relative error we have obtain through many tests is 3.936 × 2^−2p = 3.936 * 2^(-106).
# Algorithm 12 from [Joldes, Muller, Popescu 2017]
function Base.:(*)(x::Double64, y::Double64)
    hi, lo = two_prod(Hi(x), Hi(y))
    t = Lo(x) * Lo(y)
    t = fma(Hi(x), Lo(y), t)
    t = fma(Lo(x), Hi(y), t)
    t = lo + t
    hi, lo = two_hilo_sum(hi, t)
    return Double64((hi, lo))
end

biterror(::typeof(*)) = 3.936

function Base.:(*)(x::Double64, y::Float64)
    hi, lo = two_prod(Hi(x), y)
    t = Lo(x) * y
    t = t + lo
    hi, lo = two_hilo_sum(hi, t)
    return Double64((hi, lo))
end

function Base.:(*)(x::Float64, y::Double64)
    hi, lo = two_prod(x, Hi(y))
    t = x * Lo(y)
    t = lo + t
    hi, lo = two_hilo_sum(hi, t)
    return Double64((hi, lo))
end

# from DoubleFloats
# Base.:(/)(x::Double64, y::Double64) = Double64(divide_dddd_dd(Hi(x), Lo(x), Hi(y), Lo(y)))
(div0)(x::Double64, y::Double64) = Double64(divide_dddd_dd(Hi(x), Lo(x), Hi(y), Lo(y)))

@inline function divide_dddd_dd(xhi::T, xlo::T, yhi::T, ylo::T) where {T<:Float64}
    hi = xhi / yhi
    uh, ul = two_prod(hi, yhi)
    lo = ((((xhi - uh) - ul) + xlo) - hi*ylo)/yhi
    hi,lo = two_hilo_sum(hi, lo)
    return hi, lo
end

# Algorithm 17 from [Joldes, Muller, Popescu 2017]
# relative error <= 15u² + 56u³
function div17(x::Double64, y::Double64)
    xhi, xlo = HiLo(x)
    yhi, ylo = HiLo(y)
    thi = xhi / yhi
    rhi, rlo = DWTimesFP3(yhi, ylo, thi)
    dhi = xhi - rhi
    dlo = xlo - rlo
    d = dhi + dlo
    tlo = d / yhi
    return Double64(two_hilo_sum(thi, tlo))
end

function (div18)(x::Double64, y::Double64)
    yhi, ylo = HiLo(y)
    th = inv(yhi)
    rh = fma(-yhi, th, 1.0)
    rl = -ylo * th
    eh, el = two_hilo_sum(rh, rl)
    dh, dl = DWTimesFP3(eh, el, th) # (eh, el) * rh
    mh, ml = DWPlusFP(dh, dl, th) # (dh, dl) + th
    zh, zl = DWTimesDW2(Hi(x), Lo(x), mh, ml) # (xh, xl) * (mh, ml)
    return Double64((zh, zl))
end

# relative error < 9.8u², 31 FP Ops, 100.7 bits (relative)
# " 5.922 x 2^(-106) "
# Algorithm 18 from [Joldes, Muller, Popescu 2017]
function Base.:(/)(x::Double64, y::Double64)
    yhi, ylo = HiLo(y)
    th = inv(yhi)
    rh = fma(-yhi, th, 1.0)
    rl = -ylo * th
    eh, el = two_hilo_sum(rh, rl)
    dh, dl = DWTimesFP3(eh, el, th) # (eh, el) * rh
    mh, ml = DWPlusFP(dh, dl, th) # (dh, dl) + th
    zh, zl = DWTimesDW2(Hi(x), Lo(x), mh, ml) # (xh, xl) * (mh, ml)
    return Double64((zh, zl))
end

    biterror(::typeof(/)) = 5.922
biterror(::typeof(inv)) = 5.922

function Base.:(/)(x::Double64, y::Float64)
    th = inv(y)
    rh = fma(-y, th, 1.0)
    dh, dl = DWTimesFP3(rh, 0.0, th) # (eh, el) * rh
    mh, ml = DWPlusFP(dh, dl, th) # (dh, dl) + th
    zh, zl = DWTimesDW2(Hi(x), Lo(x), mh, ml) # (xh, xl) * (mh, ml)
    return Double64((zh, zl))
end

function Base.:(/)(x::Float64, y::Double64)
    yhi, ylo = HiLo(y)
    th = inv(yhi)
    rh = fma(-yhi, th, 1.0)
    rl = -ylo * th
    eh, el = two_hilo_sum(rh, rl)
    dh, dl = DWTimesFP3(eh, el, th) # (eh, el) * rh
    mh, ml = DWPlusFP(dh, dl, th) # (dh, dl) + th
    zh, zl = DWTimesDW2(x, 0.0, mh, ml) # (xh, xl) * (mh, ml)
    return Double64((zh, zl))
end

# algorithm 4 from [Joldes, Muller, Popescu 2017]
@inline function DWPlusFP(xh::Float64, xl::Float64, y::Float64)
    sh, sl = two_sum(xh, y)
    v = xl + sl
    zh, zl = two_hilo_sum(sh, v)
    return zh, zl
end
# algorithm 7 from [Joldes, Muller, Popescu 2017]
# calculation of (xh, xl) * y
@inline function DWTimesFP1(xh, xl, y)
    ch, cl1 = two_prod(xh, y)
    cl2 = xl*y
    th, tl1 = two_hilo_sum(ch, cl2)
    tl2 = tl1 + cl1
    zh, zl = two_hilo_sum(th, tl2)
    return zh, zl
end
# algorithm 9 from [Joldes, Muller, Popescu 2017]
# calculation of (xh, xl) * y
@inline function DWTimesFP3(xh, xl, y)
    ch, cl1 = two_prod(xh, y)
    cl3 = fma(xl, y, cl1)
    zh, zl = two_hilo_sum(ch, cl3)
    return zh, zl
end

# algorithm 10 from [Joldes, Muller, Popescu 2017]
# calcualtes (xh, xl) * (yh, yl)
# the largest relative error found is  4.9916 × 2^(-106)
@inline function DWTimesDW1(xh, xl, yh, yl)
    ch, cl1 = two_prod(xh, yh)
    tl1 = xh * yl
    tl2 = xl * yh
    cl2 = tl1 + tl2
    cl3 = cl1 + cl2
    zh, zl = two_hilo_sum(ch, cl3)
    return zh, zl
end

# algorithm 12 from [Joldes, Muller, Popescu 2017]
# calcualtes (xh, xl) * (yh, yl)
# the largest relative error found is  3.936 × 2^(-106)
@inline function DWTimesDW2(xh, xl, yh, yl)
    ch, cl1 = two_prod(xh, yh)
    tl0 = xl * yl
    tl1 = fma(xh, yl, tl0)
    cl2 = fma(xl, yh, tl1)
    cl3 = cl1 + cl2
    zh, zl = two_hilo_sum(ch, cl3)
    return zh, zl
end

# relative error < 9.8u², 31 FP Ops, 100.7 bits (relative)
# Algorithm 18 from [Joldes, Muller, Popescu 2017]
function Base.:(inv)(y::Double64)
    yhi, ylo = HiLo(y)
    th = inv(yhi)
    rh = fma(-yhi, th, 1.0)
    rl = -ylo * th
    eh, el = two_hilo_sum(rh, rl)
    dh, dl = DWTimesFP3(eh, el, th) # (eh, el) * rh
    mh, ml = DWPlusFP(dh, dl, th) # (dh, dl) + th
    # zh, zl = DWTimesDW2(1.0, 0.0, mh, ml) # (1, 0) * (mh, ml)
    return Double64((mh, ml))
end

Base.:(\)(x::Double64, y::Double64) = (/)(y, x)
Base.:(\)(x::Double64, y::Float64) = (/)(y, x)
Base.:(\)(x::Float64, y::Double64) = (/)(y, x)

for (T1, T2) in ((:ComplexD64, :ComplexD64), (:ComplexD64, :ComplexF64), (:ComplexF64, :ComplexD64))
  @eval begin

    function Base.:(+)(x::$T1, y::$T2)
        re = real(x) + real(y)
        im = imag(x) + imag(y)
        hi = ComplexF64(Hi(re), Hi(im))
        lo = ComplexF64(Lo(re), Lo(im))
        return ComplexD64((hi, lo))
    end

    function Base.:(-)(x::$T1, y::$T2)
        re = real(x) - real(y)
        im = imag(x) - imag(y)
        hi = ComplexF64(Hi(re), Hi(im))
        lo = ComplexF64(Lo(re), Lo(im))
        return ComplexD64((hi, lo))
    end

    function Base.:(*)(x::$T1, y::$T2)
        xr = real(x); xi = imag(x); yr = real(y); yi = imag(y)
        re = xr * yr - xi * yi
        im = xr * yi + xi * yr
        hi = ComplexF64(Hi(re), Hi(im))
        lo = ComplexF64(Lo(re), Lo(im))
        return ComplexD64((hi, lo))
    end
  end
end

#=
   This is adapted from
   "Improved Complex Division"
    by Michael Baudin (DIGITEO) and Robert Smith (Stanford University)
    Version 0.1. February 2011

    "A Robust Complex Division in Scilab
    by Michael Baudin, Robert L. Smith
    2012
=#

function Base.:(/)(x::ComplexD64, y::ComplexD64)
    a = real(x); b = imag(x); c = real(y); d = imag(y)
    if abs(d) <= abs(c)
        r = d/c
        t = inv(fma(d, r, c))
        if r == 0
            d = d * t
            e = fma( d, b/c, a) # (a + d * (b/c)) * t
            f = fma(-d, a/c, b) # (b - d * (a/c)) * t
        else
            b = b * t
            a = a * t
            e = fma( b, r, a) # (a + b * r) * t
            f = fma(-a, r, b) # (b - a * r) * t
        end
    else
        r = c/d
        t = inv(fma(c, r, d)) # 1/(c * r + d )
        if r == 0
            e = fma(c, a/d,  b) * t # (c * (a/d) + b) * t
            f = fma(c, b/d, -a) * t # (c * (b/d) - a) * t
        else
            e = fma(a, r,  b) * t # (a * r + b) * t
            f = fma(b, r, -a) * t #(b * r - a) * t
        end
    end
    hi = ComplexF64(Hi(e), Hi(f))
    lo = ComplexF64(Lo(e), Lo(f))
    return ComplexD64((hi, lo))
end

function divide(x::ComplexD64, y::ComplexD64)
    a = real(x); b = imag(x); c = real(y); d = imag(y)
    if abs(d) <= abs(c)
        r = d/c
        t = inv(muladd(d, r, c))
        if r == 0
            d = d * t
            e = muladd( d, b/c, a) # (a + d * (b/c)) * t
            f = muladd(-d, a/c, b) # (b - d * (a/c)) * t
        else
            b = b * t
            a = a * t
            e = muladd( b, r, a) # (a + b * r) * t
            f = muladd(-a, r, b) # (b - a * r) * t
        end
    else
        r = c/d
        t = inv(muladd(c, r, d)) # 1/(c * r + d )
        if r == 0
            e = muladd(c, a/d,  b) * t # (c * (a/d) + b) * t
            f = muladd(c, b/d, -a) * t # (c * (b/d) - a) * t
        else
            e = muladd(a, r,  b) * t # (a * r + b) * t
            f = muladd(b, r, -a) * t #(b * r - a) * t
        end
    end
    hi = ComplexF64(Hi(e), Hi(f))
    lo = ComplexF64(Lo(e), Lo(f))
    return ComplexD64((hi, lo))
end

#=
adapted from
On complex multiplication and division with an FMA
Jean-Michel Muller, C.-P. Jeannerod, P. Kornerup and N. Louvet
INVA 2014

function mul(ab::ComplexD64, cd::ComplexD64)
    re = Kahan(real(ab), real(cd), -imag(ab), imag(cd))
    im = Kahan(real(ab), imag(cd), imag(ab), real(cd))
    return ComplexD64((ComplexF64(Hi(re), Hi(im)), ComplexF64(Lo(re), Lo(im))))
end

function divide(ab::ComplexD64, cd::ComplexD64)
    a = real(ab); b = imag(ab); c = real(cd); d = imag(cd)
    delta = c*c + d*d
    gre = Kahan(a, b,c,d) # ac + bd
    gim = Kahan(b,-a,c,d) # bc - ad
    qre = gre / delta
    qim = gim / delta
    return qre, qim
end

function Kahan(a, b, c, d)
    w = c * d
    e = fma(c, d, -w)
    f = fma(a, b, w)
    r = f + e
    return r
end

Base.fma(x::Double64,y::Double64,z::Double64) = muladd(x,y,z)
=#
#=

function normalizedenom(n,d)
    fr,xp = frexp(d)
    n = ldexp(n, -xp)
    return n, fr
end

function normalized(n::Double64, d::Double64)
    scaleby = inv(Hi(d))
    return n*scaleby, d*scaleby
end

function Goldschmidt(n::T,d::T,k) where {T}
    n, d = normalized(n,d)
    e = 1 - d
    q = n;  qold = zero(T)
    while !iszero(k) && (qold != q)
        k -= 1
        qold = q
        q = q * (1+e)
        e = e * e
    end
    return q
end

function Goldschmidt(n::Double64,d::Double64,k=4)}
    n1, d1 = normalized(n,d)
    e = 1 - d1
    q = n1;  qold = zero(T)
    while (qold != q) && !iszero(k)
       k -= 1
       qold = q
       q = q * (1+e) # q = q + q*e
       e = e * e
    end
    delta = (n1 - q*d1)/2
    return q + delta
end

=#


text:
#' ---
#' title: Space Shuttle Reentry Trajectory
#' ---

#' **Originally Contributed by**: Henrique Ferrolho

#' This tutorial demonstrates how to compute a reentry trajectory for the
#' [Space Shuttle](https://en.wikipedia.org/wiki/Space_Shuttle),
#' by formulating and solving a nonlinear programming problem.
#' The problem was drawn from Chapter 6 of
#' ["Practical Methods for Optimal Control and Estimation Using Nonlinear Programming"](https://epubs.siam.org/doi/book/10.1137/1.9780898718577),
#' by John T. Betts.

#' The motion of the vehicle is defined by the following set of DAEs:
#'
#' $$
#' \begin{align}
#' \dot{h} & = v \sin \gamma , \\
#' \dot{\phi} & = \frac{v}{r} \cos \gamma \sin \psi / \cos \theta , \\
#' \dot{\theta} & = \frac{v}{r} \cos \gamma \cos \psi , \\
#' \dot{v} & = -\frac{D}{m} - g \sin \gamma , \\
#' \dot{\gamma} & = \frac{L}{m v} \cos(\beta) + \cos \gamma \left ( \frac{v}{r} - \frac{g}{v} \right ) , \\
#' \dot{\psi} & = \frac{1}{m v \cos \gamma} L \sin(\beta) + \frac{v}{r \cos \theta} \cos \gamma \sin \psi \sin \theta , \\
#' q & \le q_U , \\
#' \end{align}
#' $$

#' where the aerodynamic heating on the vehicle wing leading edge is $q = q_a q_r$
#' and the dynamic variables are
#'
#' $$
#' \begin{align}
#'      h & \quad \text{altitude (ft)},     \qquad & & \gamma \quad \text{flight path angle (rad)}, \\
#'   \phi & \quad \text{longitude (rad)},   \qquad & & \psi   \quad \text{azimuth (rad)},           \\
#' \theta & \quad \text{latitude (rad)},    \qquad & & \alpha \quad \text{angle of attack (rad)},   \\
#'      v & \quad \text{velocity (ft/sec)}, \qquad & & \beta  \quad \text{bank angle (rad)}.
#' \end{align}
#' $$

#' The aerodynamic and atmospheric forces on the vehicle are
#' specified by the following quantities (English units):
#'
#' $$
#' \begin{align}
#'            D & = \frac{1}{2} c_D S \rho v^2,                  \qquad & a_0 & = -0.20704, \\
#'            L & = \frac{1}{2} c_L S \rho v^2,                  \qquad & a_1 & =  0.029244, \\
#'            g & = \mu / r^2,                                   \qquad & \mu & =  0.14076539 \times 10^{17}, \\
#'            r & = R_e + h,                                     \qquad & b_0 & =  0.07854, \\
#'         \rho & = \rho_0 \exp[-h/h_r],                         \qquad & b_1 & = -0.61592  \times 10^{-2}, \\
#'       \rho_0 & = 0.002378,                                    \qquad & b_2 & =  0.621408 \times 10^{-3}, \\
#'          h_r & = 23800,                                       \qquad & q_r & =  17700 \sqrt{\rho} (0.0001 v)^{3.07}, \\
#'          c_L & = a_0 + a_1 \hat{\alpha},                      \qquad & q_a & =  c_0 + c_1 \hat{\alpha} + c_2 \hat{\alpha}^2 + c_3 \hat{\alpha}^3, \\
#'          c_D & = b_0 + b_1 \hat{\alpha} + b_2 \hat{\alpha}^2, \qquad & c_0 & =  1.0672181, \\
#' \hat{\alpha} & = 180 \alpha / \pi,                            \qquad & c_1 & = -0.19213774 \times 10^{-1}, \\
#'          R_e & = 20902900,                                    \qquad & c_2 & =  0.21286289 \times 10^{-3}, \\
#'            S & = 2690,                                        \qquad & c_3 & = -0.10117249 \times 10^{-5}.
#' \end{align}
#' $$

#' The reentry trajectory begins at an altitude where the aerodynamic forces are quite small
#' with a weight of $w = 203000$ (lb) and mass $m = w / g_0$ (slug), where $g_0 = 32.174$ (ft/sec$^2$).
#' The initial conditions are as follows:
#'
#' $$
#' \begin{align}
#'      h & = 260000 \text{ ft},  \qquad & v      & = 25600 \text{ ft/sec}, \\
#'   \phi & = 0      \text{ deg}, \qquad & \gamma & = -1    \text{ deg},    \\
#' \theta & = 0      \text{ deg}, \qquad & \psi   & = 90    \text{ deg}.
#' \end{align}
#' $$

#' The final point on the reentry trajectory occurs at the unknown (free) time $t_F$,
#' at the so-called terminal area energy management (TAEM) interface, which is defined by the conditions
#'
#' $$
#' h = 80000 \text{ ft}, \qquad v = 2500 \text{ ft/sec}, \qquad \gamma = -5 \text{ deg}.
#' $$

#' As explained in the book, our goal is to maximize the final crossrange,
#' which is equivalent to maximizing the final latitude of the vehicle,
#' i.e., $J = \theta(t_F)$.

#' ![Max crossrange shuttle reentry](https://i.imgur.com/1tUN3wM.png)

#' ### Approach
#' We will use a discretized model of time, with a fixed number of discretized points, $n$.
#' The decision variables at each point are going to be the state of the vehicle and the controls commanded to it.
#' In addition, we will also make each time step size $\Delta t$ a decision variable;
#' that way, we can either fix the time step size easily,
#' or allow the solver to fine-tune the duration between each adjacent pair of points.
#' Finally, in order to approximate the derivatives of the problem dynamics, we will use either rectangular
#' or [trapezoidal](http://en.wikipedia.org/wiki/Trapezoidal_rule) integration.

#' ### Disclaimer
#' Do not try to actually land a Space Shuttle using this notebook! 😛
#' There's no mesh refinement going on, which can lead to unrealistic trajectories having position
#' and velocity errors with orders of magnitude $10^4$ ft and $10^2$ ft/sec, respectively.

using Interpolations, Plots
using JuMP, Ipopt  # , KNITRO

#+

# Global variables
const w  = 203000.0  # weight (lb)
const g₀ = 32.174    # acceleration (ft/sec^2)
const m  = w / g₀    # mass (slug)

# Aerodynamic and atmospheric forces on the vehicle
const ρ₀ =  0.002378
const hᵣ =  23800.0
const Rₑ =  20902900.0
const μ  =  0.14076539e17
const S  =  2690.0
const a₀ = -0.20704
const a₁ =  0.029244
const b₀ =  0.07854
const b₁ = -0.61592e-2
const b₂ =  0.621408e-3
const c₀ =  1.0672181
const c₁ = -0.19213774e-1
const c₂ =  0.21286289e-3
const c₃ = -0.10117249e-5

# Initial conditions
const hₛ = 2.6          # altitude (ft) / 1e5
const ϕₛ = deg2rad(0)   # longitude (rad)
const θₛ = deg2rad(0)   # latitude (rad)
const vₛ = 2.56         # velocity (ft/sec) / 1e4
const γₛ = deg2rad(-1)  # flight path angle (rad)
const ψₛ = deg2rad(90)  # azimuth (rad)
const αₛ = deg2rad(0)   # angle of attack (rad)
const βₛ = deg2rad(0)   # banck angle (rad)
const tₛ = 1.00         # time step (sec)

# Final conditions, the so-called Terminal Area Energy Management (TAEM)
const hₜ = 0.8          # altitude (ft) / 1e5
const vₜ = 0.25         # velocity (ft/sec) / 1e4
const γₜ = deg2rad(-5)  # flight path angle (rad)

# Number of mesh points (knots) to be used
const n = 2009

# Integration scheme to be used for the dynamics
const integration_rule = "rectangular"

#+

user_options_ipopt = (
    "mu_strategy" => "monotone",
    "linear_solver" => "ma27",  # For the best results, it is advised to experiment different linear solvers.
                                # If Ipopt is not compiled with MA27/MA57, it may fallback to 'MUMPS'.
                                # In general, the linear solver MA27 is much faster than MUMPS.
)

# user_options_knitro = (
#     "algorithm" => KNITRO.KN_ALG_BAR_DIRECT,
#     "bar_murule" => KNITRO.KN_BAR_MURULE_QUALITY,
#     "linsolver" => KNITRO.KN_LINSOLVER_MA27,
# )

# Create JuMP model, using Ipopt as the solver
model = Model(optimizer_with_attributes(Ipopt.Optimizer, user_options_ipopt...))
# model = Model(optimizer_with_attributes(KNITRO.Optimizer, user_options_knitro...))

@variables(model, begin
               0 ≤ scaled_h[1:n]                # altitude (ft) / 1e5
                          ϕ[1:n]                # longitude (rad)
    deg2rad(-89) ≤        θ[1:n] ≤ deg2rad(89)  # latitude (rad)
            1e-4 ≤ scaled_v[1:n]                # velocity (ft/sec) / 1e4
    deg2rad(-89) ≤        γ[1:n] ≤ deg2rad(89)  # flight path angle (rad)
                          ψ[1:n]                # azimuth (rad)
    deg2rad(-90) ≤        α[1:n] ≤ deg2rad(90)  # angle of attack (rad)
    deg2rad(-89) ≤        β[1:n] ≤ deg2rad( 1)  # banck angle (rad)
    #        0.5 ≤       Δt[1:n] ≤ 1.5          # time step (sec)
                         Δt[1:n] == 1.0         # time step (sec)
end)

# Fix initial conditions
fix(scaled_h[1], hₛ; force=true)
fix(       ϕ[1], ϕₛ; force=true)
fix(       θ[1], θₛ; force=true)
fix(scaled_v[1], vₛ; force=true)
fix(       γ[1], γₛ; force=true)
fix(       ψ[1], ψₛ; force=true)

# Fix final conditions
fix(scaled_h[n], hₜ; force=true)
fix(scaled_v[n], vₜ; force=true)
fix(       γ[n], γₜ; force=true)

# Initial guess: linear interpolation between boundary conditions
xₛ = [hₛ, ϕₛ, θₛ, vₛ, γₛ, ψₛ, αₛ, βₛ, tₛ]
xₜ = [hₜ, ϕₛ, θₛ, vₜ, γₜ, ψₛ, αₛ, βₛ, tₛ]
interp_linear = LinearInterpolation([1, n], [xₛ, xₜ])
initial_guess = mapreduce(transpose, vcat, interp_linear.(1:n))
set_start_value.(all_variables(model), vec(initial_guess))

# Functions to restore `h` and `v` to their true scale
@NLexpression(model, h[j = 1:n], scaled_h[j] * 1e5)
@NLexpression(model, v[j = 1:n], scaled_v[j] * 1e4)

# Helper functions
# @NLexpression(model, c_L[j = 1:n], a₀ + a₁ * rad2deg(α[j]))
# @NLexpression(model, c_D[j = 1:n], b₀ + b₁ * rad2deg(α[j]) + b₂ * rad2deg(α[j])^2)
@NLexpression(model, c_L[j = 1:n], a₀ + a₁ * (180 * α[j] / π))
@NLexpression(model, c_D[j = 1:n], b₀ + b₁ * (180 * α[j] / π) + b₂ * (180 * α[j] / π)^2)
@NLexpression(model,   ρ[j = 1:n], ρ₀ * exp(-h[j] / hᵣ))
@NLexpression(model,   D[j = 1:n], 0.5 * c_D[j] * S * ρ[j] * v[j]^2)
@NLexpression(model,   L[j = 1:n], 0.5 * c_L[j] * S * ρ[j] * v[j]^2)
@NLexpression(model,   r[j = 1:n], Rₑ + h[j])
@NLexpression(model,   g[j = 1:n], μ / r[j]^2)

# Motion of the vehicle as a differential-algebraic system of equations (DAEs)
@NLexpression(model, δh[j = 1:n], v[j] * sin(γ[j]))
@NLexpression(model, δϕ[j = 1:n], (v[j] / r[j]) * cos(γ[j]) * sin(ψ[j]) / cos(θ[j]))
@NLexpression(model, δθ[j = 1:n], (v[j] / r[j]) * cos(γ[j]) * cos(ψ[j]))
@NLexpression(model, δv[j = 1:n], -(D[j] / m) - g[j] * sin(γ[j]))
@NLexpression(model, δγ[j = 1:n], (L[j] / (m * v[j])) * cos(β[j]) + cos(γ[j]) * ((v[j] / r[j]) - (g[j] / v[j])))
@NLexpression(model, δψ[j = 1:n], (1 / (m * v[j] * cos(γ[j]))) * L[j] * sin(β[j]) + (v[j] / (r[j] * cos(θ[j]))) * cos(γ[j]) * sin(ψ[j]) * sin(θ[j]))

# System dynamics
for j in 2:n
    i = j - 1  # index of previous knot

    if integration_rule == "rectangular"
        # Rectangular integration
        @NLconstraint(model, h[j] == h[i] + Δt[i] * δh[i])
        @NLconstraint(model, ϕ[j] == ϕ[i] + Δt[i] * δϕ[i])
        @NLconstraint(model, θ[j] == θ[i] + Δt[i] * δθ[i])
        @NLconstraint(model, v[j] == v[i] + Δt[i] * δv[i])
        @NLconstraint(model, γ[j] == γ[i] + Δt[i] * δγ[i])
        @NLconstraint(model, ψ[j] == ψ[i] + Δt[i] * δψ[i])
    elseif integration_rule == "trapezoidal"
        # Trapezoidal integration
        @NLconstraint(model, h[j] == h[i] + 0.5 * Δt[i] * (δh[j] + δh[i]))
        @NLconstraint(model, ϕ[j] == ϕ[i] + 0.5 * Δt[i] * (δϕ[j] + δϕ[i]))
        @NLconstraint(model, θ[j] == θ[i] + 0.5 * Δt[i] * (δθ[j] + δθ[i]))
        @NLconstraint(model, v[j] == v[i] + 0.5 * Δt[i] * (δv[j] + δv[i]))
        @NLconstraint(model, γ[j] == γ[i] + 0.5 * Δt[i] * (δγ[j] + δγ[i]))
        @NLconstraint(model, ψ[j] == ψ[i] + 0.5 * Δt[i] * (δψ[j] + δψ[i]))
    else
        @error "Unexpected integration rule '$(integration_rule)'"
    end
end

# Objective: Maximize crossrange
@objective(model, Max, θ[n])

# Solve for the control and state
status = optimize!(model)

# Show final crossrange of the solution
println("Final latitude θ = ", round(objective_value(model) |> rad2deg, digits = 2), "°")

#+

ts = cumsum([0; value.(Δt)])[1:end-1]

plt_altitude = plot(ts, value.(scaled_h), legend = nothing, title = "Altitude (100,000 ft)")
plt_longitude = plot(ts, rad2deg.(value.(ϕ)), legend = nothing, title = "Longitude (deg)")
plt_latitude = plot(ts, rad2deg.(value.(θ)), legend = nothing, title = "Latitude (deg)")
plt_velocity = plot(ts, value.(scaled_v), legend = nothing, title = "Velocity (1000 ft/sec)")
plt_flight_path = plot(ts, rad2deg.(value.(γ)), legend = nothing, title = "Flight Path (deg)")
plt_azimuth = plot(ts, rad2deg.(value.(ψ)), legend = nothing, title = "Azimuth (deg)")

plt = plot(plt_altitude,  plt_velocity, plt_longitude, plt_flight_path, plt_latitude,
           plt_azimuth, layout=grid(3, 2), linewidth=2, size=(700, 700))

#+

function q(h, v, a)
    ρ(h) = ρ₀ * exp(-h / hᵣ)
    qᵣ(h, v) = 17700 * √ρ(h) * (0.0001 * v)^3.07
    qₐ(a) = c₀ + c₁ * rad2deg(a) + c₂ * rad2deg(a)^2 + c₃ * rad2deg(a)^3
    # Aerodynamic heating on the vehicle wing leading edge
    qₐ(a) * qᵣ(h, v)
end

plt_attack_angle = plot(ts[1:end-1], rad2deg.(value.(α)[1:end-1]), legend=nothing, title="Angle of Attack (deg)")
plt_bank_angle = plot(ts[1:end-1], rad2deg.(value.(β)[1:end-1]), legend=nothing, title="Bank Angle (deg)")
plt_heating = plot(ts, q.(value.(scaled_h)*1e5, value.(scaled_v)*1e4, value.(α)), legend=nothing, title="Heating (BTU/ft/ft/sec)")

plt = plot(plt_attack_angle, plt_bank_angle, plt_heating, layout=grid(3, 1), linewidth=2, size=(700, 700))

#+

plt = plot(rad2deg.(value.(ϕ)), rad2deg.(value.(θ)), value.(scaled_h),
           linewidth=2, legend=nothing, title="Space Shuttle Reentry Trajectory",
           xlabel="Longitude (deg)", ylabel="Latitude (deg)", zlabel="Altitude (100,000 ft)")


text:
# This file was generated, do not modify it. # hide
one_hundred_models = @from_network OneHundredModels(atom=atom) <= yhat

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
###
# Archive entry

###
# Basic object manipulation

# This will be changed to `Cint` in libarchive 4.0
const _la_mode_t = Cushort
if Sys.isunix()
    const _Cdev_t = UInt64
else
    const _Cdev_t = Cuint
end

mutable struct Entry
    ptr::Ptr{Cvoid}
    function Entry(ptr::Ptr{Cvoid})
        obj = new(ptr)
        finalizer(obj, free)
        obj
    end
    function Entry(archive::Archive)
        ptr = ccall((:archive_entry_new2, libarchive), Ptr{Cvoid},
                    (Ptr{Cvoid},), archive)
        ptr == C_NULL && throw(OutOfMemoryError())
        Entry(ptr)
    end
    # Mostly for testing purpose
    function Entry()
        ptr = ccall((:archive_entry_new, libarchive), Ptr{Cvoid}, ())
        ptr == C_NULL && throw(OutOfMemoryError())
        Entry(ptr)
    end
end

function Base.deepcopy_internal(entry::Entry, stackdict::IdDict)
    ptr = ccall((:archive_entry_clone, libarchive),
                Ptr{Cvoid}, (Ptr{Cvoid},), entry)
    ptr == C_NULL && throw(OutOfMemoryError())
    new_entry = Entry(ptr)
    stackdict[entry] = new_entry
    new_entry
end

function clear(entry::Entry)
    ccall((:archive_entry_clear, libarchive), Ptr{Cvoid}, (Ptr{Cvoid},), entry)
    entry
end

function free(entry::Entry)
    ptr = entry.ptr
    ptr == C_NULL && return
    ccall((:archive_entry_free, libarchive), Cvoid, (Ptr{Cvoid},), ptr)
    entry.ptr = C_NULL
    nothing
end

function Base.cconvert(::Type{Ptr{Cvoid}}, entry::Entry)
    entry.ptr == C_NULL && error("entry already freed")
    entry
end
Base.unsafe_convert(::Type{Ptr{Cvoid}}, entry::Entry) = entry.ptr

# TODO:
# As of libarchive 3.4.3, the only way to make sure all of the readers of all the string
# properties are happy is to use either the `set_*` (local encoding) function
# or the `update_*_utf8` (utf-8 encoding) function.
# We assume the input is in UTF8 so I'm going for the `update_*_utf8` one.
# The function does some encoding conversions eagerly which is not ideal.
# After https://github.com/libarchive/libarchive/pull/1389 or an equivalent fix
# we should be able to use `set_*_utf8` instead.

# Retrieve fields from an archive_entry.
#
# There are a number of implicit conversions among these fields.  For
# example, if a regular string field is set and you read the _w wide
# character field, the entry will implicitly convert narrow-to-wide
# using the current locale.  Similarly, dev values are automatically
# updated when you write devmajor or devminor and vice versa.
#
# In addition, fields can be "set" or "unset."  Unset string fields
# return NULL, non-string fields have _is_set() functions to test
# whether they've been set.  You can "unset" a string field by
# assigning NULL; non-string fields have _unset() functions to
# unset them.
#
# Note: There is one ambiguity in the above; string fields will
# also return NULL when implicit character set conversions fail.
# This is usually what you want.

atime(entry::Entry) = ccall((:archive_entry_atime, libarchive),
                            Int, (Ptr{Cvoid},), entry)
atime_nsec(entry::Entry) = ccall((:archive_entry_atime_nsec, libarchive),
                                 Clong, (Ptr{Cvoid},), entry)
atime_is_set(entry::Entry) = ccall((:archive_entry_atime_is_set, libarchive),
                                   Cint, (Ptr{Cvoid},), entry) != 0
birthtime(entry::Entry) = ccall((:archive_entry_birthtime, libarchive),
                                Int, (Ptr{Cvoid},), entry)
birthtime_nsec(entry::Entry) =
    ccall((:archive_entry_birthtime_nsec, libarchive),
          Clong, (Ptr{Cvoid},), entry)
birthtime_is_set(entry::Entry) =
    ccall((:archive_entry_birthtime_is_set, libarchive),
          Cint, (Ptr{Cvoid},), entry) != 0
ctime(entry::Entry) = ccall((:archive_entry_ctime, libarchive),
                            Int, (Ptr{Cvoid},), entry)
ctime_nsec(entry::Entry) = ccall((:archive_entry_ctime_nsec, libarchive),
                                 Clong, (Ptr{Cvoid},), entry)
ctime_is_set(entry::Entry) = ccall((:archive_entry_ctime_is_set, libarchive),
                                   Cint, (Ptr{Cvoid},), entry) != 0
dev(entry::Entry) = UInt64(ccall((:archive_entry_dev, libarchive),
                                 _Cdev_t, (Ptr{Cvoid},), entry))
dev_is_set(entry::Entry) = ccall((:archive_entry_dev_is_set, libarchive),
                                 Cint, (Ptr{Cvoid},), entry) != 0
devmajor(entry::Entry) = UInt64(ccall((:archive_entry_devmajor, libarchive),
                                      _Cdev_t, (Ptr{Cvoid},), entry))
devminor(entry::Entry) = UInt64(ccall((:archive_entry_devminor, libarchive),
                                      _Cdev_t, (Ptr{Cvoid},), entry))
filetype(entry::Entry) = Cint(ccall((:archive_entry_filetype, libarchive),
                                    _la_mode_t, (Ptr{Cvoid},), entry))
function fflags(entry::Entry)
    set = Ref{Culong}(0)
    clear = Ref{Culong}(0)
    ccall((:archive_entry_fflags, libarchive),
          Cvoid, (Ptr{Cvoid}, Ptr{Culong}, Ptr{Culong}), entry, set, clear)
    set[], clear[]
end
fflags_text(entry::Entry) =
    unsafe_string(ccall((:archive_entry_fflags_text, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
gid(entry::Entry) =
    ccall((:archive_entry_gid, libarchive), Int64, (Ptr{Cvoid},), entry)
gname(entry::Entry) =
    unsafe_string(ccall((:archive_entry_gname_utf8, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))

hardlink(entry::Entry) =
    unsafe_string(ccall((:archive_entry_hardlink_utf8, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
ino(entry::Entry) =
    ccall((:archive_entry_ino, libarchive), Int64, (Ptr{Cvoid},), entry)
ino_is_set(entry::Entry) =
    ccall((:archive_entry_ino_is_set, libarchive),
          Cint, (Ptr{Cvoid},), entry) != 0

mode(entry::Entry) =
    Cint(ccall((:archive_entry_mode, libarchive),
               _la_mode_t, (Ptr{Cvoid},), entry))

mtime(entry::Entry) = ccall((:archive_entry_mtime, libarchive),
                            Int, (Ptr{Cvoid},), entry)
mtime_nsec(entry::Entry) = ccall((:archive_entry_mtime_nsec, libarchive),
                                 Clong, (Ptr{Cvoid},), entry)
mtime_is_set(entry::Entry) = ccall((:archive_entry_mtime_is_set, libarchive),
                                   Cint, (Ptr{Cvoid},), entry) != 0

nlink(entry::Entry) =
    ccall((:archive_entry_nlink, libarchive), Cuint, (Ptr{Cvoid},), entry)
pathname(entry::Entry) =
    unsafe_string(ccall((:archive_entry_pathname_utf8, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))

perm(entry::Entry) =
    Cint(ccall((:archive_entry_perm, libarchive),
               _la_mode_t, (Ptr{Cvoid},), entry))
rdev(entry::Entry) =
    UInt64(ccall((:archive_entry_rdev, libarchive),
                 _Cdev_t, (Ptr{Cvoid},), entry))
rdevmajor(entry::Entry) =
    UInt64(ccall((:archive_entry_rdevmajor, libarchive),
                 _Cdev_t, (Ptr{Cvoid},), entry))
rdevminor(entry::Entry) =
    UInt64(ccall((:archive_entry_rdevminor, libarchive),
                 _Cdev_t, (Ptr{Cvoid},), entry))
sourcepath(entry::Entry) =
    unsafe_string(ccall((:archive_entry_sourcepath, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
Base.size(entry::Entry) =
    ccall((:archive_entry_size, libarchive), Int64, (Ptr{Cvoid},), entry)
size_is_set(entry::Entry) =
    ccall((:archive_entry_size_is_set, libarchive),
          Cint, (Ptr{Cvoid},), entry) != 0
strmode(entry::Entry) =
    unsafe_string(ccall((:archive_entry_strmode, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
symlink(entry::Entry) =
    unsafe_string(ccall((:archive_entry_symlink_utf8, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
symlink_type(entry::Entry) =
    ccall((:archive_entry_symlink_type, libarchive), Cint, (Ptr{Cvoid},), entry)
uid(entry::Entry) =
    ccall((:archive_entry_uid, libarchive), Int64, (Ptr{Cvoid},), entry)
uname(entry::Entry) =
    unsafe_string(ccall((:archive_entry_uname_utf8, libarchive),
                        Ptr{UInt8}, (Ptr{Cvoid},), entry))
is_data_encrypted(entry::Entry) =
    ccall((:archive_entry_is_data_encrypted, libarchive), Cint, (Ptr{Cvoid},), entry) != 0
is_metadata_encrypted(entry::Entry) =
    ccall((:archive_entry_is_metadata_encrypted, libarchive), Cint, (Ptr{Cvoid},), entry) != 0
is_encrypted(entry::Entry) =
    ccall((:archive_entry_is_encrypted, libarchive), Cint, (Ptr{Cvoid},), entry) != 0

set_atime(entry::Entry, t, ns) =
    ccall((:archive_entry_set_atime, libarchive),
          Cvoid, (Ptr{Cvoid}, Int, Clong), entry, t, ns)
unset_atime(entry::Entry) =
    ccall((:archive_entry_unset_atime, libarchive), Cvoid, (Ptr{Cvoid},), entry)
set_birthtime(entry::Entry, t, ns) =
    ccall((:archive_entry_set_birthtime, libarchive),
          Cvoid, (Ptr{Cvoid}, Int, Clong), entry, t, ns)
unset_birthtime(entry::Entry) =
    ccall((:archive_entry_unset_birthtime, libarchive),
          Cvoid, (Ptr{Cvoid},), entry)
set_ctime(entry::Entry, t, ns) =
    ccall((:archive_entry_set_ctime, libarchive),
          Cvoid, (Ptr{Cvoid}, Int, Clong), entry, t, ns)
unset_ctime(entry::Entry) =
    ccall((:archive_entry_unset_ctime, libarchive), Cvoid, (Ptr{Cvoid},), entry)
set_dev(entry::Entry, dev) =
    ccall((:archive_entry_set_dev, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, dev)
set_devmajor(entry::Entry, dev) =
    ccall((:archive_entry_set_devmajor, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, dev)
set_devminor(entry::Entry, dev) =
    ccall((:archive_entry_set_devminor, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, dev)
set_filetype(entry::Entry, ftype) =
    ccall((:archive_entry_set_filetype, libarchive),
          Cvoid, (Ptr{Cvoid}, Cuint), entry, ftype)
set_fflags(entry::Entry, set, clear) =
    ccall((:archive_entry_set_fflags, libarchive),
          Cvoid, (Ptr{Cvoid}, Culong, Culong), entry, set, clear)
set_fflags(entry::Entry, fflags::AbstractString) =
    (ccall((:archive_entry_copy_fflags_text, libarchive),
           Ptr{Cvoid}, (Ptr{Cvoid}, Cstring), entry, fflags); nothing)
set_gid(entry::Entry, gid) =
    ccall((:archive_entry_set_gid, libarchive),
          Cvoid, (Ptr{Cvoid}, Int64), entry, gid)
set_gname(entry::Entry, gname::AbstractString) =
    ccall((:archive_entry_update_gname_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, gname)
set_hardlink(entry::Entry, hl::AbstractString) =
    ccall((:archive_entry_update_hardlink_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, hl)
set_ino(entry::Entry, ino) =
    ccall((:archive_entry_set_ino, libarchive),
          Cvoid, (Ptr{Cvoid}, Int64), entry, ino)
set_link(entry::Entry, link::AbstractString) =
    ccall((:archive_entry_update_link_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, link)
set_mode(entry::Entry, mode) =
    ccall((:archive_entry_set_mode, libarchive),
          Cvoid, (Ptr{Cvoid}, _la_mode_t), entry, mode)
set_mtime(entry::Entry, t, ns) =
    ccall((:archive_entry_set_mtime, libarchive),
          Cvoid, (Ptr{Cvoid}, Int, Clong), entry, t, ns)
unset_mtime(entry::Entry) =
    ccall((:archive_entry_unset_mtime, libarchive), Cvoid, (Ptr{Cvoid},), entry)
set_nlink(entry::Entry, nlink) =
    ccall((:archive_entry_set_nlink, libarchive),
          Cvoid, (Ptr{Cvoid}, Cuint), entry, nlink)
set_pathname(entry::Entry, path::AbstractString) =
    ccall((:archive_entry_update_pathname_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, path)
set_perm(entry::Entry, perm) =
    ccall((:archive_entry_set_perm, libarchive),
          Cvoid, (Ptr{Cvoid}, _la_mode_t), entry, perm)
set_rdev(entry::Entry, rdev) =
    ccall((:archive_entry_set_rdev, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, rdev)
set_rdevmajor(entry::Entry, rdev) =
    ccall((:archive_entry_set_rdevmajor, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, rdev)
set_rdevminor(entry::Entry, rdev) =
    ccall((:archive_entry_set_rdevminor, libarchive),
          Cvoid, (Ptr{Cvoid}, _Cdev_t), entry, rdev)
set_size(entry::Entry, size) =
    ccall((:archive_entry_set_size, libarchive),
          Cvoid, (Ptr{Cvoid}, Int64), entry, size)
unset_size(entry::Entry) =
    ccall((:archive_entry_unset_size, libarchive), Cvoid, (Ptr{Cvoid},), entry)
set_sourcepath(entry::Entry, path::AbstractString) =
    ccall((:archive_entry_copy_sourcepath, libarchive),
          Cvoid, (Ptr{Cvoid}, Cstring), entry, path)
set_symlink(entry::Entry, sym::AbstractString) =
    ccall((:archive_entry_update_symlink_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, sym)
set_symlink_type(entry::Entry, typ) =
    ccall((:archive_entry_set_symlink_type, libarchive),
          Cvoid, (Ptr{Cvoid}, Cint), entry, typ)
set_uid(entry::Entry, uid) =
    ccall((:archive_entry_set_uid, libarchive),
          Cvoid, (Ptr{Cvoid}, Int64), entry, uid)
set_uname(entry::Entry, uname::AbstractString) =
    ccall((:archive_entry_update_uname_utf8, libarchive),
          Cint, (Ptr{Cvoid}, Cstring), entry, uname)
set_is_data_encrypted(entry::Entry, encrypted::Bool) =
    ccall((:archive_entry_set_is_data_encrypted, libarchive),
          Cvoid, (Ptr{Cvoid}, Cint), entry, encrypted)
set_is_metadata_encrypted(entry::Entry, encrypted::Bool) =
    ccall((:archive_entry_set_is_metadata_encrypted, libarchive),
          Cvoid, (Ptr{Cvoid}, Cint), entry, encrypted)

# Storage for Mac OS-specific AppleDouble metadata information.
# Apple-format tar files store a separate binary blob containing
# encoded metadata with ACL, extended attributes, etc.
# This provides a place to store that blob.
function mac_metadata(entry::Entry)
    _sz = Ref{Csize_t}()
    ptr = ccall((:archive_entry_mac_metadata, libarchive),
                Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Csize_t}), entry, _sz)
    sz = _sz[]
    data = Vector{UInt8}(undef, sz)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), data, ptr, sz)
    data
end
set_mac_metadata(entry::Entry, data::Vector{UInt8}) =
    ccall((:archive_entry_copy_mac_metadata, libarchive),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), entry, data, sizeof(data))

# ACL routines.  This used to simply store and return text-format ACL
# strings, but that proved insufficient for a number of reasons:
#   = clients need control over uname/uid and gname/gid mappings
#   = there are many different ACL text formats
#   = would like to be able to read/convert archives containing ACLs
#     on platforms that lack ACL libraries
#
# This last point, in particular, forces me to implement a reasonably
# complete set of ACL support routines.

# Set the ACL by clearing it and adding entries one at a time.
# Unlike the POSIX.1e ACL routines, you must specify the type
# (access/default) for each entry.  Internally, the ACL data is just
# a soup of entries.  API calls here allow you to retrieve just the
# entries of interest.  This design (which goes against the spirit of
# POSIX.1e) is useful for handling archive formats that combine
# default and access information in a single ACL list.
acl_clear(entry::Entry) =
    ccall((:archive_entry_acl_clear, libarchive), Cvoid, (Ptr{Cvoid},), entry)
acl_add_entry(entry::Entry, typ, perm, tag, qual, name::AbstractString) =
    @_la_call(archive_entry_acl_add_entry,
              (Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cstring),
              entry, typ, perm, tag, qual, name)

# To retrieve the ACL, first "reset", then repeatedly ask for the
# "next" entry.  The want_type parameter allows you to request only
# certain types of entries.
acl_reset(entry::Entry, want) =
    ccall((:archive_entry_acl_reset, libarchive), Cint,
          (Ptr{Cvoid}, Cint), entry, want)
function acl_next(entry::Entry, want)
    typ = Ref{Cint}()
    perm = Ref{Cint}()
    tag = Ref{Cint}()
    qual = Ref{Cint}()
    name = Ref{Ptr{UInt8}}()
    @_la_call(archive_entry_acl_next,
              (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
               Ptr{Ptr{UInt8}}), entry, want, typ, perm, tag, qual, name)
    typ[], perm[], tag[], qual[], unsafe_string(name[])
end

"""
Construct a text-format ACL.  The flags argument is a bitmask that
can include any of the following:

* `ACL.Type.ACCESS` - Include POSIX.1e "access" entries.
* `ACL.Type.DEFAULT` - Include POSIX.1e "default" entries.
* `ACL.Type.NFS4` - Include NFS4 entries.
* `ACL.Style.EXTRA_ID` - Include extra numeric ID field in
  each ACL entry. ('star' introduced this for POSIX.1e, this flag
  also applies to NFS4.)
* `ACL.Style.MARK_DEFAULT` - Include "default:" before each
  default ACL entry, as used in old Solaris ACLs.
"""
acl_text(entry::Entry, flags) =
    unsafe_string(ccall((:archive_entry_acl_to_text, libarchive), Ptr{UInt8},
                        (Ptr{Cvoid}, Cint), entry, flags))

acl_from_text(entry::Entry, text, typ) =
    @_la_call(archive_entry_acl_from_text, (Ptr{Cvoid}, Cstring, Cint), entry, text, typ)

"Return a count of entries matching `want`"
acl_count(entry::Entry, want) =
    ccall((:archive_entry_acl_count, libarchive), Cint, (Ptr{Cvoid}, Cint),
          entry, want)

"Return bitmask of ACL types in an archive entry"
acl_types(entry::Entry) =
    ccall((:archive_entry_acl_types, libarchive), Cint, (Ptr{Cvoid},), entry)

# Return an opaque ACL object.
# There's not yet anything clients can actually do with this...
# acl(entry::Entry) =
#     ccall((:archive_entry_acl, libarchive), Ptr{archive_acl},
#           (Ptr{Cvoid},), entry)

# extended attributes
xattr_clear(entry::Entry) =
    ccall((:archive_entry_xattr_clear, libarchive), Cvoid, (Ptr{Cvoid},), entry)
xattr_add_entry(entry::Entry, name::AbstractString, value) =
    ccall((:archive_entry_xattr_add_entry, libarchive), Cvoid,
          (Ptr{Cvoid}, Cstring, Ptr{Cvoid}, Csize_t),
          entry, name, value, sizeof(value))

# To retrieve the xattr list, first "reset", then repeatedly ask for the
# "next" entry.
xattr_count(entry::Entry) =
    ccall((:archive_entry_xattr_count, libarchive), Cint, (Ptr{Cvoid},), entry)
xattr_reset(entry::Entry) =
    ccall((:archive_entry_xattr_reset, libarchive), Cint, (Ptr{Cvoid},), entry)
function xattr_next(entry::Entry)
    name = Ref{Ptr{UInt8}}()
    value = Ref{Ptr{Cvoid}}()
    len = Ref{Csize_t}()
    @_la_call(archive_entry_xattr_next,
              (Ptr{Cvoid}, Ptr{Ptr{UInt8}}, Ptr{Ptr{Cvoid}}, Ptr{Csize_t}),
              entry, name, value, len)
    buff = Vector{UInt8}(undef, len[])
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
          buff, value[], len[])
    unsafe_string(name[]), buff
end

# sparse
sparse_clear(entry::Entry) =
    ccall((:archive_entry_sparse_clear, libarchive), Cvoid, (Ptr{Cvoid},), entry)
sparse_add_entry(entry::Entry, offset, len) =
    ccall((:archive_entry_sparse_add_entry, libarchive), Cvoid,
          (Ptr{Cvoid}, Int64, Int64), entry, offset, len)

# To retrieve the xattr list, first "reset", then repeatedly ask for the
# "next" entry.
sparse_count(entry::Entry) =
    ccall((:archive_entry_sparse_count, libarchive), Cint, (Ptr{Cvoid},), entry)
sparse_reset(entry::Entry) =
    ccall((:archive_entry_sparse_reset, libarchive), Cint, (Ptr{Cvoid},), entry)
function sparse_next(entry::Entry)
    offset = Ref{Int64}()
    len = Ref{Int64}()
    @_la_call(archive_entry_sparse_next,
              (Ptr{Cvoid}, Ptr{Int64}, Ptr{Int64}), entry, offset, len)
    offset[], len[]
end

# TODO
# Utility to match up hardlinks.
#
# The 'struct archive_entry_linkresolver' is a cache of archive entries
# for files with multiple links.  Here's how to use it:
#   1. Create a lookup object with archive_entry_linkresolver_new()
#   2. Tell it the archive format you're using.
#   3. Hand each archive_entry to archive_entry_linkify().
#      That function will return 0, 1, or 2 entries that should
#      be written.
#   4. Call archive_entry_linkify(resolver, NULL) until
#      no more entries are returned.
#   5. Call archive_entry_linkresolver_free(resolver) to free resources.
#
# The entries returned have their hardlink and size fields updated
# appropriately.  If an entry is passed in that does not refer to
# a file with multiple links, it is returned unchanged.  The intention
# is that you should be able to simply filter all entries through
# this machine.
#
# To make things more efficient, be sure that each entry has a valid
# nlinks value.  The hardlink cache uses this to track when all links
# have been found.  If the nlinks value is zero, it will keep every
# name in the cache indefinitely, which can use a lot of memory.
#
# Note that archive_entry_size() is reset to zero if the file
# body should not be written to the archive.  Pay attention!
# struct archive_entry_linkresolver;

# There are three different strategies for marking hardlinks.
# The descriptions below name them after the best-known
# formats that rely on each strategy:
#
# "Old cpio" is the simplest, it always returns any entry unmodified.
#    As far as I know, only cpio formats use this.  Old cpio archives
#    store every link with the full body; the onus is on the dearchiver
#    to detect and properly link the files as they are restored.
# "tar" is also pretty simple; it caches a copy the first time it sees
#    any link.  Subsequent appearances are modified to be hardlink
#    references to the first one without any body.  Used by all tar
#    formats, although the newest tar formats permit the "old cpio" strategy
#    as well.  This strategy is very simple for the dearchiver,
#    and reasonably straightforward for the archiver.
# "new cpio" is trickier.  It stores the body only with the last
#    occurrence.  The complication is that we might not
#    see every link to a particular file in a single session, so
#    there's no easy way to know when we've seen the last occurrence.
#    The solution here is to queue one link until we see the next.
#    At the end of the session, you can enumerate any remaining
#    entries by calling archive_entry_linkify(NULL) and store those
#    bodies.  If you have a file with three links l1, l2, and l3,
#    you'll get the following behavior if you see all three links:
#           linkify(l1) => NULL   (the resolver stores l1 internally)
#           linkify(l2) => l1     (resolver stores l2, you write l1)
#           linkify(l3) => l2, l3 (all links seen, you can write both).
#    If you only see l1 and l2, you'll get this behavior:
#           linkify(l1) => NULL
#           linkify(l2) => l1
#           linkify(NULL) => l2   (at end, you retrieve remaining links)
#    As the name suggests, this strategy is used by newer cpio variants.
#    It's noticeably more complex for the archiver, slightly more complex
#    for the dearchiver than the tar strategy, but makes it straightforward
#    to restore a file using any link by simply continuing to scan until
#    you see a link that is stored with a body.  In contrast, the tar
#    strategy requires you to rescan the archive from the beginning to
#    correctly extract an arbitrary link.

# struct archive_entry_linkresolver *archive_entry_linkresolver_new(void);
# void archive_entry_linkresolver_set_strategy(
#     struct archive_entry_linkresolver *, Cint /* format_code */);
# void archive_entry_linkresolver_free(struct archive_entry_linkresolver *);
# void archive_entry_linkify(struct archive_entry_linkresolver *,
#     struct archive_entry **, struct archive_entry **);
# struct archive_entry *archive_entry_partial_links(
#     struct archive_entry_linkresolver *res, Cuint *links);

# Routines to bulk copy fields to/from a platform-native "struct
# stat."  Libarchive used to just store a struct stat inside of each
# archive_entry object, but this created issues when trying to
# manipulate archives on systems different than the ones they were
# created on.
# stat(entry::Entry) =
#     StatStruct(ccall((:archive_entry_stat, libarchive),
#                      Ptr{UInt8}, (Ptr{Cvoid},), entry))
# function set_stat(entry::Entry, stat)
#     ccall((:archive_entry_copy_stat, libarchive),
#           Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), entry, stat)
# end


text:
using ITensors
using ProfileView

# Examples of optimizing simple chained matrix multiplications,
# useful for getting an idea for the overhead.
function main(d=100)
  i = Index(d, "i")
  A = randomITensor(i', dag(i))

  #
  # 2 tensors
  #

  println("\n2 tensors")

  # Simple pair contraction
  @btime $A' * $A

  #
  # 3 tensors
  #

  println("\n3 tensors")

  println("Pick an optimal sequence explicitly")

  ITensors.disable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  @btime $A'' * $A' * $A
  # Alternative syntax
  #@btime $A'' * ($A' * $A)
  #@btime contract([$A'', $A', $A]; sequence = "right_associative")
  #@btime contract([$A'', $A', $A]; sequence = "left_associative")
  #@btime contract([$A'', $A', $A]; sequence = $([[1, 2], 3]))
  #@btime contract([$A'', $A', $A]; sequence = $([[2, 3], 1]))

  println("Pick a bad sequence explicitly")

  @btime $A'' * $A * $A'
  # Alternative syntax
  #@btime contract([$A'', $A', $A]; sequence = $([[1, 3], 2]))

  println("Let it optimize")

  ITensors.enable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  # Already starting from an optimal sequence, there
  # is still overhead to checking if there isn't a better sequence
  println("Starting from optimal sequence")
  @btime $A'' * $A' * $A
  # Find a better sequence
  println("Starting from non-optimal sequence")
  @btime $A'' * $A * $A'

  f(A, N) =
    for _ in 1:N
      A'' * A' * A
    end
  @profview f(A, 1e6)

  ITensors.disable_contraction_sequence_optimization()

  #
  # 4 tensors
  #

  println("\n4 tensors")

  println("Pick an optimal sequence explicitly")

  ITensors.disable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  @btime $A''' * $A'' * $A' * $A
  # Alternative syntax
  #@btime $A''' * ($A'' * ($A' * $A))
  #@btime contract([$A''', $A'', $A', $A]; sequence = $([[[1, 2], 3], 4]))

  println("Pick a bad sequence explicitly")

  @btime $A'' * $A * $A''' * $A'

  println("Let it optimize")

  ITensors.enable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  println("Starting from optimal sequence")
  @btime $A''' * $A'' * $A' * $A
  println("Starting from non-optimal sequence")
  @btime $A'' * $A * $A''' * $A'

  #
  # 5 tensors
  #

  println("\n5 tensors")

  println("Pick an optimal sequence explicitly")

  ITensors.disable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  @btime $A'''' * $A''' * $A'' * $A' * $A

  println("Pick a bad sequence explicitly")

  @btime $A'' * $A * $A''' * $A' * $A''''

  println("Let it optimize")

  ITensors.enable_contraction_sequence_optimization()
  @show ITensors.using_contraction_sequence_optimization()

  println("Starting from optimal sequence")
  @btime $A'''' * $A''' * $A'' * $A' * $A
  println("Starting from non-optimal sequence")
  @btime $A'' * $A'''' * $A * $A''' * $A'

  return ITensors.disable_contraction_sequence_optimization()
end