text:
# TODO: parse expressions with MacroTools

"""
    @cfunc ex

Define a "calculation" function with an "adapter" methods used to parse args from model.
"calculate" functions take arguments from system state
and calculate "rate"s determining the probability weight of the reaction be selected.

# Example

For function definition:
```julia
@cfunc @inline growth_c(r, x::Vector) = r * x # function to calculate "growth rate"
```
This macro creates two methods, an "adapter" method
```julia
growth_c(args::NamedTuple) = growth_c(args.r, args.x)
```
and the origin method
```julia
@inline growth_c(r, x::Vector) = r * x
```

!!! info
    For function definition with other macros, put those macros after this macro.

!!! warning
    The argument name `t` is reserved for time.
    Don't use those variable names for other usage
"""
macro cfunc(ex)
    fname, fargs = _parsefunc(ex, (:t,))
    adapter = :(Base.@propagate_inbounds $fname(t, args::NamedTuple) = $fname($(fargs...)))
    return esc(Expr(:block, adapter, ex))
end

"""
    @ufunc ex

Define a "update" function with an "adapter" methods used to parse args from model.
"update" functions take arguments from system state and update system state.

# Example

For function definition:
```julia
@ufunc Base.@propagate_inbounds growth_u!(ind::CartesianIndex{1}, x) = x[ind] += 1
```
This macro creates two methods, an "adapter" method
```julia
growth_u!(ind, args::NamedTuple) = growth_u(ind, args.x)
```
and the origin method
```julia
Base.@propagate_inbounds growth_u!(ind::CartesianIndex{1}, x) = x[ind] += 1
```

!!! info
    For function definition with other macros, put those macros after this macro.

!!! warning
    The argument name `t` is reserved for time,
    arguments name `rng` is reserved for random number generator,
    and the argument name `ind` is preserved for index of "reaction".
    Don't use those variable names for other usage
"""
macro ufunc(ex)
    fname, fargs = _parsefunc(ex, (:t, :ind, :rng))
    adapter = :(Base.@propagate_inbounds $fname(rng, t, ind, args::NamedTuple) =
        $fname($(fargs...)))
    return esc(Expr(:block, adapter, ex))
end

"""
    @reaction name ex

Define a [`Reaction`](@ref) with the given `name` and `ex`.

```julia
@reaction growth begin
    r * x # "calculate" expression
    begin
        i = ind[1]
        x[i] += 1 # "update" expression
    end
end
```
will create a `Reaction` with a "calculation" function:
```julia
@cfunc Base.@propagate_inbounds growth_c(r, x) = r * x
```
and an "update" function:
```julia
@ufunc Base.@propagate_inbounds growth_u!(ind, x) = begin
    i = ind[1]
    x[i] += 1
end
```
where arguments of functions were collected from given expression automatically.

!!! note
    If there are global variables, this macro may can not collect arguments correctly
    Especially, for functions which accept function and types as arguments,
    those functions and types may also be collect as an arguments.
    Thus, these variables must be marked as global variables by `global` before use them,
    even functions and types.
    Besides, type annotation expressions like `x::Vector{Int}`,
    types `Vector` and `Int` will not be collected.
    Avoid to defined your reaction with a type arguments for type annotation.

!!! warning
    The expression follow the same name preserve rule as `@cfunc` and `@ufunc`,
    don't use those variable names for other usage.

This macro cannot parse macro expression.
Thus in some cases, arguments of expression with macro may not be collected correctly.
To avoid this, define reaction with anonymous functions may helpfully:
```
@reaction reaction begin
    (A, B) -> @einsum C[i, k] := A[i, j] * B[j, k] # where C, i, j, k should not be collected
    begin
        # do something
    end
end
```
in this case, only arguments of anonymous function (`c`, `A` and `B`) will be collected.
"""
macro reaction(name::Symbol, block::Expr)
    if block.head == :block && length(block.args) == 4
        cname = Symbol("$(name)_c")
        uname = Symbol("$(name)_u!")
        cbody = block.args[2]
        ubody = block.args[4]
        cargs = collectargs(cbody)
        uargs = collectargs(ubody)
        return esc(
            Expr(
                :block,
                :(Base.@propagate_inbounds $cname(t, args::NamedTuple) =
                    $cname($((_warparg(arg, (:t,)) for arg in cargs)...))),
                :(Base.@propagate_inbounds $cname($(cargs...)) = $cbody),
                :(Base.@propagate_inbounds $uname(rng, t, ind, args::NamedTuple) =
                    $uname($((_warparg(arg, (:t, :ind, :rng)) for arg in uargs)...))),
                :(Base.@propagate_inbounds $uname($(uargs...)) = $ubody),
                Expr(:(=), name, :(Reaction($cname, $uname))),
            ),
        )
    else
        throw(ArgumentError("can't parse expression"))
    end
end

# extract name and arguments of given function expression
function _parsefunc(ex::Expr, preserve)
    head = ex.head
    if head in (:function, :(=), :where) # unpack function expression
        return _parsefunc(ex.args[1], preserve)
    elseif head == :macrocall # unpack macro expression
        return _parsefunc(ex.args[end], preserve)
    elseif head == :call # the args of :call expression is the function name and the arguments
        fname = ex.args[1]
        fargs = map(arg -> _parsearg(arg, preserve), ex.args[2:end])
        return fname, fargs
    else
        error("Unknown expr")
    end
end

# unpack argument names
_parsearg(arg::Symbol, preserve) = _warparg(arg, preserve)
_parsearg(arg::Expr, preserve) = _warparg(arg.args[1]::Symbol, preserve)
function _warparg(sym::Symbol, preserve)
    if sym in preserve
        return sym
    else
        return Expr(:., :args, QuoteNode(sym))
    end
end

# collect arguments from given expression
function collectargs(ex)
    if ex.head == :(->)
        return unbox_args!(Set{Symbol}(), ex.args[1])
    else
        return collectargs!(Set{Symbol}(), Set{Symbol}(), ex)
    end
end

function unbox_args!(args, ex::Expr)
    if ex.head == :(::)
        push!(args, ex.args[1])
    elseif ex.head == :tuple
        for arg in ex.args
            unbox_args!(args, arg)
        end
    else
        error("can't parse expression")
    end
    return args
end
unbox_args!(args, s::Symbol) = (push!(args, s), args)

# `args` is a list of symbols will be treat as arguments name
# `vars` is a list of symbols which are names of internal variable or global variables
collectargs!(args, _, _) = args # ignore when third argument is neither a symbol nor a expr
function collectargs!(args, vars, sym::Symbol) # symbols neither in syms or iv will be collected
    if !in(sym, args) && !in(sym, vars) && sym != :(:)
        push!(args, sym)
    end
    return args
end
function collectargs!(args, vars, _ex::Expr)
    ex = shortdef(_ex) # force convert function definition to short form
    head = ex.head
    if head == :(=) # assignment or function definition
        arg1 = ex.args[1]
        if (arg1 isa Symbol || (arg1 isa Expr && arg1.head == :tuple)) # assignment to one or multi new symbol
            collectargs!(args, vars, ex.args[2]) # must collect exprs right of equal firstly #2
            collectargs!(vars, args, ex.args[1]) # collect internal variables
            return args
        elseif arg1.head == :call # local function definition as short form
            push!(vars, arg1.args[1]) # local function name will be internal variable
            func_args = Set{Symbol}() # local function arguments
            collectargs!(func_args, Set{Symbol}(), arg1) # collect local function arguments
            collectargs!(args, union(vars, func_args), ex.args[2]) # collect variables used in local function
            return args
        else
            collectargs!(args, vars, ex.args[1]) # collect left side of equal
            collectargs!(args, vars, ex.args[2]) # collect right side of equal
            return args
        end
    elseif head == :-> # anonymous function
        func_args = Set{Symbol}() # anonymous function arguments
        collectargs!(func_args, Set{Symbol}(), ex.args[1]) # collect function arguments
        collectargs!(args, union(vars, func_args), ex.args[2]) # collect variables used in anonymous function
        return args
    elseif head == :. # for x.y only collect x, y will be ignore
        return collectargs!(args, vars, ex.args[1])
    elseif head == :global # global variables
        return union!(vars, ex.args)
    elseif head == :(::) # type annotation
        return collectargs!(args, vars, ex.args[1]) # collect left side of :: but ignore right side
    elseif head == :ref # ignore begin and end in reference
        collectargs!(args, vars, ex.args[1])
        foreach(
            arg -> collectargs!(args, TempUnion(vars, (:begin, :end)), arg),
            ex.args[2:end],
        )
    elseif head in (:call, :macrocall) # function and macro will not be treat as a args
        foreach(arg -> collectargs!(args, vars, arg), ex.args[2:end])
    else # in other cases, collect all arguments
        foreach(arg -> collectargs!(args, vars, arg), ex.args)
    end
    return args
end

# union two collections
# but push! to the first collection
struct TempUnion{S,T}
    s::S
    t::T
end
Base.in(item, s::TempUnion) = (item in s.s) || (item in s.t)
Base.push!(s::TempUnion, item) = push!(s.s, item)


text:
function AdamsMoultonFieldIntegrator(arg0::Field, arg1::jint, arg2::jdouble, arg3::jdouble, arg4::Vector{jdouble}, arg5::Vector{jdouble})
    return AdamsMoultonFieldIntegrator((Field, jint, jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5)
end

function AdamsMoultonFieldIntegrator(arg0::Field, arg1::jint, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return AdamsMoultonFieldIntegrator((Field, jint, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function integrate(obj::AdamsMoultonFieldIntegrator, arg0::FieldExpandableODE, arg1::FieldODEState, arg2::RealFieldElement)
    return jcall(obj, "integrate", FieldODEStateAndDerivative, (FieldExpandableODE, FieldODEState, RealFieldElement), arg0, arg1, arg2)
end

ext: jl


text:
using BinDeps

@BinDeps.setup
libnames = ["libCLBlast", "libclblast", "clblast"]
libCLBlast = library_dependency("libCLBlast", aliases = libnames)
version = "1.5.1"

if Sys.iswindows()
    if Sys.ARCH == :x86_64
        uri = URI("https://github.com/CNugteren/CLBlast/releases/download/" *
                  version * "/CLBlast-" * version * "-Windows-x64.zip")
        basedir = @__DIR__
        provides(
            Binaries, uri,
            libCLBlast, unpacked_dir = ".",
            installed_libpath = joinpath(basedir, "libCLBlast", "lib"), os = :Windows
        )
    else
        error("Only 64 bit windows supported with automatic build.")
    end
end

if Sys.islinux()
    #=if Sys.ARCH == :x86_64
        name, ext = splitext(splitext(basename(baseurl * "Linux-x64.tar.gz"))[1])
        uri = URI(baseurl * "Linux-x64.tar.gz")
        basedir = joinpath(@__DIR__, name)
        provides(
            Binaries, uri,
            libCLBlast, unpacked_dir = basedir,
            installed_libpath = joinpath(basedir, "lib"), os = :Linux
        )
    else
        error("Only 64 bit linux supported with automatic build.")
    end=#
    provides(Sources, URI("https://github.com/CNugteren/CLBlast/archive/" * version * ".tar.gz"),
             libCLBlast, unpacked_dir="CLBlast-" * version)

    builddir = joinpath(@__DIR__, "src", "CLBlast-" * version, "build")
    libpath = joinpath(builddir, "libclblast.so")
    provides(BuildProcess,
        (@build_steps begin
            GetSources(libCLBlast)
            CreateDirectory(builddir)
            FileRule(libpath, @build_steps begin
                ChangeDirectory(builddir)
                `cmake -DSAMPLES=OFF -DTESTS=OFF -DTUNERS=OFF -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ..`
                `make -j 4`
            end)
        end),
        libCLBlast, installed_libpath=libpath, os=:Linux)
end

if Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "homebrew/core/clblast", libCLBlast, os = :Darwin)
end

if Sys.islinux()
    # BinDeps.jl seems to be broken, cf. https://github.com/JuliaLang/BinDeps.jl/issues/172
    wd = pwd()
    sourcedir = joinpath(@__DIR__, "CLBlast-" * version)
    builddir = joinpath(sourcedir, "build")
    libpath = joinpath(builddir, "libclblast.so")
    if !isdir(sourcedir)
        url = "https://github.com/CNugteren/CLBlast/archive/" * version * ".tar.gz"
        download_finished = false
        try
            run(pipeline(`wget -q -O - $url`, `tar xzf -`))
            global download_finished = true
        catch e
            println(e)
        end
        if download_finished == false
            try
                run(pipeline(`curl -L $url`, `tar xzf -`))
                global download_finished = true
            catch e
                println(e)
            end
        end
    end
    isdir(builddir) || mkdir(builddir)
    cd(builddir)
    run(`cmake -DSAMPLES=OFF -DTESTS=OFF -DTUNERS=OFF -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ..`)
    run(`make -j 4`)
    cd(wd)
    open(joinpath(@__DIR__, "deps.jl"), "w") do file
        write(file,
"""
using Libdl

# Macro to load a library
macro checked_lib(libname, path)
    Libdl.dlopen(path)
    quote
        const \$(esc(libname)) = \$path
    end
end

@checked_lib libCLBlast "$libpath"
""")
    end
else
    @BinDeps.install Dict(:libCLBlast => :libCLBlast)
end


text:
"""
    KPath{D}

Stores an path in k-space.

#### Fields
- kpts::Vector{SVector{D, Float64}}
- arc::Vector{Float64}
- special::Vector{Float64}
- nkpts::Int

"""
struct KPath{D}
    kpts::Vector{SVector{D, Float64}}
    arc::Vector{Float64}
    special::Vector{Float64}
    nkpts::Int
end

"""
    KMesh{D}

Stores an path in k-space.

#### Fields
- kpts::Vector{SVector{D, Float64}}
- ws::Vector{Float64}
- nkpts::Int
"""
struct KMesh{D}
    kpts::Vector{SVector{D, Float64}}
    ws::Vector{Float64}
    nkpts::Int
end


text:
# wrappers of low-level functionality

function curandGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(curandGetProperty(CUDAapi.MAJOR_VERSION),
                          curandGetProperty(CUDAapi.MINOR_VERSION),
                          curandGetProperty(CUDAapi.PATCH_LEVEL))

macro allocates(ex)
  quote
    CuArrays.extalloc(check=err->isa(err, CURANDError) &&
                                 err.code == CURAND_STATUS_ALLOCATION_FAILED) do
      $(esc(ex))
    end
  end
end


text:
module RoughSurface

# Write your package code here.
using DelimitedFiles, JLD2, FileIO
using PyCall

export write_surf, read_surf, show_surf

```@meta
CurrentModule = RoughSurface
```

include(joinpath(@__DIR__, "RoughSurface_GUI.jl"))
include(joinpath(@__DIR__, "RoughSurface_Controller.jl"))

"""
    write_surf()

将表面数据写入文件, 默认txt文件的数据用制表符分割, csv文件的数据用``,''分割,
jld2则为直接保存二进制文件
"""
function write_surf(file_path, surf)
    # 利用splitext来得到文件类型
    file_ext = splitext(file_path)[2]
    if file_ext == ".txt"
        open(file_path, "w") do io
            writedlm(io, surf)
        end # do
    elseif file_ext == ".csv"
        open(file_path, "w") do io
            writedlm(io, surf, ',')
        end # do
    elseif file_ext == ".jld2"
        save(file_path, "surf", surf)
    else
        @warn("未知的文件类型, 保存为jld2格式")
        save(file_path*".jld2", "surf", surf)
    end

end # function

"""
    read_surf()

从文件读取表面数据, 与写入数据保持一致, 默认txt文件的数据用制表符分割, csv文件的数据用``,''分割,
jld2则为直接保存二进制文件.
"""
function read_surf(file_path)
    file_ext = splitext(file_path)[2]
    if file_ext == ".txt"
        surf = readdlm(file_path)
    elseif file_ext == ".csv"
        surf = readdlm(file_path, ',')
    elseif file_ext == ".jld2"
        surf = load(file_path,"surf")
    else
        @warn("未知的文件类型, 尝试jld2格式")
        surf = load(file_path*".jld2","surf")
    end
    return surf
end # function

"""
    show_surf()

利用python的pyvista包显示规则表面数据.
"""
show_surf(surf_mat) = py"show_surf"(surf_mat)

end  # module


text:
@symbol_func function P_ohmic(cur_reactor::AbstractReactor)
  cur_P = cur_reactor.I_P ^ 2

  cur_P *= ( 1 - f_CD(cur_reactor) )

  cur_P *= R_P(cur_reactor)

  cur_P *= 1e6

  cur_P
end


text:
#=
HackerRank problem - 10 Days of Statistics: Day 1: Quartiles
https://www.hackerrank.com/challenges/s10-quartiles/problem

Calculate quartiles from user inputs.

Date: 2020-06-26

Issues:
  - The problem does not specify the screen output formatting,
e.g. 4 or 4.0. Possible answers are generally rational number, but the tests
seem to assume integer formatting output such as '4'. The test rejects
outputs in decimal numbers.
=#

using Printf
n = parse(Int, readline(stdin))
arr = sort(parse.(Int, split(readline(stdin))))

"""
  median(arr)

Calculate a median value in the way given in the problem.

# Arguments
- `arr`: An array of numbers
"""
function median(arr)
    n = size(arr)[1]
    n % 2 == 0 ? (arr[n ÷ 2] + arr[n ÷ 2 + 1]) / 2.0 : arr[n ÷ 2 + 1]
end

println("$(@sprintf("%d", median(arr[1:(n ÷ 2)])))")
println("$(@sprintf("%d", (median(arr))))")
println("$(@sprintf("%d", median(arr[(n % 2 == 0 ? n ÷ 2 + 1 : n ÷ 2 + 2):end])))")


text:
get_module(method_instance::Core.MethodInstance) = method_instance.def.module
get_module(_) = Main

# Not sure if this is the best way to do this.
macro caller_module(i)
    quote
        st = stacktrace()
        get_module(st[$i].linfo)
    end
end


text:
# DESCRIPTION
# Rigid-body equations of motion.
#
# INPUT
# - xf: fuselage states
# - Ftot: net forces acting at the CG
# - Mtot: net moments acting on the CG
#
# OUTPUT
# - xfdot: rigid-body dynamics
#
# ------------------------------------------------------------------------------

function eqnmot(xf,Ftot,Mtot)
    # unpack forces and moments
    X=Ftot[1]
    Y=Ftot[2]
    Z=Ftot[3]
    L=Mtot[1]
    M=Mtot[2]
    N=Mtot[3]
    # unpack states
    u=xf[1]
    v=xf[2]
    w=xf[3]
    p=xf[4]
    q=xf[5]
    r=xf[6]
    phi=xf[7]
    theta=xf[8]
    psi=xf[9]
    # trig functions
    cphi=cos(phi)
    sphi=sin(phi)
    cthe=cos(theta)
    sthe=sin(theta)
    cpsi=cos(psi)
    spsi=sin(psi)
    # state derivative vector
    xf_dot=zeros(12)
    # equations of motion
    # velocities
    xfdot=zeros(NFSTATES)
    xfdot[1]=X/MASS-G*sthe-q*w+r*v
    xfdot[2]=Y/MASS+G*cthe*sphi-r*u+p*w
    xfdot[3]=Z/MASS+G*cthe*cphi-p*v+q*u
    # angular rates
    gam=IX*IZ-IXZ^2;
    xfdot[4]=(IZ*L+IXZ*N+IXZ*(IX-IY+IZ)*p*q-(IZ^2-IY*IZ+IXZ^2)*q*r)/gam
    xfdot[5]=(M+(IZ-IX)*p*r-IXZ*(p^2-r^2))/IY
    xfdot[6]=(IX*N+IXZ*L-IXZ*(IX-IY+IZ)*q*r+(IX^2-IX*IY+IXZ^2)*p*q)/gam
    # attitude
    xfdot[7]=p+q*sphi*sthe/cthe+r*cphi*sthe/cthe
    xfdot[8]=q*cphi-r*sphi
    xfdot[9]=q*sphi/cthe+r*cphi/cthe
    # position
    xfdot[10]=u*cthe*cpsi+v*(sphi*sthe*cpsi-cphi*spsi)+w*(cphi*sthe*cpsi+sphi*spsi)
    xfdot[11]=u*cthe*spsi+v*(sphi*sthe*spsi+cphi*cpsi)+w*(cphi*sthe*spsi-sphi*cpsi)
    xfdot[12]=-u*sthe+v*sphi*cthe+w*cphi*cthe
    # return state derivative vector
    return xfdot
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
abstract type AbstractPositioner end
#All theory used here came from these sources:
#https://www.piezosystem.com/piezopedia/piezotheory/
#http://www.piezo.ws/piezoelectric_actuator_tutorial/Piezo_Design_part3.php

#Linear translation piezo device (one-dimension of motion only)
#TODO: add coefficient of friction? (should be very small)

struct Piezo <: AbstractPositioner
    id::String
    vmin_in::typeof(0.0V) #min acceptable amplifier input
    vmax_in::typeof(0.0V) #max acceptable amplifier input
    vmin_out::typeof(0.0V) #min MON (sensor) output
    vmax_out::typeof(0.0V) #max MON (sensor) output
    pmin_ol::typeof(0.0μm) #min position in open loop
    pmax_ol::typeof(0.0μm)
    closed_loop_pad::typeof(0.0μm) #Amount of the motion range that becomes inaccessible in closed loop mode.  This applies equally to the low and high side, so the total loss in range is twice this amount
    max_disp::typeof(0.0μm) #nominal max displacement.  One might expect this to be the same as the open loop range, but this contradicts the
                                #force specificatinos of our piezo.  (see below and more thoughts at the end of this file)
    cap::typeof(0.0μF)
    motion_axis::Int
    effective_mass::typeof(0.0g)
    resonance_freqs::NTuple{3, typeof(0.0Hz)} #resonance frequency without load (decreases with load, avoid operating at greater than 80% of resonant frequency)
    stiffs::NTuple{3, typeof(0.0N/μm)} #actuator stiffness
end

name(p::Piezo) = p.id
vmin_in(p::Piezo) = p.vmin_in
vmax_in(p::Piezo) = p.vmax_in
vspan_in(p::Piezo) = vmax_in(p) - vmin_in(p)
mon_slope(p::Piezo, is_cl::Bool) = vspan_out(p, is_cl) / pspan(p, is_cl) #Volts MON per micron.  TODO: this as well as vmin_out and vmax_out may be inaccurate in open loop mode
closed_loop_pad(p::Piezo) = p.closed_loop_pad
vmin_out(p::Piezo, is_cl::Bool) = p.vmin_out
vmax_out(p::Piezo, is_cl::Bool) = p.vmax_out
vspan_out(p::Piezo, is_cl::Bool) = vmax_out(p, is_cl) - vmin_out(p, is_cl)

pmin(p::Piezo, is_cl::Bool) = is_cl ? p.pmin_ol + closed_loop_pad(p) : p.pmin_ol
pmax(p::Piezo, is_cl::Bool) = is_cl ? p.pmax_ol - closed_loop_pad(p) : p.pmax_ol
pspan(p::Piezo, is_cl::Bool) = pmax(p, is_cl) - pmin(p, is_cl)

blocked_pad(p::Piezo) = (max_displacement(p) - pspan(p, false))/2 #If the max nominal displacement is greater than the open loop span of the piezo then there is a region that is blocked.  This is the half-size of the total blocked distance
function closed2open(pclosed::HasLengthUnits, p::Piezo)
    if pclosed < pmin(p, false) || pclosed > pmax(p, false)
        error("The input is not within the valid range of the piezo")
    end
    return pmin(p, false) + closed_loop_pad(p) + pclosed
end

function open2closed(popen::HasLengthUnits, p::Piezo)
    clp = closed_loop_pad(p)
    if popen < (pmin(p, false) + clp) || popen > (pmax(p, false) - clp)
        error("Position $popen is not achievable in closed-loop mode")
    end
    return pmin(p, false) + clp + popen #note this is trivial (equal to popen) when the closed-loop range begins with zero (all piezosystem jena piezos)
end

max_displacement(p::Piezo) = p.max_disp
capacitance(p::Piezo) = p.cap
motion_axis(p::Piezo) = p.motion_axis
stiffness(p::Piezo, dim::Int) = p.stiffs[dim]
stiffness(p::Piezo) = stiffness(p, motion_axis(p))
stiffness_all(p::Piezo) = p.stiffs
effective_mass(p::Piezo) = p.effective_mass
resonance(p::Piezo, dim::Int) = p.resonance_freqs[dim]
resonance(p::Piezo) = resonance(p, motion_axis(p))
resonance_all(p::Piezo) = p.resonance_freqs

#clamps a given position value to be within the nominal displacement range of the piezo
function clamp(pos::HasLengthUnits, p::Piezo)
    if pos < zero(pos)
        return zero(pos)
    else
        return min(pos, max_displacement(p))
    end
end

#Several functions are defined for converting input/output voltages and positions:
#amplifier voltage <-> positioner displacement <-> monitored voltage <-> monitored position

#Converts the monitor voltage into a position reading, does NOT correspond with the nominal displacement of the piezo.  Therefore the position reading should not be used
#for making force / compression calculations.
function mon2pos(mon::HasVoltageUnits, p::Piezo, is_cl::Bool)
    if mon < vmin_out(p, is_cl) || mon > vmax_out(p, is_cl)
        error("Monitor voltage lies outside the valid output range of this piezo.")
    end
    return pmin(p, is_cl) + inv(mon_slope(p, is_cl)) * (mon - vmin_out(p, is_cl))
end

#Returns the mon voltage for a given position
function pos2mon(pos::HasLengthUnits, p::Piezo, is_cl::Bool)
    clp = closed_loop_pad(p)
    max_p = pmax(p, is_cl)
    min_p = pmin(p, is_cl)
    if pos < pmin(p, false) || pos > pmax(p, false)
        error("Position lies outside the valid position range of this piezo.")
    end
    pos_frac = (pos - min_p) / (max_p - min_p)
    return vmin_out(p, is_cl) + pos_frac * vspan_out(p, is_cl)
end

#Converts nominal displacement to a physical position, clamping the position to be within the allowed range of the positioner
#(position is only accurate without external forces)
function displacement2pos(d::HasLengthUnits, p::Piezo, is_cl::Bool)
    md = max_displacement(p)
    if d > md || d < zero(d)
        error("Displacement must be nonnegative and no greater than the maximum displacement of the piezo")
    end
    extra_pad = blocked_pad(p)
    p_ol = zero(d)
    if d <= pspan(p, false) + extra_pad
        p_ol = pmin(p, false) + max(zero(d), d - extra_pad)
    else
        p_ol = pmax(p, false)
    end
    if is_cl
        return open2closed(p_ol, p)
    else
        return p_ol
    end
end

#This version accounts for length changes due to finite piezo stiffness. Thus it uses the net force acting on the positioner (positive force corresponds with positive voltage and positive change in position).
#The force calculation should include any load, the effective mass of the positioner itself, and gravity forces
displacement2pos(d::HasLengthUnits, p::Piezo, is_cl::Bool, f::HasForceUnits) = displacement2pos(d + uconvert(unit(d), f/stiffness(p)), p, is_cl)

function check_disp_extrema(pos::HasLengthUnits, p::Piezo; only_warn = false)
    if pos <= pmin(p, false) || pos >= pmax(p, false)
        f = only_warn ? warn : error
        f("The nominal displacement can't reliably be calculated from the position when position reaches the extremes of the operating range")
    end
end

#The monitored signal doesn't give the displacement that is useful for piezo force calculations.  This function transforms the mon value to that more useful value.
#This version takes the position according to piezo sensor output and converts that to a nominal piezo displacement useful for force calculations
#Note that it's unreliable at the extremes of the piezo's range because the piezo could be generating extra force (and thus a greater nominal displacement) against a stop
#If you want to proceed anyway call the unsafe_ version
unsafe_pos2displacement(pos::HasLengthUnits, p::Piezo) = blocked_pad(p) + (pos - pmin(p, false))
function pos2displacement(pos::HasLengthUnits, p::Piezo, is_cl::Bool)
    if is_cl
        pos = closed2open(pos, p)
    end
    check_disp_extrema(pos, p)
    return unsafe_pos2displacement(pos, p)
end

#This version considers a load
#nom_displacement + f/stiffness(p) = measured_displacement
#nom_displacement = measured_disp - f/stiffness(p)
pos2displacement(pos::HasLengthUnits, p::Piezo, is_cl::Bool, f::HasForceUnits) =  pos2displacement(pos, p, is_cl) - f/stiffness(p)

#get the (nominal) displacement of the positioner for a given input voltage without any load
#like the nominal max displacement, this may exceed the allowed displacement of the piezo
#We assume there is an internal stop narrowing the allowed displacement, but for the force calculation the full displacement must be used
function input2displacement(v::HasVoltageUnits, p::Piezo)
    if v < vmin_in(p) || v > vmax_in(p)
        error("Input voltage lies outside the valid input range of this piezo.")
    end
    v_frac = (v-vmin_in(p)) / vspan_in(p)
    md = max_displacement(p)
    return v_frac * md
end

#Calculate volts required to generate a given nominal displacement
#(piezo might not actually reach that displacement due to blocking forces)
function displacement2input(d::HasLengthUnits, p::Piezo)
    @assert d <= max_displacement(p) && d >= zero(d)
    return vmin_in(p) + (d/max_displacement(p)) * vspan_in(p)
end

function input2pos(v::HasVoltageUnits, p::Piezo, is_cl::Bool, f::HasForceUnits)
    d0 = input2displacement(v, p)
    return displacement2pos(d0, p, is_cl, f)
end
input2pos(v::HasVoltageUnits, p::Piezo, is_cl::Bool) = input2pos(v, p, is_cl, 0.0N)

pos2input(pos::HasLengthUnits, p::Piezo, is_cl::Bool, f::HasForceUnits) = displacement2input(pos2displacement(pos, p, is_cl, f), p)
pos2input(pos::HasLengthUnits, p::Piezo, is_cl::Bool) = displacement2input(pos2displacement(pos, p, is_cl, 0.0N), p)

#adjusted_displacement - compression = mon_displacement
#adjusted_displacement = mon_displacement + compression
#                  = mon_displacement + f / stiffness(p)
mon2displacement(mon::HasVoltageUnits, p::Piezo, is_cl::Bool, f::HasForceUnits) = pos2displacement(mon2pos(mon, p, is_cl), p, is_cl, f)

mon2displacement(mon::HasVoltageUnits, p::Piezo, is_cl::Bool) = mon2displacement(mon, p, is_cl, 0.0N)

max_force(stif::HasStiffnessUnits, displacement::HasLengthUnits) = stif * displacement

#TODO: consult with Piezosystem Jena to make sure this is accurate.  As noted elsewhere in code comments, it's not obvious how to do this calculation with the NanoSX positioners.
function max_force(p::Piezo, cur_disp::HasLengthUnits, direction::Int)
    d = direction > 0 ? max_displacement(p) - cur_disp : -cur_disp
    return max_force(stiffness(p), d)
end

function max_force(p::Piezo, pos::HasLengthUnits, is_cl::Bool, direction::Int)
    if direction == 0
        error("Direction of motion must be positive or negative")
    end
    kt = stiffness(p)
    nd = max_displacement(p)
    if is_cl
        pos = closed2open(pos, p)
    end
    cur_disp = blocked_pad(p) + pos - pmin(p, false)
    return max_force(p, cur_disp, direction)
end

function max_force(p::Piezo, cur_pos::HasLengthUnits, new_pos::HasLengthUnits, is_cl::Bool)
    direction = sign(new_pos - cur_pos)
    return max_force(p, cur_pos, is_cl, direction)
end

#Useful for NanoSX piezos which do not specify effective mass
#Instead they specify resonance frequency without load and with one more more loads.
#This solves for effective mass using two points on the curve
function effective_mass(f0, f_load, m_load)
    #77.1917g for NanoSX800 using the zero and f_load = 300g points (note the equation doesn't fit the points they gave exactly, seems to fit better at higher loads)
    #In theory this one should be equivalent:
    #effective_mass(f0, stif) = stif/(2*pi*f0)^2
    temp = f_load^2 / f0^2
    return  -m_load * temp / (temp - 1)
end

loaded_resonance(f0::HasInverseTimeUnits, m_eff::HasMassUnits, m_load::HasMassUnits) = uconvert(Unitful.Hz, f0 * sqrt(m_eff/(m_eff+m_load)))
loaded_resonance(p::Piezo, m_load::HasMassUnits) = loaded_resonance(resonance(p), effective_mass(p), m_load)

function mod2pos(mod::HasVoltageUnits, a::SimpleAmplifier, p::Piezo, is_cl::Bool)
    @assert mod >= in_min(a) && mod <= in_max(a)
    mod_frac = (mod - in_min(a)) / in_span(a)
    return pmin(p, is_cl) + mod_frac * pspan(p, is_cl)
end

function pos2mod(pos::HasLengthUnits, a::SimpleAmplifier, p::Piezo, is_cl::Bool)
    @assert pos >= pmin(p, is_cl) && pos <= pmax(p, is_cl)
    pos_frac = (pos - pmin(p, is_cl)) / pspan(p, is_cl)
    return in_min(a) + pos_frac * in_span(a)
end

mod2displacement(mod::HasVoltageUnits, a::SimpleAmplifier, p::Piezo, is_cl::Bool) = pos2displacement(mod2pos(mod, a, p, is_cl), p, is_cl)

function displacement2mod(disp::HasLengthUnits, p::Piezo, a::SimpleAmplifier, is_cl::Bool)
    if !is_cl
        error("Don't yet trust this function in open-loop mode")
    end
    frac = disp / max_displacement(p)
    return in_min(a) + frac * (in_max(a) - in_min(a))
end

POSITIONERS = Dict{String, Any}()

function retrieve_positioner(pos_name::String)
    if !haskey(POSITIONERS, pos_name)
        error("Unknown positioner")
    end
    return POSITIONERS[pos_name]
end

function register_positioner(pos::T) where {T<:AbstractPositioner}
    if haskey(POSITIONERS, name(pos))
        error("An positioner named $(name(pos)) is already registered.")
    end
    POSITIONERS[name(pos)] = pos
end

#The NanoX / NanoSX series piezos actually use two opposing ceramic piezo stacks per positioner.
#When one pushes, the other pulls.  In principle this could allow the maximum generable force of the
#device to remain constant over the entire range of motion.  However it's not clear from spec sheets
#that this is true.
#Alternatively if the push and pull forces are additive then the maximum achievable force can be calculated
#like a conventional piezo.  The force would simply be double the force that a single piezo could generate
#for a given displacement (calculated as the absolute difference between current position and position attained without
#resisting forces
#A maximum push/pull force of 100.0N is listed, along with an open loop diplacement max
#of 900um.  If one tries to apply the basic equation for force generation:
#F = kt * L (where kt is stiffness of the piezo in N/um and L is max displacement in um)
#Then it's not possible to produce a constant 100.0N force throughout the range using this push/pull strategy.
#(the displacements of the two would have to be offset by 500um in that case, which means they would go negative at the extremes of the range)
#Instead 180.0N comes up as the max force when you use 900um as the nominal max
#Some possible complications:
#1. 900um is not the "nominal max displacement"
#2. 100.0N is not the max generable force
#Option 1 makes sense only if the real displacement is greater that 900, which would mean a max force greater than 180, not less
#Option 2 is possible if they are somehow being conservative
#OR the push/pull forces are additive, meaning that the most force can be generated at the edges of the piezo's range (only
# when that force is directed toward the center of the range)
#That's my favorite option for now.  But then it's still unclear where 100N comes from.  If the nominal max disp is really
#900um then 100N can only be achieved between 250um and 650um in open loop mode.
#also the "max load" is only 50N, suggesting it can't do 100N throughout the range.
#but these numbers suggest that 50N is only achievable from 125um to 775um open loop
#Let's let nominal max displacement be larger than specified:
#to achieve 50N throughout the open loop range then the true range must be larger than that.  How much larger?
#50N = d * 0.2 * 2
#d = 125um on either side of the range
#So the nominal displacement would be 900um + 125um*2 = 1150um


text:
using ConfParser
import Dates
using AWSCore
using AWS: @service
@service S3

# load config
conf = ConfParse("conf/config.ini")
parse_conf!(conf)
download_dir = retrieve(conf, "dirs", "download")
aws_bucket   = retrieve(conf, "data", "aws_bucket")
gfs_prefix   = retrieve(conf, "data", "gfs_prefix")
gfs_postfix  = retrieve(conf, "data", "gfs_postfix")

@doc """ Downloads all gfs files for given run from NOAA's AWS bucket
"""
function download_gfs_data(run::Int64)

    # build daily files array
    current_date = Dates.now()
    current_hour = Dates.format(current_date, "HH")
    current_date_string = "gfs." * Dates.format(current_date, "yyyymmdd")
    current_files = []

    if (length(string(run)) == 1) && run in [0, 6]
        run = "0" * string(run)
    elseif length(string(run)) == 2 && run in [12, 18]
        run = string(run)
    else
        throw(ArgumentError("run must be in [0,6,12,18]"))
    end

    for i in 0:120
        if length(string(i)) == 1
            i = lpad(string(i), 3, '0')
        elseif length(string(i)) == 2
            i = lpad(string(i), 3, '0')
        end
        # build path
        push!(current_files,joinpath(current_date_string,run,(gfs_prefix * string(run) * gfs_postfix * string(i))))
    end

    for i in 120:3:384
        push!(current_files,joinpath(current_date_string,run,(gfs_prefix * string(run) * gfs_postfix * string(i))))
    end

    aws = aws_config()

    current_dir = pwd()
    download_dir_path = joinpath(current_dir, download_dir)
    daily_dir_path = joinpath(download_dir_path, current_date_string)
    run_dir_path = joinpath(daily_dir_path, run)

    for dir in [download_dir_path, daily_dir_path, run_dir_path]
        if !isdir(dir)
            mkdir(dir)
        end
    end

    for file in current_files
        println("Downloading ", file)
        try
            file_name = file[end-24:end]
            file_path = download_dir_path * "/" * file # for any reason, joinpath does not work here

            if !isfile(file_path)
                s3_object = S3.get_object(aws_bucket, file)
                io = open(file_path,"w")
                write(io, s3_object)
                close(io)
            else
                println("File ", file, " already downloaded.")
            end

        catch
            println("Error at downloading", file)
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__

    if length(ARGS) == 0
        println("Please specify run to download with julia app.py [run::int in [0,6,12,18]]")
    elseif length(ARGS) == 1
        run = parse(Int,ARGS[1])
        if (length(string(run)) == 1) && run in [0, 6]
            download_gfs_data(run)
        elseif length(string(run)) == 2 && run in [12, 18]
            download_gfs_data(run)
        else
            throw(ArgumentError("run must be in [0,6,12,18]"))
        end
    else
        throw(ArgumentError("Too many arguments"))
    end
end


text:
include(joinpath(@__DIR__, "../../src/sketches/FindModel.jl"))
using Test

"""Define example FLS"""

#########################
# COMMUTATIVE TRIANGLE  #
#########################

f = FLSinit(@acset LabeledGraph begin
    V = 3
    E = 5
    vlabel = [:A,:B,:C]
    elabel = [:f,:h,:g, :ida, :inv]
    src    = [1,  1, 2,  1,    3]
    tgt    = [2,  3, 3,  1,    3]
end)

# f;g = h
add_diagram!(f, @acset LabeledGraph begin
    V = 3
    E = 3
    vlabel = [:A,:B,:C]
    elabel = [:f,:g,:h]
    src    = [1, 2, 1]
    tgt    = [2, 3, 3]
end)

# "ida" is an identity arrow
add_diagram!(f, @acset LabeledGraph begin
    V = 1
    E = 1
    vlabel = [:A]
    elabel = [:ida]
    src = [1]
    tgt = [1]
end);

# "inv" is an involution on C
add_diagram!(f, @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:C, :C]
    elabel = [:inv, :inv]
    src = [1, 2]
    tgt = [2, 1]
end)

#############################
# LEFT INVERSE / INVOLUTION #
#############################

leftinv = FLSinit(@acset LabeledGraph begin
    V = 2
    E = 4
    vlabel = [:A,:B]
    elabel = [:f,:g, :ida, :inv]
    src    = [1,  2,  1,   2]
    tgt    = [2,  1,  1,   2]
end);

add_diagram!(leftinv, @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:A, :B]
    elabel = [:f, :g]
    src = [1, 2]
    tgt = [2, 1]
end);

add_diagram!(leftinv, @acset LabeledGraph begin
    V = 1
    E = 1
    vlabel = [:A]
    elabel = [:ida]
    src = [1]
    tgt = [1]
end)

add_diagram!(leftinv, @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:B, :B]
    elabel = [:inv, :inv]
    src = [1, 2]
    tgt = [2, 1]
end)

#############
# EQUALIZER #
#############

fxb = Symbol("f×ᵦg")

equalizer = FLSinit(@acset LabeledGraph begin
    V = 3
    E = 4
    vlabel = [:A,:B,fxb]
    elabel = [:f,:g, :π1, :π2]
    src = [1,1,3,3]
    tgt = [2,2,1,1]
end)

add_cone!(equalizer, @acset LabeledGraph begin
    V = 4
    E = 4
    vlabel = [:A,:B,:A,fxb]
    elabel = [:f,:g, :π1, :π2]
    src = [1,3,4,4]
    tgt = [2,2,1,3]
end)

#################################
# SET PERMUTATIONS / PARTITIONS #
#################################

perm = FLSinit(@acset LabeledGraph begin
    V = 1
    E = 2
    vlabel = [:X]
    elabel = [:f, :f⁻¹]
    src = [1,1]
    tgt = [1,1]
end );

add_diagram!(perm, @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:X, :X]
    elabel = [:f, :f⁻¹]
    src = [1, 2]
    tgt = [2, 1]
end)

####################
# REFLEXIVE GRAPHS #
####################

reflg = FLSinit(@acset LabeledGraph begin
    V = 2
    E = 4
    vlabel = [:V, :E]
    elabel = [:id, :refl, :src, :tgt]
    src = [1,1,2,2]
    tgt = [1,2,1,1]
end );

add_diagram!(reflg, @acset LabeledGraph begin
    V = 2
    E = 4
    vlabel = [:V, :E]
    elabel = [:id, :refl, :src, :tgt]
    src = [1,1,2,2]
    tgt = [1,2,1,1]
end)

###########################
# CONSTANT (arrow from 1) #
###########################
s1 = Symbol("1")
const2 = FLSinit(@acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [s1, :A]
    elabel = [:f, :g]
    src = [1,1]
    tgt = [2,2]
end);

add_cone!(const2, @acset LabeledGraph begin
    V = 1
    vlabel=[s1]
end)

####################################
# Infinite sequence of a's and b's #
####################################

inflist = FLSinit(@acset LabeledGraph begin
    V = 3
    E = 4
    vlabel = [s1, :d, :l]
    elabel = [:a,:b,:head,:tail]
    src = [1,1,3,3]
    tgt = [2,2,2,3]
end )

add_cone!(inflist, @acset LabeledGraph begin
    V = 1
    vlabel = [s1]
end)

add_cone!(inflist, @acset LabeledGraph begin
    V = 3
    E = 2
    vlabel = [:d, :l, :l]
    elabel = [:head, :tail]
    src = [3,3]
    tgt = [1,2]
end)

#####################
# Cartesian Product #
#####################

pairs = FLSinit(@acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:s, :s2]
    elabel = [:p1, :p2]
    src = [2,2]
    tgt = [1,1]
end ); # s s²

add_cone!(pairs, @acset LabeledGraph begin
    V = 3
    E = 2
    vlabel = [:s, :s, :s2]
    elabel = [:p1, :p2]
    src = [3,3]
    tgt = [1,2]
end )

srch_pairs = cs -> [catelems_to_cset(pairs, r) for r in search(pairs, cs)]

trips = FLSinit(@acset LabeledGraph begin
    V = 2
    E = 3
    vlabel = [:s, :s3]
    elabel = [:p1, :p2, :p3]
    src = [2,2,2]
    tgt = [1,1,1]
end ); # s s²

add_cone!(trips, @acset LabeledGraph begin
    V = 4
    E = 3
    vlabel = [:s, :s, :s, :s3]
    elabel = [:p1, :p2, :p3]
    src = [4,4,4]
    tgt = [1,2,3]
end )

##############
# SEMIGROUPS #
##############

p1p2, p2p3, idk, kid = map(Symbol, ["π₁×π₂","π₂×π₃","id×k","k×id"])
semig = FLSinit(@acset LabeledGraph begin
    V = 3
    E = 10
    vlabel = [:s, :s2, :s3]
    elabel = [:k, :π₁, :π₂, :Π₁, :Π₂, :Π₃, p1p2, p2p3, idk, kid]
    src    = [2,  2,   2,   3,    3,   3,   3,   3,    3,   3]
    tgt    = [1,  1,   1,   1,    1,   1,   2,   2,    2,   2]
end );

# Associativity
add_diagram!(semig, @acset LabeledGraph begin
    V = 4
    E = 4
    vlabel = [:s3, :s2, :s2, :s]
    elabel = [kid, idk, :k, :k]
    src = [1,1,2,3]
    tgt = [2,3,4,4]
end)

# Define p1p2
add_diagram!(semig, @acset LabeledGraph begin
    V = 4
    E = 5
    vlabel = [:s3, :s, :s, :s2]
    elabel = [:Π₁, :Π₂, p1p2, :π₁, :π₂]
    src = [1,1,1,4,4]
    tgt = [2,3,4,2,3]
end)

# Define p2p3
add_diagram!(semig, @acset LabeledGraph begin
    V = 4
    E = 5
    vlabel = [:s3, :s, :s, :s2]
    elabel = [:Π₂, :Π₃, p2p3, :π₁, :π₂]
    src    = [1,    1,   1,   4,   4]
    tgt    = [2,    3,   4,   2,   3]
end)

# Define kid
add_diagram!(semig, @acset LabeledGraph begin
    V = 5
    E = 6
    vlabel = [:s3, :s2,  :s2, :s, :s]

    elabel = [p1p2, kid, :Π₃, :π₂, :k, :π₁]
    src    = [1,    1,    1,   3,   2,  3]
    tgt    = [2,    3,    4,   4,   5,  5]
end)

# Define idk
add_diagram!(semig, @acset LabeledGraph begin
    V = 5
    E = 6
    vlabel = [:s3, :s2,  :s2, :s, :s]

    elabel = [p2p3, idk, :Π₁, :π₁, :k, :π₂]
    src    = [1,    1,    1,   3,   2,  3]
    tgt    = [2,    3,    4,   4,   5,  5]
end)

# s2 is pair
add_cone!(semig, @acset LabeledGraph begin
    V = 3
    E = 2
    vlabel = [:s, :s, :s2]
    elabel = [:π₁, :π₂]
    src = [3,3]
    tgt = [1,2]
end )

# s3 is triple
add_cone!(semig, @acset LabeledGraph begin
    V = 4
    E = 3
    vlabel = [:s, :s, :s, :s3]
    elabel = [:Π₁, :Π₂, :Π₃]
    src = [4,4,4]
    tgt = [1,2,3]
end )

# Colimits
import DataStructures: OrderedDict
function fls_pushout(l::LabeledFLS, r::LabeledFLS; kw...)::Fls
    ns, es = [OrderedDict([k=>v for (k, v) in kw if k in l.labels[i]])
              for i in [1,2]]
    overlap = Fls()
    add_parts!(overlap, :V, length(ns))
    add_parts!(overlap, :E, length(es))
    println("kw $kw Ns $ns Es $es")
    L = CSetTransformation(overlap, l.fls; V=[vi(l, x) for x in keys(ns)],
                                           E=[ei(l, x) for x in keys(es)])
    R = CSetTransformation(overlap, r.fls; V=[vi(r, x) for x in values(ns)],
                                           E=[ei(r, x) for x in values(es)])
    i, _ = pushout(L, R)
    return codom(i)
end

po = fls_pushout(reflg, equalizer; V=:B, E=:A, s=:f, t=:g)


text:
using PhaseTypeDistributions
using SparseMatrix
using NMarkov
using Test
using Distributions
using LinearAlgebra
using SparseArrays

@testset "Test of CF1" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 10.0]
    cf1 = CF1(alpha, rate)
    @test cf1.rate == [0.4, 1.4, 10.0]
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 10.0]
    cf1 = CF1(alpha, rate)
    gph = GPH(cf1, SparseCSC)
    @test gph.T.rowind == [1, 1, 2, 2, 3]
    @test gph.T.val == [-0.4, 0.4, -1.4, 1.4, -10.0]
    gph = GPH(cf1, SparseCSR)
    @test gph.T.colind == [1, 2, 2, 3, 3]
    @test gph.T.val == [-0.4, 0.4, -1.4, 1.4, -10.0]
    gph = GPH(cf1, SparseMatrixCSC)
    @test gph.T.rowval == [1, 1, 2, 2, 3]
    @test gph.T.nzval == [-0.4, 0.4, -1.4, 1.4, -10.0]
end

@testset "Test for wsample" begin
    t = rand(10)
    w = rand(10)
    dat = WeightedSample(t, w)
    @test dat.tdat[1] == minimum(t)
    t = rand(100)
    dat = PointSample(t)
    @test dat.tdat[1] == minimum(t)
    @test all(dat.wdat .== 1.0)
end

@testset "Test of estep" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    cf1 = CF1(alpha, rate)
    ph = GPH(cf1, SparseCSC)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(cf1, dat, eres)
    @time llf = estep!(cf1, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    ph = GPH(cf1, SparseCOO)
    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)
    @time llf = estep!(cf1, dat, eres)
    @time llf = estep!(cf1, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    ph = GPH(cf1, SparseCSR)
    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)
    @time llf = estep!(cf1, dat, eres)
    @time llf = estep!(cf1, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    ph = GPH(cf1, SparseMatrixCSC)
    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)
    @time llf = estep!(cf1, dat, eres)
    @time llf = estep!(cf1, dat, eres)
    @test sum(eres.ez) ≈ sum(t)
end

@testset "Test of emstep" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    cf1 = CF1(alpha, rate)
    ph = GPH(cf1)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(cf1, dat, eres)
    @time llf = estep!(cf1, dat, eres)
    println(llf)
    @test sum(eres.ez) ≈ sum(t)

    mstep!(cf1, eres)
    @time llf2 = estep!(cf1, dat, eres)
    @test llf2 > llf
end

@testset "Test of emstep gph Matrix" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    ph = GPH(CF1(alpha, rate), Matrix)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(ph, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    mstep!(ph, eres)
    @time llf2 = estep!(ph, dat, eres)
    @test llf2 > llf
end

@testset "Test of emstep gph CSR" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    ph = GPH(CF1(alpha, rate), SparseCSR)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(ph, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    mstep!(ph, eres)
    @time llf2 = estep!(ph, dat, eres)
    @test llf2 > llf
end

@testset "Test of emstep gph CSC" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    ph = GPH(CF1(alpha, rate), SparseCSC)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(ph, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    mstep!(ph, eres)
    @time llf2 = estep!(ph, dat, eres)
    @test llf2 > llf
end

@testset "Test of emstep COO" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    ph = GPH(CF1(alpha, rate), SparseCOO)

    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time llf = estep!(ph, dat, eres)
    @test sum(eres.ez) ≈ sum(t)

    mstep!(ph, eres)
    @time llf2 = estep!(ph, dat, eres)
    @test llf2 > llf
end

@testset "Test of emstep for 10" begin
    alpha = [0.1, 0.3, 0.6]
    rate = [1.4, 0.4, 2.0]
    cf1 = CF1(alpha, rate)
    ph = GPH(cf1)
    t = rand(100)
    dat = PointSample(t)
    eres = Estep(ph)

    @time prevllf = estep!(cf1, dat, eres)
    mstep!(cf1, eres)
    for k = 1:10
        @time llf = estep!(cf1, dat, eres)
        mstep!(cf1, eres)
        @test llf > prevllf
        prevllf = llf
    end
end

@testset "Test of cf1init" begin
    t = rand(100)
    dat = PointSample(t)
    cf1 = initializePH(CF1(10), dat, verbose = true)
    # println(cf1)
    @test true
end

@testset "Test of cf1init 2" begin
    dat = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    cf1 = initializePH(CF1(10), dat, verbose = true)
    @test true
end

@testset "Test of cf1init 2" begin
    dat = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    cf1 = initializePH(CF1(10), dat, verbose = true)
    phfit!(cf1, dat, verbose = true)
    # println(cf1)
    @test true
end

@testset "Test of phfit 1" begin
    res = phfit(CF1(10), verbose=[true, true]) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    println(res)
    @test true
end

@testset "Test of pdf 1" begin
    res = phfit(CF1(100), verbose=[true, true]) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    ph = GPH(res[1])
    @test phpdf(res[1], 0.0) == sum(ph.alpha .* ph.tau)
    @test phcdf(res[1], 0.0) ≈ 0.0
    @test phccdf(res[1], 0.0) ≈ 1.0
end

@testset "Test of pdf 1" begin
    res = phfit(CF1(200), verbose=[true, true]) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    t = LinRange(0.0, 2.0, 100)
    x = phpdf(res[1], t)
    y = pdf.(Weibull(2.0, 1.0), t)
    for i = eachindex(x)
        println("$(i) $(x[i]) $(y[i])")
        @test x[i] - y[i] < 1.0e-4
    end
end

@testset "Test of phfit 1" begin
    cf1, llf0, = phfit(CF1(10), verbose=[true, true]) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    data = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    llf = phllf(cf1, data)
    @test abs((llf0 - llf) / llf) < 1.0e-6
end

@testset "Test of phfit 2" begin
    cf1, llf0, = phfit(CF1(10), verbose=[true, true]) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    data = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    println(cf1)

    deriv = Dict{Symbol,CF1{Float64}}()
    for i = 1:cf1.dim
        s = Symbol(:rate, i)
        a = zero(cf1.alpha)
        x = zero(cf1.rate)
        x[i] = 1.0
        deriv[s] = CF1(cf1.dim, a, x)

        s = Symbol(:alpha, i)
        a = zero(cf1.alpha)
        x = zero(cf1.rate)
        a[i] = 1.0
        deriv[s] = CF1(cf1.dim, a, x)
    end
    llf, llfdash = phllf(cf1, deriv, data)
    println(llfdash)

    h = 0.000001
    tol = 1e-6
    for i = 1:cf1.dim
        s = Symbol(:alpha, i)
        cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
        cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
        cf1h1.alpha[i] -= h
        cf1h2.alpha[i] += h
        llf1 = phllf(cf1h1, data)
        llf2 = phllf(cf1h2, data)
        tmp = 1.0 + abs((llf2 - llf1)) / (2*h)
        error = abs((abs(llfdash[s]) + 1.0 - tmp) / (1.0 + abs(llfdash[s])))
        println((s, abs((llf2 - llf1)) / (2*h), abs(llfdash[s])))
        @test error < tol
    end
    for i = 1:cf1.dim
        s = Symbol(:rate, i)
        cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
        cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
        cf1h1.rate[i] -= h
        cf1h2.rate[i] += h
        llf1 = phllf(cf1h1, data)
        llf2 = phllf(cf1h2, data)
        tmp = 1.0 + abs((llf2 - llf1)) / (2*h)
        error = abs((abs(llfdash[s]) + 1.0 - tmp) / (1.0 + abs(llfdash[s])))
        println((s, abs((llf2 - llf1)) / (2*h), abs(llfdash[s])))
        @test error < tol
    end
end

@testset "Test of phfit 3" begin
    cf1, llf0, = phfit(CF1(5), verbose=[true, true]) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    data = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(2.0, 1.0), x)
    end
    println(cf1)

    deriv = Dict{Symbol,CF1{Float64}}()
    deriv2 = Dict{Tuple{Symbol,Symbol},CF1{Float64}}()
    for i = 1:cf1.dim-1
        s = Symbol(:alpha, i)
        a = zeros(Float64, cf1.dim)
        x = zeros(Float64, cf1.dim)
        a[i] = 1.0
        a[cf1.dim] = -1.0
        deriv[s] = CF1(cf1.dim, a, x)
        for j = i:cf1.dim-1
            s2 = Symbol(:alpha, j)
            a = zeros(Float64, cf1.dim)
            x = zeros(Float64, cf1.dim)
            deriv2[(s,s2)] = CF1(cf1.dim, a, x)
        end
        for j = 1:cf1.dim
            s2 = Symbol(:rate, j)
            a = zeros(Float64, cf1.dim)
            x = zeros(Float64, cf1.dim)
            deriv2[(s,s2)] = CF1(cf1.dim, a, x)
        end
    end
    for i = 1:cf1.dim
        s = Symbol(:rate, i)
        a = zeros(Float64, cf1.dim)
        x = zeros(Float64, cf1.dim)
        x[i] = 1.0
        deriv[s] = CF1(cf1.dim, a, x)
        for j = i:cf1.dim
            s2 = Symbol(:rate, j)
            a = zeros(Float64, cf1.dim)
            x = zeros(Float64, cf1.dim)
            deriv2[(s,s2)] = CF1(cf1.dim, a, x)
        end
    end
    llf, llfdash, llfdashdash = phllf(cf1, deriv, deriv2, data)
    llf0, llfdash0 = phllf(cf1, deriv, data)
    @test llf ≈ llf0
    for i = eachindex(llfdash0)
        @test llfdash0[i] ≈ llfdash[i]
    end

    k = 1
    index = Dict()
    for i = 1:cf1.dim-1
        s = Symbol(:alpha, i)
        index[s] = k
        k += 1
    end
    for i = 1:cf1.dim
        s = Symbol(:rate, i)
        index[s] = k
        k += 1
    end
    IM = zeros(length(index), length(index))
    for i = eachindex(llfdashdash)
        IM[index[i[1]], index[i[2]]] = llfdashdash[i]
        IM[index[i[2]], index[i[1]]] = llfdashdash[i]
    end
    values,vectors = eigen(-IM)
    println(values)
    @test all(values .> 0)
    println(inv(-IM))
end

@testset "Test of phfit 4" begin
    cf1, llf0, = phfit(CF1(10), verbose=[true, true]) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    data = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    println(cf1)

    deriv = Dict{Symbol,CF1{Float64}}() ## １階微分
    for i = 1:cf1.dim
        s = Symbol(:alpha, i)
        a = zero(cf1.alpha)
        x = zero(cf1.rate)
        a[i] = 1.0
        deriv[s] = CF1{Float64}(cf1.dim, a, x)
    end
    for i = 1:cf1.dim
        s = Symbol(:rate, i)
        a = zero(cf1.alpha)
        x = zero(cf1.rate)
        x[i] = 1.0
        deriv[s] = CF1{Float64}(cf1.dim, a, x)
    end
    deriv2 = Dict{Tuple{Symbol,Symbol},CF1{Float64}}() ## 2階微分
    for i = 1:cf1.dim
        s1 = Symbol(:alpha, i)
        for j = i:cf1.dim
            s2 = Symbol(:alpha, j)
            a = zero(cf1.alpha)
            x = zero(cf1.rate)
            deriv2[(s1,s2)] = CF1{Float64}(cf1.dim, a, x)
        end
        for j = 1:cf1.dim
            s2 = Symbol(:rate, j)
            a = zero(cf1.alpha)
            x = zero(cf1.rate)
            deriv2[(s1,s2)] = CF1{Float64}(cf1.dim, a, x)
        end
    end
    for i = 1:cf1.dim
        s1 = Symbol(:rate, i)
        for j = i:cf1.dim
            s2 = Symbol(:rate, j)
            a = zero(cf1.alpha)
            x = zero(cf1.rate)
            deriv2[(s1,s2)] = CF1{Float64}(cf1.dim, a, x)
        end
    end

    delta = 0.00001
    tol = 1e-5
    llf, llfdash, llfdashdash = phllf(cf1, deriv, deriv2, data)
    for i = 1:cf1.dim
        for j = i:cf1.dim
            s1 = Symbol(:alpha, i)
            s2 = Symbol(:alpha, j)
            cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h3 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h4 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h1.alpha[i] -= delta
            cf1h2.alpha[i] -= delta
            cf1h1.alpha[j] -= delta
            cf1h2.alpha[j] += delta
            cf1h3.alpha[i] += delta
            cf1h4.alpha[i] += delta
            cf1h3.alpha[j] -= delta
            cf1h4.alpha[j] += delta
            lf1 = phllf(cf1h1, data)
            lf2 = phllf(cf1h2, data)
            lf3 = phllf(cf1h3, data)
            lf4 = phllf(cf1h4, data)
            tmp = abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2) + 1.0
            error = (abs(llfdashdash[s1,s2]) + 1.0 - tmp) / (abs(llfdashdash[s1,s2]) + 1.0)
            println((s1, s2, abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2), abs(llfdashdash[s1,s2])))
            if i < cf1.dim && j < cf1.dim
                @test abs(error) < tol
            end
        end
        for j = 1:cf1.dim
            s1 = Symbol(:alpha, i)
            s2 = Symbol(:rate, j)
            cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h3 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h4 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h1.alpha[i] -= delta
            cf1h2.alpha[i] -= delta
            cf1h1.rate[j] -= delta
            cf1h2.rate[j] += delta
            cf1h3.alpha[i] += delta
            cf1h4.alpha[i] += delta
            cf1h3.rate[j] -= delta
            cf1h4.rate[j] += delta
            lf1 = phllf(cf1h1, data)
            lf2 = phllf(cf1h2, data)
            lf3 = phllf(cf1h3, data)
            lf4 = phllf(cf1h4, data)
            tmp = abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2) + 1.0
            error = (abs(llfdashdash[s1,s2]) + 1.0 - tmp) / (abs(llfdashdash[s1,s2]) + 1.0)
            println((s1, s2, abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2), abs(llfdashdash[s1,s2])))
            @test abs(error) < tol
        end
    end
    for i = 1:cf1.dim
        for j = i:cf1.dim
            s1 = Symbol(:rate, i)
            s2 = Symbol(:rate, j)
            cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h3 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h4 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
            cf1h1.rate[i] -= delta
            cf1h2.rate[i] -= delta
            cf1h1.rate[j] -= delta
            cf1h2.rate[j] += delta
            cf1h3.rate[i] += delta
            cf1h4.rate[i] += delta
            cf1h3.rate[j] -= delta
            cf1h4.rate[j] += delta
            lf1 = phllf(cf1h1, data)
            lf2 = phllf(cf1h2, data)
            lf3 = phllf(cf1h3, data)
            lf4 = phllf(cf1h4, data)
            tmp = abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2) + 1.0
            error = (abs(llfdashdash[s1,s2]) + 1.0 - tmp) / (abs(llfdashdash[s1,s2]) + 1.0)
            println((s1, s2, abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2), abs(llfdashdash[s1,s2])))
            @test abs(error) < tol
        end
    end
end

@testset "Test of phfit 5" begin
    cf1, llf0, = phfit(CF1(10), verbose=[true, true]) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    data = WeightedSample((0.0, Inf64)) do x
        pdf(Weibull(1.5, 1.0), x)
    end
    println(cf1)

    deriv = Dict{Symbol,CF1{Float64}}() ## １階微分
    i = cf1.dim
    s = Symbol(:alpha, i)
    a = zero(cf1.alpha)
    x = zero(cf1.rate)
    a[i] = 1.0
    deriv[s] = CF1{Float64}(cf1.dim, a, x)
    deriv2 = Dict{Tuple{Symbol,Symbol},CF1{Float64}}() ## 2階微分
    i = cf1.dim
    j = cf1.dim
    s1 = Symbol(:alpha, i)
    s2 = Symbol(:alpha, j)
    a = zero(cf1.alpha)
    x = zero(cf1.rate)
    deriv2[(s1,s2)] = CF1{Float64}(cf1.dim, a, x)
    llf, llfdash, llfdashdash = phllf(cf1, deriv, deriv2, data)

    delta = 0.000001
    tol = 1e-5
    i = cf1.dim
    j = cf1.dim
    s1 = Symbol(:alpha, i)
    s2 = Symbol(:alpha, j)
    cf1h1 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
    cf1h2 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
    cf1h3 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
    cf1h4 = CF1(cf1.dim, copy(cf1.alpha), copy(cf1.rate))
    cf1h1.alpha[i] -= delta
    cf1h2.alpha[i] -= delta
    cf1h1.alpha[j] -= delta
    cf1h2.alpha[j] += delta
    cf1h3.alpha[i] += delta
    cf1h4.alpha[i] += delta
    cf1h3.alpha[j] -= delta
    cf1h4.alpha[j] += delta
    lf1 = phllf(cf1h1, data)
    lf2 = phllf(cf1h2, data)
    lf3 = phllf(cf1h3, data)
    lf4 = phllf(cf1h4, data)
    tmp = abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2) + 1.0
    error = (abs(llfdashdash[s1,s2]) + 1.0 - tmp) / (abs(llfdashdash[s1,s2]) + 1.0)
    println((s1, s2, abs(((lf4 - lf3) - (lf2 - lf1)) / (2*delta)^2), abs(llfdashdash[s1,s2])))
    @test abs(error) < tol
end


text:
# Alias of Model{K,Kₑ,Kₘ, K², K₊₁,K₊₂, AK₊₁,AK₊₂} for full 3D Maxwell's equations.
# The type parameters AK₊₁ and AK₊₂ specify device-specific arrays types (e.g., CuArray) and
# are user-defined in the constructor.
const ModelFull{AK₊₁,AK₊₂} = Model{3,3,3, 9, 4,5, AK₊₁,AK₊₂}

# Convenience constructor
function ModelFull(grid::Grid; Atype::Type=Array)
    cmpₛ = SInt(1,2,3)  # shapes in 3D space
    cmpₑ = SInt(1,2,3)  # E-field with x-, y-, and y-components
    cmpₘ = SInt(1,2,3)  # H-field with x-, y-, and y-components

    return ModelFull{Atype{ComplexF,4},Atype{ComplexF,5}}(;grid, cmpₛ, cmpₑ, cmpₘ)
end

# Assign the material parameters on the grid and smooth them.
function calc_matparams!(mdl::ModelFull)
    # Take necessary fields of mdl.
    grid = mdl.grid

    isbloch = grid.isbloch
    N = grid.N
    l = grid.l
    σ = grid.σ

    lg = grid.ghosted.l
    τl = grid.ghosted.τl
    ∆τ = grid.ghosted.∆τ

    boundft = mdl.boundft
    gₑ = ft2gt.(EE,boundft)
    gₘ = ft2gt.(HH,boundft)

    εarr = mdl.εarr
    µarr = mdl.μarr

    oind2shp = mdl.oind2shp
    oind2εind = mdl.oind2εind
    oind2μind = mdl.oind2μind

    cmpₑ = mdl.cmpₑ
    cmpₘ = mdl.cmpₘ

    εind2ε = sub_pind2matprm(mdl.εind2ε, cmpₑ)
    μind2μ = sub_pind2matprm(mdl.μind2μ, cmpₘ)

    # Create temporary storages.
    εxx_oind3d = create_oind_array(N)
    εyy_oind3d = create_oind_array(N)
    εzz_oind3d = create_oind_array(N)
    εoo_oind3d = create_oind_array(N)

    μxx_oind3d = create_oind_array(N)
    μyy_oind3d = create_oind_array(N)
    μzz_oind3d = create_oind_array(N)
    μoo_oind3d = create_oind_array(N)

    # Assign the material parameters.
    assign_param!(εarr, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2εind, εind2ε, gₑ, τl, isbloch)  # diagonal entries of ε tensors
    assign_param!(εarr, tuple(μoo_oind3d), oind2shp, oind2εind, εind2ε, gₑ, τl, isbloch)  # off-diagonal entries of ε tensors
    assign_param!(μarr, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2μind, μind2μ, gₘ, τl, isbloch)  # diagonal entries of μ tensors
    assign_param!(μarr, tuple(εoo_oind3d), oind2shp, oind2μind, μind2μ, gₘ, τl, isbloch)  # off-diagonal entries of μ tensors

    # Smooth the material parameters.
    smooth_param!(εarr, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2εind, εind2ε, gₑ, l, lg, σ, ∆τ, mdl.ise˔shp)  # diagonal entries of ε tensors
    smooth_param!(εarr, tuple(εoo_oind3d), oind2shp, oind2εind, εind2ε, gₑ, l, lg, σ, ∆τ, mdl.ise˔shp)  # off-diagonal entries of ε tensors
    smooth_param!(μarr, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2μind, μind2μ, gₘ, l, lg, σ, ∆τ, mdl.ish˔shp)  # diagonal entries of μ tensors
    smooth_param!(μarr, tuple(μoo_oind3d), oind2shp, oind2μind, μind2μ, gₘ, l, lg, σ, ∆τ, mdl.ish˔shp)  # off-diagonal entries of μ tensors

    return nothing
end


text:
using Ferrite
using Ferrite: getnbasefunctions, CellIterator
using TopOpt.TrussTopOptProblems.TrussTopOptProblems: compute_local_axes, getE, getA
using TopOpt.TrussTopOptProblems: truss_reinit!

"""
    get_truss_Kσs(problem::TrussProblem{xdim, TT}, u, cellvalues) where {xdim, TT}

Compute the geometric stiffness matrix for **truss elements** (axial bar element, no bending or torsion). Matrix formulation defined in eq (3) and (4) in
https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf

This matrix formulation is equivalent to the one used in
    M. Kočvara, “On the modelling and solving of the truss design problem with global stability constraints,” Struct Multidisc Optim, vol. 23, no. 3, pp. 189–203, Apr. 2002, doi: 10/br35mf.

# Arguments
`problem` = truss topopt problem struct
`u` = deformation vector (solved from first-order linear elasticity)
`cellvalues` = Ferrite cellvalues of the truss system

# Outputs
`Kσs` = a list of 2*ndim x 2*ndim element geometric stiffness matrix (in global cooridnate)
"""
function get_truss_Kσs(problem::TrussProblem{xdim,TT}, u, cellvalues) where {xdim,TT}
    Es = getE(problem)
    As = getA(problem)
    dh = problem.ch.dh

    # usually ndof_pc = xdim * n_basefuncs
    ndof_pc = ndofs_per_cell(dh)
    n_basefuncs = getnbasefunctions(cellvalues)
    @assert ndof_pc == xdim * n_basefuncs "$ndof_pc, $n_basefuncs"
    @assert n_basefuncs == 2

    global_dofs = zeros(Int, ndof_pc)
    Kσs = [zeros(TT, ndof_pc, ndof_pc) for i in 1:getncells(dh.grid)]
    Kσ_e = zeros(TT, ndof_pc, ndof_pc)

    for (cellidx, cell) in enumerate(CellIterator(dh))
        truss_reinit!(cellvalues, cell, As[cellidx])
        # get cell's dof's global dof indices, i.e. CC_a^e
        celldofs!(global_dofs, dh, cellidx)
        E = Es[cellidx]
        A = As[cellidx]
        L = norm(cell.coords[1] - cell.coords[2])
        R = compute_local_axes(cell.coords[1], cell.coords[2])
        # local axial projection operator (local axis transformation)
        γ = vcat(-R[:, 1], R[:, 1])
        u_cell = @view u[global_dofs]
        # compute axial force: first-order approx of bar force
        # better approx would be: EA/L * (u3-u1 + 1/(2*L0)*(u4-u2)^2) = EA/L * (γ'*u + 1/2*(δ'*u)^2)
        # see: https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf
        q_cell = E * A / L * (γ' * u_cell)
        Kσ_e .= 0
        for i in 2:size(R, 2)
            δ = vcat(-R[:, i], R[:, i])
            # @assert δ' * γ ≈ 0
            Kσ_e .+= δ * δ'
        end
        Kσ_e .*= q_cell / L
        Kσs[cellidx] .= Kσ_e
    end
    return Kσs
end

"""
    buckling(problem::TrussProblem{xdim, T}, ginfo, einfo, vars=ones(T, getncells(getdh(problem).grid)), xmin = T(0.0); u=undef) where {xdim, T}

Assembly global geometric stiffness matrix of the given truss problem.

# Arguments
`problem` = truss topopt problem
`ginfo` = solver.globalinfo
`einfo` = solver.elementinfo
`vars` = (Optional) design variable values, default to all ones
`xmin` = (Optional) min x value, default to 0.0
`u` = (Optional) given displacement vector, if specified as undef, a linear solve will be performed on the first-order K to get u

# Outputs
`K` = global first-order stiffness matrix (same as ginfo.K)
`Kσ` = global geometric stiffness matrix
"""
function buckling(
    problem::TrussProblem{xdim,T},
    ginfo,
    einfo,
    vars=ones(T, getncells(getdh(problem).grid)),
    xmin=T(0.0);
    u=undef,
) where {xdim,T}
    dh = problem.ch.dh
    black = problem.black
    white = problem.white
    varind = problem.varind # variable index from cell index

    if u === undef
        u = ginfo.K \ ginfo.f
    end
    Kσs = get_truss_Kσs(problem, u, einfo.cellvalues)
    # TODO replace
    # assemble_k = TopOpt.AssembleK(problem)
    # Kσ = assemble_k(Kσs)

    Kσ = deepcopy(ginfo.K)

    if Kσ isa Symmetric
        Kσ.data.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(Kσ.data, T[], Int[], Int[])
    else
        Kσ.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(Kσ, T[], Int[], Int[])
    end

    # * assemble global geometric stiffness matrix
    global_dofs = zeros(Int, ndofs_per_cell(dh))
    Kσ_e = zeros(T, size(Kσs[1]))
    celliteratortype = CellIterator{typeof(dh).parameters...}
    celliterator::celliteratortype = CellIterator(dh)
    TK = eltype(Kσs)
    for (i, cell) in enumerate(celliterator)
        celldofs!(global_dofs, dh, i)
        Kσ_e = TK isa Symmetric ? Kσs[i].data : Kσs[i]
        if black[i]
            Ferrite.assemble!(assembler, global_dofs, Kσ_e)
        elseif white[i]
            # if PENALTY_BEFORE_INTERPOLATION
            px = xmin
            # else
            #     px = penalty(xmin)
            # end
            Kσ_e = px * Kσ_e
            Ferrite.assemble!(assembler, global_dofs, Kσ_e)
        else
            px = vars[varind[i]]
            # if PENALTY_BEFORE_INTERPOLATION
            # px = density(penalty(vars[varind[i]]), xmin)
            # else
            # px = penalty(density(vars[varind[i]], xmin))
            # end
            Kσ_e = px * Kσ_e
            Ferrite.assemble!(assembler, global_dofs, Kσ_e)
        end
    end

    #* apply boundary condition
    Kσ = apply_boundary_with_zerodiag!(Kσ, problem.ch)

    return ginfo.K, Kσ
end


text:
# Julia wrapper for header: /usr/include/scip/scip_debug.h
# Automatically generated using Clang.jl wrap_c

function SCIPenableDebugSol(scip)
    ccall((:SCIPenableDebugSol, libscip), Cvoid, (Ptr{SCIP_},), scip)
end

function SCIPdisableDebugSol(scip)
    ccall((:SCIPdisableDebugSol, libscip), Cvoid, (Ptr{SCIP_},), scip)
end


text:
######################################################################
# ergodic_manager/r2.jl
#
# handles stuff needed for ergodicity
######################################################################

export ErgodicManagerR2

# kpixl is a (K+1 x bins) matrix, each entry storing cos(k*pi*x / L)
#  this assumes some discretization
"""
`ErgodicManager(d::Domain, phi::Matrix{Float64}, K::Int)`

`ErgodicManager(example_name::String; K::Int=5, bins::Int=100)`

Valid `example_name` entries are:

* "single gaussian"
* "double gaussian"
"""
mutable struct ErgodicManagerR2 <: ErgodicManager
	domain::Domain				# spatial domain
	K::Int						# number of Fourier coefficients
	phi::Matrix{Float64}		# spatial distribution
	phik::Matrix{Float64}		# distribution's Fourier coefficients

	# constant regardless of phi (depend on k1,k2)
	Lambda::Matrix{Float64}
	hk::Matrix{Float64}

	# to speed up computation
	kpixl::Matrix{Float64}
	kpiyl::Matrix{Float64}

	function ErgodicManagerR2(d::Domain, phi::MF, K::Int=5)
		em = new()
		em.domain = deepcopy(d)
		em.K = K
		em.hk = zeros(K+1,K+1)
		em.phi = deepcopy(phi)
		em.phik = zeros(K+1,K+1)
		em.Lambda = zeros(K+1,K+1)
		em.kpixl = zeros(K+1, d.cells[1])
		em.kpiyl = zeros(K+1, d.cells[2])

		Lambda!(em)
		kpixl!(em)
		hk!(em)
		decompose!(em)
		return em
	end

	function ErgodicManagerR2(d::Domain, K::Int=5)
		em = new()
		em.domain = deepcopy(d)
		em.K = K
		em.hk = zeros(K+1,K+1)
		#em.phi = deepcopy(phi)
		em.phik = zeros(K+1,K+1)
		em.Lambda = zeros(K+1,K+1)
		em.kpixl = zeros(K+1, d.cells[1])
		em.kpiyl = zeros(K+1, d.cells[2])

		Lambda!(em)
		kpixl!(em)
		hk!(em)
		#decompose!(em)
		return em
	end
end

# fills each entry Lambda[k1,k2] in the Lambda matrix
function Lambda!(em::ErgodicManagerR2)
    for k1 = 0:em.K, k2 = 0:em.K
        den = (1.0 + k1*k1 + k2*k2) ^ 1.5
        em.Lambda[k1+1, k2+1] = 1.0 / den
    end
end

# TODO: check if the row/col ordering of this is julia-efficient
# TODO: also, I could probably speed up by iterating over k first?
function kpixl!(em::ErgodicManagerR2)
	Lx = em.domain.lengths[1]
	xmin = x_min(em)
	for xi = 1:x_cells(em)
		# TODO: seems like I just cancel xmin out. Ignore it maybe?
		x = xmin + (xi-0.5)*x_size(em)
		for k = 0:em.K
			em.kpixl[k+1,xi] = cos(k*pi*(x-xmin) / Lx)
		end
	end

	Ly = em.domain.lengths[2]
	ymin = y_min(em)
	for yi = 1:y_cells(em)
		y = ymin + (yi-0.5)*y_size(em)
		for k = 0:em.K
			em.kpiyl[k+1,yi] = cos(k*pi*(y-ymin) / Ly)
		end
	end
end

# generates the hk coefficients for the ergodic manager
# these coefficients only need to be computed once
function hk!(em::ErgodicManagerR2)
    for k1 = 0:em.K, k2 = 0:em.K
        em.hk[k1+1,k2+1] = hk_ij(em, k1, k2)
    end
end

# computes the coefficients for a specific value of k1 and k2
# called by hk!
function hk_ij(em::ErgodicManagerR2, k1::Int, k2::Int)
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		cx2 = cx * cx
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			val += cx2 * cy * cy * em.domain.cell_size
		end
	end

	return sqrt(val)
end

######################################################################
# Computing Fourier coefficients
######################################################################
# update the Fourier coefficients based on some distribution
# Here, I assume it is discrete, but I should change this...
# TODO: maybe I should do some bounds checking?
function decompose!(em::ErgodicManagerR2, d::Matrix{Float64})
    for k1 = 0:em.K, k2 = 0:em.K
        em.phik[k1+1,k2+1] = phi_ij(em, k1, k2, d)
    end
    em.phi = d
end

# iterate over the state space
function phi_ij(em::ErgodicManagerR2, k1::Int, k2::Int, d::Matrix{Float64})
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			val += d[xi,yi] * cx * cy * em.domain.cell_size
		end
	end
	return val / em.hk[k1+1,k2+1]
end

"""
`decompose(em, traj::VVF)`

Decomposes a set of positions into a set of `ck` Fourier coefficients.

"""
function decompose(em::ErgodicManagerR2, traj::VVF)

	# trajectory is really of length N+1
	N = length(traj)-1

	# create matrix to hold trajectory's Fourier coefficients
	ck = zeros(em.K+1, em.K+1)

	# lengths of each dimension
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]

	# minimum values in each dimension
	xmin = x_min(em)
	ymin = y_min(em)

	for k1 = 0:em.K
		kpiL1 = k1 * pi / Lx
		for k2 = 0:em.K
			kpiL2 = k2 * pi / Ly
			hk = em.hk[k1+1, k2+1]
			fk_sum = 0.0
			# now loop over time
			for n = 0:N
				xn = traj[n + 1]
				c1 = cos(kpiL1 * (xn[1]-xmin))
				c2 = cos(kpiL2 * (xn[2]-ymin))
				fk_sum += c1*c2
			end
			ck[k1+1, k2+1] = fk_sum / (hk * (N+1))
		end
	end
	return ck
end

# reconstructs from Fourier coefficients in ck
function reconstruct(em::ErgodicManagerR2, ck::Matrix{Float64})
	# iterate over all bins
	v = em.domain.cell_size
	vals = zeros(x_cells(em), y_cells(em))

	for xi = 1:x_cells(em)
		x = x_min(em) + (xi-0.5)*x_size(em)
		for yi = 1:y_cells(em)
			y = y_min(em) + (yi-0.5)*y_size(em)
			for k1 = 0:em.K
				cx = em.kpixl[k1+1,xi]
				for k2 = 0:em.K
					cy = em.kpixl[k2+1,yi]
					vals[xi,yi] += ck[k1+1,k2+1]*cx*cy/em.hk[k1+1,k2+1]
				end
			end
		end
	end
	return vals
end


text:
using Soss, Random
Random.seed!(1);

m = @model X begin
    n = size(X,1)
    k = size(X,2)
    w ~ Normal(0,1) |> iid(k)
    Xw = X * w
    y ~ For(n) do j
        Normal(Xw[j], 0.1)
    end
end;

X = randn(5,2)
y = rand(m(X=X)).y

post = dynamicHMC(m(X=X), (y=y,));
particles(post)

pred = predict(m(X=X), post);

y_rep = particles(pred).y

y_rep .< y


text:
function __init__pubmed()
    DEPNAME = "PubMed"
    LINK = "https://github.com/kimiyoung/planetoid/raw/master/data"
    DOCS = "https://github.com/kimiyoung/planetoid"
    DATA = "ind.pubmed." .* ["x", "y", "tx", "allx", "ty", "ally", "graph", "test.index"]

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $DOCS
        """,
        map(x -> "$LINK/$x", DATA),
        "0b8bf8e80564611b540655e9cbb8c5900dd3728d4ababe0b990b6f27144bd76c",  # if checksum omitted, will be generated by DataDeps
        # post_fetch_method = unpack
    ))
end

"""
    PubMed(; dir=nothing, reverse_edges=true)

The PubMed citation network dataset from Ref. [1].
Nodes represent documents and edges represent citation links.
The dataset is designed for the node classification task.
The task is to predict the category of certain paper.
The dataset is retrieved from Ref. [2].

# References

[1]: [Deep Gaussian Embedding of Graphs: Unsupervised Inductive Learning via Ranking](https://arxiv.org/abs/1707.03815)
[2]: [Planetoid](https://github.com/kimiyoung/planetoid)

"""
struct PubMed <: AbstractDataset
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
end

function PubMed(; dir=nothing, reverse_edges=true)
    metadata, g = read_planetoid_data("PubMed", dir=dir, reverse_edges=reverse_edges)
    return PubMed(metadata, [g])
end

Base.length(d::PubMed) = length(d.graphs)
Base.getindex(d::PubMed) = d.graphs
Base.getindex(d::PubMed, i) = getindex(d.graphs, i)


text:
export ReactiveMPBackend

using TupleTools

struct ReactiveMPBackend end

function write_argument_guard(::ReactiveMPBackend, argument::Symbol)
    return :(@assert !($argument isa ReactiveMP.AbstractVariable) "It is not allowed to pass AbstractVariable objects to a model definition arguments. ConstVariables should be passed as their raw values.")
end

function write_randomvar_expression(::ReactiveMPBackend, model, varexp, arguments, kwarguments)
    return :($varexp = ReactiveMP.randomvar($model, $(GraphPPL.fquote(varexp)), $(arguments...); $(kwarguments...)))
end

function write_datavar_expression(::ReactiveMPBackend, model, varexpr, type, arguments, kwarguments)
    return :($varexpr = ReactiveMP.datavar($model, $(GraphPPL.fquote(varexpr)), ReactiveMP.PointMass{ GraphPPL.ensure_type($(type)) }, $(arguments...); $(kwarguments...)))
end

function write_constvar_expression(::ReactiveMPBackend, model, varexpr, arguments, kwarguments)
    return :($varexpr = ReactiveMP.constvar($model, $(GraphPPL.fquote(varexpr)), $(arguments...); $(kwarguments...)))
end

function write_as_variable(::ReactiveMPBackend, model, varexpr)
    return :(ReactiveMP.as_variable($model, $varexpr))
end

function write_make_node_expression(::ReactiveMPBackend, model, fform, variables, options, nodeexpr, varexpr)
    return :($nodeexpr = ReactiveMP.make_node($model, $fform, $varexpr, $(variables...); $(options...)))
end

function write_autovar_make_node_expression(::ReactiveMPBackend, model, fform, variables, options, nodeexpr, varexpr, autovarid)
    return :(($nodeexpr, $varexpr) = ReactiveMP.make_node($model, $fform, ReactiveMP.AutoVar($(GraphPPL.fquote(autovarid))), $(variables...); $(options...)))
end

function write_node_options(::ReactiveMPBackend, fform, variables, options)
    return map(options) do option

        # Factorisation constraint option
        if @capture(option, q = fconstraint_)
            return write_fconstraint_option(fform, variables, fconstraint)
        elseif @capture(option, meta = fmeta_)
            return write_meta_option(fform, fmeta)
        elseif @capture(option, pipeline = fpipeline_)
            return write_pipeline_option(fform, fpipeline)
        end

        error("Unknown option '$option' for '$fform' node")
    end
end

# Meta helper functions

function write_meta_option(fform, fmeta)
    return :(meta = $fmeta)
end

# Pipeline helper functions

function write_pipeline_option(fform, fpipeline)
    if @capture(fpipeline, +(stages__))
        return :(pipeline = +($(map(stage -> write_pipeline_stage(fform, stage), stages)...)))
    else
        return :(pipeline = $(write_pipeline_stage(fform, fpipeline)))
    end
end

function write_pipeline_stage(fform, stage)
    if @capture(stage, Default())
        return :(ReactiveMP.DefaultFunctionalDependencies())
    elseif @capture(stage, RequireEverything())
        return :(ReactiveMP.RequireEverythingFunctionalDependencies())
    elseif @capture(stage, RequireInbound(args__))

        specs = map(args) do arg
            if @capture(arg, name_Symbol)
                return (name, :nothing)
            elseif @capture(arg, name_Symbol = dist_)
                return (name, dist)
            else
                error("Invalid arg specification in node's WithInbound dependencies list: $(arg). Should be either `name` or `name = initial` expression")
            end
        end

        indices  = Expr(:tuple, map(s -> :(ReactiveMP.interface_get_index(Val{ $(GraphPPL.fquote(fform)) }, Val{ $(GraphPPL.fquote(first(s))) })), specs)...)
        initials = Expr(:tuple, map(s -> :($(last(s))), specs)...)

        return :(RequireInboundFunctionalDependencies($indices, $initials))
    else
        return stage
    end
end

# Factorisation constraint helper functions

function factorisation_replace_var_name(varnames, arg::Expr)
    index = findfirst(==(arg), varnames)
    return index === nothing ? error("Invalid factorisation argument: $arg. $arg should be available within tilde expression") : index
end

function factorisation_replace_var_name(varnames, arg::Symbol)
    index = findfirst(==(arg), varnames)
    return index === nothing ? arg : index
end

function factorisation_name_to_index(form, name)
    return :(ReactiveMP.interface_get_index(Val{ $(GraphPPL.fquote(form)) }, Val{ ReactiveMP.interface_get_name(Val{ $(GraphPPL.fquote(form)) }, Val{ $(GraphPPL.fquote(name)) }) }))
end

function check_uniqueness(t)
    return TupleTools.minimum(TupleTools.diff(TupleTools.sort(TupleTools.flatten(t)))) > 0
end

function sorted_factorisation(t)
    subfactorisations = map(TupleTools.sort, t)
    firstindices      = map(first, subfactorisations)
    staticlength      = TupleTools.StaticLength(length(firstindices))
    withindices       = ntuple(i -> (i, firstindices[i]), staticlength)
    permutation       = map(first, TupleTools.sort(withindices; by = last))
    return ntuple(i -> subfactorisations[permutation[i]], staticlength)
end

function write_fconstraint_option(form, variables, fconstraint)
    if @capture(fconstraint, (*(factors__)) | (q(names__)))
        factors = factors === nothing ? [ fconstraint ] : factors

        indexed = map(factors) do factor
            @capture(factor, q(names__)) || error("Invalid factorisation constraint: $factor")
            return map((n) -> GraphPPL.factorisation_name_to_index(form, n), map((n) -> GraphPPL.factorisation_replace_var_name(variables, n), names))
        end

        factorisation = Expr(:tuple, map(f -> Expr(:tuple, f...), indexed)...)
        errorstr = """Invalid factorisation constraint: ($fconstraint). Arguments are not unique, check node's interface names and model specification variable names."""

        return :(factorisation = GraphPPL.check_uniqueness($factorisation) ? GraphPPL.sorted_factorisation($factorisation) : error($errorstr))
    elseif @capture(fconstraint, MeanField())
        return :(factorisation = MeanField())
    elseif @capture(fconstraint, FullFactorisation())
        return :(factorisation = FullFactorisation())
    else
        error("Invalid factorisation constraint: $fconstraint")
    end
end

##

function write_randomvar_options(::ReactiveMPBackend, variable, options)
    return map(options) do option
        @capture(option, name_Symbol = value_) || error("Invalid variable options specification: $option. Should be in a form of 'name = value'")
        return option
    end
end

function write_constvar_options(::ReactiveMPBackend, variable, options)
    return map(options) do option
        @capture(option, name_Symbol = value_) || error("Invalid variable options specification: $option. Should be in a form of 'name = value'")
        return option
    end
end

function write_datavar_options(::ReactiveMPBackend, variable, options)
    return map(options) do option
        @capture(option, name_Symbol = value_) || error("Invalid variable options specification: $option. Should be in a form of 'name = value'")
        return option
    end
end


text:
"""
	grdhisteq(cmd0::String="", arg1=nothing, kwargs...)

Find the data values which divide a given grid file into patches of equal area. One common use of
grdhisteq is in a kind of histogram equalization of an image.

Full option list at [`grdhisteq`]($(GMTdoc)grdhisteq.html)

Parameters
----------

- **D** | **dump** :: [Type => Str or []]

    Dump level information to file, or standard output if no file is provided.
    ($(GMTdoc)grdhisteq.html#d)
- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdhisteq(....) form.
    ($(GMTdoc)grdhisteq.html#g)
- **N** | **gaussian** :: [Type => Number or []]

    Gaussian output.
    ($(GMTdoc)grdhisteq.html#n)
- **Q** | **quadratic** :: [Type => Bool]

    Quadratic output. Selects quadratic histogram equalization. [Default is linear].
    ($(GMTdoc)grdhisteq.html#q)
- $(GMT.opt_R)
- $(GMT.opt_V)
"""
function grdhisteq(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdhisteq", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)			# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :V_params])
	cmd  = parse_these_opts(cmd, d, [[:C :n_cells], [:D :dump], [:G :outgrid], [:N :gaussian], [:Q :quadratic]])
	common_grd(d, cmd0, cmd, "grdhisteq ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdhisteq(arg1, cmd0::String=""; kw...) = grdhisteq(cmd0, arg1; kw...)


text:
const obsnames_ising = ["Energy", "Energy^2", "Specific Heat",
                        "|Magnetization|", "Magnetization^2", "Magnetization^4",
                        "Susceptibility", "Connected Susceptibility", "Binder Ratio",
                       ]
const obsnames_clock = ["Energy", "Energy^2", "Specific Heat",
                        "|Magnetization x|", "Magnetization x^2", "Magnetization x^4",
                        "Susceptibility x", "Connected Susceptibility x", "Binder Ratio x",
                        "|Magnetization y|", "Magnetization y^2", "Magnetization y^4",
                        "Susceptibility y", "Connected Susceptibility y", "Binder Ratio y",
                        "|Magnetization|", "|Magnetization|^2", "|Magnetization|^4",
                        "Susceptibility", "Connected Susceptibility", "Binder Ratio",
                       ]
function loaddata(filename, obsnames)
    Ts = zeros(0)
    res = Parameter(n=>zeros(0) for n in obsnames)
    for line in eachline(filename)
        words = split(line)
        push!(Ts, parse(Float64, words[1]))
        for (i,n) in enumerate(obsnames)
            push!(res[n], parse(Float64, words[i+1]))
        end
    end
    return Ts, res
end

function parse_filename(filename, ::Union{Type{Ising}, Type{XY}})
    m = match(r"^J_([\d.-]*)__N_([\d.-]*).dat$", filename)
    if m == nothing
        return nothing
    end
    p = Parameter()
    p["J"] = parse(Float64, m.captures[1])
    p["N"] = parse(Int, m.captures[2])
    return p
end

function parse_filename(filename, ::Union{Type{Potts}, Type{Clock}})
    m = match(r"^Q_(\d*)__J_([\d.-]*)__N_(\d*).dat$", filename)
    if m == nothing
        return nothing
    end
    p = Parameter()
    p["Q"] = parse(Int, m.captures[1])
    p["J"] = parse(Float64, m.captures[2])
    p["N"] = parse(Int, m.captures[3])
    return p
end

@testset "$modelstr" for (modelstr, pnames, obsnames) in [("Ising", ("J", "N"), obsnames_ising),
                                                          ("Potts", ("Q", "J", "N"), obsnames_ising),
                                                          ("Clock", ("Q", "J", "N"), obsnames_clock),
                                                          ("XY", ("J", "N"), obsnames_clock),
                                                         ]
    model = eval(Symbol(modelstr))
    for filename in readdir(joinpath("ref", modelstr))
        p = parse_filename(filename, model)
        if p == nothing
            continue
        end
        testname = ""
        for pname in pnames
            testname *= "$(pname)=$(p[pname]) "
        end
        @testset "$testname" begin
            Ts, exacts = loaddata(joinpath("ref", modelstr, filename),obsnames)
            nT = length(Ts)
            p["Model"] = model
            p["Lattice"] = fully_connected_lattice
            p["J"] = p["J"] / p["N"]
            p["MCS"] = MCS
            p["Thermalization"] = Therm
            @testset "$(upstr)" for upstr in ("local_update!",
                                              "SW_update!",
                                              "Wolff_update!",
                                             )
                p["Update Method"] = eval(Symbol(upstr))
                res1 = []
                res2 = []
                for i in 1:nT
                    p["T"] = Ts[i]
                    p["Seed"] = SEED
                    push!(res1, runMC(p))
                    p["Seed"] = SEED2
                    push!(res2, runMC(p))
                end
                @testset "$n" for n in obsnames
                    for i in 1:nT
                        T = Ts[i]
                        exact = exacts[n][i]
                        r1 = res1[i]
                        r2 = res2[i]
                        ## single MC test may fail.
                        mc1 = r1[n]
                        mc2 = r2[n]
                        ex = exact
                        if !(p_value(mc1, exact) > alpha/nT || p_value(mc2, exact) > alpha/nT)
                            @show T
                            @show exact
                            @show mc1, p_value(mc1, exact)
                            @show mc2, p_value(mc2, exact)
                        end
                        @test p_value(mc1, exact) > alpha/nT || p_value(mc2, exact) > alpha/nT
                    end
                end
            end
        end
    end
end


text:
using PtFEM

data = Dict(
# Plane(ndim, nst, nxe, nye, nip, direction=(:x|:y), finite_element(nod, nodof), axisymmetric=false)
# Plane(ndim, nst, nxe, nye, nip, direction=:z, finite_element(nod, nodof), axisymmetric=true)
# Plane(ndim, nst, nre, nze, nip, direction=:r, finite_element(nod, nodof), axisymmetric=false)
  :struc_el => Plane(2, 6, 1, 5, 4, :r, Quadrilateral(8, 3), false),
  :properties => [1.0e5 0.3;],
  :r_coords => [0.0, 0.5],
  :z_coords => [10.0, 8.0, 6.0, 4.0, 2.0, 0.0],
  :support => [
    ( 1, [1 0 1]),
    ( 4, [1 0 1]),
    ( 6, [1 0 1]),
    ( 9, [1 0 1]),
    (11, [1 0 1]),
    (14, [1 0 1]),
    (16, [1 0 1]),
    (19, [1 0 1]),
    (21, [1 0 1]),
    (24, [1 0 1]),
    (26, [0 0 0]),
    (27, [0 0 0]),
    (28, [0 0 0])
    ],
  :loaded_nodes => [
    (3, [1/pi 0.0 0.0])
    ],
  :lth => 1,
  :iflag => 1,
  :chi => 0.0
)

data |> display
println()

@time m, dis_df, fm_df = p52(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:")
fm_df |> display
println()


text:
include("maximin.jl")
include("cc_maximin.jl")
include("minimax.jl")
include("custom.jl")
include("no_obj.jl")


text:
module AWSCLIProvider

using ...SimpleTasks.Types

import AWS
import SimpleTasks.Services.Bucket
import SimpleTasks.Services.CLIBucket

export Details

const AWS_KEY_PATH = "$(homedir())/.aws/"
const AWS_KEY_FILE = "$AWS_KEY_PATH/config"
const AWS_PREFIX = `s3:/`
const AWS_COMMAND = `aws s3`

Details(env::AWS.AWSEnv) =
    create_access!(env.aws_id, env.aws_seckey) &&
    CLIBucket.Provider(AWS_PREFIX, AWS_COMMAND)

Details(aws_access_key_id::AbstractString,
    aws_secret_access_key::AbstractString) =
    create_access!(aws_access_key_id, aws_secret_access_key) &&
    CLIBucket.Provider(AWS_PREFIX, AWS_COMMAND)

function create_access!(aws_access_key_id::AbstractString,
        aws_secret_access_key::AbstractString)
    if stat(AWS_KEY_FILE).size == 0
        mkpath(AWS_KEY_PATH)
        println("No access file found for aws cli, creating one!")
        file = open(AWS_KEY_FILE, "w")
        write(file, "[default]\n")
        write(file, "aws_access_key_id = $(aws_access_key_id)\n")
        write(file, "aws_secret_access_key = $(aws_secret_access_key)\n")
        close(file)
    end
    return true
end

end # module AWSCLIProvider

max_forks_repo_name: seung-lab/SimpleTasks


text:
# Download and install RADMC-3D
binary_file = "radmc-3d_v0.39_17.03.15.zip"

# Change to the src directory
cd("src")
println("Current working directory ", pwd())

# Download and extract the RADMC-3D executable
download("http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/$binary_file", binary_file)
println("Current working directory ", pwd())

# Unzipping the file
run(`unzip $binary_file`)

cd("radmc-3d/version_0.39/src")
run(`make`)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using JSON3
using HTTP

json_string = read("jsonschema_examples/beerjson-1.0.1/json/beer.json", String)
json_schema = JSON3.read(json_string)

@show keys(json_schema)
for key in (Symbol("\$schema"), Symbol("\$id"))
    if key in keys(json_schema)
        println("$key: $(json_schema[key])")
    end
end

uri = "https://raw.githubusercontent.com/beerjson/beerjson/master/json/beer.json"
@assert isvalid(HTTP.URI(uri)) "Invalid URI"
json_string = String(HTTP.get(uri).body)
json_schema = JSON3.read(json_string)

@show keys(json_schema)
for key in (Symbol("\$schema"), Symbol("\$id"))
    if key in keys(json_schema)
        println("$key: $(json_schema[key])")
    end
end


text:
import NLopt

type MLE2
reps
iterations
lb
ub
maxtime
maxeval
jitter
method
end
MLE2(lb, ub;reps=5, iterations=50, maxtime=10, maxeval=1000, jitter=1, method=:LN_COBYLA) = MLE2(reps, iterations, lb, ub, maxtime, maxeval, jitter, method)

type MLE2_DIRECT end

function wrap_model(m::GP)
	function wrapper_nll(hyp, grad)
		if any(isnan(hyp))
			return Inf
		end
		try
			m.k, m.lik = inflate_hyp(m, hyp)
			return nll(m)
		catch e
			return Inf
		end
	end
	return wrapper_nll
end

function infer!(o::MLE2, m::GP)
	init = flatten_hyp(m)
	init_best = (nll(m), init)
	best = (nll(m), init, :ORIGINAL)
	opt = NLopt.Opt(o.method, length(init))
	NLopt.lower_bounds!(opt, o.lb)
	NLopt.upper_bounds!(opt, o.ub)
	# NLopt.xtol_rel!(opt,1e-4)
	NLopt.ftol_rel!(opt, 1e-12)
	NLopt.maxtime!(opt, o.maxtime)
	NLopt.maxeval!(opt, o.maxeval)
	NLopt.min_objective!(opt, wrap_model(m))
	for r = 1:o.reps
		this_init = init + randn(size(init)) * o.jitter
		this_init = min(init, o.ub)
		this_init = max(init, o.lb)
		candidate = NLopt.optimize(opt, this_init)
		if candidate[1] < best[1]
			best = candidate
		end
	end
	println("Hyperparameter inf: $(best[3])")
	m.k, m.lik = inflate_hyp(m, best[2])
end

function infer!(o::MLE2, m::FITCModel)
	init = flatten_hyp(m)
	init_best = (nll(m), init)
	best = (nll(m), init, :ORIGINAL)
	opt = NLopt.Opt(o.method, length(init))
	NLopt.lower_bounds!(opt, o.lb[1:length(init)])
	NLopt.upper_bounds!(opt, o.ub[1:length(init)])
	NLopt.xtol_rel!(opt,1e-4)
	NLopt.ftol_rel!(opt, 1e-12)
	NLopt.maxtime!(opt, o.maxtime)
	NLopt.maxeval!(opt, o.maxeval)
	NLopt.min_objective!(opt, wrap_model(m))
	for r = 1:o.reps
		this_init = init + randn(size(init)) * o.jitter
		this_init = min(init, o.ub[1:length(init)])
		this_init = max(init, o.lb[1:length(init)])
		candidate = NLopt.optimize(opt, this_init)
		if candidate[1] < best[1]
			best = candidate
		end
	end
	m.k, m.lik = inflate_hyp(m, best[2][1:length(flatten_hyp(m))])
end


text:
using Gtk

btn = GtkButton("Start")
sp = GtkSpinner()
ent = GtkEntry()

grid = GtkGrid()
grid[1,1] = btn
grid[2,1] = sp
grid[1:2,2] = ent

id = addprocs(1)[1]

signal_connect(btn, "clicked") do widget
 start(sp)
 @Gtk.sigatom begin
   @async begin
    s = @fetchfrom id begin
      sleep(4)
      return "I am back"
    end
    @Gtk.sigatom begin
      stop(sp)
      set_gtk_property!(ent,:text,s)
    end
  end
 end
end

win = GtkWindow(grid, "Progress Bar", 200, 200)
showall(win)

max_stars_repo_name: ashwani-rathee/Chatbots

max_issues_repo_name: ashwani-rathee/Chatbots

max_forks_repo_name: ashwani-rathee/Chatbots


text:
using Discord
include("base.jl")

function fin(client::Client,m::Message,damage::Int)
    global  management_dict

    channel_id_str = string(m.channel_id)

    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"予約確認板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    boss_number = string(management_dict[channel_id_str]["boss_number"])
    reserve_list =  management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"]

    dict_number = 0

    #予約ずみかどうか事前に凸宣言をしていたかどうかをチェックする
    for i = 1:length(reserve_list)
        if reserve_list[i]["id"] == m.author.id
            if !reserve_list[i]["attack"]
                return reply(client,m,"凸宣言がされておりません．")
            end
            dict_number = i
            break
        end
    end

    #予約していなかっtら怒る
    if dict_number == 0
        return reply(client,m,"凸予約がされていません.")
    end

    remain_hp = management_dict[channel_id_str]["reserve"][boss_number]["remain_hp"] - damage

    if remain_hp <0
        return reply(client,m,"ボスを討伐した場合には`.la`コマンドを使用して下さい")
    end

    management_dict[channel_id_str]["reserve"][boss_number]["remain_hp"] = remain_hp
    deleteat!(management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"],dict_number)
    management_dict[channel_id_str]["reserve"][boss_number]["plan_remain_hp"] = calc_plan_remain_hp(management_dict[channel_id_str]["reserve"][boss_number])

    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    #残凸状況を修正
    if management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["remain"] > 0
        management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["remain"] -= 1
        management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["over"] = false
        management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["over_content"]  = ""
        remain_attack_message_content = make_remain_attack_message_content(management_dict[channel_id_str]["remain_attack"],management_dict[channel_id_str]["today"])
        edit_message(client,management_dict[channel_id_str]["remain_attack_channel"],management_dict[channel_id_str]["remain_attack_message"],content=remain_attack_message_content)
    end

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')

end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle