text:
# Specify model parameters and/or initial values to optimize
function get_search_index()::Tuple{Array{Int64,1},Array{Int64,1}}
    # parameters
    search_idx_params::Vector{Int} = [
        C.V1
        C.Km1
        C.V5
        C.Km5
        C.V10
        C.Km10
        C.n10
        C.p11
        C.p12
        C.p13
        C.V14
        C.Km14
        C.V15
        C.Km15
        C.KimDUSP
        C.KexDUSP
        C.V20
        C.Km20
        C.V21
        C.Km21
        C.V24
        C.Km24
        C.V25
        C.Km25
        C.KimRSK
        C.KexRSK
        C.V27
        C.Km27
        C.V28
        C.Km28
        C.V29
        C.Km29
        C.V30
        C.Km30
        C.V31
        C.Km31
        C.n31
        C.p32
        C.p33
        C.p34
        C.V35
        C.Km35
        C.V36
        C.Km36
        C.V37
        C.Km37
        C.KimFOS
        C.KexFOS
        C.V42
        C.Km42
        C.V43
        C.Km43
        C.V44
        C.Km44
        C.p47
        C.m47
        C.p48
        C.p49
        C.m49
        C.p50
        C.p51
        C.m51
        C.V57
        C.Km57
        C.n57
        C.p58
        C.p59
        C.p60
        C.p61
        C.KimF
        C.KexF
        C.p63
        C.KF31
        C.nF31
        C.a
    ]

    # initial values
    search_idx_initials::Vector{Int} = [
        # V.(variableName)
    ]

    return search_idx_params, search_idx_initials
end

function get_search_region()::Matrix{Float64}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
    search_param::Vector{Float64} = init_search_param(search_idx, p, u0)

    search_rgn::Matrix{Float64} = zeros(2, length(p)+length(u0))
    #=
    # Default: 0.1 ~ 10x
    for (i,j) in enumerate(search_idx[1])
        search_rgn[1,j] = search_param[i]*0.1  # lower bound
        search_rgn[2,j] = search_param[i]*10.0  # upper bound
    end

    # Default: 0.5 ~ 2x
    for (i,j) in enumerate(search_idx[2])
        search_rgn[1,j+length(p)] = search_param[i+length(search_idx[1])]*0.5  # lower bound
        search_rgn[2,j+length(p)] = search_param[i+length(search_idx[1])]*2.0  # upper bound
    end
    =#

    # search_rgn[:, C.param_name] = [lower_bound, upper_bound]
    # search_rgn[:, V.var_name+length(p)] = [lower_bound, upper_bound]

    search_rgn[:, C.V1] = [7.33e-2, 6.60e-01]
    search_rgn[:, C.Km1] = [1.83e+2, 8.50e+2]
    search_rgn[:, C.V5] = [6.48e-3, 7.20e+1]
    search_rgn[:, C.Km5] = [6.00e-1, 1.60e+04]
    search_rgn[:, C.V10] = [exp(-10), exp(10)]
    search_rgn[:, C.Km10] = [exp(-10), exp(10)]
    search_rgn[:, C.n10] = [1.00, 4.00]
    search_rgn[:, C.p11] = [8.30e-13, 1.44e-2]
    search_rgn[:, C.p12] = [8.00e-8, 5.17e-2]
    search_rgn[:, C.p13] = [1.38e-7, 4.84e-1]
    search_rgn[:, C.V14] = [4.77e-3, 4.77e+1]
    search_rgn[:, C.Km14] = [2.00e+2, 2.00e+6]
    search_rgn[:, C.V15] = [exp(-10), exp(10)]
    search_rgn[:, C.Km15] = [exp(-10), exp(10)]
    search_rgn[:, C.KimDUSP] = [2.20e-4, 5.50e-1]
    search_rgn[:, C.KexDUSP] = [2.60e-4, 6.50e-1]
    search_rgn[:, C.V20] = [4.77e-3, 4.77e+1]
    search_rgn[:, C.Km20] = [2.00e+2, 2.00e+6]
    search_rgn[:, C.V21] = [exp(-10), exp(10)]
    search_rgn[:, C.Km21] = [exp(-10), exp(10)]
    search_rgn[:, C.V24] = [4.77e-2, 4.77e+0]
    search_rgn[:, C.Km24] = [2.00e+3, 2.00e+5]
    search_rgn[:, C.V25] = [exp(-10), exp(10)]
    search_rgn[:, C.Km25] = [exp(-10), exp(10)]
    search_rgn[:, C.KimRSK] = [2.20e-4, 5.50e-1]
    search_rgn[:, C.KexRSK] = [2.60e-4, 6.50e-1]
    search_rgn[:, C.V27] = [exp(-10), exp(10)]
    search_rgn[:, C.Km27] = [1.00e+2, 1.00e+4]
    search_rgn[:, C.V28] = [exp(-10), exp(10)]
    search_rgn[:, C.Km28] = [exp(-10), exp(10)]
    search_rgn[:, C.V29] = [4.77e-2, 4.77e+0]
    search_rgn[:, C.Km29] = [2.93e+3, 2.93e+5]
    search_rgn[:, C.V30] = [exp(-10), exp(10)]
    search_rgn[:, C.Km30] = [exp(-10), exp(10)]
    search_rgn[:, C.V31] = [exp(-10), exp(10)]
    search_rgn[:, C.Km31] = [exp(-10), exp(10)]
    search_rgn[:, C.n31] = [1.00, 4.00]
    search_rgn[:, C.p32] = [8.30e-13, 1.44e-2]
    search_rgn[:, C.p33] = [8.00e-8, 5.17e-2]
    search_rgn[:, C.p34] = [1.38e-7, 4.84e-1]
    search_rgn[:, C.V35] = [4.77e-3, 4.77e+1]
    search_rgn[:, C.Km35] = [2.00e+2, 2.00e+6]
    search_rgn[:, C.V36] = [exp(-10), exp(10)]
    search_rgn[:, C.Km36] = [1.00e+2, 1.00e+4]
    search_rgn[:, C.V37] = [exp(-10), exp(10)]
    search_rgn[:, C.Km37] = [exp(-10), exp(10)]
    search_rgn[:, C.KimFOS] = [2.20e-4, 5.50e-1]
    search_rgn[:, C.KexFOS] = [2.60e-4, 6.50e-1]
    search_rgn[:, C.V42] = [4.77e-3, 4.77e+1]
    search_rgn[:, C.Km42] = [2.00e+2, 2.00e+6]
    search_rgn[:, C.V43] = [exp(-10), exp(10)]
    search_rgn[:, C.Km43] = [1.00e+2, 1.00e+4]
    search_rgn[:, C.V44] = [exp(-10), exp(10)]
    search_rgn[:, C.Km44] = [exp(-10), exp(10)]
    search_rgn[:, C.p47] = [1.45e-4, 1.45e+0]
    search_rgn[:, C.m47] = [6.00e-3, 6.00e+1]
    search_rgn[:, C.p48] = [2.70e-3, 2.70e+1]
    search_rgn[:, C.p49] = [5.00e-5, 5.00e-1]
    search_rgn[:, C.m49] = [5.00e-3, 5.00e+1]
    search_rgn[:, C.p50] = [3.00e-3, 3.00e+1]
    search_rgn[:, C.p51] = [exp(-10), exp(10)]
    search_rgn[:, C.m51] = [exp(-10), exp(10)]
    search_rgn[:, C.V57] = [exp(-10), exp(10)]
    search_rgn[:, C.Km57] = [exp(-10), exp(10)]
    search_rgn[:, C.n57] = [1.00, 4.00]
    search_rgn[:, C.p58] = [8.30e-13, 1.44e-2]
    search_rgn[:, C.p59] = [8.00e-8, 5.17e-2]
    search_rgn[:, C.p60] = [1.38e-7, 4.84e-1]
    search_rgn[:, C.p61] = [exp(-10), exp(10)]
    search_rgn[:, C.KimF] = [2.20e-4, 5.50e-1]
    search_rgn[:, C.KexF] = [2.60e-4, 6.50e-1]
    search_rgn[:, C.p63] = [exp(-10), exp(10)]
    search_rgn[:, C.KF31] = [exp(-10), exp(10)]
    search_rgn[:, C.nF31] = [1.00, 4.00]
    search_rgn[:, C.a] = [1.00e+2, 5.00e+2]

    search_rgn = conv_lin2log!(search_rgn, search_idx)

    return search_rgn
end

function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Float64,1}}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()

    for (i,j) in enumerate(search_idx[1])
        @inbounds p[j] = indiv[i]
    end
    for (i,j) in enumerate(search_idx[2])
        @inbounds u0[j] = indiv[i+length(search_idx[1])]
    end

    # constraints --------------------------------------------------------------
    p[C.V6] = p[C.V5]
    p[C.Km6] = p[C.Km5]
    p[C.KimpDUSP] = p[C.KimDUSP]
    p[C.KexpDUSP] = p[C.KexDUSP]
    p[C.KimpcFOS] = p[C.KimFOS]
    p[C.KexpcFOS] = p[C.KexFOS]
    p[C.p52] = p[C.p47]
    p[C.m52] = p[C.m47]
    p[C.p53] = p[C.p48]
    p[C.p54] = p[C.p49]
    p[C.m54] = p[C.m49]
    p[C.p55] = p[C.p50]
    p[C.p56] = p[C.p51]
    p[C.m56] = p[C.m51]
    # --------------------------------------------------------------------------

    return p, u0
end

function decode_gene2val(indiv_gene::Vector{Float64})::Vector{Float64}
    search_rgn::Matrix{Float64} = get_search_region()
    indiv::Vector{Float64} = zeros(length(indiv_gene))

    for i in eachindex(indiv_gene)
        indiv[i] = 10^(
            indiv_gene[i] * (
                search_rgn[2,i] - search_rgn[1,i]
            ) + search_rgn[1,i]
        )
    end

    return round.(indiv,sigdigits=7)
end

function encode_val2gene(indiv::Vector{Float64})
    search_rgn::Matrix{Float64} = get_search_region()
    indiv_gene::Vector{Float64} = zeros(length(indiv))

    for i in eachindex(indiv)
        indiv_gene[i] = (
            log10(indiv[i]) - search_rgn[1,i]
        ) / (
            search_rgn[2,i] - search_rgn[1,i]
        )
    end

    return indiv_gene
end

function encode_bestIndivVal2randGene(
        idx::Int64,
        best_indiv::Vector{Float64},
        p0_bounds::Vector{Float64})::Float64
    search_rgn::Matrix{Float64} = get_search_region()
    rand_gene::Float64 = (
        log10(
            best_indiv[idx]*10^(
                rand() * log10(p0_bounds[2]/p0_bounds[1]) + log10(p0_bounds[1])
            )
        ) - search_rgn[1,idx]
    ) / (
        search_rgn[2,idx] - search_rgn[1,idx]
    )
    return rand_gene
end

function init_search_param(
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}},
        p::Vector{Float64},
        u0::Vector{Float64})::Vector{Float64}
    duplicate::Vector{String} = []
    if length(search_idx[1]) != length(unique(search_idx[1]))
        for idx in findall([count(x->x==i,search_idx[1])
                            for i in unique(search_idx[1])] .!= 1)
            push!(duplicate, C.NAMES[search_idx[1][idx]])
        end
        error(
            "Duplicate parameters (C.): $duplicate"
        )
    elseif length(search_idx[2]) != length(unique(search_idx[2]))
        for idx in findall([count(x->x==i,search_idx[2])
                            for i in unique(search_idx[2])] .!= 1)
            push!(duplicate, V.NAMES[search_idx[2][idx]])
        end
        error(
            "Duplicate species (V.): $duplicate"
        )
    end
    search_param = zeros(
        length(search_idx[1]) + length(search_idx[2])
    )
    for (i,j) in enumerate(search_idx[1])
        @inbounds search_param[i] = p[j]
    end
    for (i,j) in enumerate(search_idx[2])
        @inbounds search_param[i+length(search_idx[1])] = u0[j]
    end

    if any(x -> x == 0.0, search_param)
        msg::String = "search_param must not contain zero."
        for idx in search_idx[1]
            if p[idx] == 0.0
                error(
                    @sprintf(
                        "`C.%s` in search_idx_params: ", C.NAMES[idx]
                    ) * msg
                )
            end
        end
        for idx in search_idx[2]
            if u0[idx] == 0.0
                error(
                    @sprintf(
                        "`V.%s` in search_idx_initials: ", V.NAMES[idx]
                    ) * msg
                )
            end
        end
    end

    return search_param
end

function conv_lin2log!(
        search_rgn::Matrix{Float64},
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}})::Matrix{Float64}
    for i=1:size(search_rgn,2)
        if minimum(search_rgn[:,i]) < 0.0
            msg = "search_rgn[lower_bound,upper_bound] must be positive.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
            end
        elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0
            msg = "lower_bound must be larger than 0.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
            end
        elseif search_rgn[2,i] - search_rgn[1,i] < 0.0
            msg = "lower_bound must be smaller than upper_bound.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
            end
        end
    end

    nonzero_idx::Vector{Int} = []
    for i=1:size(search_rgn,2)
        if search_rgn[:,i] != [0.0,0.0]
            push!(nonzero_idx,i)
        end
    end
    difference::Vector{Int} = collect(
        symdiff(
            Set(nonzero_idx),
            Set(append!(search_idx[1], C.NUM .+ search_idx[2]))
        )
    )
    if length(difference) > 0
        for idx in difference
            if j <= C.NUM
                println(@sprintf("`C.%s`", C.NAMES[Int(idx)]))
            else
                println(@sprintf("`V.%s`", V.NAMES[Int(idx)-C.NUM]))
            end
        end
        error(
            "Set these search_params in both search_idx and search_rgn."
        )
    end

    search_rgn = search_rgn[:,nonzero_idx]

    return log10.(search_rgn)
end

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim


text:
module LinuxPerf

using Printf
using PrettyTables
using Formatting

export @measure, @measured, @pstats
export make_bench, enable!, disable!, reset!, reasonable_defaults, counters

import Base: show, length

macro measure(expr, args...)
    esc(quote
        local bench
        _, bench = $LinuxPerf.@measured($expr, $(args...))
        $counters(bench)
    end)
end
macro measured(expr, events = reasonable_defaults)
    quote
        local bench = $make_bench($events);
        local v
        try
            $enable!(bench)
            v = $(esc(expr))
        finally
            $disable!(bench)
        end
        (v, bench)
    end
end

mutable struct perf_event_attr
    typ::UInt32
    size::UInt32
    config::UInt64
    sample_period_or_freq::UInt64
    sample_type::UInt64
    read_format::UInt64
    flags::UInt64
    wakeup_events_or_watermark::UInt32
    bp_type::UInt32
    bp_addr_or_config1::UInt64
    bp_len_or_config2::UInt64
    branch_sample_type::UInt64

    sample_regs_user::UInt64
    sample_stack_user::UInt32
    clockid::Int32
    sample_regs_intr::UInt64
    aux_watermark::UInt32
    __reserved_2::UInt32

end

perf_event_attr() = perf_event_attr(ntuple(x->0, fieldcount(perf_event_attr))...)

const PERF_TYPE_HARDWARE = 0
const PERF_TYPE_SOFTWARE = 1
const PERF_TYPE_TRACEPOINT = 2
const PERF_TYPE_HW_CACHE = 3
const PERF_TYPE_RAW = 4
const PERF_TYPE_BREAKPOINT = 5

const EVENT_TYPES =
    [
     (:hw, PERF_TYPE_HARDWARE, # PERF_TYPE_HARDWARE
      [(:cycles, 0), # PERF_COUNT_HW_CPU_CYCLES
       (:instructions, 1), # PERF_COUNT_HW_INSTRUCTIONS
       (:cache_access, 2), # PERF_COUNT_HW_CACHE_REFERENCES
       (:cache_misses, 3), # PERF_COUNT_HW_CACHE_MISSES
       (:branches, 4), # PERF_COUNT_HW_BRANCH_INSTRUCTIONS
       (:branch_mispredicts, 5), # PERF_COUNT_HW_BRANCH_MISSES
       (:bus_cycles, 6), # PERF_COUNT_HW_BUS_CYCLES
       (:stalled_cycles_frontend, 7), # PERF_COUNT_HW_STALLED_CYCLES_FRONTEND
       (:stalled_cycles_backend, 8), # PERF_COUNT_HW_STALLED_CYCLES_BACKEND
       (:scaled_cycles, 9) # PERF_COUNT_HW_REF_CPU_CYCLES
       ]),
     (:sw, PERF_TYPE_SOFTWARE,
      [(:cpu_clock, 0), # PERF_COUNT_SW_CPU_CLOCK
       (:task_clock, 1), # PEF_COUNT_SW_TASK_CLOCK
       (:page_faults, 2), # PERF_COUNT_SW_PAGE_FAULTS
       (:ctx_switches, 3), # PERF_COUNT_SW_CONTEXT_SWITCHES
       (:cpu_migrations, 4), # PERF_COUNT_SW_CPU_MIGRATIONS
       (:minor_page_faults, 5), # PERF_COUNT_SW_PAGE_FAULTS_MIN
       (:major_page_faults, 6), # PERF_COUNT_SW_PAGE_FAULTS_MAJ
       (:alignment_faults, 7), # PERF_COUNT_SW_ALIGNMENT_FAULTS
       (:emulation_faults, 8), # PERF_COUNT_SW_EMULATION_FAULTS
       (:dummy, 9), # PERF_COUNT_SW_DUMMY
       (:bpf_output, 10), # PERF_COUNT_SW_BPF_OUTPUT
       ])
     ]

# cache events have special encoding, PERF_TYPE_HW_CACHE
const CACHE_TYPES =
    [(:L1_data, 0),
     (:L1_insn, 1),
     (:LLC,   2),
     (:TLB_data, 3),
     (:TLB_insn, 4),
     (:BPU, 5),
     (:NODE, 6)]
const CACHE_OPS =
    [(:read, 0),
     (:write, 1),
     (:prefetch, 2)]
const CACHE_EVENTS =
    [(:access, 0),
     (:miss, 1)]

const PERF_FORMAT_TOTAL_TIME_ENABLED = 1 << 0
const PERF_FORMAT_TOTAL_TIME_RUNNING = 1 << 1
const PERF_FORMAT_GROUP = 1 << 3

struct EventType
    category::UInt32
    event::UInt64
end

function all_events()
    evts = EventType[]
    for (cat_name, cat_id, events) in EVENT_TYPES
        for (type_name, type_id) in events
            push!(evts, EventType(cat_id, type_id))
        end
    end
    evts
end

function Base.show(io::IO, e::EventType)
    if e.category == PERF_TYPE_HW_CACHE
        print(io, "cache:")
        cache = e.event & 0xff
        idx = findfirst(k -> k[2] == cache, CACHE_TYPES)
        print(io, idx == 0 ? cache : CACHE_TYPES[idx][1], ":")
        cache_op = (e.event & 0xff00) >> 8
        idx = findfirst(k -> k[2] == cache_op, CACHE_OPS)
        print(io, idx == 0 ? cache : CACHE_OPS[idx][1], ":")
        cache_event = (e.event & 0xff0000) >> 16
        idx = findfirst(k -> k[2] == cache_event, CACHE_EVENTS)
        print(io, idx == 0 ? cache : CACHE_EVENTS[idx][1])
    else
        for (cat_name, cat_id, events) in EVENT_TYPES
            cat_id == e.category || continue
            print(io, cat_name, ":")
            for (type_name, type_id) in events
                type_id == e.event || continue
                print(io, type_name)
                return
            end
            print(io, e.event)
            return
        end
        print(io, "event(", e.category, ":", e.event, ")")
    end
end

const SYS_perf_event_open = 298

"""
    perf_event_open(attr::perf_event_attr, pid, cpu, fd, flags)
"""
function perf_event_open(attr::perf_event_attr, pid, cpu, leader_fd, flags)
    r_attr = Ref(attr)
    GC.@preserve r_attr begin
        # Have to do a manual conversion, since the ABI is a vararg call
        ptr = Base.unsafe_convert(Ptr{Cvoid}, Base.cconvert(Ptr{Cvoid}, r_attr))
        fd = ccall(:syscall, Cint, (Clong, Clong...), SYS_perf_event_open,
                   ptr, pid, cpu, leader_fd, flags)
    end
    return fd
end

function EventType(cat::Symbol, event::Symbol)
    cat !== :cache || error("cache events needs 3 arguments")
    for (cat_name, cat_id, events) in EVENT_TYPES
        cat_name === cat || continue
        for (type_name, type_id) in events
            type_name === event || continue
            return EventType(cat_id, type_id)
        end
        error("event $event not found in $cat")
    end
    error("category $cat not found")
end

function EventType(cat::Symbol, cache::Symbol, op::Symbol, evt::Symbol)
    cat === :cache || error("only cache events takes 3 arguments")
    idx = findfirst(x -> x[1] === cache, CACHE_TYPES)
    idx != 0 || error("cache not found $cache")
    cache_id = CACHE_TYPES[idx][2]
    idx = findfirst(x -> x[1] === op, CACHE_OPS)
    idx != 0 || error("op not found $op")
    op_id = CACHE_OPS[idx][2]
    idx = findfirst(x -> x[1] === evt, CACHE_EVENTS)
    idx != 0 || error("cache event not found $evt")
    evt_id = CACHE_EVENTS[idx][2]
    return EventType(PERF_TYPE_HW_CACHE,
                     cache_id | (op_id << 8) | (evt_id << 16))
end

const EXCLUDE_NONE       = UInt(0)
const EXCLUDE_USER       = UInt(1) << 0
const EXCLUDE_KERNEL     = UInt(1) << 1
const EXCLUDE_HYPERVISOR = UInt(1) << 2

function exclude_flags(u, k, h)
    exclude = EXCLUDE_NONE
    u || (exclude |= EXCLUDE_USER)
    k || (exclude |= EXCLUDE_KERNEL)
    h || (exclude |= EXCLUDE_HYPERVISOR)
    return exclude
end

struct EventTypeExt
    event::EventType
    modified::Bool
    exclude::UInt  # bit flags
end

mutable struct EventGroup
    leader_fd::Cint
    fds::Vector{Cint}
    event_types::Vector{EventType}
    leader_io::IOStream

    function EventGroup(types::Vector{<:Union{EventType,EventTypeExt}};
                        warn_unsupported = true,
                        userspace_only = true,
                        pinned = false,
                        exclusive = false,
                        pid = Cint(0),
                        )
        my_types = EventType[]
        group = new(-1, Cint[], EventType[])

        for (i, evt_type) in enumerate(types)
            attr = perf_event_attr()
            attr.size = sizeof(perf_event_attr)
            if evt_type isa EventTypeExt
                attr.typ = evt_type.event.category
                attr.config = evt_type.event.event
            else
                attr.typ = evt_type.category
                attr.config = evt_type.event
            end
            attr.sample_period_or_freq = 0
            attr.flags = 0
            # first attribute becomes group leader
            if group.leader_fd == -1
                attr.flags |= (1 << 0) # start disabled
            end
            if pinned
                attr.flags |= (1 << 2)
            end
            if exclusive
                attr.flags |= (1 << 3)
            end
            # (1 << 4) exclude_user
            if userspace_only
                attr.flags |= (1 << 5) # exclude kernel
            end
            # (1 << 6) exclude hypervisor
            # (1 << 7) exclude idle

            if evt_type isa EventTypeExt
                if evt_type.exclude & EXCLUDE_USER != 0
                    attr.flags |= (1 << 4)
                end
                if evt_type.exclude & EXCLUDE_KERNEL != 0
                    attr.flags |= (1 << 5)
                end
                if evt_type.exclude & EXCLUDE_HYPERVISOR != 0
                    attr.flags |= (1 << 6)
                end
            end

            attr.read_format =
                PERF_FORMAT_GROUP |
                PERF_FORMAT_TOTAL_TIME_ENABLED |
                PERF_FORMAT_TOTAL_TIME_RUNNING

            fd = perf_event_open(attr, pid, -1, group.leader_fd, 0)
            if fd < 0
                errno = Libc.errno()
                if errno in (Libc.EINVAL,Libc.ENOENT)
                    if warn_unsupported
                        @warn("$evt_type not supported, skipping")
                    end
                    continue
                else
                    if errno == Libc.EACCES && !userspace_only
                        @warn("try to adjust /proc/sys/kernel/perf_event_paranoid to a value <= 1 or use user-space only events")
                    end
                    error("perf_event_open error : $(Libc.strerror(errno))")
                end
            end
            if evt_type isa EventTypeExt
                push!(group.event_types, evt_type.event)
            else
                push!(group.event_types, evt_type)
            end
            push!(group.fds, fd)
            if group.leader_fd == -1
                group.leader_fd = fd
                group.leader_io = fdio(fd)
            end
        end
        reset!(group)
        group
    end
end

Base.length(g::EventGroup) = length(g.event_types)

function Base.show(io::IO, g::EventGroup)
    println(io, "EventGroup(")
    for e in g.event_types[1:end-1]
        println(io, "\t", e, ",")
    end
    print(io, "\t", g.event_types[end], ")")
end

const PERF_EVENT_IOC_ENABLE =  UInt64(0x2400)
const PERF_EVENT_IOC_DISABLE = UInt64(0x2401)
const PERF_EVENT_IOC_RESET =   UInt64(0x2403)

function ioctl(group::EventGroup, x)
    res = ccall(:ioctl, Cint, (Cint, Clong, Clong), group.leader_fd, x, 1)
    Base.systemerror(:ioctl, res < 0)
    return nothing
end

enable!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_ENABLE)
disable!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_DISABLE)
reset!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_RESET)

function Base.close(g::EventGroup)
    for fd in g.fds
        fd == g.leader_fd && continue # close leader_fd last
        ccall(:close, Cint, (Cint,), fd)
    end
    ccall(:close, Cint, (Cint,), g.leader_fd)
end

mutable struct PerfBench
    pid::Cint
    groups::Vector{EventGroup}
end

struct Counter
    event::EventType
    value::UInt64
    enabled::UInt64
    running::UInt64
end

struct Counters
    counters::Vector{Counter}
end

function Base.show(io::IO, c::Counters)
    events = map(x -> x.event, c.counters)
    stats  = mapreduce(vcat, c.counters) do c
        c.enabled == 0 ? ["never enabled" "0 %"] :
            c.running == 0 ? ["did not run" "0 %"] :
                [format(Int64(c.value), commas=true) @sprintf("%.1f %%", 100*(c.running/c.enabled))]
    end
    return pretty_table(io, stats, ["Events", "Active Time"], row_names=events, alignment=:l, crop=:none, body_hlines=collect(axes(stats, 1)))
end

enable!(b::PerfBench) = foreach(enable!, b.groups)
disable!(b::PerfBench) = foreach(disable!, b.groups)
reset!(b::PerfBench) = foreach(reset!, b.groups)

Base.close(b::PerfBench) = foreach(close, b.groups)

function counters(b::PerfBench)
    c = Counter[]
    for g in b.groups
        values = Vector{UInt64}(undef, length(g)+1+2)
        read!(g.leader_io, values)
        #?Ref@assert(length(g) == values[1])
        enabled, running = values[2], values[3]
        for i = 1:length(g)
            push!(c, Counter(g.event_types[i], values[3+i],
                             enabled, running))
        end
    end
    Counters(c)
end

const reasonable_defaults =
    [EventType(:hw, :cycles),
     [EventType(:hw, :cache_access),
      EventType(:hw, :cache_misses)],
     [EventType(:hw, :branches),
      EventType(:hw, :branch_mispredicts),
      EventType(:hw, :instructions)],
     [EventType(:sw, :ctx_switches),
      EventType(:sw, :page_faults),
      EventType(:sw, :minor_page_faults),
      EventType(:sw, :major_page_faults),
      EventType(:sw, :cpu_migrations)],
#=     [EventType(:cache, :L1_data, :read, :access),
      EventType(:cache, :L1_data, :read, :miss)],
     [EventType(:cache, :L1_data, :write, :access),
      EventType(:cache, :L1_data, :write, :miss)]=#]

function make_bench(x)
    groups = EventGroup[]
    for y in x
        if isa(y, EventType)
            push!(groups, EventGroup([y]))
        else
            push!(groups, EventGroup(y))
        end
    end
    PerfBench(0, groups)
end

make_bench() = make_bench(reasonable_defaults)

struct PerfBenchThreaded
    data::Vector{PerfBench}
end

enable!(b::PerfBenchThreaded) = foreach(enable!, b.data)
disable!(b::PerfBenchThreaded) = foreach(disable!, b.data)
reset!(b::PerfBenchThreaded) = foreach(reset!, b.data)

Base.close(b::PerfBenchThreaded) = foreach(close, b.data)

function make_bench_threaded(groups; threads = true)
    data = PerfBench[]
    for tid in (threads ? alltids() : zero(getpid()))
        push!(data, PerfBench(tid, [EventGroup(g, pid = tid, userspace_only = false) for g in groups]))
    end
    return PerfBenchThreaded(data)
end

alltids(pid = getpid()) = parse.(typeof(pid), readdir("/proc/$(pid)/task"))

# Event names are taken from the perf command.
const NAME_TO_EVENT = Dict(
    # hardware events
    "branch-instructions" => EventType(:hw, :branches),
    "branch-misses" => EventType(:hw, :branch_mispredicts),
    "cache-misses" => EventType(:hw, :cache_misses),
    "cache-references" => EventType(:hw, :cache_access),
    "cpu-cycles" => EventType(:hw, :cycles),
    "instructions" => EventType(:hw, :instructions),
    "stalled-cycles-backend" => EventType(:hw, :stalled_cycles_backend),
    "stalled-cycles-frontend" => EventType(:hw, :stalled_cycles_frontend),

    # software events
    "alignment-faults" => EventType(:sw, :alignment_faults),
    "bpf-output" => EventType(:sw, :bpf_output),
    "context-switches" => EventType(:sw, :ctx_switches),
    "cpu-clock" => EventType(:sw, :cpu_clock),
    "cpu-migrations" => EventType(:sw, :cpu_migrations),
    "dummy" => EventType(:sw, :dummy),
    "emulation-faults" => EventType(:sw, :emulation_faults),
    "major-faults" => EventType(:sw, :major_page_faults),
    "minor-faults" => EventType(:sw, :minor_page_faults),
    "page-faults" => EventType(:sw, :page_faults),
    "task-clock" => EventType(:sw, :task_clock),

    # hardware cache events
    "L1-dcache-load-misses" => EventType(:cache, :L1_data, :read, :miss),
    "L1-dcache-loads" => EventType(:cache, :L1_data, :read, :access),
    "L1-icache-load-misses" => EventType(:cache, :L1_insn, :read, :miss),
    "L1-icache-loads" => EventType(:cache, :L1_insn, :read, :access),
    "dTLB-load-misses" => EventType(:cache, :TLB_data, :read, :miss),
    "dTLB-loads" => EventType(:cache, :TLB_data, :read, :access),
    "iTLB-load-misses" => EventType(:cache, :TLB_insn, :read, :miss),
    "iTLB-loads" => EventType(:cache, :TLB_insn, :read, :access),
)
const EVENT_TO_NAME = Dict(event => name for (name, event) in NAME_TO_EVENT)

function is_supported(event::EventType; space::Symbol)
    attr = perf_event_attr()
    attr.typ = event.category
    attr.size = sizeof(perf_event_attr)
    attr.config = event.event
    if space == :user
        attr.flags |= (1 << 5)
        attr.flags |= (1 << 6)
    elseif space == :kernel
        attr.flags |= (1 << 4)
        attr.flags |= (1 << 6)
    elseif space == :hypervisor
        attr.flags |= (1 << 4)
        attr.flags |= (1 << 5)
    else
        throw(ArgumentError("unknown space name: $(space)"))
    end
    fd = perf_event_open(attr, 0, -1, -1, 0)
    if fd ≥ 0
        ret = ccall(:close, Cint, (Cint,), fd)
        if ret != 0
            @warn "failed to close file descriptor for some reason"
        end
        return true
    end
    return false
end

is_supported(name::AbstractString; kwargs...) = haskey(NAME_TO_EVENT, name) && is_supported(NAME_TO_EVENT[name]; kwargs...)

function list()
    for t in [PERF_TYPE_HARDWARE, PERF_TYPE_SOFTWARE, PERF_TYPE_HW_CACHE]
        events = collect(filter(x -> x[2].category == t, NAME_TO_EVENT))
        sort!(events, by = x -> x[1])  # sort events by name
        if t == PERF_TYPE_HARDWARE
            println("hardware:")
        elseif t == PERF_TYPE_SOFTWARE
            println("software:")
        elseif t == PERF_TYPE_HW_CACHE
            println("cache:")
        else
            @assert false
        end
        for (name, event) in events
            spaces = String[]
            if is_supported(event, space = :user)
                push!(spaces, "user")
            end
            if is_supported(event, space = :kernel)
                push!(spaces, "kernel")
            end
            if is_supported(event, space = :hypervisor)
                push!(spaces, "hypervisor")
            end
            if isempty(spaces)
                msg = "not supported"
            else
                msg = join(spaces, ", ")
            end
            @printf "  %-25s%s" name msg
            println()
        end
        t != PERF_TYPE_HW_CACHE && println()
    end
end

function parse_pstats_options(opts)
    # default events
    events = :(parse_groups("
        (cpu-cycles, stalled-cycles-frontend, stalled-cycles-backend),
        (instructions, branch-instructions, branch-misses),
        (task-clock, context-switches, cpu-migrations, page-faults)
    "))
    # default spaces
    user = true
    kernel = hypervisor = false
    # default threads
    threads = true
    for (i, opt) in enumerate(opts)
        if i == 1 && !(opt isa Expr && opt.head == :(=))
            events = :(parse_groups($(esc(opt))))
        elseif opt isa Expr && opt.head == :(=)
            key, val = opt.args
            val = esc(val)
            if key == :user
                user = val
            elseif key == :kernel
                kernel = val
            elseif key == :hypervisor
                hypervisor = val
            elseif key == :threads
                threads = val
            else
                error("unknown key: $(key)")
            end
        else
            error("unknown option: $(opt)")
        end
    end
    return (events = events, spaces = :($(user), $(kernel), $(hypervisor)), threads = threads,)
end

# syntax: groups = (group ',')* group
function parse_groups(str)
    groups = Vector{EventTypeExt}[]
    i = firstindex(str)
    next = iterate(str, i)
    while next !== nothing
        i = skipws(str, i)
        group, i = parse_group(str, i)
        push!(groups, group)
        i = skipws(str, i)
        next = iterate(str, i)
        if next === nothing
            continue
        end
        c, i = next
        if c == ','
            # ok
        else
            error("unknown character: $(repr(c))")
        end
    end
    return groups
end

# syntax: group = event | '(' (event ',')* event ')' modifiers?
function parse_group(str, i)
    group = EventTypeExt[]
    next = iterate(str, i)
    if next === nothing
        error("no events")
    elseif next[1] == '('
        # group
        i = next[2]
        while true
            i = skipws(str, i)
            event, i = parse_event(str, i)
            push!(group, event)
            i = skipws(str, i)
            next = iterate(str, i)
            if next === nothing
                error("unpaired '('")
            end
            c, i = next
            if c == ','
                # ok
            elseif c == ')'
                break
            else
                error("unknown character: $(repr(c))")
            end
        end
        i = skipws(str, i)

        # parse group-level modifiers (if any)
        next = iterate(str, i)
        if next !== nothing && next[1] == ':'
            (u, k, h), i = parse_modifiers(str, i)
            group = map(group) do event
                event.modified && return event
                exclude = exclude_flags(u, k, h)
                return EventTypeExt(event.event, true, exclude)
            end
        end
    else
        # singleton group
        i = skipws(str, i)
        event, i = parse_event(str, i)
        push!(group, event)
    end
    return group, i
end

# syntax: event = [A-Za-z0-9-]+ modifiers?
function parse_event(str, i)
    # parse event name
    isevchar(c) = 'A' ≤ c ≤ 'Z' || 'a' ≤ c ≤ 'z' || '0' ≤ c ≤ '9' || c == '-'
    start = i
    next = iterate(str, start)
    while next !== nothing && isevchar(next[1])
        i = next[2]
        next = iterate(str, i)
    end
    stop = prevind(str, i)
    if start > stop
        error("empty event name")
    end
    name = str[start:stop]
    if !haskey(NAME_TO_EVENT, name)
        error("unknown event name: $(name)")
    end
    event = NAME_TO_EVENT[name]

    # parse event-level modifiers (if any)
    modified = false
    exclude = EXCLUDE_NONE
    i = skipws(str, i)
    next = iterate(str, i)
    if next !== nothing && next[1] == ':'
        (u, k, h), i = parse_modifiers(str, i)
        modified = true
        exclude = exclude_flags(u, k, h)
    end

    return EventTypeExt(event, modified, exclude), i
end

# syntax: modifiers = ':' [ukh]*
function parse_modifiers(str, i)
    next = iterate(str, i)
    @assert next[1] == ':'
    ismodchar(c) = 'A' ≤ c ≤ 'Z' || 'a' ≤ c ≤ 'z'
    # u: user, k: kernel, h: hypervisor
    u = k = h = false  # exclude all
    i = skipws(str, next[2])
    next = iterate(str, i)
    while next !== nothing && ismodchar(next[1])
        c, i = next
        if c == 'u'
            u = true
        elseif c == 'k'
            k = true
        elseif c == 'h'
            h = true
        else
            error("unsupported modifier: $(repr(c))")
        end
        next = iterate(str, i)
    end
    return (u, k, h), i
end

# skip whitespace if any
function skipws(str, i)
    @label head
    next = iterate(str, i)
    if next !== nothing && isspace(next[1])
        i = next[2]
        @goto head
    end
    return i
end

struct ThreadStats
    pid::Cint
    groups::Vector{Vector{Counter}}
end

function ThreadStats(b::PerfBench)
    groups = Vector{Counter}[]
    for g in b.groups
        values = Vector{UInt64}(undef, length(g)+1+2)
        read!(g.leader_io, values)
        #?Ref@assert(length(g) == values[1])
        enabled, running = values[2], values[3]
        push!(groups, [Counter(g.event_types[i], values[3+i], enabled, running) for i in 1:length(g)])
    end
    return ThreadStats(b.pid, groups)
end

function Base.haskey(stats::ThreadStats, name::AbstractString)
    event = NAME_TO_EVENT[name]
    return any(counter.event == event for group in stats.groups for counter in group)
end

function Base.getindex(stats::ThreadStats, name::AbstractString)
    event = NAME_TO_EVENT[name]
    for group in stats.groups, counter in group
        if counter.event == event
            return counter
        end
    end
    throw(KeyError(name))
end

function Base.show(io::IO, stats::ThreadStats)
    println(io, stats.pid)
    printcounters(io, stats.groups)
end

isenabled(counter::Counter) = counter.enabled > 0
isrun(counter::Counter) = counter.running > 0
fillrate(counter::Counter) = counter.running / counter.enabled
scaledcount(counter::Counter) = counter.value * (counter.enabled / counter.running)

struct Stats
    threads::Vector{ThreadStats}
end

Stats(b::PerfBenchThreaded) = Stats(map(ThreadStats, b.data))

Base.show(io::IO, stats::Stats) = printsummary(io, stats)

printsummary(stats::Stats; kwargs...) = printsummary(stdout, stats; kwargs...)

"""
    printsummary([io,] stats::Stats; expandthreads = false, skipinactive = true)

Print summary of event statistics.

If `expandthreads` is `true`, the statistics of each thread are printed with
its thread ID (TID). If `skipinactive` is `true`, the statistics from
unmeasured (inactive) threads are ignored.
"""
function printsummary(io::IO, stats::Stats; expandthreads::Bool = false, skipinactive::Bool = true)
    printsep(io, '━')
    println(io)
    if isempty(stats.threads)
        print(io, "no threads")
        return
    end

    # aggregate all counts
    n_aggregated = 0
    counters = [[Counter(c.event, 0, 0, 0) for c in g] for g in stats.threads[1].groups]
    for t in stats.threads
        if skipinactive && !any(isenabled, c for g in t.groups for c in g)
            continue
        end
        n_aggregated += 1
        if expandthreads
            println(io, "Thread #$(n_aggregated) (TID = $(t.pid))")  # label
            printcounters(io, t.groups)
            printsep(io, '┄')
            println(io)
        end
        for (j, g) in enumerate(t.groups)
            for (k, c) in enumerate(g)
                c′ = counters[j][k]
                @assert c′.event == c.event
                counters[j][k] = Counter(
                    c.event,
                    c.value   + c′.value,
                    c.enabled + c′.enabled,
                    c.running + c′.running
                )
            end
        end
    end

    for g in counters, c in g
        if !isrun(c)
            @warn "Some events are not measured"
            break
        end
    end

    expandthreads && n_aggregated > 1 && println(io, "Aggregated")  # label
    printcounters(io, counters)
    if n_aggregated > 1
        println(io, lpad("aggregated from $(n_aggregated) threads", TABLE_WIDTH))
    end
    printsep(io, '━')
end

const TABLE_WIDTH = 2 + 23 + 18
printsep(io::IO, c::Char) = print(io, c^TABLE_WIDTH)

function printcounters(io::IO, groups::Vector{Vector{Counter}})
    for group in groups
        function findcount(name)
            event = NAME_TO_EVENT[name]
            # try to find within the same group
            for c in group
                c.event == event && return c
            end
            # fall back to other groups
            for g in groups, c in g
                c.event == event && return c
            end
            return nothing
        end
        for (i, counter) in enumerate(group)
            # grouping character
            c = length(group) == 1 ? '╶' :
                i == 1             ? '┌' :
                i == length(group) ? '└' : '│'
            event = counter.event
            name = EVENT_TO_NAME[event]
            @printf io "%-2s%-23s" c name
            if !isenabled(counter)
                @printf(io, "%18s", "not enabled")
            elseif !isrun(counter)
                @printf(io, "%10s%7.1f%%", "NA", 0.0)
            else
                @printf(io, "%10.2e%7.1f%%", scaledcount(counter), fillrate(counter) * 100)
            end
            if isrun(counter)
                # show a comment
                if name == "cpu-cycles"
                    @printf(io, "  # %4.1f cycles per ns", counter.value / counter.running)
                elseif name == "instructions" && (cycles = findcount("cpu-cycles")) !== nothing
                    @printf(io, "  # %4.1f insns per cycle", scaledcount(counter) / scaledcount(cycles))
                elseif name == "cpu-clock" || name == "task-clock"
                    clk = float(scaledcount(counter))
                    if clk ≥ 1e9
                        clk /= 1e9
                        unit = "s"
                    elseif clk ≥ 1e6
                        clk /= 1e6
                        unit = "ms"
                    elseif clk ≥ 1e3
                        clk /= 1e3
                        unit = "μs"
                    else
                        unit = "ns"
                    end
                    @printf(io, "  # %4.1f %s", clk, unit)
                else
                    for (num, den, label) in [
                            ("stalled-cycles-frontend", "cpu-cycles", "cycles"),
                            ("stalled-cycles-backend", "cpu-cycles", "cycles"),
                            ("branch-instructions", "instructions", "instructions"),
                            ("branch-misses", "branch-instructions", "branch instructions"),
                            ("cache-misses", "cache-references", "cache references"),
                            ("L1-dcache-load-misses", "L1-dcache-loads", "dcache loads"),
                            ("L1-icache-load-misses", "L1-icache-loads", "icache loads"),
                            ("dTLB-load-misses", "dTLB-loads", "dTLB loads"),
                            ("iTLB-load-misses", "iTLB-loads", "iTLB loads"),
                        ]
                        if name == num && (d = findcount(den)) !== nothing
                            @printf(io, "  # %4.1f%% of %s", scaledcount(counter) / scaledcount(d) * 100, label)
                            break
                        end
                    end
                end
            end
            println(io)
        end
    end
end

function set_default_spaces(groups, (u, k, h))
    map(groups) do group
        map(group) do event
            if event.modified
                return event
            end
            exclude = EXCLUDE_NONE
            u || (exclude |= EXCLUDE_USER)
            k || (exclude |= EXCLUDE_KERNEL)
            h || (exclude |= EXCLUDE_HYPERVISOR)
            return EventTypeExt(event.event, event.modified, exclude)
        end
    end
end

# for debug
function dump_groups(groups)
    buf = IOBuffer()
    println(buf, "Event groups")
    for (i, group) in enumerate(groups)
        println(buf, "group #$(i)")
        for (j, event) in enumerate(group)
            if j < length(group)
                print(buf, "  ├─ ")
            else
                print(buf, "  └─ ")
            end
            print(buf, event.event)
            event.modified && print(buf, " → modified")
            exclude = event.exclude
            exclude != 0 && print(buf, ", exclude ")
            exclude & EXCLUDE_USER       != 0 && print(buf, 'u')
            exclude & EXCLUDE_KERNEL     != 0 && print(buf, 'k')
            exclude & EXCLUDE_HYPERVISOR != 0 && print(buf, 'h')
            println(buf)
        end
    end
    String(take!(buf))
end

"""
    @pstats [options] expr

Run `expr` and gather its performance statistics.

This macro basically measures the number of occurrences of events such as CPU
cycles, branch prediction misses, page faults, and so on. The list of
supported events can be shown by calling the `LinuxPerf.list` function.

Due to the resource limitation of performance measuring units (PMUs)
installed in a CPU core, all events may not be measured simultaneously,
resulting in multiplexing several groups of events in a single measurement.
If the running time is extremely short, some event groups may not be measured
at all.

The result is shown in a table. Each row consists of four columns: an event
group indicator, an event name, a scaled count and a running rate. A comment
may follow these columns after a hash (#) character.
1. The event group indicated by a bracket is a set of events that are
   measured simultaneously so that their count statistics can be meaningfully
   compared.
2. The event name is a conventional name of the measured event.
3. The scaled count is the number of occurrences of the event, scaled by the
   reciprocal of the running rate.
4. The running rate is the ratio of the time of running and enabled.

The macro can take some options. If a string object is passed, it is a
comma-separated list of event names to measure. A group of events is
surrounded by a pair of parentheses. Modifiers can be added to confine
measured events to specific space. Currently, three space modifiers are
supported: user (`u`), kernel (`k`), and hypervisor (`h`) space. The
modifiers follow an event name separated by a colon. For example,
`cpu-cycles:u` ignores all CPU cycles except in user space (which is the
default). It is also possible to pass `user`, `kernel`, and `hypervisor`
parameters to the macro, which affect events without modifiers. Only user
space is activated by default (i.e., `user` is `true` but `kernel` and
`hypervisor` are `false`). To measure kernel events, for example, add the `k`
modifier to events you are interested in or pass `kernel=true` to the macro,
which globally activates events in kernel space.

All threads are measured and event counts are aggregated by default. Passing
`threads=false` to the macro disables this feature and only measures events
that occurred in the current thread invoking the macro.

For more details, see perf_event_open(2)'s manual page.

# Examples

```
julia> xs = randn(1_000_000);

julia> sort(xs[1:9]);  # compile

julia> @pstats sort(xs)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
┌ cpu-cycles               2.60e+08   49.7%  #  3.4 cycles per ns
│ stalled-cycles-frontend  1.09e+07   49.7%  #  4.2% of cycles
└ stalled-cycles-backend   7.07e+06   49.7%  #  2.7% of cycles
┌ instructions             1.96e+08   50.3%  #  0.8 insns per cycle
│ branch-instructions      4.02e+07   50.3%  # 20.5% of instructions
└ branch-misses            8.15e+06   50.3%  # 20.3% of branch instructions
┌ task-clock               7.61e+07  100.0%  # 76.1 ms
│ context-switches         7.00e+00  100.0%
│ cpu-migrations           0.00e+00  100.0%
└ page-faults              1.95e+03  100.0%
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

julia> @pstats "(cpu-cycles,instructions,branch-instructions,branch-misses),page-faults" sort(xs)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
┌ cpu-cycles               2.64e+08  100.0%  #  3.5 cycles per ns
│ instructions             1.86e+08  100.0%  #  0.7 insns per cycle
│ branch-instructions      3.74e+07  100.0%  # 20.1% of instructions
└ branch-misses            8.21e+06  100.0%  # 21.9% of branch instructions
╶ page-faults              1.95e+03  100.0%
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
```
"""
macro pstats(args...)
    if isempty(args)
        error("@pstats requires at least one argument")
    end
    opts, expr = parse_pstats_options(args[1:end-1]), args[end]
    quote
        (function ()
            groups = set_default_spaces($(opts.events), $(opts.spaces))
            @debug dump_groups(groups)
            bench = make_bench_threaded(groups, threads = $(opts.threads))
            try
                enable!(bench)
                val = $(esc(expr))
                disable!(bench)
                # trick the compiler not to eliminate the code
                stats = rand() < 0 ? val : Stats(bench)
                return stats::Stats
            catch
                rethrow()
            finally
                close(bench)
            end
        end)()
    end
end

end

ext: jl


text:
using Base.Iterators: Iterators, product
using ..Datasets: getsavefolder, Cube

"""
    savecube(cube,name::String)

Save a [`YAXArray`](@ref) to the folder `name` in the YAXArray working directory.
"""
function savecube(
    c,
    name::AbstractString;
    chunksize = Dict(),
    max_cache = 1e8,
    backend = :zarr,
    backendargs...,
)
    allax = caxes(c)
    if !isa(chunksize, AbstractDict) && !isa(chunksize, NamedTuple)
        @warn "Chunksize must be provided as a Dict mapping axis names to chunk size in the future"
        chunksize = OrderedDict(axname(i[1]) => i[2] for i in zip(allax, chunksize))
    end
    firstaxes = sort!([findAxis(String(k), allax) for k in keys(chunksize)])
    lastaxes = setdiff(1:length(allax), firstaxes)
    allax = allax[[firstaxes; lastaxes]]
    dl = cumprod(length.(caxes(c)))
    isplit = findfirst(i -> i > max_cache / sizeof(eltype(c)), dl)
    isplit isa Nothing && (isplit = length(dl) + 1)
    forcesingle = (isplit + 1) < length(firstaxes)
    axin = allax[1:isplit-1]
    axn = axname.(axin)
    indims = InDims(axn...)
    outdims = OutDims(
        axn...,
        backend = backend,
        chunksize = chunksize,
        path = name;
        backendargs...,
    )
    function cop(xout, xin)
        xout .= xin
    end
    if forcesingle
        nprocs() > 1 && println("Forcing single core processing because of bad chunk size")
        o = mapCube(
            cop,
            c,
            indims = indims,
            outdims = outdims,
            ispar = false,
            max_cache = max_cache,
        )
    else
        o = mapCube(cop, c, indims = indims, outdims = outdims, max_cache = max_cache)
    end
end

function loadcube(s)
    Cube(getsavefolder(s, true))
end

function rmcube(s)
    p = getsavefolder(s, true)
    if isfile(p) || isdir(p)
        rm(p, recursive = true)
    end
end

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
module MultiBasisFESpaces
##
using Gridap.DOFBases: LagrangianDOFBasis

using Gridap, Test

using Gridap.Helpers

using Base.Cartesian

import Gridap: ∇
# 1D reffe

include("ScaledNodesLagrangianRefFEs.jl")
include("MultiBasisFESpaces.jl")
F = Gridap.ConformingFESpaces

D = 2
n = 1
# ℓ = 1.0
β = 400.0
# h = 1.0 / n
# h0 = 1.0 / (n-1+β)
# r = h/h0
# e = h/(ℓ-h0*(n-1))
# ℓ0 = n*h
# geomap = ( x -> ( x[1] <= (n-1)*h ) ? x[1]/r : (n-1)*h0+(x[1]-(n-1)*h)/e )
# geomap = ( x -> x[1] <= (ℓ-h) ? x[1] : (x[1]-ℓ+h)*β+ℓ-h)
# newl = geomap(1.0)
# typeof(newl)
geomap = ( x -> x )
model = CartesianDiscreteModel(partition=(n,n),map=geomap)

# Construct the discrete model
ufun(x) = x[1]*(x[1]-2.0)
ufun_grad(x) = VectorValue(2*x[1]-2.0,0.0)
bfun(x) = -2.0
∇(::typeof(ufun)) = ufun_grad

function solve_problem(β,model,ufun)

  # diri = [7]
  diri = [1,2,3,4,5,6,7,8]
  diri = "boundary"
  p = Polytope(1,1)
  order = 2

  trian = Triangulation(model)
  reffe = LagrangianRefFE(Float64,p,order)
  sreffe = ScaledNodesLagrangianRefFE(Float64,p,order,1.0/β)
  values = LagrangianRefFE{2,Float64}[reffe,sreffe]
  pointer = [ i == ncells(trian) ? 2 : 1 for i in 1:ncells(trian)]
  reffes = CompressedCellValue(values,pointer)
  typeof(reffes)

  graph = GridGraph(model)
  labels = FaceLabels(model)
  mydiri_tags = Gridap.BoundaryGrids._setup_tags(labels,diri)

  fespace = MultiShapeConformingFESpace(reffes,trian,graph,labels,mydiri_tags)

  V = TestFESpace(fespace)
  U = TrialFESpace(fespace,ufun)
  trian = Triangulation(model)
  quad = CellQuadrature(trian,degree=order*2)

  a(v,u) = inner(∇(v), ∇(u))
  b(v) = inner(v,bfun)

  assem = SparseMatrixAssembler(V,U)
  op = LinearFEOperator(a,b,V,U,assem,trian,quad)
  ls = BackslashSolver()
  solver = LinearFESolver(ls)
  uh = solve(solver,op)

  # Define exact solution and error
  u = CellField(trian,ufun)
  e = u - uh

  # Define norms to measure the error
  l2(u) = inner(u,u)
  h1(u) = a(u,u) + l2(u)
  # Compute errors
  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
  eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

  @show el2 #< 1.e-8
  @show eh1 #< 1.e-8

  return op
end

op = solve_problem(1,model,ufun)
op.mat

op2 = solve_problem(100,model,ufun)
op2.mat

# writevtk(trian,"onedres",cellfields=["uh"=>uh,"e"=>e])
##
end # module


text:
using Test

@testset "Learning ConstrainedRootSolvers" begin
    @test true;
end

max_stars_repo_name: Plant-Form-Function/juliant-template

max_issues_repo_name: Plant-Form-Function/juliant-template

max_forks_repo_name: Plant-Form-Function/juliant-template


text:
using PrettyTables

Base.show(io::IO, model::LocalDephasingModel) = begin
    type = typeof(model.Jx)
    m = Base.format_bytes(Base.summarysize(model))
    println(io, "Continuous monitoring model with $(model.params.Nj)-spin Dicke state\n")

    println(io, "Operator type: $(typeof(model.M0))")
    sup_size = model.second_term.size
    nvals = length(model.second_term.values)
    density = nvals / (*(sup_size...))^2
    println(io, "Superoperator: $sup_size with $nvals stored entries (density $density)")

    println(io, "Total memory: $m\n")

    show(io, model.params)
end

Base.show(io::IO, model::CollectiveDephasingModel) = begin
    type = typeof(model.Jx)
    m = Base.format_bytes(Base.summarysize(model))
    println(io, "Continuous monitoring model with $(model.params.Nj)-spin Dicke state\n")

    println(io, "Operator type: $(typeof(model.M0))")
    sup_size = size(model.second_term)
    nvals = length(model.second_term.nzval)
    density = nvals / (*(sup_size...))
    println(io, "Superoperator: $sup_size with $nvals stored entries (density $density)")

    println(io, "Total memory: $m\n")

    show(io, model.params)
end

Base.show(io::IO, state::BlockDiagonalState) = begin
    n = size(state.ρ, 1)
    type = typeof(state.ρ)
    m = Base.format_bytes(Base.summarysize(state))
    println(io, "Dicke state for $(nspins(n)) spins. $n × $n matrix of type $type. Total memory: $m")
end

Base.show(io::IO, state::FixedjState) = begin
    n = size(state.ρ, 1)
    type = typeof(state.ρ)
    m = Base.format_bytes(Base.summarysize(state))
    println(io, "Dicke state: $n × $n matrix of type $type. Total memory: $m")
end

Base.show(io::IO, x::ModelParameters) = begin
    println(io, "Model Parameters")
    data = ["Nj"    x.Nj    "Number of spins";
            "κ_i"   x.kind  "Independent noise rate";
            "κ_c"   x.Gamma "Collective noise rate";
            "ω"     x.omega "Frequency";
            "η"     x.eta   "Measurement efficiency"]
    pretty_table(io, data, ["Parameter", "Value", "Description"], alignment=[:l, :r, :l])

    println(io, "\nSimulation Parameters")
    data = ["Tfinal" x.Tfinal "Final time";
            "dt"     x.dt     "Time step";
            "Ntime"  x.Ntime  "N. of simulation points";
            "outpoints" x.outpoints "N. of output points"]
    pretty_table(io, data, ["Parameter", "Value", "Description"], alignment=[:l, :r, :l])
end


text:
export AbstractSpace, AbstractContinuousSpace, AbstractDiscreteSpace,
       sample

import StatsBase:sample
import Base:in, size, ==, length, eltype

"""
    AbstractSpace

Supertype of [`AbstractContinuousSpace`](@ref) and [`AbstractDiscreteSpace`](@ref).
"""
abstract type AbstractSpace end

"""
    AbstractContinuousSpace <: AbstractSpace

Supertype of different kinds of continuous spaces.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `sample(space)` | Get a random sample from the space |
| `Base.in(space, x)` | Test whether `x` is in the `space` |
| **Optional Methods** | |
| `Base.eltype(space)` | Return the type of the sample in a space |

See also: [`AbstractDiscreteSpace`](@ref)
"""
abstract type AbstractContinuousSpace <: AbstractSpace end

"""
    AbstractDiscreteSpace <: AbstractSpace

Supertype of different kinds of discrete spaces.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `sample(space)` | Get a random sample from the space |
| `Base.in(space, x)` | Test whether `x` is in the `space` |
| `Base.size(space)` | Return the size of the space in all dimensions |
| **Optional Methods** | |
| `Base.length(space)` | Return the number of elements in that space. By default it will be `*(Base.size(space))`. |
| `Base.eltype(space)` | Return the type of the sample in a space |

See also: [`AbstractContinuousSpace`](@ref)
"""
abstract type AbstractDiscreteSpace <: AbstractSpace end

length(s::AbstractDiscreteSpace) = *(size(s)...)

# Tuple Support
sample(s::Tuple{Vararg{<:AbstractSpace}}) = map(sample, s)
in(a::Tuple{Vararg{T, N} where T}, b::Tuple{Vararg{<:AbstractSpace, N}}) where N = all(map((x, y) -> in(x, y), a, b))

# Dict Support
sample(s::Dict{String}) = Dict(map((k, v) -> (k, sample(v)), s))
in(a::Dict{String}, b::Dict{String}) = length(a) == length(b) &&
    all(p -> haskey(a, p.first) ? in(a[p.first], p.second) : false,
        b)


text:
immutable EnvironmentRegion{N, T}
    position::SimpleHRepresentation{N, T}
    force::SimpleHRepresentation{N, T}
end
dimension{N, T}(::Type{EnvironmentRegion{N, T}}) = N

function isfree(region::EnvironmentRegion)
    vr = vrep(polyhedron(region.force))
    !hasrays(vr) && all(all(p .== 0) for p in points(vr))
end

function lphrep_with_bounds(vr::VRepresentation)
    hr = hrep(polyhedron(vr, CDDLibrary()))
    lphr = LPHRepresentation(hr)
    l, u = axis_aligned_bounding_box(vr)
    LPHRepresentation(lphr.A, l, u, lphr.lb, lphr.ub)
end

# Forces
function friction_cone_constraints(n, μ, fnmax)
    n⟂ = perp(n)
    β1 = normalize(n + μ * n⟂) * fnmax
    β2 = normalize(n - μ * n⟂) * fnmax
    vr = SimpleVRepresentation([zeros(1, length(n⟂)); β1'; β2'])
    SimpleHRepresentation(hrep(polyhedron(vr, CDDLibrary())))
end

function no_force_constraints(n)
    SimpleHRepresentation(eye(2), zeros(2), IntSet([1, 2]))
end

# Positions
function line_segment_constraints(p1, p2)
    vr = SimpleVRepresentation([p1'; p2'])
    SimpleHRepresentation(hrep(polyhedron(vr, CDDLibrary())))
end

function axis_aligned_free_box(p1, p2)
    l = min.(p1, p2)
    u = max.(p1, p2)
    n = length(l)
    SimpleHRepresentation([eye(n); -eye(n)], [u; -l])
end

# Position × Force environment regions
function contact_region(p1, p2, μ, fnmax)
    n = normalize(perp(p2 - p1))
    position = line_segment_constraints(p1, p2)
    force = friction_cone_constraints(n, μ, fnmax)
    EnvironmentRegion(position, force)
end

function axis_aligned_free_box_region(p1, p2)
    position = axis_aligned_free_box(p1, p2)
    force = no_force_constraints(length(p1))
    EnvironmentRegion(position, force)
end

function force_bounds{T<:EnvironmentRegion}(environment::AbstractVector{T})
    # TODO: could do this on a contact-by-contact basis if certain contact-region assignments are not allowed
    n = dimension(T)
    fmin = zeros(n)
    fmax = zeros(n)
    for region in environment
        fl, fu = axis_aligned_bounding_box(region.force)
        fmin .= min.(fmin, fl)
        fmax .= max.(fmax, fu)
    end
    fmin, fmax
end


text:
using ImageEdgeDetection
using Test

@testset "ImageEdgeDetection.jl" begin
    # Write your own tests here.
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
#
# Some examples of how to use the Chaboche material model.

using Parameters
using ForwardDiff
using DelimitedFiles, Test

using Materials

function simple_integration_test()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    chabmat = Chaboche(parameters=parameters, ddrivers=ddrivers)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
end
simple_integration_test()

function test_chaboche()
    path = joinpath(@__DIR__, "one_elem_disp_chaboche", "unitelement_results.rpt")
    data = readdlm(path, Float64; skipstart=4)
    ts = data[:,1]
    s11_ = data[:,2]
    s12_ = data[:,3]
    s13_ = data[:,4]
    s22_ = data[:,5]
    s23_ = data[:,6]
    s33_ = data[:,7]
    e11_ = data[:,8]
    e12_ = data[:,9]
    e13_ = data[:,10]
    e22_ = data[:,11]
    e23_ = data[:,12]
    e33_ = data[:,13]

    strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)]

    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)
    chabmat = Chaboche(parameters=parameters)
    s33s = [chabmat.variables.stress[3,3]]
    for i=2:length(ts)
        dtime = ts[i]-ts[i-1]
        dstrain = fromvoigt(Symm2{Float64}, strains[i]-strains[i-1]; offdiagscale=2.0)
        chabmat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain)
        integrate_material!(chabmat)
        update_material!(chabmat)
        push!(s33s, chabmat.variables.stress[3,3])
    end
    @test isapprox(s33s, s33_; rtol=0.01)
end
test_chaboche()

# Profile.clear_malloc_data()
# test_chaboche()
# using BenchmarkTools
# @btime test_chaboche()

function simple_integration_test_fd_tangent()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    chabmat = Chaboche(parameters=parameters, ddrivers=ddrivers)

    function get_stress(dstrain::Symm2)
        chabmat.ddrivers.strain = dstrain
        integrate_material!(chabmat)
        return chabmat.variables_new.stress
    end

    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    # https://kristofferc.github.io/Tensors.jl/stable/man/automatic_differentiation/
    # TODO: doesn't work, a Nothing ends up in the type for some reason?
    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D_mat = $(tovoigt(chabmat.variables_new.jacobian))"
    @info "D = $(tovoigt(D))"

    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
    chabmat.variables_new = typeof(chabmat.variables_new)()
    chabmat.ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    # stress = get_stress(0.25*dstrain_dtime)
    # @info "stress = $stress"

    D, dstress = Tensors.gradient(get_stress, 0.25*dstrain_dtime, :all)
    @info "D = $(tovoigt(D))"
end
simple_integration_test_fd_tangent()

function simple_integration_test_fd_tangent2()
    parameters = ChabocheParameterState(E = 200.0e3,
                                        nu = 0.3,
                                        R0 = 100.0,
                                        Kn = 100.0,
                                        nn = 10.0,
                                        C1 = 10000.0,
                                        D1 = 100.0,
                                        C2 = 50000.0,
                                        D2 = 1000.0,
                                        Q = 50.0,
                                        b = 0.1)

    dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
    ddrivers = ChabocheDriverState(time=0.25, strain=0.25*dstrain_dtime)
    chabmat = Chaboche(parameters=parameters, ddrivers=ddrivers)
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    g! = Materials.ChabocheModule.create_nonlinear_system_of_equations(chabmat)

    x0 = [tovoigt(chabmat.variables_new.stress); chabmat.variables_new.R; tovoigt(chabmat.variables_new.X1); tovoigt(chabmat.variables_new.X2)]
    drdx = ForwardDiff.jacobian(debang(g!), x0)
    @info "size(drdx) = $(size(drdx))"
    @info "drdx = $drdx"
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = parameters
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))
    jacobian = isotropic_elasticity_tensor(lambda, mu)
    drde = zeros((19,6))
    drde[1:6, 1:6] = -tovoigt(jacobian)
    @info "drde = $drde"
    @info "size(drde) = $(size(drde))"
    jacobian2 = (drdx\drde)[1:6, 1:6]
    @info "jacobian = $(tovoigt(jacobian))"
    @info "jacobian2 = $jacobian2"
    jacobian3 = (drdx[1:6, 1:6] +  drdx[1:6,7:end]*(drdx[7:end,7:end]\-drdx[7:end, 1:6]))\drde[1:6, 1:6]
    @info "jacobian3 = $jacobian3"
    @info "jacobian4 = $(tovoigt(chabmat.variables_new.jacobian))"
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"

    chabmat.ddrivers = ddrivers
    integrate_material!(chabmat)
    g! = Materials.ChabocheModule.create_nonlinear_system_of_equations(chabmat)

    x0 = [tovoigt(chabmat.variables_new.stress); chabmat.variables_new.R; tovoigt(chabmat.variables_new.X1); tovoigt(chabmat.variables_new.X2)]
    drdx = ForwardDiff.jacobian(debang(g!), x0)
    @info "size(drdx) = $(size(drdx))"
    @info "drdx = $drdx"
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = parameters
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))
    jacobian = isotropic_elasticity_tensor(lambda, mu)
    drde = zeros((19,6))
    drde[1:6, 1:6] = -tovoigt(jacobian)
    @info "drde = $drde"
    @info "size(drde) = $(size(drde))"
    jacobian2 = (drdx\drde)[1:6, 1:6]
    @info "jacobian = $(tovoigt(jacobian))"
    @info "jacobian2 = $jacobian2"
    jacobian3 = (drdx[1:6, 1:6] +  drdx[1:6,7:end]*(drdx[7:end,7:end]\-drdx[7:end, 1:6]))\drde[1:6, 1:6]
    @info "jacobian3 = $jacobian3"
    @info "jacobian4 = $(tovoigt(chabmat.variables_new.jacobian))"
    update_material!(chabmat)
    @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)"
end
simple_integration_test_fd_tangent2()


text:
# # Synthetic models
#
# ## Packages

using CSV
using CairoMakie
using CalibrationErrors
using CalibrationErrorsDistributions
using CalibrationTests
using DataFrames
using Distributions
using FillArrays
using ProgressLogging
using Query
using Showoff
using StatsBase

using LinearAlgebra
using Printf
using Random

using CairoMakie.AbstractPlotting.ColorSchemes: Dark2_8

using Logging: with_logger
using TerminalLoggers: TerminalLogger

## set random seed
Random.seed!(1234)

## create path before saving
function wsavefig(file, fig=current_figure())
    mkpath(dirname(file))
    return save(file, fig)
end

## define progress logging frontend
const PROGRESSLOGGER = TerminalLogger()

## define non-intrusive plotting style
set_theme!(
    Theme(;
        Axis=(
            rightspinevisible=false,
            topspinevisible=false,
            xgridvisible=false,
            ygridvisible=false,
        ),
        Legend=(framevisible=false,),
    ),
)

# ## Synthetic models
#
# We study two setups with $d$-dimensional targets $Y$ and normal distributions $P_X$
# of the form $\mathcal{N}(c \mathbf{1}_d, 0.1^2 \mathbf{I}_d)$ as predictions,
# where $c \sim \mathrm{U}(0, 1)$.
# Since calibration analysis is only based on the targets and predicted distributions,
# we neglect features $X$ in these experiments and specify only the distributions of
# $Y$ and $P_X$.
#
# ### Calibrated setup
#
# In the first setup we simulate a calibrated model. We achieve this by sampling
# targets from the predicted distributions, i.e., by defining the conditional distribution
# of $Y$ given $P_X$ as
# ```math
# Y \,|\, P_X = \mathcal{N}(\mu, \Sigma) \sim \mathcal{N}(\mu, \Sigma).
# ```

function calibrated_model(dim::Int, nsamples::Int)
    ## sample predictions
    predictions = [MvNormal(Fill(rand(), dim), 0.1) for _ in 1:nsamples]

    ## sample targets
    targets = map(rand, predictions)

    return predictions, targets
end

# ### Uncalibrated setup
#
# In the second setup we simulate an uncalibrated model of the form
# ```math
# Y \,|\, P_X = \mathcal{N}(\mu, \Sigma) \sim \mathcal{N}([0.1, \mu_2, \ldots, \mu_d], \Sigma).
# ```

function uncalibrated_model(dim::Int, nsamples::Int)
    ## sample predictions
    predictions = [MvNormal(Fill(rand(), dim), 0.1) for _ in 1:nsamples]

    ## sample targets
    targets = map(rand, predictions)
    altdist = Normal(0.1, 0.1)
    for t in targets
        t[1] = rand(altdist)
    end

    return predictions, targets
end

# ## Convergence and computation time of estimators
#
# We perform an evaluation of the convergence and computation time of the biased estimator
# $\widehat{\mathrm{SKCE}}_k$, the unbiased estimator $\widehat{\mathrm{SKCE}}_{k,B}$ with
# blocks of size $B \in \{2, \sqrt{n}, n\}$. We use the tensor product kernel
# ```math
# \begin{aligned}
# k\big((p, y), (p', y')\big) &= \exp{\big(- W_2(p, p')\big)} \exp{\big(-(y - y')^2/2\big)} \\
# &= \exp{\big(-\sqrt{(m_p - m_{p'})^2 + (\sigma_p - \sigma_{p'})^2}\big)} \exp{\big( - (y - y')^2/2\big)},
# \end{aligned}
# ```
# where $W_2$ is the 2-Wasserstein distance and $m_p, m_{p'}$ and $\sigma_p, \sigma_{p'}$
# denote the mean and the standard deviation of the normal distributions $p$ and $p'$.

# ### Ground truth
#
# For both models, we have to "evaluate" the true calibration error. Generally, the error
# depends on the model (and hence also dimension $d$) and the kernel. If the model is
# calibrated, we know that the calibration error is zero. For the uncalibrated model, we
# estimate the ground truth with the minimum-variance unbiased estimator as the mean of
# SKCE estimates for 1000 randomly sampled datasets with 1000 data points.

true_SKCE(::typeof(calibrated_model), kernel; dim::Int) = 0.0
function true_SKCE(model::typeof(uncalibrated_model), kernel; dim::Int)
    estimator = UnbiasedSKCE(kernel)
    return mean(calibrationerror(estimator, model(dim, 1_000)...) for _ in 1:1_000)
end

# ### Benchmarking
#
# The following two functions implement the benchmarking. We sample 500 datasets of
# 4, 16, 64, 256, and 1024 data points each for the models of dimensions $d=1$ and $d=10$.
# For each of the datasets, we evaluate the different SKCE estimators. We compute the
# mean absolute error, the variance, and the minimum computation time for the estimates,
# grouped by the dimension of the model and the number of samples in the dataset.

function benchmark_estimator(estimator, model; dim::Int, nsamples::Int, groundtruth)
    ## compute the estimator (potentially depending on number of samples)
    _estimator = estimator(nsamples)

    ## cache for calibration error estimates
    estimates = Vector{Float64}(undef, 500)

    mintime = Inf

    name = @sprintf("benchmarking (dim = %2d, nsamples = %4d)", dim, nsamples)
    @progress name = name for i in eachindex(estimates)
        ## sample predictions and targets
        predictions, targets = model(dim, nsamples)

        ## define benchmark function
        benchmark_f =
            let estimator = _estimator, predictions = predictions, targets = targets
                () -> @timed calibrationerror(estimator, predictions, targets)
            end

        ## precompile function
        benchmark_f()

        ## compute calibration error and obtain elapsed time
        val, t = benchmark_f()

        ## only keep minimum execution time
        mintime = min(mintime, t)

        ## save error estimate
        estimates[i] = val
    end

    ## save the mean absolute deviation and the variance of the estimates
    meanerror = mean(abs(x - groundtruth) for x in estimates)
    variance = var(estimates)

    return (; dim, nsamples, meanerror, variance, mintime)
end

function benchmark_estimators(model)
    ## output file
    filename = joinpath("data", "synthetic", "errors_$(model).csv")

    ## check if results exist
    isfile(filename) && return DataFrame(CSV.File(filename))

    ## define kernel
    kernel = WassersteinExponentialKernel() ⊗ SqExponentialKernel()

    ## define estimators
    estimators = (
        "SKCE" => _ -> BiasedSKCE(kernel),
        "SKCE (B = 2)" => _ -> BlockUnbiasedSKCE(kernel, 2),
        "SKCE (B = √n)" => n -> BlockUnbiasedSKCE(kernel, max(2, Int(floor(sqrt(n))))),
        "SKCE (B = n)" => _ -> UnbiasedSKCE(kernel),
    )

    ## define number of samples
    nsamples = 2 .^ (2:2:10)

    ## ensure that output directory exists and open file for writing
    mkpath(dirname(filename))
    open(filename, "w") do file
        ## write headers
        println(file, "estimator,dim,nsamples,meanerror,variance,mintime")

        ## for dimensions ``d=1`` and ``d=10``
        for d in (1, 10)
            ## compute/estimate ground truth
            groundtruth = true_SKCE(model, kernel; dim=d)

            for (i, (name, estimator)) in enumerate(estimators)
                ## benchmark estimator
                @info "benchmarking estimator: $(name)"

                for n in nsamples
                    stats = benchmark_estimator(
                        estimator, model; dim=d, nsamples=n, groundtruth=groundtruth
                    )

                    ## save statistics
                    print(file, name, ",")
                    join(file, stats, ",")
                    println(file)
                end
            end
        end
    end

    ## load results
    return DataFrame(CSV.File(filename))
end

# We benchmark the estimators with the calibrated model.

Random.seed!(100)
with_logger(PROGRESSLOGGER) do
    benchmark_estimators(calibrated_model)
end

# We repeat the benchmark with the uncalibrated model.

Random.seed!(100)
with_logger(PROGRESSLOGGER) do
    benchmark_estimators(uncalibrated_model)
end

# ### Visualization
#
# We show a visualization of the results below.

function logtickformat(base::Int)
    function format(values)
        return map(Base.Fix2(logformat, base), showoff(values))
    end
    return format
end

function logformat(digits::String, base::Int)
    buf = IOBuffer()
    print(buf, base)
    for c in digits
        if '0' ≤ c ≤ '9'
            print(buf, Showoff.superscript_numerals[c - '0' + 1])
        elseif c == '-'
            print(buf, '⁻')
        elseif c == '.'
            print(buf, '·')
        end
    end
    return String(take!(buf))
end

function plot_benchmark_estimators(model; dim::Int)
    ## load and preprocess data
    filename = joinpath("data", "synthetic", "errors_$(model).csv")
    groups = @from i in DataFrame(CSV.File(filename)) begin
        @where i.dim == dim
        @orderby i.nsamples
        @select {
            i.estimator,
            log2_nsamples = log2(i.nsamples),
            log10_meanerror = log10(i.meanerror),
            log10_variance = log10(i.variance),
            log10_mintime = log10(i.mintime),
        }
        @collect DataFrame
    end

    ## create figure
    fig = Figure(; resolution=(960, 800))

    ## create axes to plot mean error and variance vs number of samples
    ax1 = Axis(
        fig[1, 1];
        xlabel="# samples",
        ylabel="mean error",
        xticks=2:2:10,
        xtickformat=logtickformat(2),
        ytickformat=logtickformat(10),
    )
    ax2 = Axis(
        fig[2, 1];
        xlabel="# samples",
        ylabel="variance",
        xticks=2:2:10,
        xtickformat=logtickformat(2),
        ytickformat=logtickformat(10),
    )

    ## create axes to plot mean error and variance vs timings
    ax3 = Axis(
        fig[1, 2];
        xlabel="time [s]",
        ylabel="mean error",
        xtickformat=logtickformat(10),
        ytickformat=logtickformat(10),
    )
    ax4 = Axis(
        fig[2, 2];
        xlabel="time [s]",
        ylabel="variance",
        xtickformat=logtickformat(10),
        ytickformat=logtickformat(10),
    )

    ## plot benchmark results
    estimators = ["SKCE", "SKCE (B = 2)", "SKCE (B = √n)", "SKCE (B = n)"]
    markers = ['●', '■', '▲', '◆']
    for (i, (estimator, marker)) in enumerate(zip(estimators, markers))
        group = filter(:estimator => ==(estimator), groups)
        color = Dark2_8[i]

        ## plot mean error vs samples
        scatterlines!(
            ax1,
            group.log2_nsamples,
            group.log10_meanerror;
            color=color,
            linewidth=2,
            marker=marker,
            markercolor=color,
        )

        ## plot variance vs samples
        scatterlines!(
            ax2,
            group.log2_nsamples,
            group.log10_variance;
            color=color,
            linewidth=2,
            marker=marker,
            markercolor=color,
        )

        ## plot mean error vs time
        scatterlines!(
            ax3,
            group.log10_mintime,
            group.log10_meanerror;
            color=color,
            linewidth=2,
            marker=marker,
            markercolor=color,
        )

        ## plot variance vs time
        scatterlines!(
            ax4,
            group.log10_mintime,
            group.log10_variance;
            color=color,
            linewidth=2,
            marker=marker,
            markercolor=color,
        )
    end

    ## link axes and hide decorations
    linkxaxes!(ax1, ax2)
    hidexdecorations!(ax1)
    linkxaxes!(ax3, ax4)
    hidexdecorations!(ax3)
    linkyaxes!(ax1, ax3)
    hideydecorations!(ax3)
    linkyaxes!(ax2, ax4)
    hideydecorations!(ax4)

    ## add legend
    elems = map(1:length(estimators)) do i
        [
            LineElement(; color=Dark2_8[i], linestyle=nothing, linewidth=2),
            MarkerElement(; color=Dark2_8[i], marker=markers[i], strokecolor=:black),
        ]
    end
    Legend(fig[end + 1, :], elems, estimators; orientation=:horizontal, tellheight=true)

    return fig
end

# We obtain the following plots:

plot_benchmark_estimators(calibrated_model; dim=1)
wsavefig("figures/synthetic/estimators_calibrated_model_dim=1.pdf") #jl
#!jl wsavefig("figures/synthetic/estimators_calibrated_model_dim=1.svg");
#!jl # ![](figures/synthetic/estimators_calibrated_model_dim=1.svg)

plot_benchmark_estimators(calibrated_model; dim=10)
wsavefig("figures/synthetic/estimators_calibrated_model_dim=10.pdf") #jl
#!jl wsavefig("figures/synthetic/estimators_calibrated_model_dim=10.svg");
#!jl # ![](figures/synthetic/estimators_calibrated_model_dim=10.svg)

plot_benchmark_estimators(uncalibrated_model; dim=1)
wsavefig("figures/synthetic/estimators_uncalibrated_model_dim=1.pdf") #jl
#!jl wsavefig("figures/synthetic/estimators_uncalibrated_model_dim=1.svg");
#!jl # ![](figures/synthetic/estimators_uncalibrated_model_dim=1.svg)

plot_benchmark_estimators(uncalibrated_model; dim=10)
wsavefig("figures/synthetic/estimators_uncalibrated_model_dim=10.pdf") #jl
#!jl wsavefig("figures/synthetic/estimators_uncalibrated_model_dim=10.svg");
#!jl # ![](figures/synthetic/estimators_uncalibrated_model_dim=10.svg)

# ## Test errors and computation time of calibration tests
#
# We fix the significance level $\alpha = 0.05$.
# Test predictions are sampled from the same distribution as $P_X$, and test targets are
# sampled independently from $\mathcal{N}(0, 0.1^2 \mathbf{I}_d)$.
#
# ### Benchmarking

iscalibrated(::typeof(calibrated_model)) = true
iscalibrated(::typeof(uncalibrated_model)) = false

function benchmark_test(test, model; dim::Int, nsamples::Int)
    ## number of simulations
    nrepeat = 500

    ## initial values
    ntesterrors = 0
    mintime = Inf

    name = @sprintf("benchmarking (dim = %2d, nsamples = %4d)", dim, nsamples)
    @progress name = name for _ in 1:nrepeat
        ## sample predictions and targets
        predictions, targets = model(dim, nsamples)

        ## define benchmark function
        benchmark_f = let test = test, predictions = predictions, targets = targets
            () -> @timed pvalue(test(predictions, targets))
        end

        ## precompile function
        benchmark_f()

        ## compute calibration error and obtain elapsed time
        val, t = benchmark_f()

        ## only keep minimum execution time
        mintime = min(mintime, t)

        ## update number of empirical test errors for
        ## significance level ``\alpha = 0.05``
        ntesterrors += iscalibrated(model) ⊻ (val ≥ 0.05)
    end

    ## compute empirical test error rate
    testerror = ntesterrors / nrepeat

    return (; dim, nsamples, testerror, mintime)
end

function benchmark_tests(model)
    ## output file
    filename = joinpath("data", "synthetic", "tests_$(model).csv")

    ## check if results exist
    isfile(filename) && return DataFrame(CSV.File(filename))

    ## define kernel
    kernel = WassersteinExponentialKernel() ⊗ SqExponentialKernel()

    ## define number of samples
    nsamples = 2 .^ (2:2:10)

    ## ensure that output directory exists and open file for writing
    mkpath(dirname(filename))
    open(filename, "w") do file
        ## write headers
        println(file, "test,dim,nsamples,testerror,mintime")

        ## for dimensions ``d=1`` and ``d=10``
        for d in (1, 10)
            ## define tests
            testpredictions = [MvNormal(rand(d), 0.1) for _ in 1:10]
            testtargets = [rand(MvNormal(d, 0.1)) for _ in 1:10]
            tests = (
                "SKCE (B = 2)" =>
                    (predictions, targets) -> AsymptoticBlockSKCETest(
                        BlockUnbiasedSKCE(kernel, 2), predictions, targets
                    ),
                "SKCE (B = √n)" =>
                    (predictions, targets) -> AsymptoticBlockSKCETest(
                        BlockUnbiasedSKCE(kernel, Int(floor(sqrt(length(predictions))))),
                        predictions,
                        targets,
                    ),
                "SKCE (B = n)" =>
                    (predictions, targets) ->
                        AsymptoticSKCETest(kernel, predictions, targets),
                "CME" =>
                    (predictions, targets) -> AsymptoticCMETest(
                        UCME(kernel, testpredictions, testtargets), predictions, targets
                    ),
            )

            for (i, (name, test)) in enumerate(tests)
                ## benchmark estimator
                @info "benchmarking test: $(name)"

                for n in nsamples
                    stats = benchmark_test(test, model; dim=d, nsamples=n)

                    ## save statistics
                    print(file, name, ",")
                    join(file, stats, ",")
                    println(file)
                end
            end
        end
    end

    ## load results
    return DataFrame(CSV.File(filename))
end

# First we benchmark the calibrated model.

Random.seed!(100)
with_logger(PROGRESSLOGGER) do
    benchmark_tests(calibrated_model)
end

# We repeat the analysis with the uncalibrated model.

Random.seed!(100)
with_logger(PROGRESSLOGGER) do
    benchmark_tests(uncalibrated_model)
end

# ### Visualization
#
# Again we visualize the results of our benchmarks. However, this time we
# compare the results for the calibrated and the uncalibrated model in the
# same plot.

function plot_benchmark_tests(; dim::Int)
    ## load and preprocess data
    df = mapreduce(vcat, (calibrated_model, uncalibrated_model)) do model
        filename = joinpath("data", "synthetic", "tests_$(model).csv")
        df = DataFrame(CSV.File(filename))
        df[!, :model] .= string(model)
        return df
    end
    groups = @from i in df begin
        @where i.dim == dim
        @orderby i.nsamples
        @select {
            i.test,
            i.model,
            log2_nsamples = log2(i.nsamples),
            i.testerror,
            log10_mintime = log10(i.mintime),
        }
        @collect DataFrame
    end

    ## create figure
    fig = Figure(; resolution=(960, 400))

    ## add labels
    Label(fig[1:2, 1], "empirical test error"; rotation=π / 2, tellheight=false)
    Label(fig[1, 2:3, Top()], "calibrated model"; padding=(0, 0, 10, 0))
    Label(fig[2, 2:3, Top()], "uncalibrated model"; padding=(0, 0, 10, 0))

    ## create axes to plot test error vs number of samples
    ax1 = Axis(
        fig[1, 2];
        ylabel="type I error",
        xticks=2:2:10,
        xtickformat=logtickformat(2),
        xticklabelsize=12,
        yticklabelsize=12,
    )
    ax2 = Axis(
        fig[2, 2];
        xlabel="# samples",
        ylabel="type II error",
        xticks=2:2:10,
        xtickformat=logtickformat(2),
        xticklabelsize=12,
        yticklabelsize=12,
    )

    ## create axes to plot test error vs timings
    ax3 = Axis(
        fig[1, 3]; xtickformat=logtickformat(10), xticklabelsize=12, yticklabelsize=12
    )
    ax4 = Axis(
        fig[2, 3];
        xlabel="time [s]",
        xtickformat=logtickformat(10),
        xticklabelsize=12,
        yticklabelsize=12,
    )

    ## plot benchmark results
    tests = ["SKCE (B = 2)", "SKCE (B = √n)", "SKCE (B = n)", "CME"]
    markers = ['●', '■', '▲', '◆']
    for (i, (test, marker)) in enumerate(zip(tests, markers))
        color = Dark2_8[i]

        ## for both calibrated and uncalibrated model
        for (axes, model) in
            zip(((ax1, ax3), (ax2, ax4)), (calibrated_model, uncalibrated_model))
            group = filter(x -> x.test == test && x.model == string(model), groups)

            ## plot test error vs samples
            scatterlines!(
                axes[1],
                group.log2_nsamples,
                group.testerror;
                color=color,
                linewidth=2,
                marker=marker,
                markercolor=color,
            )

            ## plot test error vs timings
            scatterlines!(
                axes[2],
                group.log10_mintime,
                group.testerror;
                color=color,
                linewidth=2,
                marker=marker,
                markercolor=color,
            )
        end
    end

    ## plot horizontal lines for significance level
    for axis in (ax1, ax3)
        hlines!(axis, 0.05; color=:black, linestyle=:dash, linewidth=2)
    end

    ## link axes and hide decorations
    linkxaxes!(ax1, ax2)
    hidexdecorations!(ax1)
    linkxaxes!(ax3, ax4)
    hidexdecorations!(ax3)
    linkyaxes!(ax1, ax3)
    hideydecorations!(ax3)
    linkyaxes!(ax2, ax4)
    hideydecorations!(ax4)

    ## add legend
    elems = map(1:length(tests)) do i
        [
            LineElement(; color=Dark2_8[i], linestyle=nothing, linewidth=2),
            MarkerElement(; color=Dark2_8[i], marker=markers[i], strokecolor=:black),
        ]
    end
    push!(elems, [LineElement(; color=:black, linestyle=:dash, linewidth=2)])
    Legend(
        fig[1:2, end + 1],
        elems,
        vcat(tests, "significance level");
        tellwidth=true,
        gridshalign=:left,
    )

    return fig
end

plot_benchmark_tests(; dim=1)
wsavefig("figures/synthetic/tests_dim=1.pdf") #jl
#!jl wsavefig("figures/synthetic/tests_dim=1.svg");
#!jl # ![](figures/synthetic/tests_dim=1.svg)

plot_benchmark_tests(; dim=10)
wsavefig("figures/synthetic/tests_dim=10.pdf") #jl
#!jl wsavefig("figures/synthetic/tests_dim=10.svg");
#!jl # ![](figures/synthetic/tests_dim=10.svg)


text:
###############################################################################
#
# Solar angle type
#
###############################################################################
"""
    struct SolarAngles{FT}

Struct for observation and solar angles

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SolarAngles{FT}
    "Hill zenith angle `[°]`, hill slope angle"
    hza::FT = 0
    "Hill azimuth angle `[°]`, 0 for North, 180 for south"
    haa::FT = 180
    "Solar azimuth angle `[°]`, a function of time"
    saa::FT = 0
    "Solar zenith angle `[°]`, a function of lat and time"
    sza::FT = 30
    "Viewing azimuth angle `[°]`"
    vaa::FT = 0
    "Viewing zenith angle `[°]`"
    vza::FT = 0
    "Relative azimuth angle `[°]`, difference between saa and vaa"
    raa::FT = 0
end

max_stars_repo_name: climate-machine/Land


text:
# Wrap ‘llvm.prefetch‘ Intrinsic
# https://llvm.org/docs/LangRef.html#llvm-prefetch-intrinsic

# address is the address to be prefetched, rw is the specifier determining if the fetch
# should be for a read (0) or write (1), and locality is a temporal locality specifier
# ranging from (0) - no locality, to (3) - extremely local keep in cache. The cache type
# specifies whether the prefetch is performed on the data (1) or instruction (0) cache.
# The rw, locality and cache type arguments must be constant integers.
@inline function _prefetch(address::Ptr{T}, rw::Integer, locality::Integer, cachetype::Integer) where T
    Base.llvmcall((""" declare void @llvm.prefetch(i8*, i32, i32, i32) """,
                   """
                   %ptr = inttoptr i64 %0 to i8*
                   call void @llvm.prefetch(i8* %ptr, i32 %1, i32 %2, i32 %3)
                   ret void
                   """),
                  Nothing, Tuple{UInt64, Int32, Int32, Int32},
                  UInt64(address), Int32(rw), Int32(locality), Int32(cachetype))
end

@inline _prefetch_r(ptr::Ptr{T}) where T = _prefetch(ptr, 0, 3, 1)
@inline _prefetch_w(ptr::Ptr{T}) where T = _prefetch(ptr, 1, 3, 1)

@inline function prefetch_r(::Type{Val{M}}, ::Type{Val{N}}, ::Type{Val{Rem}}, ::Type{Val{Shift}}, ptr::Ptr{T}, ld) where {M,N,Rem,Shift,T}
    for n in 0:N-1
        for m in 0:(M÷64 + (M % 64 >= Rem) - 1)
            _prefetch_r(ptr + M * 64 + Shift + ld * N)
        end
    end
end


text:
println("===============================================")
println("Testing maxwell.jl")
println("===============================================")

function Fract_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = relaxmod(Fract_Maxwell, t, [1.0, 1.0, 1.0, 0.0])
    maxwell = relaxmod(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(t))
end
@test Fract_Maxwell_G_reduce()

function FractS_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = relaxmod(FractS_Maxwell, t, [1.0, 1.0, 1.0])
    maxwell = relaxmod(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(t))
end
@test FractS_Maxwell_G_reduce()

function FractD_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = relaxmod(FractD_Maxwell, t, [1.0, 1.0, 0.0])
    maxwell = relaxmod(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(t))
end
@test FractD_Maxwell_G_reduce()

function Fract_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = creepcomp(Fract_Maxwell, t, [1.0, 1.0, 1.0, 0.0])
    maxwell = creepcomp(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(t))
end
@test Fract_Maxwell_J_reduce()

function FractS_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = creepcomp(FractS_Maxwell, t, [1.0, 1.0, 1.0])
    maxwell = creepcomp(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(t))
end
@test FractS_Maxwell_J_reduce()

function FractD_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = creepcomp(FractD_Maxwell, t, [1.0, 1.0, 0.0])
    maxwell = creepcomp(Maxwell, t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(t))
end
@test FractD_Maxwell_J_reduce()

function Fract_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = storagemod(Fract_Maxwell, chirp, [1.0, 1.0, 1.0, 0.0])
    maxwell = storagemod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(chirp))
end
@test Fract_Maxwell_Gp_reduce()

function FractS_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = storagemod(FractS_Maxwell, chirp, [1.0, 1.0, 1.0])
    maxwell = storagemod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(chirp))
end
@test FractS_Maxwell_Gp_reduce()

function FractD_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = storagemod(FractD_Maxwell, chirp, [1.0, 1.0, 0.0])
    maxwell = storagemod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(chirp))
end
@test FractD_Maxwell_Gp_reduce()

function Fract_Maxwell_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = lossmod(Fract_Maxwell, chirp, [1.0, 1.0, 1.0, 0.0])
    maxwell = lossmod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(chirp))
end
@test Fract_Maxwell_Gpp_reduce()

function FractS_Maxwell_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = lossmod(FractS_Maxwell, chirp, [1.0, 1.0, 1.0])
    maxwell = lossmod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(chirp))
end
@test FractS_Maxwell_Gpp_reduce()

function FractD_Maxwell_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = lossmod( FractD_Maxwell, chirp, [1.0, 1.0, 0.0])
    maxwell = lossmod(Maxwell, chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(chirp))
end
@test FractD_Maxwell_Gpp_reduce()


text:
module ManySum
using DataStructures
export pairwiseSum!
export Manysum

function pairwiseSum!(a::Array{T})::T where T<:AbstractFloat
  i=1
  n=length(a)
  while i<n
    for j=1:2*i:n-i
      a[j]+=a[j+i]
    end
    i*=2
  end
  if n>0
    a[1]
  else
    0
  end
end

struct Manysum{T <: AbstractFloat}
  bucket ::SortedDict{Int,T}
  function Manysum()::T
    msum=new{T}()
    bucket=SortedDict{Int,T}()
    msum
  end
end

end # module


text:
using Printf
using Plots
using Statistics
using LinearAlgebra
using MDBM
using PyPlot
pygui(true);
using Interpolations

#interpolate
function it(A)
    matrdim=size(A,2)-1
    scaleitp=real(A[:,1])
    knots=(scaleitp,)
    ARe=real(A)
    AIm=imag(A)

    imRe=[interpolate(knots, ARe[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imRetemp=[interpolate(knots, real(ARe[:,j]), Gridded(Linear()))]
        imRe=vcat(imRe,imRetemp)
    end
    imIm=[interpolate(knots, AIm[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imImtemp=[interpolate(knots, real(AIm[:,j]), Gridded(Linear()))]
        imIm=vcat(imIm,imImtemp)
    end
    return(hcat(imRe,imIm))
end

B=transpose([0 0.5 1 1.5 2; 3*im 2+4*im -5+1*im 6*im 5*im-4; 7*im -8+3*im -5+2*im 7*im 3*im-2])
real(B)
typeof(real(B))

iterB=it(B)

interpolate((real(B[:,1]),), real(B[:,2]),Gridded(Linear()))
#substitution in interpolating function
function sub(it,t)
    subdim=size(it,1)
    out=zeros(ComplexF64,subdim)
    for j = 1:subdim
        out[j]=it[j,1](t)+it[j,2](t)*im
    end
    return(out)
end

sub(iterB,1.25)


text:
abstract type InitialConditionQuantity end
struct DevicePower <: InitialConditionQuantity end
struct DeviceStatus <: InitialConditionQuantity end
struct TimeDurationON <: InitialConditionQuantity end
struct TimeDurationOFF <: InitialConditionQuantity end
struct DeviceEnergy <: InitialConditionQuantity end

mutable struct InitialCondition{T<:Union{PJ.ParameterRef, Float64}}
    device::PSY.Device
    update_ref::UpdateRef
    value::T
    cache::Union{Nothing, Type{<:AbstractCache}}
end

function InitialCondition(device::PSY.Device, access_ref::UpdateRef, value::T) where {T<:Union{PJ.ParameterRef, Float64}}
    return InitialCondition(device, access_ref, value, nothing)
end

struct ICKey{IC<:InitialConditionQuantity, D<:PSY.Device}
    quantity::Type{IC}
    device_type::Type{D}
end

const DICKDA = Dict{ICKey, Array{InitialCondition}}

function value(p::InitialCondition{Float64})
    return p.value
end

function value(p::InitialCondition{PJ.ParameterRef})
    return PJ.value(p.value)
end

get_condition(ic::InitialCondition) = ic.value

device_name(ini_cond::InitialCondition) = PSY.get_name(ini_cond.device)


text:
using ObjCWrap
using Test

@classes NSString, NSNumber, NSMutableDictionary, NSArray, NSMutableArray

@testset "All Tests" begin
    include("messaging.jl")
    include("syntax.jl")
    include("collections.jl")
end


text:
"""
    QuantumGlue

Interface between quantum science and technology packages.

Yao: quantum information/computing
QuantumOptics: ...

Spin up in `QuantumOptics.SpinBasis(1//2)` corresponds to `0` in the default basis of `Yao`.

`QuantumOptics` is aliased to `QO`.
"""
module QuantumGlue
import Yao
import QuantumOptics
const QO = QuantumOptics

export Ket, register

"""
    qubitbasis(n::Integer)

Create an `n`-qubit basis of type QO.CompositeBasis. The basis is the tensor
product of built-in `Spin(1//2)` bases.
"""
function qubitbasis(n::Integer)
    return QO.CompositeBasis((QO.SpinBasis(1//2) for i in 1:n)...)
end

"""
    Ket(reg::Yao.AbstractRegister)

Convert `reg` to a instance of type `QuantumOptics.Ket`. The input is not checked.
This method is an inverse of `register` acting on some objects of type `QuantumOptics.Ket`.
"""
function QuantumOptics.Ket(reg::Yao.AbstractRegister)
    return QuantumOptics.Ket(reg |> Yao.nqubits |> qubitbasis, reg |> Yao.state |> Yao.Intrinsics.matvec)
end

"""
    register(s::QO.Ket)

Convert `s` to an instance of `Yao.DefaultRegister`. The data is not copied, but
is shared between the input and output. The input is not checked. `s` must be a
`QuantumOptics.CompositeBasis` of `SpinBasis(1//2)`.  This method is an inverse
of `spins`.
"""
function Yao.register(s::QO.Ket)
    return Yao.register(s.data)
end

end # module QuantumGlue